Phytoplankton transcriptomic and physiological responses to fixed nitrogen in the California current system

Marine phytoplankton are responsible for approximately half of photosynthesis on Earth. However, their ability to drive ocean productivity depends on critical nutrients, especially bioavailable nitrogen (N) which is scarce over vast areas of the ocean. Phytoplankton differ in their preferences for N substrates as well as uptake efficiencies and minimal N requirements relative to other critical nutrients, including iron (Fe) and phosphorus. In this study, we used the MicroTOOLs high-resolution environmental microarray to examine transcriptomic responses of phytoplankton communities in the California Current System (CCS) transition zone to added urea, ammonium, nitrate, and also Fe in the late summer when N depletion is common. Transcript level changes of photosynthetic, carbon fixation, and nutrient stress genes indicated relief of N limitation in many strains of Prochlorococcus, Synechococcus, and eukaryotic phytoplankton. The transcriptomic responses helped explain shifts in physiological and growth responses observed later. All three phytoplankton groups had increased transcript levels of photosynthesis and/or carbon fixation genes in response to all N substrates. However, only Prochlorococcus had decreased transcript levels of N stress genes and grew substantially, specifically after urea and ammonium additions, suggesting that Prochlorococcus outcompeted other community members in these treatments. Diatom transcript levels of carbon fixation genes increased in response to Fe but not to Fe with N which might have favored phytoplankton that were co-limited by N and Fe. Moreover, transcription patterns of closely related strains indicated variability in N utilization, including nitrate utilization by some high-light adapted Prochlorococcus. Finally, up-regulation of urea transporter genes by both Prochlorococcus and Synechococcus in response to filtered deep water suggested a regulatory mechanism other than classic control via the global N regulator NtcA. This study indicated that co-existing phytoplankton strains experience distinct nutrient stresses in the transition zone of the CCS, an understudied region where oligotrophic and coastal communities naturally mix.


