Response of Spring Diatoms to CO2 Availability in the Western North Pacific as Determined by Next-Generation Sequencing

Next-generation sequencing (NGS) technologies have enabled us to determine phytoplankton community compositions at high resolution. However, few studies have adopted this approach to assess the responses of natural phytoplankton communities to environmental change. Here, we report the impact of different CO2 levels on spring diatoms in the Oyashio region of the western North Pacific as estimated by NGS of the diatom-specific rbcL gene (DNA), which encodes the large subunit of RubisCO. We also examined the abundance and composition of rbcL transcripts (cDNA) in diatoms to assess their physiological responses to changing CO2 levels. A short-term (3-day) incubation experiment was carried out on-deck using surface Oyashio waters under different pCO2 levels (180, 350, 750, and 1000 μatm) in May 2011. During the incubation, the transcript abundance of the diatom-specific rbcL gene decreased with an increase in seawater pCO2 levels. These results suggest that CO2 fixation capacity of diatoms decreased rapidly under elevated CO2 levels. In the high CO2 treatments (750 and 1000 μatm), diversity of diatom-specific rbcL gene and its transcripts decreased relative to the control treatment (350 μatm), as well as contributions of Chaetocerataceae, Thalassiosiraceae, and Fragilariaceae to the total population, but the contributions of Bacillariaceae increased. In the low CO2 treatment, contributions of Bacillariaceae also increased together with other eukaryotes. These suggest that changes in CO2 levels can alter the community composition of spring diatoms in the Oyashio region. Overall, the NGS technology provided us a deeper understanding of the response of diatoms to changes in CO2 levels in terms of their community composition, diversity, and photosynthetic physiology.


