Diel Transcriptional Oscillations of a Plastid Antiporter Reflect Increased Resilience of Thalassiosira pseudonana in Elevated CO2

Acidification of the ocean due to high atmospheric CO2 levels may increase the resilience of diatoms causing dramatic shifts in abiotic and biotic cycles with lasting implications on marine ecosystems. Here, we report a potential bioindicator of a shift in the resilience of a coastal and centric model diatom Thalassiosira pseudonana under elevated CO2. Specifically, we have discovered, through EGFP-tagging, a plastid membrane localized putative Na+(K+)/H+ antiporter that is significantly upregulated at >800 ppm CO2, with a potentially important role in maintaining pH homeostasis. Notably, transcript abundance of this antiporter gene was relatively low and constant over the diel cycle under contemporary CO2 conditions. In future acidified oceanic conditions, dramatic oscillation with >10-fold change between nighttime (high) and daytime (low) transcript abundances of the antiporter was associated with increased resilience of T. pseudonana. By analyzing metatranscriptomic data from the Tara Oceans project, we demonstrate that phylogenetically diverse diatoms express homologs of this antiporter across the globe. We propose that the differential between night- and daytime transcript levels of the antiporter could serve as a bioindicator of a shift in the resilience of diatoms in response to high CO2 conditions in marine environments.


INTRODUCTION
Diatoms (Bacillariophyta) represent a diverse group of marine phytoplankton widely distributed across the globe, from coastal habitats to the open oceans (Malviya et al., 2016). They account for 40% of the marine primary productivity (20% globally) (Nelson and Gordon, 1982;Falkowski et al., 1998;Field et al., 1998) and facilitate carbon export to the deep ocean through the biological carbon pump Smetacek, 1999). Diatoms can form large seasonal blooms when conditions are favorable, generating vast amounts of biomass for marine food webs (Platt et al., 2003). Diatoms have brought stability to key ecosystems (Armbrust, 2009), adjusting to changes in nutrient availability, such as nitrogen (N) (Levitan et al., 2015), phosphorus (P) (Brembu et al., 2017), silicon (Si), iron (Fe) (Marchetti et al., 2009;Smith et al., 2016) and trace elements, as well as physiochemical conditions, including light (Wilhelm et al., 2014;Smith et al., 2016;Lepetit et al., 2017) levels (photosynthetically active radiation; PAR) (Aguirre et al., 2018), temperature (Armbrust, 2009), salinity, pH, and CO 2 Raven et al., 2019), over daily, seasonal, and decadal cycles. The robustness of diatoms to combinatorial environmental changes has been investigated through controlled photobioreactor studies, demonstrating how the model diatom Thalassiosira pseudonana regulates >5,400 genes (>40% of all genes in its genome) to adopt distinct physiologic states matched to relevant environmental conditions (light, dark, nutrient replete, and nutrient deplete) (Ashworth et al., 2013). Hundreds of genes encoding sensory, signaling, and regulatory functions execute environment-specific gene regulatory programs to drive transitions of T. pseudonana across these physiological states (Ashworth et al., 2013). For example, previous work implicated the photoreceptor AUREO1c in regulation of photosynthesis and sugar metabolism during dawn and exponential growth; while the light-responsive transcription factor bZIP7a/PAS was implicated in regulating other transcription factors. Furthermore, bZIP24b was implicated in regulating nucleosome and chromatin organization at high CO 2 conditions. Under increasing CO 2 concentrations, T. pseudonana remodels chromatin and transcriptionally downregulates photosynthesis, respiration, and carbon-concentrating mechanisms (CCMs) (Ashworth et al., 2013;Hennon et al., 2015), which are utilized by many diatoms (Badger et al., 1998;Crawfurd et al., 2011) to grow in the CO 2 -limited oceans of today (<10 -30 µmol CO 2 ) (Hopkinson et al., 2011;Reinfelder, 2011).
The reduced need to concentrate carbon at an elevated CO 2 level may alleviate the need for CCMs (Hopkinson et al., 2011;Matsuda et al., 2017), which might allow diatoms to reallocate resources to manage stressful environments. Indeed, a study designed to quantify resilience of T. pseudonana discovered that at higher CO 2 levels, the diatom was able to withstand and recover from incrementally larger amounts of stress (i.e., ultraviolet radiation or UVR), while transitioning through ecologically-relevant physiologic states (i.e., light, dark, nutrientreplete, and nutrient-limited) (Valenzuela et al., 2018). Thus, elevated CO 2 conditions of ocean acidification may stabilize diatom populations under NO 3 limiting conditions and may lead to the expansion of their geographic distribution and ecological niches. However, predicting how diatoms and other phytoplankton will fare in future oceans is complicated by nonlinear consequences of combined increases in atmospheric CO 2 , global temperatures, and nutrient availabilities (Feely et al., 2004). For instance, while ocean acidification may improve diatom resilience when they are not Fe limited (Valenzuela et al., 2018), diatom productivity will suffer in environments where Fe bioavailability may decline with increasing temperature and acidification (Shi et al., 2010). The combined effects of light and CO 2 on marine productivity can also vary depending on the season, latitude, depth, and the photo-acclimation capabilities of phytoplankton (Jones, 1998;Ihnken et al., 2011;Heiden et al., 2016;Heydarizadeh et al., 2017), with light intensity setting the upper limit of productivity (Jones, 1998;Platt et al., 2003). Hence, understanding how diatom resilience and productivity will change in response to shifts in the complex interplay of biotic and abiotic factors due to ocean acidification is essential to predict species succession or niche expansions in seasonally acidified environments, such as upwelling systems and nutrientlimited zones (Bruland et al., 2001;Cohen et al., 2017;Godhe and Rynearson, 2017;Valenzuela et al., 2018).
In order to investigate how varying levels of light will affect diatoms in high and low carbon conditions, we have generated and analyzed the transcriptome responses of T. pseudonana to combinatorial changes in light intensity (75 µmol photons·m −2 ·s −1 ; hereafter low light or "LL" and 300 µmol photons·m −2 ·s −1 ; hereafter high light or "HL") and CO 2 level (400 ppm; hereafter low carbon or "LC" and 800 ppm; hereafter high carbon or "HC"). In addition to recapitulating the large-scale physiological state-transitions of T. pseudonana, the transcriptome analysis discovered the transcriptional modulation of a putative antiporter that might play an important role in pH homeostasis during the transition between day and night under elevated CO 2 . In particular, we have discovered a dramatic switch in the dynamical day/night expression patterns of a putative Na + (K + )/H + antiporter, which indicates a shift in diatom response to high CO 2 conditions regardless of changes in other environmental conditions (e.g., PAR, UVR, nutrient limitation). Notably, this dynamic diel oscillation of the antiporter transcript levels in response to high CO 2 was reproduced in a mesocosm experiment, demonstrating ecological relevance of this pattern in a complex ocean environment. Through EGFP-tagging, we determined that the putative antiporter is localized to the plastid membrane, where it might function to help maintain an optimal pH for the activity of enzymes related to carbon fixation, photosynthesis, transport, and other metabolic processes (Launay et al., 2020). This is consistent with the expectation that seawater acidification will alter electrochemical gradients that govern pH gradients across membranes (Taylor et al., 2012). While diatoms typically maintain a neutral pH in the cytoplasm (Taylor et al., 2012;Goldman et al., 2017), the pH within the plastid is higher (pH ∼8) during the light period, and decreases during the dark (pH ∼7) (Launay et al., 2020). This pH homeostasis is essential for plastid function, and therefore, the expression dynamics of the antiporter captures the non-linear consequences of combined changes in multiple environmental factors that influence the overall physiological state of diatoms, potentially making the antiporter a useful composite bioindicator of diatom resilience. By data-mining metatranscriptome sequences from the Tara Oceans project (Carradec et al., 2018;Villar et al., 2018), we determined that phylogenetically diverse diatoms transcribe homologs of this putative antiporter, making it a candidate bioindicator. Importantly, we postulate that it is the differential between night-and daytime transcript expression that will ultimately reveal the shift in diatom resilience in response to elevated CO 2 associated ocean acidification.