Introduction
Marine phytoplankton are responsible for about half of photosynthesis on Earth [1]. The growth and productivity of phytoplankton are constrained by the availability of critical nutrients, primarily nitrogen (N), phosphorus (P), and iron (Fe) [2][3][4][5][6][7]. Over wide areas of the ocean, N limits phytoplankton growth [7,8] because most phytoplankton cannot use the gaseous form, dinitrogen (N 2 ). Instead they use a variety of other chemical forms of N, including organic forms such as urea, as well as inorganic forms such as nitrite (NO 2 -) and nitrate (NO 3 -), with preferences differing among phytoplankton [9][10][11][12]. Physical transport processes, remineralization, and other factors, often seasonal and region-specific, affect which N forms are available [9, [13][14][15][16], which in turn affects species composition and trophic dynamics [17,18]. Much of our knowledge about phytoplankton communities comes from either oligotrophic ocean gyres or coastal regions, the latter of which typically have higher levels of available N. Less is known about boundary current systems where the two communities naturally mix, such as the transition zone of the California Current System (CCS). The present study examines CCS phytoplankton transcriptomic responses to added N substrates and links those responses to observed physiological responses. The CCS is a highly dynamic environment where N availability varies with location and season and, together with other factors, controls phytoplankton community composition and productivity [13,16,19,20]. Along the coast, the upwelling of cold nutrient rich waters by Ekman transport leads to high NO 3 concentrations [13,21]. About 200 km offshore, the CCS system is bounded by the warm oligotrophic California Current (CC) which flows toward the equator [13]. The transition zone (TZ) between the coastal mesotrophic region and the oligotrophic CC is characterized by high-energy eddies and cross-stream jets that drive mesoscale variability in nutrients and phytoplankton productivity [22][23][24]. In the late spring and summer, upwelled nutrient-rich water travels offshore across the TZ and, along the way, becomes depleted of NO 3 due to biological activities and physical forces [13,14]. In the late summer and early fall, weaker nearshore upwelling followed by mixing of TZ and oligotrophic CC waters can make NO 3 scarce and limit phytoplankton primary production [13]. Moreover, iron (Fe) can be depleted faster than NO 3 leading to Fe limitation or N and Fe co-limitation for phytoplankton in the TZ [3,14].
In the coastal regions of the CCS, diatoms and other large photosynthetic eukaryotes drive primary production, whereas in the TZ, diverse and abundant photosynthetic picoeukaryotes and picocyanobacteria are major contributors [25][26][27]. Picocyanobacteria of the genera Synechococcus and Prochlorococcus are major contributors to primary production due to their abundance throughout the world ocean, including the CCS [28][29][30]. Both genera consist of several phylogenetic clades and ecological types (ecotypes) that occupy different niches based on temperature, nutrient, and light availability [7,[30][31][32][33][34][35][36][37][38]. Multiple clades of Synechococcus and Prochlorococcus have been observed across the CCS from coastal to transition zones including novel clades or those that lack close reference genomes [39,40]. Generally, Prochlorococcus genomes are smaller and more GC-rich, and thus require less N compared to Synechococcus [34, [41][42][43][44]. Under N limitation, some Prochlorococcus conserve N by using alternative transcription start sites to produce shorter proteins [45], whereas some Synechococcus make N available by degrading photosynthetic pigments [46]. The two genera vary in their abilities to utilize different N species. For example, while nearly all Synechococcus strains can assimilate NO 3 -and NO 2 -, only some Prochlorococcus strains can [34,47]. Moreover, intrastrain variations in environmental sensing and nutrient assimilation capabilities have been observed [47][48][49] which makes it challenging to infer functional diversity from 16S rRNA gene sequences alone.
Deep sequencing approaches have become widely used in marine meta-omics and are effective for studying abundant community members, genes, and transcripts [50][51][52][53][54]. Alternatively, high-resolution microarrays have advantages for detecting rare community members and for differentiating among closely related strains [55][56][57]. The MicroTOOLs microarray targets functional genes in abundant and rare members of oligotrophic and coastal surface marine microbial communities, including picocyanobacteria (Prochlorococcus and Synechococcus), N 2 -fixing cyanobacteria, photosynthetic eukaryotic phytoplankton, as well as a small number of viruses [55]. Probes on the array can distinguish among closely related strains known from culture and environmental samples [55], and thus can elucidate the physiological state of individual microbial populations within a complex community [56].
In August 2014, we conducted experiments to determine the effects of different N substrates on phytoplankton communities at two stations, one in the North Pacific Ocean and one in the transition zone of the California Current (Stn. TZ), originally described in Shilova et al.
[10]. Surface water samples were incubated with NO 3 -, ammonium (NH 4 + ), urea, or filtered deep water (FDW) for 48 hours (T48). Two treatments had added Fe 3+ , either alone ("Fe") or with a mix of N substrates ("N+Fe"), to determine the effects of Fe on the utilization of N substrates. After 48 hours, all treatments resulted in changes in phytoplankton cell abundances and photosynthetic activity at both locations, with differences between phytoplankton groups. Prochlorococcus had large increases in biomass in response to NH 4 + and urea, while both eukaryotic phytoplankton and Synechococcus had highest biomass increases in response to FDW, and to N+Fe for Synechococcus. Moreover, distinct physiological responses were observed within sub-populations of Prochlorococcus and Synechococcus. In order to better understand the variable responses to N substrates among phytoplankton groups and sub-populations in the CCS transition zone, the present work used the MicroTOOLs microarray to examine transcriptomic changes that occurred 24 hours (T24) after the substrates were added. We hypothesized that transcript level changes at T24 would indicate which phytoplankton taxa were N-limited, and thus help explain changes in cell abundances for individual phytoplankton groups observed at T48. Furthermore, we hypothesized that the diversity in physiological responses within Prochlorococcus and Synechococcus would be evident in the transcriptomic responses measured at sub-population resolution.

RNA extraction and processing
Total RNA was extracted from the samples using the Ambion 1 RiboPure RNA purification kit (ThermoFisher Scientific) with the addition of a bead-beating step during TRI Reagent extraction as described in Shilova et al. [55]. DNA was removed from the total RNA extracts in a solution using RNase-Free DNase Kit (Qiagen, Germantown, MD, USA), and RNA was purified again with RNA Clean & Concentrator™-25 (Zymo Research, Irvine, CA, USA) according to the manufacturers' protocols. The RNA quality and quantity were evaluated using the Agilent BioAnalyzer RNA Nano Kit (Agilent, Santa Clara, CA, USA) and Qiagen Qubit. T0 and T24 RNA samples with an RNA Integrity Number greater than 9 were processed for microarray analyses. Three hundred ng of total RNA was used for synthesis of double-stranded cDNA using the TransPlex Whole Transcriptome Amplification 2 Kit (Sigma-Aldrich, St. Louis, MO, USA) following the manufacturer's instructions. For amplification control, 0.5 μl of 1:100 dilution of the mRNA Ambion ERCC mix 1 (ThermoFisher Scientific) was spiked into each RNA sample prior to cDNA synthesis. Amplified transcriptome samples were purified with the GenElute PCR CleanUp Kit (Sigma Aldrich). cDNA was labeled at the Roy J. Carver Center for Genomics (The University of Iowa, USA) using the Agilent SureTag DNA Labeling Kit (Cat# 5190-3400) following a protocol described in "Agilent Oligonucleotide ArrayBased CGH for Genomic DNA Analysis: Enzymatic Labeling for Blood, Cells, or Tissues (Version 7.3 March 2014)" (S1 File). Cy3-labeled cDNA was hybridized using a Gene Expression Hybridization Kit (Cat# 5188-5242) and following the protocol described in "One-Color Microarray-Based Gene Expression Analysis: Low Input Quick Amp Labeling (Version 6.7, September 2014)" (S1 File). Microarray platform GPL24371 Agilent-073391 MicroTOOLs_171K_oligo_v2.0 was used in this study. Microarrays were scanned at the Roy J. Carver Center for Genomics using an Agilent SureScan Microarray Scanner G2600D (Serial #: SG13134301) and the Agilent scanning protocol GE1_1200_Jun14 (Feature Extractor software version 11.5.1.1). The microarray data were submitted to The National Center for Biotechnology Information (NCBI) Gene Expression Omnibus (GEO) under accession GSE130464.

Microarray analyses
Microarray analyses were done using the MicroTOOLs R package (https://www.jzehrlab.com/ microtools) using the same approaches and parameters described in Robidart et al. [56] (robust multi-array averaging of probes and quantile normalization across samples), except that differentially expressed (DE) genes were identified based on a fold change significantly greater than 1.2 (Benjamini-Hochberg adjusted p-value <0.05; [58]). A total of 23 RNA samples were analyzed in this study, including three replicates for controls taken at T0 and T24 and three replicates for each treatment collected at T24, with the exception of the NH 4 + treatment which had two replicates (the third NH 4 + replicate did not pass microarray hybridization quality control and was excluded from the analysis). Multiple approaches were applied to examine responses at the phylogroup and individual strain levels and also due to variability among some treatment replicates (Table 1). We used an Ensemble of Gene Set Enrichment Analyses (EGSEA; [62]) to identify sets of genes that collectively were significantly differentially expressed based on a consensus of twelve GSEA algorithms (Benjamini-Hochberg adjusted p-value < 0.01, unadjusted p-value calculated using Wilkinson's method [64] to combine p-values from the GSEA algorithms). We defined each gene set to contain MicroTOOLs gene targets from a specific phylogroup and physiological response, e.g. genes from high-light (HL) adapted Prochlorococcus that typically increase during nitrogen limitation ( Table 2). Each gene set was analyzed for differential expression in each treatment relative to controls at T24.
The same volume of seawater was processed for metatranscriptomic analysis of each T24 sample, and neither the cell abundances for major phylogroups nor the microbial community composition, as measured with the 16S rRNA gene V4 sequencing, changed significantly by T24 [10]. Thus, the transcriptomic changes we report are not artifacts of community composition differences among treatments. However, the same amount of cDNA was hybridized to each microarray, therefore transcript level increases from one phylogroup could result in transcript level decreases for other phylogroups. We found that this artifact was too small to account for differentially expressed gene sets [described in S1 File and S1 Fig]. For each gene set g from phylogroup p analyzed for differential expression, the average fold change of member genes (a g,p,t ) in treatment t vs. the control at T24 was compared to the average fold change in total transcripts from the phylogroup across replicates of the treatment (a tot,p,t shown in S1 Fig for each phylogroup and treatment). Only differentially expressed gene sets with a g,p,t / a tot, p,t > 1 if an increase or a g,p,t / a tot,p,t < 1 if a decrease were further analyzed. Finally, we used weighted correlation network analysis (WGCNA; [63]) to identify genes with highly correlated expression patterns across all 20 of the T24 samples, mainly for genes that had strong patterns across taxa (e.g. urtA). Strains were considered to be present in a sample based on the detection of at least five of their target genes in that sample, or at least half of their targets if they had fewer than five targets (mainly photosynthetic eukaryotes).

Transcripts were detected from oligotrophic and coastal microbial taxa
The surface microbial community at Stn. TZ 24 h after N addition was diverse and transcriptionally active. A total of 9760 genes had detectable transcripts in one or more of the 23 total Table 1. Transcriptomic response analyses.

Analysis
Identifies How applied Ref.
Sample clusters compared within this study and to a NPSG study [56]. [59]

Single-gene DE
Individual genes that are significantly differentially expressed between two conditions. FDW, NO 3 -, and urea treatments vs. controls at T24. Controls at T24 vs. at T0. [60,61] EGSEA Sets of genes that collectively are significantly differentially expressed between two conditions. All treatments vs. controls at T24. Controls at T24 vs. at T0. [62] WGCNA Modules of genes with highly correlated expression patterns across samples. Gene expression profiles across T24 samples were correlated. [63] Transcription patterns were analyzed using three approaches. The single-gene DE analysis is the traditional approach for identifying differentially expressed genes. In the main text "DE" is used only for results from the single-gene analysis. NMDS-non-metric multidimensional scaling; DE-differentially expressed; EGSEA-Ensemble of Gene Set Enrichment Analyses; WGCNA-weighted correlation network analysis.
https://doi.org/10.1371/journal.pone.0231771.t001 samples (Materials and methods). Based on detected genes, samples had on average 575±24 distinct strains of the 924 included on the MicroTOOLs array. All taxa identified in control samples at T0 were also detected at T24 and represented major phylogroups found in the open ocean and coastal/transitional environments: picocyanobacteria (Prochlorococcus and Synechococcus), alpha-, gamma-and beta-proteobacteria, and eukaryotic phytoplankton including stramenopiles (e.g. diatoms) and haptophytes ( Table 3). Some of the highest transcript relative abundances, hereafter referred to as "transcript levels," were detected from the dominant picocyanobacteria that Shilova et al.
[10] identified The EGSEA analysis identified collectively significant changes from the genes in each set. When EGSEA was applied to a phylogroup, all of the genes in each set had MicroTOOLs targets from multiple strains. When EGSEA was applied to a specific strain, all available genes for the strain were included. Some strains lacked targets for some genes, but usually each strain had multiple targets for every gene listed. Note that each stress gene set has member genes that all increase or decrease together when the stress changes. https://doi.org/10.1371/journal.pone.0231771.t002 using 16S rRNA gene V4 sequencing and oligotyping analysis, including high-light adapted (HL) Prochlorococcus strains MED4 and MIT9515 (Table 4, S1 Table). However, other strains detected in the present study were rare or not observed by Shilova et al.
[10], and had even higher transcript levels than some of the "dominant" strains from the same ecotype. For example, HL Prochlorococcus strain MIT9301 comprised just 0.33% of the total Prochlorococcus 16S rRNA gene sequences, but in 18 of the 23 samples it had a higher total transcript intensity (normalized to detected genes) than MIT9515. Similar to Prochlorococcus, we detected transcripts from two dominant Synechococcus strains previously identified at Stn. TZ using 16S rRNA gene sequencing [10], CC9605 (clade II) and CC9902 (clade IV), as well as other clades; however, transcripts from CC9902 were rare (Table 4, S1 Table). Given the rarity of Synechococcus at Stn. TZ (relative abundancẽ 0.8% based on 16S rRNA gene abundances, or 3.9±0.7 × 10 3 cells mL -1 ; [10]), the diversity of detected Synechococcus strains demonstrated the sensitivity of MicroTOOLs for studying rare but transcriptionally active community members.
Altogether, we observed numerous strains of PE, Prochlorococcus (43 types), and Synechococcus (57 types), based on transcripts detected from multiple gene targets for the picocyanobacteria (Materials and methods). The numbers of Prochlorococcus and Synechococcus strains were higher than reported by Shilova et al.
[10] using 16S rRNA gene oligotyping (11 Prochlorococcus and 31 Synechococcus). The complex ensemble of detected open ocean and coastal strains of phytoplankton may reflect mesoscale processes in the CCS, such as the circulation of baroclinic jets and eddies [23], as well as late summer physical forces that mix open ocean and transition zone waters [69]. However, the stability of and dynamics within such mixed communities are unknown. Physiologically, they respond differently to nutrients than open ocean

Nutrient additions resulted in distinct transcription patterns
The transcription pattern for each sample, its metatranscriptome, was defined by the transcript relative abundances for the 9760 total detected genes. Metatranscriptomic changes that occurred within 24 hours of nutrient additions (Fig 1 and S2 Fig)   The MicroTOOLs design includes far more probes that target Prochlorococcus and Synechococcus than PE and heterotrophic bacteria. Thus, the similarity of NO 3 and FDW metatranscriptomes to T0 metatranscriptomes, together with the significant increases in Chl a concentrations and primary productivity rates, suggest that other bacteria and eukaryotes had advantages over Prochlorococcus for responding to higher concentrations of NO 3 and other nutrients in FDW [70,71]. Furthermore, the responses to FDW were similar to those in another NPSG study that used MicroTOOLs (see NMDS in [56]; in situ samples in that study are analogous to controls at T0 here). Both studies have captured metatranscriptomic changes from surface microbial communities responding to influxes of nutrients, such as those expected as a result of anticyclonic eddies. Fe was likely not a limiting nutrient for the microbial community. Metatranscriptomes from treatments with Fe or N+Fe were more similar to controls at T24 and also more variable compared to treatments with only N substrates (except for NH 4 + ; Fig 1). This suggests that the community, as measured with the microarray, responded more to N than to Fe additions. Consequently, cell abundances were not significantly higher in N+Fe treatments compared to NO 3 alone [10].
In a separate NMDS analysis, the metatranscriptomes from this CCS study clustered more tightly than metatranscriptomes from a NPSG study that investigated the effects of deep water mixing on the surface microbial community ( [56]; S3 Fig, mean within-study Euclidean distances, before NMDS, were 45.5 for CCS versus 139.6 for NPSG). This suggests that the CCS surface community metatranscriptomes were less perturbed over a 24 h period by the addition of nutrients than were the NPSG metatranscriptomes by the addition of deep water (S1 File). The small relative abundance changes and lack of changes in photosynthetic efficiency by T24 in the CCS samples relative to NPSG samples from the same study [10] support this conclusion.

Prochlorococcus had similar responses to NH 4 + and urea
Prochlorococcus showed signs of alleviation of N stress 24 h after the NH 4 + or urea addition as evidenced by transcript decreases for N stress associated genes and increases for photosynthesis (PS) and carbon fixation genes (rbcL ; Figs 2A, 3A, 4A, 4B, 5A and 5B). The Ensemble of Gene Set Enrichment Analyses (EGSEA) showed that N stress gene transcripts decreased for HL Prochlorococcus strains overall (Fig 2A), as well as for the dominant strains MED4 and MIT9515 (S1 Table). Although N stress genes across LL strains did not change significantly in the NH 4 + or urea treatments, their carbon fixation genes were up-regulated in both treatments.
PS genes were up-regulated in the NH 4 + treatment across the low-light adapted (LL) Prochlorococcus strains (Fig 2A) and for multiple LL strains in the urea treatment (NATL2A, MIT9211, CCMP1375; S1 Table). When cyanobacteria experience low NH 4 + conditions, the global nitrogen transcriptional regulator NtcA responds to high intracellular C-N ratios by activating transcription of genes associated with N acquisition and/or N stress, including ntcA itself [72][73][74][75]. Although we did not observe significantly differentially expressed (DE) ntcA from Prochlorococcus in the NH 4 + or urea treatments, transcripts from most N transport and metabolism genes decreased (non-DE) in the urea treatments (Fig 4A and 4B). These genes included the ammonium transporter  amt, carbamoyl-phosphate synthase carA, and glutamine synthetase glnA. Moreover, EGSEA indicated that collectively all of these N transport and metabolism genes (Materials and methods Table 2) had significant decreases for HL Prochlorococcus in both treatments, as described above (Fig 2A). Analysis of relative transcript levels for rbcL versus ntcA indicated that treatments with NH 4 + or urea substrates provided enough N to shift the internal C-N balance in Prochlorococcus MED4 within 24 hours (S1 File, S4 Fig). The up-regulation of PS and rbcL genes and down-regulation of N stress genes in response to urea or NH 4 + are consistent with the highest increases of Prochlorococcus cell abundances observed in these treatments at T48 [10] (Fig 2A) and further support the ability of natural populations of Prochlorococcus to assimilate either substrate [12,43,76].   Table), suggesting that it did not utilize NO 3 which is consistent with nongrowth of MED4 on NO 3 in culture [78]. However, 8 of the 11 detected MED4-like ntcA genes had decreased transcript levels suggesting the presence of MED4-like subpopulations at Stn. TZ that utilized NO 3 -or reduced N available from other cells utilizing NO 3 - (Fig 4A; S1 Table). Subpopulations with different abilities to utilize NO 3 were also suspected for the second most abundant Prochlorococcus population, MIT9515-like, as well as for the rare Prochlorococcus MIT9301-like population, based on mixed PS responses within each strain (Fig 5A; S1 Table). The FDW treatment had a similar NO 3 concentration as the treatment with NO 3 alone and also resulted in diverse responses from Prochlorococcus PS and N-stress genes (Figs 4A, 4B, 5A and 5B; S1 Table). Altogether, the results showed that specific Prochlorococcus subpopulations in the CCS were utilizing NO 3 and that utilization varied within closely related populations. A mixed population might explain why Prochlorococcus cell abundance increases by T48 were less (but significant) in the treatments that added NO 3 than in the treatments that added urea or NH 4 + [10,43] (Fig 2A). It is also possible that the Prochlorococcus subpopulations that responded positively to NO 3 were either minor members of the Prochlorococcus community and/or growing more slowly on NO 3 than on reduced nitrogen sources [78]. Mixed populations of Prochlorococcus strains with regards to their NO 3 utilization capabilities have been observed at Pacific and Atlantic sampling sites [79] and specifically within MED4-like  subpopulations in the CCS [80]. The high variability in NO 3 uptake rates was also reported recently among cells within phytoplankton groups including Prochlorococcus at both CCS and NPSG locations [12] suggesting that intra-population heterogeneity in NO 3 utilization is likely widespread.

Synechococcus was N-limited but had dissimilar response patterns to Prochlorococcus
Transcriptomic changes indicated that Synechococcus likely assimilated the added N substrates. For each of the treatments with N substrates alone (urea, NO 3 -, or NH 4 + ) or with FDW, transcript levels increased from rbcL and PS genes from Synechococcus strains collectively (Fig 2B), as well as from individual strains (Figs 3B and 5C; S1 Table). For example, both the abundant Synechococcus CC9605 and CC9902 and the rare CC9311 had increased transcript levels of rbcL (all DE) and PS genes (all DE except for CC9605) in the NO 3 treatment. These results were in line with the observed Synechococcus abundance increases by 48h in all treatments except for the NH 4 + treatment (Fig 2B; [10]), and consistent with N assimilation by Synechococcus in culture studies [81][82][83].
In contrast to the strong up-regulation observed for rbcL and PS genes across Synechococcus, a mixed response was observed for N stress genes though mainly from non-dominant strains. In the NO 3 treatment, transcript levels decreased for one ntcA gene from strain

PLOS ONE
BL107, and for the nitrate/nitrite permease (nrtP) gene from strains WH7805 and RCC307 (Fig 4C). In the urea treatment, amt transcript levels decreased for RS9917. However, the EGSEA analysis did not show a collective decrease for Synechococcus N stress genes in response to N substrates alone, unlike for HL Prochlorococcus (Fig 2A and 2B). Despite the low abundance of Synechococcus, its ntcA transcript levels were on average 3× higher than those of (abundant) Prochlorococcus (0.45±0.02 vs. 0.16±0.01% of total transcripts from all taxa, respectively). The high ntcA transcript levels, low variability in rbcL to ntcA transcript levels (S1 File, S4 Fig), and lack of changes for N stress genes from Synechococcus in N treatments might be due to its different strategy for adapting to low reduced N environments in comparison to Prochlorococcus. Some Synechococcus strains, such as open ocean strain WH8103, maintain the ability to transport and utilize oxidized N forms even if NH 4 + is present [83,84].
Moreover, Synechococcus WH7803 expresses ntcA at a low level in the presence of NH 4 + [85].
Maintaining the ability to utilize oxidized N forms might be an energetically favorable strategy for Synechococcus in a high light and low reduced N environment such as the open ocean. The lack of an overall response from Synechococcus N stress genes might also indicate that some strains were N and Fe co-limited as supported by physiological results [10]. For example, cell abundance increases for a strain of KORDI-100 were nearly twice as much in response to the N+Fe treatment (significant) compared to the Fe addition alone (not significant) [10]. In contrast, treatments with N+Fe and Fe alone resulted in similar cell abundance increases for dominant strains within CC9605 and CC9902 [10]. Transcriptional changes from Fe stress genes differed by strain (S5 Fig; S1 Table), suggesting strain-specific Fe requirements in Synechococcus and in line with our hypothesis of N and Fe co-limitation. For example, in response to FDW or NO 3 -, transcript levels from ferric transcriptional regulator fur genes decreased for CC9902 (DE for one fur target), suggesting reduced Fe stress, while fur transcript levels mainly increased for other strains (non-DE; S5 Fig; S1 Table). The differences in Fe requirements [34, 86,87] may dictate which Synechococcus strains dominate CCS waters under Fe limitation [6,88]. Consistent with strain-specific Fe requirements, the collective responses from Synechococcus to Fe (alone or with N) was indistinct. In the Fe treatment, genes associated with Fe stress as a whole did not decrease, nor did rbcL or PS genes increase as was observed in treatments that added N alone, or FDW which produced the largest Synechococcus cell increases (Fig 2B). Thus, the physiological and transcriptional changes together suggest that Fe benefited some strains but that N was the primary limiting nutrient for the majority of Synechococcus strains.
In summary, the transcriptomic results suggest that most Synechococcus spp. at Stn. TZ were N-limited. The addition of N substrates led to increases in PS and rbcL transcripts, and the results support findings that natural populations of Synechococcus utilize NH 4 + , NO 3 -, and urea [12]. However, there was a lack of response from Synechococcus N stress genes, in comparison to Prochlorococcus, which may have resulted from differences in Fe requirements or N metabolic strategies. The more abundant community members, such as Prochlorococcus and high nucleic acid (HNA) heterotrophs, likely outcompeted Synechococcus for N substrates, especially in urea and NH 4 + , partly due to diffusion advantages afforded by their larger surface to volume ratios [89][90][91]. These factors likely prevented substantial growth of Synechococcus in comparison to Prochlorococcus and HNA in N treatments by T48 [10].

Upregulation of urea uptake genes by picocyanobacteria in FDW
In the FDW treatment, the up-regulation of picocyanobacterial urea transporter genes urtA, but not other N genes under NtcA regulation, suggested additional regulatory controls for urtA besides classic NtcA regulation. Specifically, transcripts for nearly all Prochlorococcus urtA genes increased while transcripts for other genes associated with N stress and under NtcA control, such as ntcA and amt, decreased (Figs 2A, 4A and 4B; S1 Table). Synechococcus urtA also was up-regulated in the FDW treatment (Fig 4C), and WGCNA analysis indicated that urtA transcription patterns were similar across picocyanobacteria but distinct from nearly all other detected N stress genes (S1 File). Similar to other N stress genes that are under NtcA control, known urtA promoters in Prochlorococcus and Synechococcus have a putative binding site for NtcA [92], and urtA is up-regulated under N limitation in Prochlorococcus HL and LL strains [93]. However, it seems unlikely that the increases of picocyanobacterial urtA transcripts observed in the FDW treatment were due to classic NtcA control. Diel expression of cyanobacteria urtA [94] also seems an unlikely explanation because urtA decreases would have been expected in our pre-dawn samples (S1 File). Possibly, the urtA transcript level increases were specific to the availability of urea, which would benefit picocyanobacteria as both a N and C source [95]. A recent NanoSIMS study showed that both Prochlorococcus and Synechococcus assimilated urea during short (~5 h) incubations [12]. Urea might have been generated by other community members in the FDW treatment. Heterotrophic bacteria have been shown to produce substantial amounts of urea in the Southern California Bight [96], and in Shilova et al.
[10], heterotrophic bacteria cell density (HNA cells mL -1 ) more than doubled (by T48) in response to FDW. Eukaryotic phytoplankton, which responded rapidly to FDW in our study, might also have produced urea [97][98][99]. The delay between FDW addition and urea production might explain why picocyanobacteria urtA was up-regulated at T24 in FDW but down-regulated at T24 in the urea treatment. If urtA is a fast responding gene, 24 hours might have been too late to detect its up-regulation in the urea addition treatment.
Alternatively, chemical or biological interactions after mixing of surface waters with deep waters could have triggered the up-regulation of Prochlorococcus and Synechococcus urtA. Robidart et al. [56] also observed urtA increases from Prochlorococcus 24 h after the addition of subeuphotic zone (130 m) FDW to the surface community in the NPSG. However, other studies did not observe up-regulation of urtA in response to the addition of high-molecular weight dissolved organic matter [52] or the addition of deep seawater from 700 m to surface NPSG microbial communities [100]. Thus, the urtA response observed in the present study and [56] might be specific to conditions or picocyanobacterial populations. Future culture-based and in situ experiments will help determine whether urea transporter responses reflect the regulation of one or multiple copies of urtA, by NtcA and possibly other transcription factors [34].

Photosynthetic eukaryote RuBisCO transcript levels increased in N treatments
All nutrient additions elicited rbcL transcriptional responses from photosynthetic eukaryotes (PE) and PE cell abundance increases by T48. The PE rbcL response to NO 3 was robust, with significant (DE) transcript level increases for 283 of the 296 detected PE rbcL targets, including all detected haptophyte targets and most detected stramenopile targets (mainly for diatoms on MicroTOOLs; Fig 3C). Transcript levels also increased for PE nitrate reductase (Nr) targets in the NO 3 treatment (3 DE for diatoms from Amphora, Phaeodactylum, and an undetermined genus; S1 Table), consistent with the activation of Nr transcription in the presence of nitrate in PE [101,102]. Diatoms use nitrate opportunistically [70,71] rather than as a preferred N form [12], and haptophytes can utilize nitrate as well [103,104]

PLOS ONE
Phytoplankton responses to fixed nitrogen in the CCS and stramenopiles (88 rbcL, and 1 Nr from Phaeodactylum tricornutum; Fig 3C; S1 Table). The urea treatment resulted in no DE genes from PE (however, 294 rbcL genes had non-DE transcript increases). Like the single-gene DE analysis described above, the EGSEA analysis found significant rbcL transcript level increases for PE in the NO 3 and FDW treatments, but it also identified large (>1.2×), significant increases in the urea, NH 4 + , and Fe treatments (Fig 2C). The rbcL increases were followed by significant PE cell abundance increases at T48 in all treatments with the highest increase in FDW (Fig 2C) [10].
The concentration of Fe at Stn. TZ was below the level of detection (0.058 nmol L -1 ; [10]). The Fe treatment added Fe 3+ to a final concentration of 2 nmol L -1 [10], whereas [6] estimated that Fe concentrations >0.5 nmol L -1 are enough for even the largest coastal CCS diatoms to grow if other nutrients are not limiting. The Fe treatment resulted in rbcL transcript increases from stramenopiles and haptophytes and modest PE abundance increases by T48 [10] (Fig 2C). Thus, for PE the Fe addition was sufficient for growth and the low initial concentrations of NO 3 and NH 4 + (above) were not limiting. It was therefore surprising that the N+Fe treatment, which also added Fe 3+ to 2 nmol L -1 , resulted in rbcL transcript decreases from stramenopiles ( Fig 2C). One potential explanation is that the added N substrates were used by N and Fe co-limited microorganisms which then outcompeted those detected by MicroTOOLs. Competition might also explain why haptophyte rbcL transcript increases were smaller in the N+Fe treatment ( Fig 2C). Moreover, the differences in rbcL responses between stramenopiles and haptophytes to N+Fe indicate that haptophytes have lower Fe requirements than stramenopiles at Stn. TZ [105]; Fig 2C).

Conclusions
This study and the previous by Shilova et al.
[10] together provide a high-resolution snapshot of phytoplankton biodiversity in the CCS transition zone along three dimensions: taxonomic, transcriptional, and functional. The differential transcriptomic and physiological responses to N forms revealed in these two studies indicated that N was the primary limiting nutrient, but the responses differed with N substrate among Prochlorococcus, Synechococcus, and PE. Transcriptomes of some Synechococcus and PE populations indicated possible N and Fe co-limitation, in line with highest cell abundance increases after the addition of N+Fe for Synechococcus or FDW for both phylogroups [10]. Diverse transcriptomic responses were observed among closely related strains and sub-populations within Prochlorococcus and Synechococcus, indicative of different N assimilation capabilities and/or degrees of N limitation. For example, we observed heterogeneous populations of Prochlorococcus in their capacity to utilize NO 3 -, supporting previous single-cell nutrient uptake rate findings [12]. In the treatment with FDW, the unexpected increase of picocyanobacterial urtA while other N stress genes decreased highlights our incomplete understanding of urea utilization by marine microorganisms. The differences among natural phytoplankton populations in transcriptional and physiological responses were likely due to many factors including genetics, competition, and prior environmental conditions. Gene transcripts were detected from a mixed community of open ocean, transitional, and coastal strains reflecting the dynamic but poorly understood physical-biological interactions in the CCS transition zone. In fact, heterogeneity in responses within this mixed CCS microbial community might be a reason for the less pronounced whole transcriptome responses to added N observed in this study in comparison to the responses of the oligotrophic gyre community observed in a previous MicroTOOLs study [56]. Future studies along the dimensions of biodiversity at multiple locations and seasons will provide a more complete picture of how available N forms impact phytoplankton communities in this dynamic and productive part of the Pacific Ocean.
Supporting information S1 File. Additional details on methods and results. (DOCX) S1 Table. Detected strains and genes. Spreadsheet containing normalized log 2 transcript levels for all 9760 detected genes and whether they were DE, as well as their taxonomic and functional annotation. (CSV) S2 Table. Distances between the sample metatranscriptomes. Spreadsheet containing the Euclidean distances between all 23 samples. Each sample is represented by the normalized log 2 transcript levels for all 9760 detected genes. The NMDS analysis (Fig 1) was performed on these distances. Treatments in the spreadsheet are ordered by increasing distance to T0 controls, calculated as the median of pairwise distances between treatment and T0 control replicates. (CSV) S1 Fig. Transcript proportions from major phylogroups. For each major phylogroup and treatment, the mean proportion of transcripts are shown. Transcripts are for all detected genes within the phylogroup and averaged over replicates within the treatment. For each phylogroup the proportions did not vary much across treatments. Consequently, the differential expression we report for a phylogroup mainly reflects changes in how its transcripts were distributed across its gene targets on MicroTOOLs, i.e. up-and down-regulation of its genes. Note that the large proportion of transcripts from Synechococcus, which was rare by relative abundance, is because Synechococcus has many probes on MicroTOOLs. PS = photosynthetic. (PDF) S2 Fig. Heat map of differentially expressed genes. A total of 3805 significantly differentially expressed (DE) genes were identified by comparing controls at T0 vs. T24, or treatments with NO3 -, FDW, or urea at T24 vs. controls at T24. Genes (rows) cluster mainly by phylogroup, evident by the opposite patterns for Prochlorococcus and Synechococcus genes, mainly for photosynthesis and carbon metabolism. However, also apparent are sub-clusters by Process or Element (e.g., several blocks of N metabolism genes) and by Gene (e.g., rbcL across eukaryotic phytoplankton, mostly denoted as "Other"). Sample clusters (columns) were robust (solid discs indicate �70% support out of 1000 bootstraps). With the exception of the Fe and NH4 + treatments, at least two of the three replicates for each treatment cluster tightly as in Fig 1. White cells indicate that the strain for that gene was not detected in the sample (Materials and methods). (PDF)

S3 Fig. NMDS of CCS and NPSG samples.
NMDS of metatranscriptomes from the present CCS study and from an NPSG study [56]. In comparison to the CCS community, the NPSG community shows larger shifts by 24 h in response to added FDW (130 m) that had high concentrations of nitrite. The CCS community appeared to be N-limited (see main text) while the NSPG community was N-starved. Except for unfiltered deep water (UDW), all metatranscriptomes are from surface samples. The in situ samples were collected at dawn and dusk from 14-16 September 2011 with an Environmental Sample Processor. In order to compare metatranscriptomes from the CCS and NPSG in a single NMDS analysis, microarray data for samples from both experiments were processed together (including normalization; Materials and methods). The stress was 0.09. HL Prochlorococcus Ni-containing superoxide dismutase (NiSOD) encoding genes. The heat map shows the mean transcript levels (normalized as described in Materials and methods) of all 19 detected NiSOD-encoding genes ("NiSOD targets", rows) from HL Prochlorococcus strains in each treatment. Transcript levels are not centered or scaled. Therefore, deeply red NiSOD targets (e.g. #5590 for MED4) represent subpopulations that were more actively transcribing or were more abundant. For each target gene, transcript levels are ranked across treatments from highest (cells with 1) to lowest (cells with 8). Treatments are sorted by the mean of the ranks in each column: NiSOD targets usually had their highest transcript levels in the T0 controls, second highest levels in the urea treatment, and lowest levels in FDW. Annotation at the left of the heat map indicates for each target gene whether it was DE, the HL strain, and the WGCNA module to which the target was assigned.
For targets that were DE in multiple conditions, only one is shown with preference given to treatments over controls. All DE conditions are in S1 Table. (PDF) (PDF) experiment set-up and sampling. The authors also appreciate the efforts of the anonymous reviewers who critically read and helped to clarify the manuscript.