Introduction
Progressive increases in the seawater partial pressure of CO 2 (pCO 2 ) and decreases in pH (i.e., ocean acidification, [1]) caused by industrial CO 2 emissions could affect biological processes in the ocean [2]. Because CO 2 is the primary substrate for photosynthesis, ocean acidification can enhance CO 2 availability for marine phytoplankton. Accordingly, ocean acidification could play key roles in controlling productivity and organic matter production [3][4][5], thereby potentially affecting the biogeochemical and ecological processes in the ocean.
Studies on the impacts of increased CO 2 levels have been conducted on a variety of phytoplankton species [3,6,7,8]. In the context of the global carbon cycle and feedbacks to climate change, diatoms are an important class of phytoplankton because they are responsible for approximately 40% of total primary production in the global ocean [9]. Comparative studies suggested that the growth of diatoms tended to increase with elevated CO 2 levels [10,11]. Indeed, previous field studies have provided some evidence that elevated CO 2 levels could enhance the growth of diatoms in the Equatorial Pacific [12], the North Atlantic [13], and the Southern Ocean [14]. However, no significant or negative effects of elevated CO 2 on diatoms have been reported from the bloom-inducing field incubation experiments in the Raunefjord [15], the Western Subarctic Gyre (WSG) of the North Pacific [16], and the Bering Sea basin [17,18]. These discrepancies could be caused by the species-specific differences in CO 2 response [19], although experimental conditions and other environmental factors might also affect the outcomes. Hence, detailed taxonomic information on diatoms would be indispensable for understanding their responses to changes in CO 2 levels in seawater.
In addition, marine diatoms can experience a decrease in CO 2 availability in some oceanic regions. In the Oyashio region of the western North Pacific, the pCO 2 in surface seawater decrease significantly from winter (ca. 400 μatm) to spring (< 200 μatm) as a consequence of intense diatom blooms [20]. Previous studies have suggested that diatoms can overcome low CO 2 availability by using biophysical CO 2 -concentrating mechanisms (CCMs), which provide a high CO 2 concentration at the site of carboxylation [21,22]. However, some experiments demonstrated that the growth and productivity of diatoms decreased under low CO 2 (100-220 μatm) conditions [7,23]. At present, it is largely unknown how diatom assemblages respond to decrease in CO 2 availability, which occurs during spring blooms in the Oyashio region.
Recently, molecular biological techniques have become important tools for understanding the community composition and biodiversity of natural microbial assemblages including phytoplankton based on DNA sequence information [24][25][26]. It is reported that the abundance of particular genes can be an indicator of phytoplankton pigment biomass [27]. These techniques can also be applied to assess the responses of phytoplankton assemblages under different CO 2 levels [28][29][30]. For example, Hopkinson et al. [28] used 18S rRNA and nitrate reductase (NR) gene fragments to estimate the community compositions of haptophytes and diatoms, respectively, in different CO 2 levels. In addition, the rbcL gene, which encodes the large subunit of the enzyme ribulose bisphosphate carboxylase/oxygenase (RubisCO), was also used to estimate the abundance and community composition of phytoplankton assemblages during the CO 2 -controlled mesocosm experiment conducted off the Norwegian coast [29]. These authors demonstrated a rapid shift in the community composition of pico-sized prasinophytes (i.e. Bathycoccus and Micromonas), which were difficult to identify by light microscopy, in response to increased CO 2 levels.
In addition, mRNA (cDNA) transcripts of functional genes could be a proxy for physiological responses of phytoplankton assemblages to environmental change. For example, the transcript levels of the rbcL gene could be an indicator of productivity because CO 2 is assimilated via RubisCO in the Calvin-Benson-Bassham (CBB) cycle [31,32]. Changes in gene expression and protein synthesis levels of RubisCO occur within a few hours in response to environmental change such as light and nutrient availability [33][34][35]. Therefore, it can be used to infer the potential effect of different CO 2 levels on the natural phytoplankton assemblages even from short-term incubation experiment. Endo et al. [27] showed that diatom-specific rbcL transcripts decreased in response to elevated CO 2 levels after 2 or 3 days of incubation in the Bering Sea. They also showed that a shift occurred in the community composition of photosynthetically active diatoms with an increase in seawater CO 2 levels. However, the conventional molecular cloning method with Sanger sequencing [36], which was used in previous studies, is sometimes inadequate for the extraction of comprehensive information from environmental samples, due to limited throughput, and may therefore underestimate the taxonomic richness of microbial assemblages [37,38].
Recent advances in next-generation sequencing (NGS) technologies can overcome this limitation because these deep sequencing technologies can now generate several hundred thousand reads per sample [39]. Practically, metagenomic and amplicon sequencing using NGS platforms have been used to reveal the community composition and/or diversity of bacteria [40], phytoplankton [41], and zooplankton [42] in marine environments. However, to the best of our knowledge, these new technologies have not been used to estimate the effects of CO 2 availability on marine phytoplankton assemblages.
Here, we report the use of NGS technologies in combination with real-time PCR (qPCR) and HPLC pigment analysis to improve our understanding of the effects of CO 2 availability on the community structure and photosynthetic physiology of diatoms in the Oyashio region. Because of the massive diatom blooms in spring, this region has one of the greatest capacities for seasonal biological drawdown of pCO 2 in surface waters among the world's oceans [43]. Field observations revealed that the blooms were mainly composed of centric diatoms, such as Chaetoceros spp. and Thalassiosira spp. [44][45][46][47]. Taniguchi [48] noted that these diatoms contribute to the efficient energy transfer to higher trophic levels in the Oyashio region from spring to summer. Consequently, the spring diatom blooms significantly contribute to the formation of the highly productive fishing grounds in Oyashio and its surrounding waters [48,49]. However, despite the crucial roles of diatom blooms in the ecosystems and biogeochemical processes, there are no reports on the impacts of different CO 2 availability on phytoplankton assemblages in the Oyashio region. We therefore conducted an on-deck incubation experiment using spring Oyashio waters under different pCO 2 levels.

Materials and Methods
Seawater samples used in this study were not collected in a protected area. No specific permissions were required for sampling of seawater in the locality. The samples taken for our study did not involve endangered or protected species.

Experimental setup and sampling
The study was carried out aboard the R/V Tansei Maru (JAMSTEC) during the KT-11-7 cruise in May 2011. Water samples were collected from 10-m depth at a station (41°30' N, 144°00' E, 1668 m depth, Fig 1) in the Oyashio region of the western North Pacific on 9 May with Niskin-X bottles attached to a CTD-CMS system. A total of 200 L of seawater was poured into four 50 L polypropylene carboys through silicon tubing with 197 μm mesh Teflon nets to remove large particles. Prior to incubation, FeCl 3 solution (5 nmol L −1 in final concentration) was added to the carboys to promote the development of phytoplankton blooms because the growth of bloom-forming diatoms can be limited by low iron (Fe) availability in this region [45,46]. resulting value was used as the control benchmark. To control pCO 2 in the incubation bottles, air mixtures containing 180, 350, 750, and 1000 μatm CO 2 were bubbled into the incubation bottles. The flow rate of the gases was set at 100 mL min -1 for the first 24 h, and thereafter maintained at 50 mL min -1 . Incubation was conducted on deck in a tank with running surface seawater for 3 days, maintaining in situ temperature of 5°C and 50% surface irradiance, which was adjusted using a neutral density screen. The temperature of the incubation tanks ranged between 4.7 and 6.0°C throughout the incubation. All of the samples, except on day 0, were collected between 5:00 and 6:00 a.m.