Experimental Design
Thalassiosira pseudonana (CCMP 1335) was grown at two CO 2 levels (400 and 800 ppm) and two light levels, undersaturated (75 µmol photons·m −2 ·s −1 ) and saturated for growth (300 µmol photons·m −2 ·s −1 ) (Sobrino and Neale, 2007;Torres et al., 2013). The cells were pre-acclimated to the respective CO 2 and light conditions for approximately 4 days (∼4 generations) prior to bottle reactor inoculation. RNA samples were collected during the exponential and late-exponential phases of growth. Each experiment was conducted in triplicate, with three serial replicates per combination of CO 2 and light levels. In total, 24 microarrays were collected and analyzed (8 conditions × 3 replicates each = 24 microarrays).

Batch Culture Growth, Monitoring, Sampling, and Analysis
Batch cultures of T. pseudonana were grown in triplicate in 5 L glass bottles containing enriched artificial seawater (ESAW) media with a nitrate concentration of ∼100 µM. All cultures were grown under a continuous light regime. During all growth experiments, pH and cell counts were monitored twice daily (∼ 12 h apart). The pH was measured using a spectrophotometrically calibrated pH probe and cell counts were quantified immediately after sampling using a hemocytometer. The cultures were equilibrated to 400 ppm or elevated CO 2 (800 ppm) by bubbling and stirring (50 rpm). The two CO 2 levels were generated by first stripping the water from the laboratory air with DU-CAL (Drierite company, Xenia, OH, United States) and CO 2 with Sodasorb (Divers Supply Inc., Gretna, LA, United States) and then using mass-flow controllers (model GFC-17, Aalborg, Orangeburg, NY, United States) to mix with 99.99% pure CO 2 at exact ratios (Praxair, Danbury, CT, United States). The CO 2 concentration of the gas was measured with a CO 2 analyzer (model S151, Qubit Systems, Kingston, ON, Canada). The gas flow (0.4 LPM: liters per minute) on each bottle was controlled with a flow meter (model 6A0101BV-AB, Dakota Instruments Inc., Orangeburn, NY, United States). The CO 2 -air mixture was filtered through a 0.2 µm Millipore filter before flowing into the culturing bottles to maintain axenic cultures. Total dissolved inorganic carbon (DIC) was measured using an Apollo SciTech (DIC analyzer) Model AS-C3 and Li-COR LI-7000 CO 2 /H 2 O analyzer from filtered (0.2 µm) samples. The flasks were inoculated with 5 × 10 4 cells/mL of acclimated axenic T. pseudonana cells. Photosynthetic efficiency (maximal PSII quantum yield, F v /F m ) was obtained from the maximal fluorescence (F m ) and variable fluorescence (F v ) using the Phyto-PAM (Pulse Amplitude Modulated) Phytoplankton Analyzer (Walz). Variable fluorescence was calculated from F m to F o , where F o is the fluorescence yield when cells are dark acclimated. Nutrients (NO 3 , PO 4 , and Si) were sampled twice daily and measured at the Chemical Services Lab (University of Washington).

