Representative Diatom and Coccolithophore Species Exhibit Divergent Responses throughout Simulated Upwelling Cycles

C. decipiens and E. huxleyi displayed clear physiological responses to the different UCBC conditions under iron-replete conditions (Fig. 2, black bars). Both reduced chlorophyll a under stationary growth phases (T2 and T3) (Fig. 2A and B; see also Fig. S1 in the supplemental material). Photosynthetic efficiency (F_v:F_m) was also lower at the stationary-phase time points (T2 and T3) than at the initial exponential growth time point (T1) (Fig. 2C and D). In C. decipiens, these declines were significant (P < 0.001) (see Data Set S1 at https://doi.org/10.6084/m9.figshare.c.5315135.v1). Following the dark period, both species increased photosynthetic efficiency when the cells returned to exponential growth, which is consistent with previous field observations (7).

Cellular C:N ratios for both species also fluctuated throughout the cycle (Fig. 2E and F). Both species were growing at near-Redfield-predicted values (6.6:1) during the initial exponential growth phase and increased their C:N ratios during stationary phase following N depletion (see Fig. S2). As with chlorophyll a and F_v:F_m in C. decipiens, these differences were highly significant (P < 0.0001) (see Data Set S1 at https://doi.org/10.6084/m9.figshare.c.5315135.v1). Previous field observations of particulate C:N ratios from phytoplankton communities at depth have also shown values far exceeding the Redfield ratio, reaching as high as 60:1, although uncertainty remained as to whether these results were due to cellular modifications or simply C-rich detrital material (7, 8, 10). With a maximum C:N ratio of 20.2 ± 2.1 for C. decipiens at T3 and 29.96 ± 8.4 for E. huxleyi at T2 under the iron-replete conditions, it can be inferred that a large proportion of the high carbon relative to that of nitrogen can be attributed to cells suffering from N and/or light limitation, although in E. huxleyi, these data include inorganic carbon that is incorporated into coccoliths (Fig. S2). Similarly, under N limitation, several Thalassiosira species were shown to have C:N ratios of >20, which can be attributed to an accumulation of carbohydrates (24). Previous laboratory experiments that mimicked sinking cells with the diatom Pseudo-nitzschia australis grown in excess N relative to Si did not show as dramatic of an increase in the C:N ratio, suggesting that Si limitation, as can often occur in conjunction with iron-limited diatoms in the California coastal upwelling region, may influence these elemental stoichiometries (25, 26).

To examine gene expression profiles in both organisms throughout the UCBC, transcriptome-wide coexpression was analyzed with weighted gene coexpression network analysis (WGCNA). WGCNA resulted in 23 modules (clusters of interconnected genes) for C. decipiens and 58 modules for E. huxleyi. Significant positive correlations (Pearson; P < 0.05) between modules and all time points apart from T6 for C. decipiens were found, indicating that there is distinct and coordinated expression of certain genes within each UCBC stage for both organisms (Fig. 3 and Fig. S3). Several modules were significantly associated with the time point pairs that represented similar environmental conditions (T1 and T5; T2 and T3). Eighteen modules in C. decipiens and 40 modules in E. huxleyi had significant positive correlations with each time point. Importantly, these significant associates were largely different between time points, indicating that both organisms exhibited a high degree of transcriptional responsiveness to the different UCBC conditions.

Genes associated with each UCBC stage were identified by first examining modules that correlated with each time point and then by determining which genes within those modules were also correlated with each time point (Fig. 3; see also Data Sets S2 and S3 at https://doi.org/10.6084/m9.figshare.c.5315135.v1). By comparing KEGG orthologs (KOs) and protein family (Pfam) annotations of contigs associated at each time point, similarities or differences in expression for both organisms throughout the UCBC were examined. Overall, C. decipiens and E. huxleyi showed highly divergent responses throughout the UCBC. When solely comparing KOs and Pfams detected in both transcriptomes, 513 unique KOs and 518 unique Pfams were significantly associated at any time point in both species; however, most of these were not associated with the same time points in both organisms. Furthermore, more KOs or Pfams were significantly associated with one of the organisms rather than both at each time point (Table 1).

Broadly, this divergence in responses represents differing metabolic investments and potential niche partitioning between diatoms and coccolithophores; however, some similarities can be observed by examining specific KEGG orthologs, KEGG modules, and Pfams associated with each organism at each time point (see Data Sets S2 and S3 at https://doi.org/10.6084/m9.figshare.c.5315135.v1). In both organisms, many photosynthesis-related genes were associated with time points when light was available. Expression of glycolysis and fatty acid biosynthesis genes during exponential growth was also common. At 12 h after returning to the light, both organisms highly expressed genes associated with translation and amino acid biosynthesis, as they were in transition to resume growth.

In examining KOs found in both organisms that were correlated with different time points, differences in the timing of transcription-related gene expression were observed (Fig. 4; see Data Sets S2 and S3 at the above-mentioned URL). These include RNA polymerases and spliceosome genes to generate mature RNA as well as TRAMP and exosome complexes that can interact to degrade abnormal rRNAs and spliced-out introns (27). C. decipiens displayed elevated expression of these gene sets during stationary growth, particularly at T3, compared to that during its exponential growth phases. In an opposite pattern, E. huxleyi decreased expression of these genes during stationary phase and elevated expression during exponential growth and at T4 in response to a return to light and nutrients. WGCNA indicated that these gene sets were significantly associated with T3 in C. decipiens but with T4 in E. huxleyi (see Data Set S2 and S3 at the URL mentioned above).

These opposing patterns indicate relatively high cellular investments in transcription during differing growth phases. In C. decipiens, this higher expression corresponds to periods when conditions are not optimal for growth, i.e., a lack of light and nutrients, whereas higher expression corresponds to periods leading to or within exponential growth in E. huxleyi. This difference in timing of expression may drive transcriptional differences observed in the field, where diatoms were found to highly express, or frontload, nitrogen assimilation genes prior to upwelling conditions in contrast to other groups, including haptophytes (7). An increase in transcriptional machinery may be linked to the frontloading response by enabling an accumulation of nitrogen assimilation transcripts that ultimately assists the diatoms in outcompeting other groups by being transcriptionally proactive rather than reactive in response to the relative rapid and frequent changes in environmental conditions (28).

Frontloading of nitrogen assimilation genes and an acceleration of nitrate uptake, or the shift-up response, is hypothesized to be linked to diatom success during the upwelling portion of the UCBC (2, 5, 8). As a result, the upwelling portion of the UCBC simulation was sampled after 12 h (T4) to investigate the physiological capability for relatively rapid shift-up. By comparing T3 through T5, a view from pre-upwelling, 12 h after upwelling, and a return to exponential growth can be examined (Fig. 1).

A rapid shift-up response was observed in C. decipiens compared to that in E. huxleyi. C. decipiens returned to exponential growth within 24 h of returning to light and nutrients, whereas E. huxleyi did not respond appreciably for approximately 48 h (see Fig. S4). UCBC simulations with various dark phases suggest that the capability to respond in E. huxleyi is impacted by the length of the dark phase, with an inability to reinstate growth after 20 days in the dark (Fig. S5). Concurrently, the C:N ratio in C. decipiens rapidly declined from 20.18 ± 2.13 to 11.89 ± 0.70 within the first 12 h of returning to light and nutrients, suggesting that C. decipiens rapidly assimilated nitrate once returned to optimal growth conditions (Fig. 2E and F). After 36 h in the light, cellular C:N ratios in C. decipiens returned to the Redfield ratio. In contrast, E. huxleyi maintained a high C:N ratio (22.10 ± 4.91) after returning to light for 12 h, but the C:N ratio dramatically declined to 1.80 ± 0.40 at 3 days after returning to the light when the cells reached exponential growth (Fig. 2F and Fig. S4). This initial lack of decline indicates that E. huxleyi has a slower response to the newly available nitrate than C. decipiens; however, once the cells reach exponential growth, E. huxleyi may exhibit luxury uptake and store excess nitrate prior to an increase in cell division, as evidenced by the later decline in their C:N ratio in conjunction with the highest cellular nitrogen quotas (see Fig. S6). Altogether, these results align with the relatively rapid shift-up response of diatoms compared to that of other phytoplankton groups and show that nitrate uptake rates are likely rapidly increasing to result in the observed changes.

Results from a field-based simulated upwelling experiment suggest that frontloading of nitrogen assimilation genes is an underlying molecular mechanism behind these responses (7). Diatoms were observed to highly express nitrate assimilation genes (nitrate transporter [NRT], nitrate reductase [NR], and nitrite reductases [NiR]) prior to upwelling events in order to be transcriptionally proactive, or frontload, to rapidly assimilate nitrate and outcompete other phytoplankton groups once conditions are optimal for growth.

In comparing the relative and summed expression of nitrogen metabolism genes in C. decipiens and E. huxleyi from the UCBC simulation, C. decipiens exhibited this distinct and proactive frontloading response, whereas E. huxleyi did not (Fig. 5A and B). C. decipiens significantly upregulated the aforementioned nitrogen assimilation genes during stationary phase (P < 0.05), with highest expression at 10 days in the dark, i.e., the simulated pre-upwelling condition. WGCNA also placed these genes in the blue module, which was significantly associated with the dark stationary-phase time points (Fig. 3; Data Set S2 at https://doi.org/10.6084/m9.figshare.c.5315135.v1). Upon returning to high light and nutrients, expression of these genes significantly declined and returned to previous levels observed during exponential growth (P < 0.01) (Fig. 5B). Diatoms also possess two nitrite reductases (Nir), both of which are localized to the chloroplast; however, one utilizes ferredoxin (Fd) and the other utilizes NAD(P)H (9). Under exponential growth, Fd-Nir transcript abundance was higher than NAD(P)H-Nir; however, under stationary growth in the light and dark, NAD(P)H-Nir transcript abundance was higher, suggesting potential metabolic advantages of each nitrite reductase under these different growth conditions (Fig. 5A).

In contrast, E. huxleyi displayed an opposing expression pattern for their nitrate assimilation genes (Fig. 5A). During stationary growth, E. huxleyi significantly downregulated transcripts for nitrate transporters and nitrate reductase (P < 0.01); nitrite reductase expression remained similar across time points (Fig. 5B). Once returned to the light and nutrients (T4), E. huxleyi significantly upregulated transcripts for nitrate transporters (P < 0.0001) in response to nitrate availability and optimal growth conditions. Nitrate reductase and nitrite reductase expression also remained low compared to that of other E. huxleyi transcripts (16th to 40th percentiles) throughout the experiment; therefore, in conjunction with low cellular C:N ratios, these data support that E. huxleyi may store nitrate for some time post-upwelling.

Unlike the nitrate assimilation genes, ammonium transporters (AMT) and urea transporters (UT) displayed a similar transcriptional pattern in both organisms; however, while nitrate transporters were among the most highly expressed nitrogen metabolism genes in C. decipiens (87th to 100th percentiles), E. huxleyi highly expressed ammonium transporters (94th to 100th percentiles) despite ammonium not being added to the seawater medium. For AMT, these patterns were driven by the mostly highly expressed copies in both organisms, with some copies displaying different gene expression patterns (see Fig. S7). High expression of AMT and UT during stationary phase suggests that both organisms are tuned to exploit released ammonium and urea as nitrate was depleted, although urease expression remained low in both (4th to 21st percentiles) (Fig. 5B and C and Fig. S7). E. huxleyi has also previously been shown to prefer ammonium over nitrate (29). The high total level of expression of AMT relative to the expression of other genes in E. huxleyi and compared to that of AMT expression in C. decipiens corroborates this preference for ammonium and implies that the two species exhibit different nitrogen preferences.

The divergent responses between the two species with frontloading observed in the diatom is increasingly apparent when examining the nitrogen metabolism genes downstream of the assimilation genes, the glutamine synthetase (GS)-glutamate synthase (GOGAT) cycle and the ornithine-urea cycle (OUC) (Fig. 5B and C and Fig. S8). The GS-GOGAT cycle catalyzes ammonium into organic nitrogen. The ornithine-urea cycle (OUC) has been shown to aid the model diatom Phaeodactylum tricornutum in returning from N-limiting conditions by serving as a hub for nitrogen redistribution and pathway to satisfy arginine demand (9, 30). E. huxleyi also possesses an OUC which is hypothesized to serve a similar role (31).

In C. decipiens, significant differences in expression between time points were nearly absent in both the GS-GOGAT cycle under the iron-replete conditions and OUC genes under both iron treatments (Fig. 5B). These results correspond to previous field observations of constitutive and relatively high expression in diatoms pre- and postupwelling such that the cell already possesses the cellular machinery to channel assimilated nitrogen toward growth (7). Like that for the nitrate assimilation genes, E. huxleyi exhibits a different transcriptional response for the GS-GOGAT and OUC genes in response to the changing UCBC conditions. OUC genes were generally downregulated under stationary growth and upregulated in response to returning to light at T4 (P < 0.05). In response to returning to light, E. huxleyi strongly upregulated GS-GOGAT genes (P < 0.001) (Fig. 5C). Both C. decipiens and E. huxleyi displayed low and constitutive expression of arginase and urease (3rd to 31st percentiles), the final steps of the OUC pathway which similarly were not shown to exhibit a frontloading response in the same degree in the field (see Fig. S9) (7). The OUC thus likely serves to satisfy demand for arginine in both organisms, as predicted through metabolic flux models in P. tricornutum (9). Much like that for the nitrate assimilation genes, C. decipiens appears to frontload the GS-GOGAT and OUC genes, whereas E. huxleyi exhibits a more responsive transcriptional pattern.

N limitation has previously been studied extensively in several diatoms and E. huxleyi, providing points of comparison. In the model pennate diatom P. tricornutum, genes encoding nitrate transporters, nitrate reductases, and nitrite reductases were upregulated under N limitation (32, 33), although a study examining these genes in a short-term response associated them with nitrogen-replete conditions (9). Upregulation of these genes under N starvation in most cases was further shown in Fragilariopsis cylindrus, Pseudo-nitzschia multiseries, and several Skeletonema species (33, 34); therefore, most studies are consistent with our findings here of upregulation under N limitation, suggesting that this response is shared among diatoms.

For GS-GOGAT and OUC-related genes, a pattern in gene expression is less apparent. Most studies reported changes in expression of these genes under N-replete or -deplete conditions in contrast to the lack of change observed here in C. decipiens, although the directionality of the expression among studies is not consistent (9, 32, 33). However, in F. cylindrus, P. tricornutum, and P. multiseries, stable expression of OUC genes was shown between N conditions apart from that for the first step, carbamoyl-phosphate synthase, similar to observations in this study (33). Therefore, unlike nitrate assimilation genes, transcriptional regulation of these genes appears to be variable among diatoms as previously suggested (33). The coordinated frontloading response of GS-GOGAT and OUC genes with the nitrate assimilation genes may be a strategy more exclusively employed by bloom-forming diatoms residing in upwelling environments.

Aligning with the results presented here, E. huxleyi has previously been shown to downregulate or show reduced activity of NRT, NR, and NiR while upregulating AMT under N limitation, although these results are also not entirely consistent across studies (29, 31, 35–37). E. huxleyi also was shown to downregulate the initial OUC genes under N limitation, the one exception here being the gene for arginosuccinate synthase at T2 (31). GS-GOGAT gene expression was similarly lower in previous studies under N limitation albeit not significantly different here (31).

Iron bioavailability is an important control on phytoplankton growth in upwelling areas characterized by narrow continental shelves (1). Biogeochemical modeling of the California Current ecosystem indicates that iron limitation is highly prevalent, although only occurring at small scales and short durations relative to N limitation (38). This iron limitation occurs in both the surface mixed layer and at depth as shown with subsurface chlorophyll maximum layers (13, 39). Furthermore, projected increases in upwelled nitrate unmatched by iron coupled to ocean acidification is projected to expand iron limitation in these regions (1, 17, 18). Here, UCBC simulations were performed under iron-limiting conditions to examine the effects of iron status on responses throughout the cycle.

Within the first exponential phase, the growth rate in C. decipiens declined from 0.86 ± 0.05 day⁻¹ in iron-replete cells to 0.36 ± 0.08 day⁻¹ in iron-limited cells (Fig. S1). The reduction in growth rates was much smaller for E. huxleyi (0.68 ± 0.05 day⁻¹ to 0.63 ± 0.08 day⁻¹). Reductions in F_v:F_m were also observed in both organisms, indicating iron stress in the low-iron treatments. With F_v:F_m already impacted by iron stress, the variation as a result of light and/or nitrogen stress was comparatively less than when iron was replete. Previous experiments have also suggested that E. huxleyi is well suited to grow at low-Fe concentrations (40). Higher growth rates may suggest that E. huxleyi has a competitive advantage under Fe limitation, yet during the UCBC, iron-limited E. huxleyi cells were unable to reinstate exponential growth upon returning to light in two of three replicates, indicating that E. huxleyi can be severely impacted by the combination of both prolonged darkness and Fe limitation (Fig. S4). Importantly, in the one E. huxleyi replicate that was able to grow, it did not return to exponential growth for 48 h, while C. decipiens appeared to immediately grow, similar to the observations of the iron-replete cells. This observation is in line with diatoms being able to outcompete coccolithophores under dynamic conditions (41).

Iron limitation can influence nitrate assimilation, as the required enzymes, nitrate and nitrite reductase, can account for greater than 15% of the cellular iron requirement (42, 43). Iron limitation has also been shown to influence the expression of the genes that encode these proteins (44, 45). In C. decipiens, C:N ratios followed a similar pattern when comparing the iron-replete to iron-stressed conditions (Fig. 2). In early stationary phase (T2 and T6), cellular C:N ratios were significantly reduced relative to those of iron-replete cells at the same time points but reached similar values after 10 days in the dark. Importantly, C:N ratios in iron-limited cells still quickly approached the Redfield ratio upon returning to the light (T4), indicating that the shift-up response in N assimilation may not be significantly impacted by low-iron conditions in C. decipiens.

Iron limitation altered the expression levels of certain N-related genes at specific time points in C. decipiens (see Fig. S10). The nitrate assimilation pathway was downregulated at T3 and T6. This decreased expression at T3 made expression levels between T3 and T4 similar, which may explain similar responses observed pre- and post-upwelling in the field-based incubations, where iron stress in the initial community may have occurred (7).

Genes encoding ammonium and urea transporters as well as urease were upregulated under low-iron conditions, as expected with nitrate assimilation being impacted. Under low-iron conditions, variability in expression between time points increased for mitochondrial GOGAT genes compared to the fairly stable expression in the iron-replete treatment, with significant upregulation under low-iron conditions at T3 (Fig. 5B and S10). The mitochondrial GOGAT genes are specialized for the recovery of recycled nitrogen and/or assimilation of urea (9); therefore, in conjunction with the ammonium and urea transporters being upregulated under low-iron conditions, this transcriptional response is consistent with a greater emphasis on utilizing recycled forms of N. As with the lack of change in OUC genes throughout the UCBC, there were no significant differences in OUC transcript abundances between iron treatments at all time points.

Other differentially expressed genes show how C. decipiens responds to iron limitation during the UCBC (Fig. 6A; see Data set S4 at https://doi.org/10.6084/m9.figshare.c.5315135.v1). Similar to other diatoms, several transcriptional indicators of iron stress were significantly upregulated in the iron-limited treatment. Transcripts for the inorganic iron uptake gene phytotransferrin (pTF; previously known as ISIP2a) were upregulated at all time points in iron-limited cells (P < 0.0001) (17). Flavodoxin (FLDA) is an iron-free replacement for the photosynthetic electron transfer protein ferredoxin (46), and diatoms can possess multiple copies of FLDA, with some that are not iron responsive (47). Two copies of FLDA were detected, one of which was significantly upregulated at all time points in iron-limited cells (P < 0.0001) and the other was not. A gene encoding a metacaspase similar to one identified in Thalassiosira pseudonana (TpMC4) was significantly upregulated at all time points (P < 0.0001); this protein has been implicated in programmed cell death, suggesting a triggering of this process in iron-limited cells (48).

Some genes displayed differential expression only at certain time points throughout the UCBC. ISIP3, a gene of unknown function normally highly expressed under iron limitation in diatoms, was only upregulated at T1 (P = 0.06) and T5 (P = 0.005), indicating that this gene may not be a reliable marker for iron stress in Chaetoceros under serial limitation or colimitation scenarios, unlike that for Thalassiosira oceanica (49–51). A ferric reductase in the cytochrome b₅ family was differentially expressed only in the light (P < 0.001). This contig is homologous to FRE3 and FRE4 in P. tricornutum which are speculated to be involved in heme-based iron uptake (44). Two class I fructose 1,6-bisphosphate aldolase (FBA) contigs were upregulated under low-iron conditions, although one was only upregulated in the light and other in the dark (T3) (P < 0.05). These enzymes are important for cellular C metabolism and are commonly found to be upregulated in iron-limited diatoms (45, 52).

Three contigs encoding aureochromes were differentially expressed as a function of iron status at certain time points. Aureochromes are light-induced transcription factors specific to stramenopiles, and all three of these in C. decipiens contain the characteristic bZIP domain for DNA binding and the light-oxygen-voltage (LOV) domain for blue-light sensing (53, 54). These proteins have also previously been implicated in diel regulation of gene expression (55, 56). None were differentially expressed between iron treatments during exponential growth. Two (12603_2 and 14552) increased in expression under iron limitation when in stationary phase in the dark, suggesting that these genes have an increased role in regulating gene expression when both nitrogen and iron are limiting, a colimitation scenario which may be commonplace in the ocean (57). The third aureochrome (12603_1) was upregulated in the high-iron treatment under stationary growth, suggesting it has a more pronounced role in gene regulation when iron is replete but nitrogen is limiting.

Despite its prevalence and ability to grow under low-iron conditions, differential expression in E. huxleyi between varied iron treatments had not been analyzed previous to this study. This analysis was restricted to T1 and T2, as only these time points had replicates available (Fig. 6B; see Data Set S5 at https://doi.org/10.6084/m9.figshare.c.5315135.v1). More genes were significantly differentially expressed at T1 than at T2 (Fig. 6B). Only two contigs showed significant opposing regulation between T1 and T2, that is, one was upregulated at one time point but downregulated at the other.

Many of the genes upregulated under iron limitation were similar to those that have been frequently observed in diatoms. These include the previously mentioned genes ISIP3, FBA class I, and FLDA. High expression of FBP1 and ZIP family divalent metal transporters were also detected. FBP1 is involved in siderophore-based Fe uptake in diatoms, suggesting a similar functional role in E. huxleyi (58). ZIP family transporters may also permit transport of ferrous iron (59). Other highly expressed genes have no annotated function (e.g., E. huxleyi CCMP1516 identifiers [IDs] 103255, 110131, and 112797) suggesting that E. huxleyi possesses or employs unique mechanisms for coping with iron limitation. These genes and the proteins they encode are clear targets for future study and possible markers of iron stress for E. huxleyi in the natural environment.

Phytoplankton in upwelling regions are exposed to highly dynamic conditions as a result of upwelling cycles. In these upwelling cycle simulations, both a representative diatom, C. decipiens, and the coccolithophore, E. huxleyi, displayed physiological and transcriptional plasticity to these cycles. Part of this physiological response included transitioning to high cellular C:N ratios, which has previously been observed in natural populations and influences the biogeochemical cycling of these elements. Although both phytoplankton species were transcriptionally responsive, their gene expression patterns at each time point were highly divergent. In particular, C. decipiens exhibited a pattern of frontloading transcriptional and nitrate assimilation genes, which may explain certain diatoms’ ability to rapidly respond and outcompete other groups once upwelling events occur. As observed in the diatom Cylindrotheca fusiformis, nitrate reductase transcripts accumulate under nitrogen starvation to enable extremely rapid translation into protein upon addition of nitrate (60); therefore, the diatom frontloading response observed here may be restricted to high transcription and rely on posttranscriptional regulation, as the synthesis of mature proteins is expensive.

As a highly diverse group, diatoms likely occupy a wide range of niches and growth strategies (20). The response characterized here may be relatively exclusive to bloom-forming diatoms in upwelling regions, although transcriptional responses in nitrate assimilation genes appear fairly consistent, as indicated from other diatoms in laboratory studies. Furthermore, it is likely even among bloom-forming diatoms that there are further divergences in the transcriptional responses, as observed in a previous field simulation (7). Additional study of other phytoplankton species, particularly other bloom-forming diatoms such as Pseudo-nitzschia and Thalassiosira, will be necessary to more fully understand the spectrum of responses to UCBC conditions.

Both C. decipiens and E. huxleyi employed mechanisms in response to iron stress that are common to responses previously observed in diatoms; however, certain genes such as those encoding aureochromes appear to only be differentially expressed when nitrogen and/or light is simultaneously limiting, indicating a unique response under colimitation scenarios. As Si and Fe limitation are often linked in diatoms (25, 61), examination of Si and Fe colimitation responses throughout the UCBC would provide further environmentally relevant insights.

In E. huxleyi, many of the genes that were differentially expressed under iron limitation are not functionally characterized. Importantly, E. huxleyi was unable to respond to prolonged darkness under iron-limiting conditions, whereas C. decipiens’s shift-up response once returned to the light was unaffected by iron status. Some projections suggest increased iron limitation in upwelling environments as a result of ocean acidification and increases in upwelled nitrate relative to iron (1, 17); however, these findings suggest that diatoms will continue to exhibit the shift-up response and outcompete other phytoplankton groups under both high- and low-iron upwelling conditions.

Clearly, there are a number of differences between the simulation experiments presented here and the conditions found in the natural environment (Fig. 1). Temperature and pressure would not be constant as they were in the simulations. Upwelled waters are also more acidic, and pH adjustments were not made (62). Although nutrient concentrations would increase as the cells sink to depth, nitrate, the depleted nutrient, was not supplemented until the cells returned to light following the dark phase. In addition, iron in nature is bound to a wide variety of organic ligands and supply is dynamic, whereas our simulations utilized steady-state buffered conditions with iron bound to a single chelator, EDTA (63). The laboratory cultures also contained significantly higher cell densities than those found in the natural environment, and the phytoplankton strains studied individually here, albeit nonaxenically, prevented these responses from being examined with direct interactions and competition among a larger number of species. Such complexities will require further field-based observations and experimentation to more fully understand phytoplankton responses to the UCBC.

Two phytoplankton species isolated from the California Upwelling zone were used in this study: a diatom, Chaetoceros decipiens (UNC1416), and a coccolithophore, Emiliania huxleyi (UNC1419). Species isolation and identification were described by Lampe et al. (7). Cultures were grown in artificial Aquil* medium following trace metal clean (TMC) techniques at 12°C and continuous light at approximately 115 μmol photons m⁻² s⁻¹ (64–66). Starting macronutrient concentrations were modified from standard Aquil* amounts such that nitrate rather than silicate would be drawn down and depleted first by the phytoplankton isolates (50 μmol liter⁻¹ NO₃, 10 μmol liter⁻¹ PO₄, 200 μmol liter⁻¹ H₄SiO₄). The culture medium was treated with Chelex in a TMC room, microwave sterilized, allowed to cool, and then supplemented with filter-sterilized (0.2-μm Acrodisc) EDTA-trace metals (minus iron) and vitamins (B₁₂, thiamine, and biotin). Trace metals were buffered using 100 μmol liter⁻¹ EDTA and added at the concentrations listed by Sunda et al. (66) except Cd, which was added at a total concentration of 121 nmol liter⁻¹. Premixed Fe-EDTA (1:1) was added separately at a total concentration of 1,370 nmol liter⁻¹ for the high-iron treatments or 3.1 nmol liter⁻¹ for the low-iron treatments. The resulting concentrations of iron not complexed to EDTA (Fe′) were estimated as 2,730 pmol liter⁻¹ (high iron) and 6 pmol liter⁻¹ (low iron) (66), resulting in iron-replete and iron-limited growth, respectively. Vitamins were added at concentrations for f/2 medium. Media were allowed to equilibrate overnight before use.

Cultures were acclimated to media at the two iron concentrations using the semicontinuous batch culture technique in acid-cleaned 28-ml polycarbonate centrifuge tubes until growth rates were consistent between transfers. Cultures were then grown in a 1-liter polycarbonate bottle, and when in exponential phase, 5 ml (high iron) or 20 ml (low iron) of each culture was then transferred to triplicate clean 2-liter polycarbonate bottles fitted with Teflon tubing for subsampling and bubbling. Seawater was continuously stirred and bubbled with air passed through 1.2 mol liter⁻¹ HCl and then through TMC Milli-Q water before entering the culture bottles to remove any potential contaminants in the air. Cultures were not axenic, although sterile techniques were used for all culture work to minimize contamination.

The upwelling conveyer belt cycle (UCBC) was simulated by growing or transitioning the cultures to various light and nutrient conditions (Fig. 1B; see also Fig. S1 in the supplemental material). First, to simulate bloom and termination phases, cultures were grown in exponential phase until nutrient exhaustion, at which time, cells entered stationary phase. After 3 days in stationary phase, bubbling and stirring were halted, and the bottles were moved to a dark incubator (0 μmol photons m⁻² s⁻¹, 12°C) to simulate sinking out of the euphotic zone (2, 3, 7). Although temperature, pressure, and nutrient concentrations would change in deeper waters in the natural environment, temperature and pressure were kept constant and cells were not supplied with additional nutrients during this stage. After 10 days, 500 ml of culture was transferred into 1.5 liters of fresh medium and returned to the light with stirring and bubbling resumed. Cultures were subsequently grown until stationary growth was measured for 2 days, thus completing the cycle. Subsampling was always performed at the same time of day (approximately 09:00 Eastern time) except T4, which occurred 12 h after T3 (21:00).

The 10-day time frame was selected based on our previous field-based simulated upwelling experiments, where satellite data indicated that upwelling had not occurred in the region for approximately 13 days, although intervals are certainly variable in nature (7). Additionally, preliminary iron-replete UCBC experiments suggested that E. huxleyi is incapable of reinstating growth after 20 days of darkness and thus unable to complete the UCBC simulation with that duration of darkness (Fig. S5). This provided an upper limit for the length of the dark period, further leading to the selection of 10 days.

Subsampling was performed by first dispensing into sterilized bottles under a laminar flow hood in a TMC room and then immediately aliquoting from this subsample. Relative abundances and growth throughout the experiment were assessed by regular measurements of blank-corrected raw fluorescence units (RFUs) with a Turner 10-AU fluorometer. Specific growth rates were calculated during exponential growth phases from the natural log of RFUs versus time (67).

Experimental trials were conducted with both phytoplankton isolates to determine the effect of different lengths of darkness on growth during UCBC simulations (Fig. S5). The two isolates were grown in acid-cleaned 28-ml polycarbonate centrifuge tubes with synthetic Aquil* medium that contained the same macronutrient and iron concentrations as the iron-replete treatment. The dark phase was varied with 5-day intervals between 5 and 20 days in duration. After the dark phase was complete, 7 ml of culture was transferred into 21 ml of fresh medium. Changes in phytoplankton biomass and calculation of growth rates were assessed with near-daily measurements of blank-corrected raw fluorescence units (RFUs) using a Turner 10-AU fluorometer (67).

Five-milliliter samples were preserved in 2% Lugol’s iodine in glass vials. Cells were then enumerated on an Olympus CKX41 inverted microscope using a 1-ml Sedgwick-Rafter counting chamber after allowing the cells to settle for 5 min; a minimum of 300 cells were counted within 10 to 30 fields of view.

Fifty milliliters of culture was filtered through a 0.45-μm mixed cellulose ester filter (25 mm; EMD Millipore Corporation, Burlington, MA, USA) under gentle vacuum pressure (<100 mm Hg) and immediately frozen at −80°C until analysis. Chlorophyll a extraction was performed using 90% acetone at −20°C for 24 h and measured via in vitro fluorometry on a Turner Designs 10-AU fluorometer using the acidification method (68).

Particulate carbon (PC) and nitrogen (PN) were obtained by gentle vacuum filtration of 50 ml of culture onto a precombusted (450°C for 6 h) GF/F filter. Filters were immediately stored in petri dishes at −20°C. Prior to analysis, filters were dried at 65°C for 24 h and then encapsulated in tin. Total nitrogen and carbon were quantified with a Costech 4010 CHNOS elemental combustion system according to U.S. Environmental Protection Agency method 440.0 (69). Three blanks were run alongside the samples and were all below the detection limits (0.005 mg N and 0.071 mg C). Samples were not acidified and included particulate inorganic carbon, including that from the coccolith in E. huxleyi. Carbon and nitrogen on a per cell basis were calculated by dividing the particulate carbon or nitrogen concentration by cell concentration.

The maximum photochemical yield of Photosystem II (F_v:F_m) was measured using a Satlantic fluorescence induction and relaxation system (FIRe) (70, 71). Samples were acclimated to low light for 20 min prior to measuring the minimum (F_o) and maximum (F_m) fluorescence yields. Data were blank corrected using microwave-sterilized Aquil medium. The resulting F_v:F_m was derived from the induction profile using a saturating pulse (single turnover flash; 20,000 μmol photons m² s⁻¹) for a duration of 100 μs. The gain was optimized for each sample (400, 600, or 800), and the average from 50 iterations was obtained.

Filtrate from the 0.45-μm filters used for RNA was transferred to polypropylene tubes and immediately frozen at −80°C. Dissolved nitrate plus nitrite (NO₃ plus NO₂) concentrations were quantified with an OI Analytical Flow Solutions IV autoanalyzer according to EPA method 353.4 (72). The detection limit for NO₃ plus NO₂ was 0.2 μmol liter⁻¹.

Approximately 300 ml of culture was filtered onto 0.45-μm Pall Supor polyethersulfone filters (47 mm) using gentle vacuum pressure and immediately frozen and stored at −80°C. For the Chaetoceros experiments, total RNA was extracted using the RNAqueous-4PCR Total RNA isolation kit (Ambion, Foster City, CA, USA) according to the manufacturer’s protocol with an initial bead beating step to disrupt cells. For the E. huxleyi experiments, total RNA was extracted using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer’s protocol except for an initial bead beating step and two chloroform steps instead of one to separate proteins and DNA. Trace DNA contamination was removed by DNase 1 (Ambion) digestion at 37°C for 45 min, and RNA was purified with the RNeasy MinElute Cleanup kit (Qiagen, Germantown, MD, USA). Quality and quantity were checked on 1% agarose gels and with a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA).

Phytoplankton cultures were grown to late exponential phase and extracted as described under “RNA extraction.” RNA libraries were created with either the Illumina TruSeq stranded mRNA library preparation kit for C. decipiens or the KAPA stranded mRNA-Seq kit for Illumina platforms for E. huxleyi. Samples were sequenced on an Illumina MiSeq (300 bp, paired-end reads) and an Illumina HiSeq 2500, with one lane in high-output mode (100 bp, paired-end reads) and another lane in rapid-run mode (150 bp, paired-end reads) (see Table S1 at https://doi.org/10.6084/m9.figshare.c.5315135.v1).

Raw reads were trimmed for quality with Trimmomatic v0.36 (73) and then assembled de novo with Trinity v2.5.1 with the default parameters for paired reads and a minimum contig length of 90 bp (74). Contigs were clustered based on 99% similarity using CD-HIT-EST v4.7 (75), and then protein sequences were predicted with GeneMark S-T (76). Protein sequences were annotated by best homology (lowest E value) with the KEGG (release 86.0) (77), UniProt (release 2018_03) (78), and PhyloDB (v1.076) (79) databases via BLASTP v2.7.1 (E value ≤ 10⁻⁵) and with Pfam 31.0 (80) via HMMER v3.1b2 (see Table S2 at the URL mentioned above). KEGG ortholog (KO) annotations were assigned from the top hit with a KO annotation from the top 10 hits (https://github.com/ctberthiaume/keggannot).

RNA libraries for the UCBC simulations were constructed using a custom protocol for 3′ poly(A)-directed mRNA sequencing (mRNA-seq; also known as TagSeq) based on the report by Meyer et al. (81) and adapted for Illumina HiSeq based on the studies by Lohman et al. (82) and Strader et al. (83). For most samples, 1 μg but as low as 250 ng of total RNA was fragmented by incubating at 95°C for 15 min.

First-strand cDNA was synthesized with SMARTScribe reverse transcriptase (TaKaRa Bio, Mountain View, CA, USA), an oligo(dT) primer, and template switching to attach known sequences to each end of the poly(A) mRNA fragments. cDNA was then amplified with a PCR consisting of 32 μl sterile H₂O, 5 μl deoxynucleoside triphosphates (dNTPs; 2.5 mM each), 5 μl 10× Titanium Taq buffer (TaKaRa Bio), 1 μl of each primer (10 μM) designed to amplify the sequences attached during cDNA synthesis, and 1 μl of Titanium Taq DNA polymerase (TaKaRa Bio). The PCR was run with an initial denaturing step at 95°C for 5 min and then 17 cycles consisting of 95°C for 1 min, 63°C for 2 min, and 72°C for 2 min. cDNA amplification was then verified on a 1% agarose gel and purified with the QIAquick PCR purification kit (Qiagen). cDNA was quantified with the Quant-iT double-stranded DNA (dsDNA) high-sensitivity assay (Invitrogen), and cDNA concentrations were then normalized to the same volume.

Samples were then barcoded with a PCR consisting of 11 μl sterile H₂O, 3 μl dNTPs (2.5 mM each), 3 μl 10× Titanium Taq buffer (TaKaRa Bio), 0.6 μl TruSeq universal primer (10 μM), 6 μl barcoding primer (1 μM), Titanium Taq polymerase (TaKaRa Bio), and 6 μl of purified cDNA with an initial denaturing step at 95°C for 5 min and then 5 cycles consisting of 95°C for 40 s, 63°C for 2 min, and 72°C for 1 min. Samples were barcoded from both ends using unique combinations of 6-bp sequences on both the TruSeq universal primer and barcoding primers. Products were then again confirmed on a 2% agarose gel and combined into small pools of 6 to 8 samples. The 400- to 500-bp region of these pools was extracted with the QIAquick gel extraction kit (Qiagen), quantified with the Quant-iT dsDNA high-sensitivity assay, and then mixed in equal proportions. The library was sequenced at the University of Texas at Austin Genomic Sequencing and Analysis Facility on an Illumina HiSeq 2500 (three lanes, single-end 50-bp reads) with a 15% PhiX spike-in to target approximately 8 million reads per sample (see Table S2 at https://doi.org/10.6084/m9.figshare.c.5315135.v1) (81).

PCR duplicates were removed based on a matching degenerate leader sequence incorporated during cDNA synthesis and a match of the first 20 bases as described by Dixon et al. (84). Reads were trimmed for adapter removal and quality with a minimum length of 20 bases with Trimmomatic v0.38 (73). Transcripts were quantified with Salmon v0.9.1 with the GC bias correction option against predicted proteins from de novo transcriptome assemblies generated via paired-end Illumina sequencing of the same strains (85). Normalized counts, fold change values, and associated Benjamini-Hochberg adjusted P values were calculated with DESeq2 v1.22.2 (86). Normalized counts are based on the median of the ratios of the observed counts (87). Reference transcriptome and experimental results are further detailed in Tables S1 and S2 and Data Set S6 at https://doi.org/10.6084/m9.figshare.c.5315135.v1.

Correlation network analysis of gene expression was performed using WGCNA using variance-stabilizing transformed counts generated by DESeq2 as input (88). All available samples from both the iron-replete and iron-limited treatments were used for this analysis (see Data Set S6 at the URL mentioned above). Contigs with fewer than 10 normalized counts in more than 90% of the samples were removed before analysis as suggested by the WGCNA manual. A signed adjacency matrix was constructed for each organism using a soft-thresholding power of 14 for the C. decipiens samples and 22 for the E. huxleyi samples (see Fig. S11 at the URL listed above). The soft-thresholding power was selected as the lowest power that satisfies the approximate scale-free topology criteria as described by Langfelder and Horvath (88). The adjacency matrix was then transformed into a topological overlap matrix (TOM). Subnetworks, or modules, were created by hierarchical clustering of the TOM with the following parameters: minModuleSize, 30; deepSplit, 2. The gray module is reserved for genes that are not assigned to a module. The first principal component of each module, i.e., eigengene, was correlated (Pearson) with each time point, and associated P values were calculated. Contigs in modules that were significant for a specific time point were then also correlated (Pearson) to the same time points to evaluate significance on an individual contig basis (see Data Sets S2 and S3 at the URL listed above). The associated P values were corrected for multiple testing using the Benjamini-Hochberg false-discovery rate controlling procedure (89); significant correlations were those with an adjusted P value of less than 0.05. Parameters for all software used are listed in Table S3 at the URL mentioned above.

One- and two-way analyses of variance (ANOVAs) followed by Tukey’s multiple-comparison tests were performed on the biological and chemical properties of the seawater media in GraphPad PRISM v8.4.3.

Supplementary data sets and tables were deposited in Figshare (https://doi.org/10.6084/m9.figshare.c.5315135.v1). Raw reads for the C. decipiens and E. huxleyi experiments were deposited in the National Center for Biotechnology Information (NCBI) Sequence Read Archive (SRA) under the accession numbers (no.) SRP176464 and SRP178167 (BioProject accession no. PRJNA513340 and PRJNA514092), respectively. Raw reads for the reference transcriptomes were also deposited in SRA under the accession no. SRP234548 and SRP234650 (BioProject accession no. PRJNA593314 and PRJNA593538). Assemblies and annotations for each reference transcriptome were deposited in Zenodo (90, 91). Example code for the DESeq2 and WGCNA analyses are available at https://github.com/rhlampe/simulated_upwelling_exps.