Carbonate chemistry, nutrients, and Chl a analyses
Samples for total alkalinity (TA), dissolved inorganic carbon (DIC), nutrients, and size-fractionated chlorophyll a (Chl a) were collected on days 0, 0.6, 1.6, and 2.6 (hereafter days 0, 1, 2, and 3, respectively) from the incubation bottles. Samples for TA and DIC were collected in airtight glass vials and poisoned with HgCl 2 prior to their storage at 4°C for analysis on shore. Concentrations of TA and DIC were measured with a total alkalinity analyzer using the potentiometric Gran Plot method (ATT-05, Kimoto Electric) following Edmond [50]. The values of pCO 2 and pH were calculated from the TA and DIC using the CO2SYS program [51]. Nutrient samples were poured into 10 mL plastic tubes and stored at −20°C until they were analyzed on shore. Concentrations of nitrate plus nitrite, nitrite, phosphate, and silicic acid were measured using a QuAATro-2 continuous-flow analyzer (Gran+Luebbe). Size-fractionated Chl a samples were collected onto a 10 μm pore size polycarbonate filter without vacuum or onto Whatman GF/F filters under a gentle vacuum (< 0.013 MPa). Pigments were extracted in N,N-dimethylformamide (DMF; [52]), and Chl a concentrations were then determined using a Turner Design fluorometer (model 10-AU) with the non-acidification method of Welschmeyer [53].

FIRe fluorometry
To obtain the maximum photochemical quantum efficiency (F v /F m ) of photosystem II (PSII) for phytoplankton in each bottle, a Fluorescence Induction and Relaxation (FIRe) fluorometer (Satlantic) was used. The water samples were collected daily from the incubation bottles and were measured following Sugie et al. [18].

HPLC and CHEMTAX analyses
Pigment samples for high-performance liquid chromatography (HPLC) were collected on days 0, 2, and 3. The water samples (600-1000 mL) were filtered onto Whatman GF/F filters under a gentle vacuum (< 0.013 MPa), flash-frozen in liquid nitrogen and then stored in a deep freezer (−80°C) until analysis. HPLC pigment analysis was conducted following the method of Endo et al. [16]. The net growth rate (μ, day −1 ) of each pigment was calculated from the following equation: where C is the pigment concentration (μg L −1 ) and t is incubation time (i.e., 2.6 days). Based on the chemotaxonomic pigment concentrations, phytoplankton community structure was estimated using the CHEMTAX program [54] following Endo et al. [16]. The initial and final matrices of pigment:Chl a ratios are shown in S1 and S2 Tables, respectively.

qPCR and qRT-PCR
Samples for DNA and RNA analyses were collected on days 0, 2 and 3. DNA samples (300-400 mL) were collected onto 0.2 μm pore size polycarbonate Nuclepore membrane filters (Whatman) with gentle vacuum (< 0.013 MPa), flash frozen in liquid nitrogen and then stored in a deep freezer at −80°C until analysis. Seawater samples (300-400 mL) for RNA analysis were filtered onto 0.2 μm pore size polycarbonate Nuclepore filters (Whatman) with gentle vacuum (< 0.013 MPa) and then stored in 1.5-mL cryotubes previously filled with 0.2 g of muffled 0.1-mm glass beads and 600 μL of RLT buffer (Qiagen) with 10 μL mL −1 β-mercaptoethanol (Sigma-Aldrich). The RNA samples were flash-frozen in liquid nitrogen and stored in a deep freezer at −80°C until analysis. The detailed methodologies of total DNA and RNA extractions were described in Endo et al. [16] and Endo et al. [27], respectively. The extracted RNA was reverse-transcribed into cDNA using the PrimeScript RT Reagent Kit with gDNA Eraser (Takara). To quantify the diatom-specific rbcL gene (DNA) and its transcripts (cDNA), qPCR and qRT-PCR were performed according to the method of Endo et al. [27].