RNA Extraction and Expression Microarrays
RNA extraction and expression microarrays were performed according to Ashworth et al. (2013). Total RNA was harvested from a total of 3 × 10 7 cells from each reactor during (a) exponential and (b) late-exponential growth, using the mirVANA kit from Invitrogen. A total of 24 RNA samples were collected, based on the following experimental matrix: {400 ppm, 800 ppm} × {unsaturated light, saturated light} × {exponential, late-exponential} × {3 parallel replicates} × {3 serial replicates}. RNA in each sample was reverse transcribed and labeled using the Agilent Quick Amp 56 Labeling Kit. The resulting Cy3-labeled cRNA was hybridized along with a uniform Cy5-labeled reference sample (as an internal standard) to custom gene-specific 8 × 60 k oligonucleotide Agilent arrays (Agilent design ID: 037886; GEO platform ID: GPL18682). Arrays were scanned using an Agilent two-color array scanner.

Microarray Data Processing and Expression Analysis
Two-color microarray data were processed and normalized using the limma package in R (Ritchie et al., 2015). Relative log 2 expression ratios were calculated in relation to the Cy5-labeled universal standard reference RNA sample. The effects of three dichotomous experimental factors (CO 2 level, light level, and growth phase) on the expression changes of all genes during the experiments were assessed using a three-way ANOVA in R, with multiple hypothesis correction for significance assessment using the Benjamini-Hochberg method. Functional enrichment analysis was performed using the g:Profiler toolset (Raudvere et al., 2019). These microarray expression datasets are available in the GEO database under the accession number(s): GSE57737 Series C (n = 24).

RNA-seq Processing and Expression Analysis From Diel Cycle Photobioreactor Experiments
Transcriptional analysis was performed as described in Valenzuela et al. (2018). In that study, T. pseudonana was cultured in photobioreactors under 12:12-h (light:dark) cycles with nitrate limitation at HC (1,000 ppm CO 2 ) and LC (300 ppm CO 2 ) to force the transition across four physiological states (light, dark, early, and late-phase exponential phase). Each growth cycle represented a "stage" and at the end of the first stage an aliquot was transferred to fresh growth medium to initiate the next stage. Starting in stage 2, at midday, 0.5 mW/cm 2 UVR was applied for 1 h and incrementally increased by +0.5 mW/cm 2 in subsequent stages. During each stage, transcriptomes from both HC and LC conditions across the four principal states were analyzed for signatures of resilience. A total of 56 samples were processed for RNA-sequencing. Read files were cleaned with Trimmomatic (Bolger et al., 2014) and mapped with aligner STAR (Dobin et al., 2013) to T. pseudonana reference genome (Armbrust et al., 2004;Bowler et al., 2008). Gene expression levels were quantified as fragments per kilobase of transcript per million mapped reads (FPKM) using Cufflinks (Trapnell et al., 2012). These RNA datasets are publicly available from the National Center for Biotechnology (NCBI) Sequence Read Archive (SRA), accession code PRJNA38016.

RNA-seq Processing and Expression Analysis From Diel Cycle Mesocosm Experiments
Indoor mesocosm experiments were conducted at Ocean Acidification Environmental Laboratory (OAEL), Friday Harbor Laboratories (FHL), to study the effects of CO 2induced acidification over the diel cycle on T. pseudonana in seawater pumped directly from the adjacent Puget Sound. Four independent mesocosm tanks (9 L) were set up with continuous flow (10-12 mL/min) of filtered (0.2 µm) and UV sterilized seawater combined with a CO 2 blending system to control pH and simulate mid-century (pH 7.9) and acidified oceanic conditions (pH 7.6) in duplicate. T. pseudonana was acclimated to FHL seawater supplemented with F/2 nutrients for 24 h prior to inoculation into mesocosm reservoirs, which were outfitted with custom enclosures to simulate a 12:12 light:dark diel cycle (275 µmol photons·m −2 ·s −1 ) and supplemented with additional nutrients to support eutrophic growth over 48 h. Mesocosm water samples were collected twice daily during the diel cycle (n = 16) and then processed for transcriptional analysis. Read files were cleaned with Trimmomatic (Bolger et al., 2014) followed by pseudoalignment to a T. pseudonana reference genome (Armbrust et al., 2004;Bowler et al., 2008) and read quantification using Kallisto v.46.2 (Bray et al., 2016). Differential gene expression analysis was performed using DESeq2 package v1.28.1 (Love et al., 2014) in R with shrinking LFC estimates applied (Zhu et al., 2018). The mesocosm RNA datasets are publicly available from the NCBI GEO database under the accession number: GSE168812 (n = 16).

Vector Construction and Thalassiosira pseudonana Transformation
The expression vector used for protein localization was constructed as described in Shrestha and Hildebrand (2015), using MultiSite Gateway Technology (Life Technologies) (Shrestha and Hildebrand, 2015). The CDS for gene 262258 was amplified with PCR primers flanked with attB sites to insert the fragment into the pDONR221 entry vector. The entry clone facilitated the insertion of the 262258 fragment upstream and in-frame with enhanced green fluorescent protein (EGFP) in the final destination vector pMHL_079, for constitutive transcription under the control of a fucoxanthin chlorophyll a/c binding protein promoter (FCPp). A list of primers used is listed in Supplementary Table 1. The FCPp-262258-EGFP-FCPt expression vector was co-transformed with pMHL_009 expressing the nat1 gene under the control of the acetyl-coenzyme A carboxylase promoter, which confers resistance to the antibiotic nourseothricin (NAT). The protocol for T. pseudonana transformation was performed as described in Shrestha and Hildebrand (2015). Following transformation, 250 mL of culture was added to 50 ml of ASW medium containing NAT at a final concentration of 100 µg/mL. The liquid culture was grown to exponential phase (5 × 10 5 cells/mL) and the top 5% of EGFP-expressing cells (10,000 cells total) were sorted by FACS into fresh ASW medium using the BD Influx Cell Sorter. Exponential-phase cells from the sorted culture were used for fluorescence microscopy.

Confocal Fluorescence Microscopy
Thalassiosira pseudonana cells expressing the EGFP-tagged transmembrane antiporter protein (262258) were imaged with a Zeiss Axio Observer Z1 inverted microscope equipped with an ApoTome and a Zeiss AxioCam MRm camera (Carl Zeiss Microimaging, Inc., Thornwood, NY, United States). The Zeiss #16 filter set was used to image chlorophyll autofluorescence (Excitation BP 485/20 nm, Dichromatic mirror FT 510 nm, Emission LP 515 nm) and the Zeiss #38HE filter set was used to image EGFP fluorescence (Excitation BP 470/40 nm, Dichromatic mirror FT 495 nm, Emission BP 525/50 nm). Nonfluorescent images were taken using differential interference contrast (DIC). Z-stack images through the cells were acquired with 63x/1.4 oil immersion Plan-Apochromat objective and a 1.6x optovar module. Images were processed using Axiovision 4.7.2 software. The average composition of three Z-stack images was done using ImageJ software (Schneider et al., 2012).

Diatom Antiporter Sequences in the Environment
The Ocean Gene Atlas web server 1 (Villar et al., 2018) was used to query the Marine Atlas of Tara Oceans Unigenes version 1 Metatranscriptomes (MATOUv1 + T) (Carradec et al., 2018) against the amino acid sequence of the putative antiporter (262258). Biological samples from the planktonic eukaryote communities of the MATOUv1 + T were fractionated in four size ranges (0.8-5, 5-20, 20-180, and 180-2,000 µm). We removed any hits that did not align to the taxonomic phylum Bacillariophyta and had a bitscore below 300 (largest non-zero e-value 7.69e −107 ) to ensure sequences were diatom specific with high confidence.