Ion Torrent Personal Genome Machine (PGM) sequencing
Gene fragments of diatom-specific rbcL sequences were amplified from the extracted DNA or cDNA samples with barcoded fusion primer pairs, which included the primer set designed by John et al. [55]: forward primer, 5'-GATGATGARAAYATTAACTC-3'; reverse primer, 5'-TA WGAACCTTTWACTTCWCC-3'. The forward primer included the A-adaptor (5'-CCA TCTCATCCCTGCGTGTCTCCGAC-3'), key (5'-TCAG-3') and multiplex identifier (MID) sequences set by the manufacturer (Life Technologies), whereas the reverse primer included the truncated Pi-adapter (trP1: 5'-CCTCTCTATGGGCAGTCGGTGAT-3') sequence. Triplicate PCR amplifications were conducted for each samples using the TaKaRa Ex Taq Hot Start Version (Takara). The thermal cycling conditions were 60 s at 94°C for initial denaturation, then 30 cycles consisting of 10 s at 98°C, 30 s at 52°C, and 60 s at 72°C, followed by 10 min at 72°C for final extension. The PCR amplification was checked with 1.5% agarose gel electrophoresis. Amplicons were purified using AMPure beads (Beckman Coulter) and then quantified by an Agilent 2100 Bioanalyzer using the DNA 1000 Kit (Agilent Technologies) following the manufacturer's instructions. Then, the PCR templates were diluted to the final concentration of 26 pmol L −1 , and these were mixed equally in terms of concentration. Emulsion PCR was performed using an Ion One Touch 2 system with Ion PGM Template OT2 200 kit (Life Technologies) according to the manufacturer's protocol. The emulsion PCR products were enriched using an Ion One Touch ES (Life Technologies) and then loaded onto an Ion 318 v2 chip (Life Technologies). Sequencing of the amplicon libraries was performed using an Ion Torrent PGM system with the Ion PGM sequencing 200 kit v2 (Life Technologies) according to the manufacturer's protocol.