RESULTS
In order to dissect the combinatorial effects of light and CO 2 , we cultured T. pseudonana in four conditions representing combinations of low (75 µmol photons·m −2 ·s −1 , LL) and high (300 µmol photons·m −2 ·s −1 , HL) light levels, and current (400 ppm, LC) and projected (800 ppm, HC) CO 2 concentrations. We observed similar growth rates across all combinations of light and CO 2 conditions (HL:HC, HL:LC, LL:HC, and LL:LC; Figure 1A). Cultures in HC achieved a significantly higher carrying capacity (p-value = 0.0027, t-test), as expected (Valenzuela et al., 2018). Nutrient (N, P, and Si) drawdown was similar across all four treatments, with complete utilization of phosphate within 24 h (Supplementary Figure 1). Even after depletion of phosphate in the growth medium, the continued growth of diatoms is likely attributable to intracellular polyphosphate reserves (Dyhrman et al., 2012;Valenzuela et al., 2012Valenzuela et al., , 2013Supplementary Figure 1B). The complete utilization of all macronutrients by ∼60 h (i.e., when RNA was sampled for transcriptomic analysis) coincided with a decline in growth rate, a decrease in F v /F m [i.e., maximal fluorescence (F m ) and variable fluorescence (F v )], and transition of all cultures into late-exponential phase. Consistent with a previous report (Li et al., 2014), F v /F m declined faster during HL growth regardless of CO 2 level ( Figure 1B). DIC decreased proportionally with an increase in biomass, and DIC levels were restored during transition to stationary phase when nutrients were depleted, and photosynthetic efficiency (F v /F m ) decreased ( Figure 1C).
Uninterrupted photosynthesis due to culturing in continuous light conditions did not allow for external pH recovery in the media, which typically occurs during nighttime (i.e., no photosynthesis) (Ashworth et al., 2013;Valenzuela et al., 2018). As a consequence, pH continued to increase in all cultures through late-exponential phase. In general, the level and rate of pH change was higher (∼0.52/day in both LC and HC) in HL conditions compared to LL (∼0.39/day in LC and ∼0.23/day in HC), irrespective of CO 2 levels. HC cultures maintained a relatively lower pH throughout the experiment. The highest pH (∼8.75) was observed in the HL:LC condition; by contrast, LL:HC reached a maximum pH of ∼8.2. Notably, the pH decreased as the cultures transitioned into late-exponential phase, which happened earlier in HL (∼36 h), relative to LL (∼48 h) grown cultures, and coincided with a decrease in PE. As expected, changes in DIC and pH of the culture media were inversely correlated (Figures 1C,D). In other words, the carbon cycling dynamics during the exponential phase of growth in these continuous light experiments resembled that of phytoplankton during a daytime bloom, where drawdown of nutrients can be attributed to growth of diatoms, simultaneously causing increases in daily pH (>1 pH unit) during the light periods (Wallace et al., 2014;Raven et al., 2019).
We performed whole transcriptome analysis on cells harvested during the mid-(∼24 h) and late-exponential (∼60 h) phases of growth to investigate molecular and systems-level changes associated with the response of diatoms to combinatorial changes in light, and CO 2 . Using three-way analysis of variation (ANOVA), we identified over 1,400 genes with greater than a twofold change in transcript abundance (False Discovery Rate: FDR ≤ 0.001; Table 1 and Supplementary Data File 1). This analysis also allowed attribution of differential regulation of genes due to changes in each of the three variables in this experiment (growth phase, light intensity, and CO 2 level). As expected, transition from mid-to late-exponential growth phase accounted for the majority of transcriptional variation, samples could also be distinguished by carbon condition ( Supplementary  Figures 2A,B). The 705 genes that were upregulated in midexponential phase to support growth and biomass production were subsequently downregulated in late-exponential phase when nutrients were depleted. Functions enriched within these  (7) Significant transcriptome changes were determined with three-way ANOVA (CO 2 level, light level, and growth phase). Number of gene models with transcriptional changes that were significantly associated with CO 2 level, light level, and growth phase. Bold numbers in parentheses indicate the number genes with significant changes greater than twofold. 1 Benjamini-Hochberg was applied to all ANOVA p-values with n = 15,059.
genes included photosynthesis (e.g., light-harvesting complexes and photosystems I and II reaction center proteins) and carbon fixation [e.g., RuBisCO large subunit (bd2088), glyceraldehyde 3-phosphate dehydrogenase (31383), and ribose-5-phosphate isomerase (32252)] (Supplementary Figures 3A,B). By contrast, 676 genes were downregulated at mid-exponential phase (upregulated at late-exponential phase) and associated with functions related to ribosome metabolism and nutrient transport (Supplementary Figures 3C,D). In particular, the upregulation of nitrate, phosphate, and silicic acid transporters (269274, 27414, 261414, and 268895) was consistent with the depletion of these micronutrients in stationary phase (Mock et al., 2008;Valenzuela et al., 2012;Ashworth et al., 2013;Lomana et al., 2015). The coordinated upregulation of multiple RNA helicases, a ubiquitin-conjugating enzyme (36045), and an endonuclease (bd410) suggests that T. pseudonana may utilize these genes to recycle N and P from nucleotides (Bourgeois et al., 2016;Brembu et al., 2017;Alipanah et al., 2018) during nutrient starvation-induced growth arrest. In sum, transition of cells into late-exponential phase was associated with downregulation of energy production, upregulation of RNA metabolism, and nutrient transport, suggesting a shift toward energy conservation and nutrient scavenging. Interestingly, only 15 genes had a greater than a twofold differential regulation (FDR ≤ 0.001) as a direct consequence of the change in light intensity (6 upregulated in HL and 9 in LL) (Supplementary Data File 1). Similarly, differential regulation of just 31 genes could be attributed to CO 2 level change, with 14 genes upregulated in HC (Figure 2A) and 17 genes upregulated in LC (Table 1 and Supplementary Data File 1). While there was no significant enrichment of known functions among the differentially regulated genes responsive to CO 2 , the analysis did recapitulate high CO 2 -responsive downregulation of three genes [6528, 10360 (a putative transcriptional regulator: Tp_bZIP24a) and 233 (a putative carbonic hydrase)], that were previously identified as part of a CCM and photorespiration sub-cluster Hennon et al., 2015; Supplementary  Table 2 and Supplementary Data File 1). Among the five most upregulated genes in HC, three had unknown functions, and two had putative transporter functional annotations (Supplementary Table 2). We investigated the robustness of CO 2 -responsive upregulation of the five genes by analyzing an independent RNA-seq dataset from a study designed to investigate the resilience of T. pseudonana at two CO 2 levels (300 and 1,000 ppm) over the diel cycle (Valenzuela et al., 2018). In brief, T. pseudonana cultures were grown under a 12:12-h light:dark cycle in nitrogen-limited growth medium, and an aliquot from late-exponential phase was inoculated into fresh medium -each growth cycle representing a "stage." Each stage represents a batch culture of T. pseudonana growing in 12:12 light:dark diel cycle over different growth phases. In this setup, the T. pseudonana culture transitions across four ecologically relevant physiological states associated with daytime, nighttime, nutrient replete (early growth phase) and nutrient limited (late-exponential growth phase) conditions. Starting with stage 2, the cultures were subjected to 1 h of an ecologically relevant dose of UVR (0.5 mW/cm 2 ) at midday; and the UVR dosage was incrementally increased by +0.5 mW/cm 2 in each subsequent stage (i.e., 1.0 mW/cm 2 in stage 3 and 1.5 mW/cm 2 in stage 4), until cultures collapsed. This study demonstrated that T. pseudonana cultures were more resilient in HC conditions, which manifested in their ability to recover from greater number of exposures to incrementally higher doses of UVR. Analysis of transcript level changes and absolute expression (FPKM: Fragments Per Kilobase of transcript per Million mapped reads) of the five genes in RNAseq data from all three stages of this experiment demonstrated that only the putative antiporter (262258) was significantly differentially expressed with a relatively high absolute transcript abundance. The other four genes had relatively low transcript abundance (Supplementary Figure 4) and only gene 5647 was significantly upregulated in HC conditions. The transcript abundance of the putative antiporter was higher at 1,000 ppm CO 2 by up to 52-fold (mean = 129.3, max = 251.4 FPKM) relative to its transcript level in 300 ppm CO 2 (mean = 2.49, max = 8.29 FPKM) (Figure 2B and Supplementary Data File 2). The difference in the transcriptional patterns of the other four transcripts across the two studies suggests that the regulation of these genes may be sensitive to constant light or diel conditions, making them unsuitable for use as bioindicators of the CO 2 -responsive change in the physiological state and resilience of the diatom. However, the comparative analysis of T. pseudonana transcriptome responses to high and low CO 2 in different contexts [i.e., with constant light at two different intensities (microarray data generated in this study)] and over diel conditions with intermittent UVR stress [RNA-seq data from Valenzuela et al. (2018)], demonstrated that diel changes in abundance of the antiporter transcript might serve as a robust predictor of a shift in the resilience of diatoms in the natural environment at elevated CO 2 , irrespective of changes in other factors (i.e., light intensity, diel cycle, nutrients, or UVR stress).  Valenzuela et al., 2018). (C) Putative antiporter sequences are widely distributed across the globe, particularly at coastal regions. The area of each circle reflects the abundance values for the sequence at that site (larger area = higher abundance). The color scheme reflects the pH of the environment the sample was collected. (D) Transcript abundance of the antiporter homolog sequences from the query of the MATOUv1 + T dataset were observed at high and low pH conditions, which was negatively correlated with CO 2 conditions. Log scales were used for the y-axis for C and D.
We further investigated whether the putative antiporter and its homologs in other diatoms were also expressed in the natural environment by querying the Marine Atlas of Tara Oceans Unigenes Eukaryotic Metatranscriptomic database (MATOUv1 + T) using the Ocean Gene Atlas webserver (Villar et al., 2018). We ascertained diatom-specificity of antiporter homologs by imposing stringent filters to eliminate potential homologs that did not align to the taxonomic phylum Bacillariophyta or had a bitscore below 300 (equated to the largest non-zero e-value 7.69e −107 with sequence identity between 62.8 and 82 percent). Altogether 66 diatom-specific (Phylum Bacillariophyta) antiporter sequences were identified with high confidence, within these homologs, some were assigned to the lower taxonomical levels (class: Coscinodiscophyceae = 7, order: Chaetocerotales = 6), including six homologs to the genus Chaetoceros, the most abundant diatom genus in the world (Rines and Theriot, 2003;Malviya et al., 2016). The diatom antiporter transcripts were discovered in ocean samples across the globe, exclusively in the surface water layer, and mostly in coastal regions (Figure 2C), which is typical of planktonic diatoms that inhabit nutrient-rich environments (Smetacek, 2012). We report the relative abundance of each putative antiporter transcript at a specific location as the percentage of all mapped reads in the associated samples (RPKM: Reads Per Kilobase covered per Million of mapped reads) (Villar et al., 2018). The wide variation in abundance of transcripts in samples of similar pH, including samples from the same location, demonstrated that absolute transcript level of the antiporter "alone" may have no predictive value, vis-à-vis the status of CO 2 -responsiveness of diatoms at a given location ( Figure 2D). Therefore, using the same diel dataset from Valenzuela et al. (2018), described above, we investigated whether dynamic changes in the expression level of the antiporter transcript (262258) over the diel cycle was diagnostic of a shift in the resilience of T. pseudonana in FIGURE 3 | Expression patterns of the putative antiporter transcript over the diel cycle in LC and HC, and its plastid localization. (A) Line plots of absolute expression demonstrate the dynamic daytime (white bar; low expression) and nighttime (gray bar; high expression) oscillations in HC (red) but not LC (blue) conditions, with maximal expression in nighttime and late-exponential phase of growth, when nutrients were depleted. Each stage represents a growth cycle (early to late-exponential phase), under a 12:12-h diel cycle, thus cells transition between four principal states (i.e., the light and dark phase of both early and late-exponential growth). UVR exposures are represented by the purple vertical line plots above the line plots of transcript abundance. During stage 1 there was no UVR, however, in stage 2 the cultures received 1 h of an ecologically relevant dose of UVR (0.5 mW/cm 2 ) at midday to stress the cultures. UVR doses were incrementally increased by 0.5 mW/cm 2 in each subsequent stage. future oceanic conditions. The antiporter expression level did not change significantly over the diel cycle at 300 ppm CO 2 and had low absolute expression (<8.3 FPKM). Strikingly, at 1,000 ppm CO 2 , the transcript level of the permease oscillated dramatically with low daytime levels (19.2 FPKM) and high nighttime levels (251.4 FPKM) ( Figure 3A). These findings suggest that the increased resilience of T. pseudonana in laboratory simulated future oceanic conditions with elevated CO 2 is reflected in the differential in daytime and nighttime expression levels of the antiporter.
To confirm whether an elevated CO 2 level that could drive seawater below pH 7.9 would trigger dynamic transcriptional oscillations of the antiporter in the natural environment, we designed a controlled mesocosm study at Ocean Acidification Environmental Laboratory (OAEL) in Friday Harbor Laboratories (FHL). Briefly, four independent mesocosms were set up with continuous flow (10-12 mL/min) of filtered seawater from the Puget Sound, two in mid-century (pH 7.9) conditions and two in future acidified oceanic (pH 7.6) conditions. Mesocosm reservoirs were supplemented with nutrients and inoculated with T. pseudonana acclimated in FHL seawater. Mesocosms were outfitted with custom enclosures to simulate a 12:12 light:dark diel cycle. Water samples collected from the mesocosms were then processed for transcriptional analysis as described in section "Methods". Altogether, 14 genes were significantly differentially expressed (FC ≥ 2 and p-value ≤ 0.001) between HC (n = 8) and LC (n = 8) mesocosms. Strikingly, the putative antiporter was the most significantly upregulated gene in HC (∼10.1 Fold Change, p-value 6.5 × 10 −16 ) (Supplementary Data File  3), corroborating that while the antiporter was expressed in low abundance over the diel cycle in LC, its overall abundance increased dramatically in HC, especially at nighttime, reproducing the oscillation pattern observed in the laboratory experiment ( Figure 3B). Furthermore, we sought to understand the putative role of the antiporter that may explain how its dynamic expression patterns over the diel cycle is associated with the increased resilience of T. pseudonana in HC. Secondary structure, organelle localization (TargetP) (Emanuelsson et al., 2000), and functional domain classifications (GO0005215: transporter activity; KOG2639: inorganic ion transport), suggested that the putative antiporter (262258) is a plastid-localized (Supplementary Table 3), 12transmembrane Na + /H + antiporter (Wang et al., 2018a,b;Supplementary Figure 5). In order to confirm its localization to the plastid membrane, we chromosomally integrated an engineered DNA construct for constitutive expression of a C-terminal EGFP-tagged antiporter (see section "Methods" for details) to interrogate its intracellular distribution using confocal fluorescence microscopy ( Figure 3C and Supplementary Figure 6). Analysis of the fluorescence and confocal images confirmed the predicted localization of the transmembrane antiporter to the plastid membrane. Notably, we discovered that the ortholog of the Na + /H + antiporter in the model pennate diatom, Phaeodactylum tricornutum (Phatr3_EG02645) is also a plastid localized transmembrane protein with putative inorganic ion transport activity. Using an RNA-seq dataset from an independent study investigating the diel cycling dynamics of P. tricornutum under Fe limitation, we confirmed that the antiporter homolog was also upregulated in the dark with low expression in light (Smith et al., 2016; Supplementary  Table 3 and Supplementary Figure 7). Together, the ∼5-fold upregulation (FDR < 1.1×10 −4 ) and dynamic day/night oscillations in HC, localization to the plastid membrane, and putative function of the transmembrane protein suggest the Na + /H + antiporter may potentially play an important role in maintaining electrochemical gradients across the plastid membrane in the context of diel cycles and elevated CO 2 levels (Hennon et al., 2015). If so, then the function of the antiporter might become even more essential in an acidified ocean, further supporting its potential utility as a bioindicator of a shift in diatom resilience in the natural environment.

DISCUSSION
In this study, we observed that > 1,400 genes were differentially regulated (FDR < 0.001) during the growth of T. pseudonana in combinations of perturbations in light (HL and LL), and CO 2 levels (HC and LC), which simulated daytime bloomlike dynamics of rapid growth in nutrient-replete conditions followed by a transition to a more quiescent state in nutrientdeplete stationary phase. The discovery that relatively few genes were differentially expressed exclusively because of changes in individual factors demonstrated the degree of combinatorial regulation of diatom response by light, CO 2 , and nutrient availability. It is noteworthy that the differential expression of the putative antiporter was directly reflective of the response to CO 2 , across different intensities of light, UV stress, and nutrient availability. We postulate that this is because the gene is potentially an intracellular pH-sensitive Na + (K + )/H + antiporter with a direct role in maintaining pH homeostasis across the plastid membrane, similar to the chx23 Na + (K + )/H + antiporter in Arabidopsis thaliana, which by itself has been shown to maintain pH homeostasis in the chloroplast stroma (Song et al., 2004). In fact, A. thaliana also uses a K + /H + antiporter to modify the proton motive force in the thylakoid lumen during the light phase (Kunz et al., 2014;Finazzi et al., 2015).
During dark to light transitions, in many plants as well as diatoms, the internal pH of the plastid can increase from ∼7 to ∼8 (Colman and Rotatore, 1995), which is the optimal pH for the activity of most plastid enzymes, including RuBisCO and Calvin cycle enzymes (Goldman et al., 2017). Given that the cytosolic pH of diatoms in contemporary CO 2 levels is in the range of 7.0-7.5 (Taylor et al., 2012;Höhner et al., 2016;Goldman et al., 2017;Launay et al., 2020), active efflux of H + or counter exchange of ions by means of antiporters may be necessary to achieve a pH of 8 inside the plastid. Studies on T. weissflogii have demonstrated that acidification of seawater proportionally decreases intracellular pH causing internal acidosis in diatoms (Goldman et al., 2017). While the increased passive diffusion of CO 2 will be sufficient to support photosynthesis without a need for CCMs, which are downregulated in inverse proportion to rise in CO 2 levels (Hennon et al., 2015), this will also result in increased intracellular acidosis. In fact, by 2100, increased atmospheric CO 2 is expected to drive down the pH of seawater from 8.1 to 7.6 (Caldeira and Wickett, 2005), dramatically altering intracellular electrochemical gradients in diatoms (Taylor et al., 2012;Goldman et al., 2017) and elevating the need for active mechanisms to restore cytosolic/plastid pH homeostasis (Launay et al., 2020), particularly at nighttime. Our findings suggest that it is in this context that the proton efflux activity of the antiporter will be important to restore pH homeostasis within the plastid, especially at nighttime when the cells experience greater intracellular acidosis. This hypothesis that the role of the putative Na + (K + )/H + antiporter might become even more important in regulating pH homeostasis during diel transitions in an acidified ocean, is demonstrated by a switch in the expression pattern of the antiporter from relatively constant, low abundance across the day/night cycle at 300 ppm CO 2 to dramatic oscillations at 1,000 ppm CO 2 , with considerably higher abundance at nighttime.
Absolute transcript level of the antiporter at a given location is proportional to the relative abundance of diatoms in that location and their response to their environmental context (e.g., nutrient availability, CO 2 level, light condition, etc.). Therefore, absolute transcript level alone at any given time is by itself not indicative of the CO 2 -response of the diatom. Instead, we postulate that the differential between nighttime and daytime expression of the antiporter is a better proxy for the overall physiologic response of the diatom to non-linear consequences of changes in multiple factors associated with ocean acidification as we have demonstrated through laboratory photobioreactor and mesocosm experiments. Importantly, the switch in expression pattern of the putative antiporter from constitutively low levels to dynamic oscillations over the diel cycle may inform when a diatom has or is likely to transition to a high CO 2 physiological state associated with increased resilience, potentially driving a shift in ecosystem dynamics of a marine habitat.

DATA AVAILABILITY STATEMENT
All microarray expression datasets are available in the GEO database under the accession number(s): GSE57737 Series C (n = 24). In addition, all RNA-seq datasets (n = 56) are publicly available from the National Center for Biotechnology (NCBI) Sequence Read Archive (SRA), accession code PRJNA38016 and previously published in Valenzuela et al., 2018. RNA-seq datasets for mesocosm datasets are publicly available from NCBI's GEO database under accession number: GSE168812 (n = 16).

AUTHOR CONTRIBUTIONS
JA, AC, EA, MO, and NB designed the experiments, which were performed at EA's laboratory. JA, AC, and MO performed all growth experiments, assays, RNA extraction, and preparation. JA and JV performed microarray and transcriptome analyses. JV performed meta-transcriptomic and additional RNA-seq analyses. MH provided the laboratory and resources for RA to construct all vectors and perform fluorescence microscopy. JV, JA, AC, RA, EA, MO, and NB contributed to the discussion of results. JV, MO, JA, and NB wrote the manuscript. All authors contributed to the article and approved the submitted version.

ACKNOWLEDGMENTS
We would like to thank MH who unfortunately passed away during this study. MH touched many lives both personally and scientifically, he was a compassionate leader and gracious colleague.