Sequence analyses
2.7.1. Quality filtering. Sequences with low quality, polyclonal sequences, and sequences that did not match against the A-adapter were initially filtered out with the Torrent Suite Software (Life Technologies), and the data so obtained were exported as FASTQ files. The complete run files in each sample were deposited in the DDBJ Sequence Read Archive (DRA) with the accession number DRA003722. Additional quality controls were performed using the software FASTX-Toolkit (http://hannonlab.cshl.edu/fastx_toolkit/). Sequencing reads were excluded if they met any of the following criteria: (i) reads not containing the trP1 adapter sequence; (ii) reads not matching the reverse primer sequence; (iii) reads shorter than 112 bp; and (iv) reads with average quality score below 25. Sequence quality statistics of raw and quality-filtered libraries are shown in S3 Table. 2.7.2. Taxonomic assignment. For taxonomic classification, contig assembly was achieved using the software SeqMan NGen (DNASTAR). Fifty thousand reads for each library were assembled into contigs with 97% sequence identity. The representative sequences of each contig containing with 100 reads ( 0.2% for total reads) were compared with rbcL sequences deposited in GenBank database (http://www.ncbi.nlm.nih.gov) using the BLAST query engine. Contigs that have 93% sequence similarity with known rbcL sequences were classified into the following groups: Chaetocerotaceae, Coscinodiscaceae, Cymatosiraceae, Rhizosoleniaceae, Stephanodiscaceae, Thalassiosiraceae, Achananthaceae, Bacillariaceae, Naviculaceae, Fragilariaceae, unidentified diatoms, and other eukaryotes. Contigs related to two or more diatom families with the same sequence similarity were assigned as unidentified diatoms. Other eukaryotes consisted of the contigs that were closely related to the sequences from organisms other than diatoms deposited in GenBank. A phylogenetic tree constructed with rbcL reference sequences is shown in S4 Fig. 2.7.3. Sequence-based statistical analyses. Principal Component Analysis (PCA) was used to recognize the relationships among pCO 2 levels in the DNA or cDNA libraries by reducing the multivariate information using the software R (http://www.r-project.org/). The ordination was performed based on the matrix of relative compositions at the family level in the libraries. In addition, ten thousand reads in each library were extracted for diversity analyses by the software mothur v. 1.25.0 [56]. All libraries were merged and clustered into operational taxonomic units (OTUs) with 97% sequence similarity cutoff. Singleton OTUs were then removed from the libraries. Genetic diversity was assessed based on the number of OTUs, Shannon-Wiener index (H', [57]), and Simpson's index (1-D, [58]).

Significance test
Statistical analyses were performed with the software R (http://www.r-project.org). To evaluate the statistically significant differences among CO 2 treatments, Kruskal-Wallis one-way analysis of variance (ANOVA) was used. Holm's test for multiple comparisons was also used to identify the source of variance. For all tests, p < 0.05 was considered significant.

Carbonate chemistry and nutrients
The initial concentrations of TA and DIC in the seawater were 2239.4 and 2074.7 μmol kg -1 , respectively (Table 1). Initial values of calculated pH and pCO 2 were 8.1 and 333.1 μatm, respectively (Table 1). Our incubation system successfully created significant gradients in pCO 2 and pH (S1 Fig) in seawater without any variation in TA between CO 2 treatments. Significant gradients in pCO 2 were achieved within 0.6 days, and it took at least 1.5 days to reach the intended pCO 2 values (S1A Fig). The mean pCO 2 values after day 2 in the 180, 350, 750, and 1000 μatm CO 2 treatments were 232, 328, 628, and 822 μatm, respectively. The initial concentrations of nitrate, phosphate, and silicic acid in seawater were 13.72, 0.99, and 11.76 μmol L -1 , respectively, and the macronutrients remained until the end of the incubation (Table 1).

Phytoplankton pigments
In the initial seawater, the concentrations of fucoxanthin (Fuco) and alloxanthin (Allo), which are biomarkers of diatoms in the Oyashio region during spring [47] and of cryptophytes [59], respectively, were relatively high among the carotenoids detected in this study. During incubation, the net growth rates of Fuco were the highest among all of the chemotaxonomic pigments in all of the incubation bottles (Fig 2). Significant differences in the net growth rates of Fuco were found among the CO 2 treatments (Kruskal-Wallis ANOVA, p < 0.05), and the value for the control treatment was higher than that for the other treatments (Holm's test, p < 0.05) (Fig  2). The net growth rates of Chl a and Allo were not significantly different among the CO 2 treatments (Kruskal-Wallis ANOVA, p > 0.05).
The initial phytoplankton community was mainly dominated by cryptophytes and diatoms (38.1% and 34.8% contributions to the Chl a biomass, respectively) (Fig 3). Contributions of diatoms to the Chl a biomass increased (> 47% in all CO 2 treatments on day 3), whereas those of cryptophytes decreased (< 30%). The diatom contributions in the high-CO 2 treatments are significantly lower than those in the low or control CO 2 treatments on day 3 (Kruskal-Wallis ANOVA, Holm's test, 180 and 350 μatm > 750 and 1000 μatm, p < 0.05). On the other hand, contributions of cryptophytes, haptophytes, and pelagophytes to Chl a biomass were

Copy numbers of rbcL DNA and cDNA of diatoms
A significant positive correlation between Fuco concentration and the diatom-specific rbcL gene copy count was found during the incubation (R 2 = 0.910, p < 0.001, n = 25) (Fig 4). In addition, a significant negative correlation was found between seawater pCO 2 and diatom-specific rbcL cDNA copy count on day 3 (R 2 = 0.506, p < 0.01, n = 12) (Fig 5A). However, no significant correlation existed between the seawater pCO 2 and the diatom-specific rbcL cDNA abundance normalized to the DNA abundance (cDNA/DNA) on day 3 (R 2 = 0.035, p > 0.05, n = 9) (Fig 5B).

Molecular OTU-based analyses of diatom-specific rbcL DNA and cDNA
The sequence data ranged from 340,809 to 1,194,481 reads per sample were generated from the Torrent Suite software, and approximately 80% of the raw sequences were removed by the subsequent procedures (S3 Table). As a result of the data filtering, average quality score of our sequence data increased from ca. 27.7 to ca. 31.4, indicating that the average error rate of our libraries was below 0.1%. The rarefaction curves plateaued when 10,000 sequences were  Table 3). The values of H' and 1-D were significantly higher in the 180 and 350 μatm CO 2 treatments than in the 750 and 1000 μatm CO 2 treatments on day 3 (t-test, p < 0.05). For the DNA and cDNA libraries, more than 85% of the sequences were closely related to the known diatom sequences with 93% sequence similarity. In the initial seawater, the dominant diatom groups in the rbcL DNA were Bacillariaceae and Fragilariaceae (52.0% and 18.3%, respectively), followed by unidentified diatoms, Chaetocerotaceae, Thalassiosiraceae, and Coscinodiscaceae (5.6%, 3.7%, 3.4%, and 1.9%, respectively) ( Fig 6A). Initial rbcL cDNA consisted mainly of Bacillariaceae and Chaetocerotaceae (43.8% and 15.0%, respectively), followed by unidentified diatoms, Fragilariaceae, Coscinodiscaseae and Thalassiosiraceae (12.4%, 11.0%, 3.2% and 2.9%, respectively) (Fig 6B). At the end of the incubation, the contributions of Thalassiosiraceae and Fragilariaceae to the total population increased relative to the initial contribution for all of the CO 2 treatments in the DNA and cDNA libraries (Fig 6). The contributions of Bacillariaceae to the total became higher in the 180, 750 and 1000 μatm CO 2 treatments compared with the control (Fig 7A and 7B). On the other hand, the contributions of Fragilariaceae, Thalassiosiraceae, and Chaetocerotaceae to the total decreased in the 180, 750 and 1000 μatm CO 2 treatments relative to the control on day 3. Percent differences in Bacillariaceae, Fragilariaceae, Thalassiosiraceae, and Chaetocerotaceae were higher in the cDNA than in the DNA. The contributions of Coscinodiscaceae to the total population increased in the low CO 2 treatment relative to the control in both the DNA and cDNA libraries on day 3. In the PCA plot (Fig 8), libraries were divided into DNA and cDNA by the first axis (PC1), and these were further divided to high (750 and 1000 μatm) and other (180 and 350 μatm) CO 2 treatments by the second axis (PC2), except for the 750 μatm CO 2 treatment in the DNA library. In the PCA analysis, the first and second axes from the ordination explained 50.5% and 27.6% of the total variability, respectively (cumulative contribution of 78.1%).

Chemical and biological characteristics of the incubation seawater
In our experiment, initial seawater pCO 2 was intermediate value between pre-and post-bloom phases in the Oyashio region reported by Midorikawa et al. [20]. In addition, in the initial Response of Diatoms to CO 2 Availability seawater, macronutrients were not depleted, and the Chl a concentration was lower than the values previously reported during the spring diatom blooms in this region (> 5 μg L -1 , [20,47,60]). These results indicate that the intense diatom bloom has not yet occurred in the seawater used for the incubation experiment. This speculation is supported by our CHEMTAX analysis, which revealed that the initial phytoplankton assemblages were dominated by a mixture of diatoms and cryptophytes (Fig 3) as previously observed during the pre-bloom period [47,60]. The initial value of F v /F m was also similar to that in the pre-bloom phase in spring Oyashio water measured with fast repetition rate fluorometry (FRRF) by Yoshie et al. [61]. Concentrations of Chl a increased exponentially during the incubation (S2 Fig) and it was accompanied by the increases in diatom contribution to the Chl a biomass (Fig 3), suggesting that our incubation could simulate the early growth phase of the diatom bloom.

Effects of CO 2 availability on the phytoplankton pigments and F v /F m
The net growth rates of Fuco were less in the low (180 μatm) and high (750 and 1000 μatm) CO 2 treatments relative to the control (350 μatm), while no significant CO 2 difference was  observed in the other chemotaxonomic pigments (Fig 2), suggesting that the growth of diatoms could be selectively diminished by changes in CO 2 availability. This assumption is supported by size-fractionated Chl a analysis, which showed that elevated CO 2 levels selectively decreased large-sized (10 μm) phytoplankton cells, that were most likely composed of diatoms. Our results were inconsistent with the previous laboratory culture experiments, which indicated that diatoms tended to increase their growth under elevated CO 2 levels [10,11]. In addition, Wu et al. [62] demonstrated that the CO 2 enrichment stimulated the growth of larger diatom species rather than that of smaller cells. However, field studies using natural phytoplankton assemblages showed positive [12,13,28] and negative [17,27,63] effects of increased CO 2 availability on diatom growth. This pattern may have been caused by differences in diatom species [4,19,62], an effect of CO 2 on zooplankton grazing [64], and/or the bioavailability of iron in response to the change in carbonate chemistry [18]. Our findings further support that the sensitivity of phytoplankton assemblages to elevated CO 2 differs geographically, according to the differences in species composition and environmental conditions [5]. In the present study, decreased contributions of diatoms resulted in the relative increases in other phytoplankton taxa, such as cryptophytes and haptophytes (Fig 3). Our results suggest that CO 2 enrichment may decrease the diatom contribution to the phytoplankton assemblages during the growth phase of spring blooms in the Oyashio region. However, it should be noted that our incubation period was short (i.e., 1-2 days under fully manipulated CO 2 conditions) and may not have allowed for full acclimation of phytoplankton community to changes in CO 2 . Therefore, the biological responses observed in our study could be transient, and further work with a longer incubation period is also required to refine the results observed in our study. Differences in F v /F m values among treatments were not significant throughout the incubation (S3 Fig), indicating that CO 2 availability had a minimal effect on the photochemical quantum efficiency of PSII for phytoplankton assemblages. Similar results were also obtained from Fe-fertilized field experiments [16,18].

Effects of CO 2 availability on the diatom-specific rbcL gene and its transcripts
The copy number of rbcL gene targeting diatoms exhibited a good correlation with Fuco concentration in this study (Fig 4). This indicates a close coupling between rbcL gene and the chemotaxonomic pigment in diatoms. A similar result was also observed during the CO 2 - controlled incubation experiment in the oceanic Bering Sea [27]. Given that Fuco can be a strong indicator of diatom carbon biomass in the Oyashio region [47], our result suggests that the diatom-specific rbcL gene could also serve as a potential indicator for diatom carbon biomass in the area.
The abundance of diatom rbcL gene transcripts decreased with an increase in pCO 2 ( Fig  5A), indicating that RubisCO expression could decrease in high CO 2 environments. However, DNA-normalized rbcL transcription did not vary as a function of pCO 2 level, indicating that the transcriptional activity was unaffected by CO 2 availability. Pigment analysis by HPLC also indicated the decrease in diatom biomass under high CO 2 conditions (Figs 2 and 3). Therefore, the decrease in rbcL transcript abundance in diatom community might be caused primarily by the decrease in diatom biomass rather than the decrease in transcriptional activity in diatom cells. A laboratory culture experiment using the diatom Thalassiosira weissflogii also showed a decrease in RubisCO concentration with an increase in CO 2 concentration over a range of 182 to 750 ppm [65]. Decreases in rbcL transcripts or RubisCO concentrations under elevated CO 2 conditions were also observed from the field incubation experiments conducted off the coast of California [65] and in the Bering Sea [27]. However, there is little evidence on CO 2 -induced transcriptional regulation of rbcL in diatom cells. Moreover, Losh et al. [65,66] demonstrated that the growth and productivity of phytoplankton assemblages were not affected by increased CO 2 level, although the RubisCO expression decreased. Their study suggested that the decrease in RubisCO content could be compensated by an increase in CO 2 concentration around RubisCO at higher CO 2 levels. In our study, therefore, the decrease in diatom growth rate under high CO 2 levels might be unrelated to the transcriptional response of rbcL gene. One possible mechanism for the decrease in diatom biomass is that the increase in photoinhibition under elevated CO 2 levels [62]. The shift in diatom community composition could also affect the abundance of diatom-specific rbcL and its transcripts (Fig 7).
Although our results cannot extrapolate to longer-term behaviors in the study area, shortterm incubation may be suitable for determining physiological responses in the natural phytoplankton assemblages, because that can reduce artifacts derived from incubation (i.e., so-called bottle effects; [67,68]). In addition, environmental parameters other than CO 2 (e.g., nutrient availability) would be differentiated among CO 2 treatments after longer incubation period [16,69], and that complicated the source of the algal responses.
Our rarefaction curves reached plateaus with 10,000 reads for all samples (S5 Fig), suggesting that the sequencing efforts were sufficient to address the comparative analyses among libraries. Taxonomic analysis based on NGS revealed that community composition of rbcL gene and its transcripts differed largely among CO 2 treatments on day 3. The result suggests that the diatom community structure shifted rapidly in response to changes in CO 2 levels, whereas transcriptional activity of total diatom community was little affected by CO 2 availability. In our NGS libraries, CO 2 -induced changes in the taxonomic composition of diatoms exhibited similar trends between DNA and cDNA. For example, the contributions of Chaetocerotaceae, Thalassiosiraceae, and Fragilariaceae consistently decreased in both the DNA and cDNA libraries in response to changes in CO 2 levels compared to the control, whereas those of the Bacillariaceae increased (Fig 7). A possible explanation for the consistency is that the decrease and increase in rbcL transcription could be caused by the decrease and increase in rbcL gene abundance, respectively. Our results indicate that the effects of CO 2 availability could differ among diatom taxa, resulting in changes in the community composition.
Elevated CO 2 availability decreased the OTU richness and diversity of diatom-specific rbcL genes and their transcripts in our experiment (Table 3). Although the decreases in diversity indices were rather small, the diatom community lost 5-10% in terms of the rbcL OTU richness during the incubation. The PCA ordination revealed that the rbcL libraries were clearly separated by high and low CO 2 treatments in both DNA and cDNA libraries (Fig 8), suggesting that the decreases in OTU diversity were accompanied with shifts in the diatom-specific rbcL composition.
In contrast to elevated CO 2 treatments, the diversity, taxonomic composition, and transcription activity of diatom-specific rbcL were weakly affected by decreases in CO 2 availability (Table 3, Figs 5 and 8), indicating that the diatom community in the study area was rather stable to low CO 2 conditions compared with higher CO 2 conditions. Previous field observations also showed that the diversity and productivity of diatom assemblages were maintained during the spring blooms in the Oyashio region [47,61]. Our results suggest that the dominant diatom species during the spring blooms in this region could be acclimated to the bloom-induced low CO 2 conditions every spring [20], whereas they may not acclimate well to a high CO 2 environment that had not been experienced before. The acclimation would be accomplished by the use of CCMs mediated by the enzyme carbonic anhydrase (CA; [21,22]). Burkhardt et al. [70] reported that the activity of CA increased in the low CO 2 (36-180 μatm) treatments compared to the ambient (360 μatm) treatment. Furthermore, it has been proposed that bloom-forming diatom species possess highly efficient and tightly regulated CCMs to maintain their growth even under low CO 2 environments [71]. However, decreases in the growth rate of Fuco and the number of rbcL OTUs observed in the low CO 2 treatment, implying that some diatom clades could not respond to abrupt decrease in CO 2 availability. Because there is inter-specific variation in the activity and plasticity in CCMs among diatoms [72,73], one possible explanation of our results is that diatom species with low CCM capacity were eliminated in the low CO 2 treatments.

Conclusions
In the present study, we focused on the effects of different CO 2 availability on the growth, community composition, and photosynthetic capacity of diatoms in the spring Oyashio waters of the western North Pacific. The comprehensive NGS analyses for the diatom-specific rbcL gene and its transcripts enabled us, for the first time, to detect significant changes in taxonomic composition and diversity of diatoms among CO 2 treatments even under the short-term incubation. Consequently, we revealed that the elevated CO 2 availability decreased the OTU richness and diversity of diatom-specific rbcL gene and their transcripts. Interestingly, the increase in seawater CO 2 levels reduced the growth rate of diatoms and their relative contribution to the Chl a biomass as estimated from photosynthetic pigment signatures. In addition, the transcript abundance of diatom-specific rbcL gene also decreased under elevated CO 2 levels. According to Chiba et al. [74], a decrease in springtime Chl a concentration accompanied the changes in the predominant diatoms from centrics to pennates during the 1960s and 1990s in the Oyashio region. Additionally, the bloom dynamics could be controlled by physical factors such as the Pacific Decadal Oscillation (PDO, [75]). Therefore, it is crucial to investigate the relationship between environmental changes in Oyashio waters and the dynamics and mechanisms of its spring diatom blooms. We believe that molecular biological techniques can help us to gain a deeper understanding of the changes in the composition and physiology of spring Oyashio diatoms.  Table. Quality statistics of pre-and post-quality control sequence data. Sequence data were obtained from DNA and cDNA samples collected on days 0 and 3. (DOCX)