The S. costatum strain was kindly provided by the Culture Collection Center of Marine Algae, Xiamen University, China. S. costatum cells were maintained in K-medium with 48 μM silicon at 20°C under a 14 h:10 h light:dark photoperiod at a light intensity of approximately 100 μmol m^-2 s^-1 provided by fluorescent lamps (Keller et al., 1987). Before the experiment, a mixture of antibiotics containing penicillin G (1 g L^-1), neomycin (250 mg L^-1), gentamicin (1 g L^-1), and kanamycin (0.5 g L^-1) was added to the culture media to eliminate bacterial contamination, and the culture was checked periodically for contamination with 4′,6-diamidino-2-phenylindole stain by microscopic inspection (Cottrell and Suttle, 1993; Ou et al., 2008).

The experiment included four treatments: inorganic P-replete, P-deficient, DIP-resupplied, and DOP-resupplied. Each group had triplicate biological repeats. S. costatum cells in the exponential growth phase were collected using centrifugation (2,500 g for 15 min at 20°C), then washed twice with sterile seawater, and finally cultured in K-medium without P for 2 days to exhaust intracellular P. Then the cultures were inoculated into 12 bottles each with 5 L of culture medium to yield an initial density of 9.0 × 10³ cells mL^-1. In P-replete cultures, 10 μM Na₂HPO₄ was added to the culture media at the beginning of the experiment. For P-deficient cultures, 0.2 μM was added to the culture media to maintain cell activity. At day 4, three P-deficient cultures were resupplied with phosphate to the final concentration of 10 μM as the DIP-resupplied group, and another three P-deficient cultures were resupplied with glucose-6-phosphate (G-6-P) to the final concentration of 10 μM as the DOP-resupplied group. The remaining three bottles were maintained as the P-deficient group.

Physiological parameters including cell density, the photochemical efficiency of photosystem II (Fv/Fm), DIP, DOP, particulate P (PP), and bulk alkaline phosphatase activity (APA) were monitored daily. In addition, samples were collected 4 h and 28 h after the resupply of DIP and DOP.

Three 1 mL aliquots of each bottle were collected every day and fixed in 2% Lugol’s solution for subsequent cell counting under a light microscope (Zhang et al., 2014). The specific growth rate of S. costatum was calculated using the following equation: μ = (ln N₂ – ln N₁)/(t₂ – t₁), where N₁ and N₂ were the cell densities at time t₁ and t₂ (Ou et al., 2014). Three 5 mL aliquots of each bottle collected at 11:00 am each day were dark acclimated for 15 min and the Fv/Fm was measured using PHYTO-PAM (Pulse Amplitude Modulation, ED, Walz, Effeltrich, Germany; Zhang et al., 2007). The bulk APA of S. costatum was measured using the method reported by Ou et al. (2015). Cells were collected on pre-combusted GF/F filters (450°C, 2 h) for PP analysis using the method of Solórzano and Sharp (1980), which employs magnesium sulfate (MgSO₄) and digested for 2 h. The filtrates were used for DIP and DOP analyses; the concentration of DIP was measured using the molybdenum blue method described by Murhpy and Riley (1962) and the concentration of DOP was analyzed using the methods reported by Jeffries et al. (1979). DOP was calculated from the difference between DIP and total dissolved phosphorus (TDP) which employs acid potassium persulfate (K₂S₂O₈; Ou et al., 2015).

Skeletonema costatum cells were collected for transcriptomic analyses during exponential growth phase in the P-replete cultures (day 3). The P-deficient cells were harvested on day 4 as well as 4 h and 28 h after the resupply of DIP and DOP. For each sample, cells were collected onto the polycarbonate membrane filter (pore-size 3.0 μm, Millipore) and resuspended in 1 mL Trizol Reagent (Invitrogen, Carlsbad, CA, USA), immediately frozen in liquid nitrogen and subsequently stored at -80°C for RNA isolation (Zhang et al., 2014).

Total RNA was isolated using TRI-Reagent (MRC, Cincinnati, OH, USA) and dissolved in RNase-free water as previously reported by Zhang et al. (2014). For each sample, equal amounts of total RNA from three replicates were mixed together for transcriptome sequencing (Xia et al., 2011; Yu et al., 2016). RNA-Seq libraries were constructed using the Illumina TruSeq^TM RNA Sample Preparation Kit (Illumina, San Diego, CA, USA) following the Illumina TruSeq RNA-Seq library protocol. Poly(A) mRNA was enriched using poly-T oligo-attached magnetic beads. First-strand cDNA was synthesized using random oligonucleotides and SuperScript II (Life Technologies, Carlsbad, CA, USA). Second-strand cDNA synthesis was subsequently performed using DNA polymerase 1 (DNAP I; New England BioLabs) and RNase H (Invitrogen; Zhang et al., 2014). In total, one individual paired-end cDNA library was constructed for each set of samples. The cDNA library quality was checked using an Agilent high sensitivity DNA assay on an Agilent Bioanalyzer 2100 system (Santa Clara, CA, USA). Subsequently, the libraries were sequenced on the Illumina Hiseq^TM 2000 platform (Expression Analysis Inc., San Diego, CA, USA).

Before assembling, raw reads with adapter contamination, low quality reads, and reads with unknown nucleotides (>5%) were removed. On average, the Q20 percentage of clean reads from all six samples was ∼98. All the downstream analyses were based on clean reads. Next, all the clean reads were assembled as contigs based on the overlap of short reads using the software Trinity de novo assembler (Release-20130225¹) with the min_kmer_cov_set to 1 and all other parameters set to default (Grabherr et al., 2011; Cheng et al., 2014; Zhang et al., 2014).

The unigenes extracted with in-house Perl scripts were aligned against the NCBI-NR database (Release-20130408), NCBI-NT database (Release-20130408), COG (Release-20090331), Swiss-Prot database (Release-63.0), and KEGG database (Release-63.0) using BLAST (v2.2.26 + x64-linux) with a threshold e-value ≤ 1e^-5 as in Zhang et al. (2014).

The RNA-Seq clean reads of each sample were mapped back to the transcript reference database assembled with Trinity, and the expression levels of unigenes were calculated using the number of mapped reads as EdgeR inputs², followed by normalization of read count number to fragments per kilobase of transcript per million mapped reads (Mortazavi et al., 2008; Cheng et al., 2014; Zhang et al., 2014). Differentially expressed genes (DEGs) were identified through 13 pair-wise comparisons with different physiological/ecological relevance (Supplementary Table S1). According to the comparison method developed by the Beijing Genomics Institute (BGI, China), the probability of one gene being expressed equally between two samples was judged according to the p-value corresponding to the differential gene expression test and false discovery rate (FDR; Audic and Claverie, 1997; Benjamini and Yekutieli, 2001; Li et al., 2014). A FDR of 1% or less and a fold change ≥4 were set as the threshold for significant differential expression. DEGs were identified as enriched in GO terms (p ≤ 0.05) and metabolic pathways (q ≤ 0.05) by searching against GO and KEGG, respectively (Kanehisa et al., 2006, 2008; Sun et al., 2014).

The total RNA of each sample used for quantitative RT-PCR (qRT-PCR) analysis was extracted as described above. The expression levels of alkaline phosphatase (scoap), inorganic phosphate transporter (PiT), phospholipase D (PLD), cryptochrome 1 (cry 1), and phytochrome B (phy B) were examined using qRT-PCR. First-strand cDNA was synthesized using the FastQuant cDNA RT Kit with gDNase (TIANGEN, Beijing, China). qRT-PCR was performed on an ABI 7500 System (Applied Biosystems) using a SuperReal PreMix Plus (SYBR Green) Kit (TIANGEN, Beijing, China). Thermocycling was conducted as follows: 50°C for 2 min, 95°C for 10 min, 40 temperature cycles at 95°C for 30 s, and 60°C for 60 s. The primers designed for qRT-PCR in this study are listed in Supplementary Table S2. Calmodulin (calm), a commonly used housekeeping gene (Shi et al., 2013; Zhang et al., 2014), was chosen as the reference gene to normalize the expression of scoap, PiT, PLD, cry 1, and phy B. Relative expressions of DEGs were calculated based on the 2^-ΔΔCt relative response method (Zhang et al., 2014).

A Student’s t-test was performed to compare the differences between the control (P-replete) and each treatment group for cell density, Fv/Fm, DIP, DOP, PP, and bulk APA; a P-value < 0.05 was regarded as a significant difference. Prior to analysis, data were tested for the equality of variances. All tests were performed using SPSS 17.0 (SPSS Inc., Chicago, IL, USA).

Physiological responses of S. costatum to P-depletion and -resupply are shown in Figure 1. For the P-replete cultures, cell density increased rapidly and reached a peak (6.5 × 10⁵ cells mL^-1) at day 4, and then entered into stationary phase (Figure 1A). In the other three groups with 0.2 μM P, the growth trends were similar to the P-replete group in the first 3 days with low cell density. From day 4, cells in the P-deficient group maintained a stable level with a cell density ∼4.0 × 10⁵ cells mL^-1. However, growth of P-deficiency cultures recovered rapidly after P resupply, and cell densities reached high levels of 9.0 × 10⁵ cells mL^-1 at day 6 (Figure 1A).

In the P-replete group, Fv/Fm reached a peak (∼0.5) at day 3 and then declined gradually to a low value (less than 0.1). Fv/Fm of the P-deficient group decreased from day 3 to the end of the experiment. For the P-resupplied group, Fv/Fm began to recover after DIP or DOP resupply at day 4. However, Fv/Fm of both groups did not reach the maximum value (∼0.5) as that obtained in the P-replete group (Figure 1B).

The concentration of soluble reactive P decreased rapidly in all cultures. In the P-replete group, P concentration decreased to 6.0 μM at day 1, and was undetectable at day 2 (Figure 1C). However, bulk APA varied a little from day 1 to 6 with a slight increase in the last 2 days (Figure 1D). In the other three treatment groups, APA increased from day 1 to day 4 before the resupply of DIP or DOP (Figures 1C,D). After resupplement of P, both DIP and DOP concentrations decreased rapidly and DIP was exhausted after 28 h of P-resupply. APA decreased rapidly in both DIP and DOP-resupplied groups until the end of the experiment (Figures 1C,D).

Variations of PP contents are shown in Figure 1E. PP contents decreased rapidly in the first 3 days in all groups, and then maintained stable levels in both the P-replete and P-deficient groups at day 4, but PP contents in P-deficient cells were much less than those in P-replete cells (Figure 1E). After resupply of P, PP contents increased rapidly in P-resupplied-4 h cells but decreased in P-resupplied-28 h cells until the end of the experiment (Figure 1E).

There were no significant differences between DIP- and DOP-resupplied groups, indicating that S. costatum could utilize both inorganic and organic P as the P source for cell growth.

In this study, six samples were sequenced with an average read length of 90 bp (Table 1), and all clean reads have been deposited in the SRA database of GenBank³ with the BioProject accession number PRJNA313486. After removing low quality reads and trimming, poly-N and adapter sequences, the Q20 of the clean reads data ranged from 98.05 to 98.42% and the GC content of either left end or right end read was constant at approximately 47% for each sample. Using Trinity software, these clean reads were de novo assembled to 34,942 unique unigenes (transcripts): 31,959 unigenes were obtained from P-replete; 34,590 from P-deficient; 31,530 from DIP-resupplied-4 h; 31,562 from DOP-resupplied-4 h; 34,679 from DIP-resupplied-28 h; and 34,253 from DOP-resupplied-28 h (Table 1). The average length of unigenes was 1,430 bp with an N50 length of 2,157 bp (Supplementary Figure S1).

Overall, 26,374 unigenes were successfully annotated in at least one database (Table 2). A total of 26,102 unigenes presented significant similarity to known proteins in the NCBI-NR database. Among them, 44.6% of the annotated unigenes were matched with sequences from the diatom Thalassiosira pseudonana CCMP 1335 and 38.6% with sequences from the diatom T. oceanica (Supplementary Figure S2).

Functional classification of all unigenes was conducted using the COG and GO databases. 11,063 unigenes were annotated against the COG database and classified into 25 categories based on sequence homology (Figure 2A). Among these categories, the largest group was the most common and the non-specific category of general function prediction only (41.5%), followed by transcription (19.6%), cell wall/membrane/envelope biogenesis (17.7%), and carbohydrate transport and metabolism (17.6%).

Of the 26,102 unigenes, 8,289 unigenes were classified into three ontologies and 49 sub-categories (Figure 2B). Among the ontology of biological processes, metabolic process (54.4%) and cellular process (52.2%) were the two dominant groups, followed by single-organism process (25.7%), response to stimulus (13.6%), localization (12.1%), and establishment of localization (11.9 %).

KEGG annotation revealed that 12,090 unigenes were annotated against the KEGG database (Table 2). In the second level, all annotated unigenes were classified into 18 categories (Supplementary Figure S3B), and most of the unigenes were associated with translation (16.5%), carbohydrate metabolism (14.4%), and lipid metabolism (12.6%), followed by nucleotide metabolism (6.3%), environmental adaption (3.1%), and membrane transport (1.0%). The top 25 abundant biochemical pathways with numbers of assigned unigenes are shown in Supplementary Figure S3A.

Differentially expressed gene were identified through 13 pair-wise comparisons and approximately 20.8% of the unigenes presented significant differential expressions (Figure 3). Among them, 16.8% were significantly up-regulated and 4.0% were significantly down-regulated in P-deficient cells compared with P-replete cells. Compared with the P-deficient cells, 14.0% (∼5.0% up-regulated and ∼9.0% down-regulated) and 9.8% (∼4.5% up-regulated and ∼5.3% down-regulated) of the unigenes exhibited significant differences in the DIP-resupplied-4 h and DOP-resupplied-4 h cells while 19.1% of DEGs shared between these two groups and 8.46% of DEGs shared among P-replete, DIP-resupplied-4 h, and DOP-resupplied-4 h (Figure 3A). After P-resupply for 28 h, 14.1% of the unigenes (∼3.4% up-regulated and ∼10.7% down-regulated) in the DIP-resupplied-28h cells and 20.4% of the unigenes (4.4% up-regulated and 16.0% down-regulated) in the DOP-resupplied-28 h cells presented significant differences compared with the P-deficient cells, and 12.5% (∼7.9% up-regulated and ∼4.6% down-regulated) and 17.2% (∼7.5% up-regulated and ∼9.7% down-regulated) of unigenes showed significant difference in expression compared with the P-replete cells, indicating that cells had not recovered to normal homeostasis. Compared with the P-deficient cells, 21.36% of DEGs were shared among P-replete, DIP-resupplied-28 h, and DOP-resupplied-28 h, and 38.64% of DEGs were shared between DIP-resupplied-28 h and DOP-resupplied-28 h cells (Figure 3A). Comparison of the DIP-resupplied-28 h and DOP-resupplied-28 h cells showed that only a small number of the unigenes presented significant differential expressions (Figure 3B).

KEGG pathway analysis revealed that all DEGs were assigned to 124 specific pathways. Among the 15 most frequently represented pathways, the majority of the DEGs were enriched in those pathways which needed P or which were related to cell growth, such as RNA transport, glycerophospholipid metabolism, ether lipid metabolism, nucleotide metabolism, and ribosome biogenesis (Supplementary Table S3; Supplementary Figure S4).

Cells can increase the competitive advantage for phosphate by inducing higher affinity transporters and/or synthesizing more transporters under P-deficient conditions (Dyhrman et al., 2012). In the diatom T. pseudonana, a phosphate transporter (PID: 24435) and a high-affinity phosphate transporter are up-regulated at transcription and/or protein level in response to P deficiency (Dyhrman et al., 2012; Fu et al., 2013). The phosphate transporter gene expression pattern of P-deficient S. costatum is also consistent with its homolog in Skeletonema spp. and T. rotula to P-limitation in the field (Alexander et al., 2015). High and low transport systems have also been identified in S. costatum (Ou, 2006). In our study, transcripts of phosphate transporters PiT (major facilitator super-family), the SLC20 family and SLC25A3 were significantly up-regulated in P-deficient cells (Supplementary Table S4), and their expression patterns were consistent with the homologs in Skeletonema spp. and T. rotula in the field (Alexander et al., 2015). As a cell surface transporter, the SLC20 family plays a fundamental housekeeping role in phosphate transport, such as absorbing phosphate from interstitial fluid (Ravera et al., 2007). In plants and fungi, essential transporters for phosphate uptake are proton-coupled transporters of the major facilitator super-family, which also function as transceptors to signal external phosphate concentration (Popova et al., 2010; Nussaume et al., 2011; Pedersen et al., 2013). However, we did not detect any genes encoding high-affinity phosphate transporters in our study. Thus, S. costatum might produce different transporters in response to ambient P deficiency, and PiT might play an important role in extracellular P sensing.

In the ocean, polyphosphate is thought to be the product of luxury uptake and storage of phosphate in phytoplankton (Diaz et al., 2008; Dyhrman et al., 2012; Martin et al., 2014). Proteins containing the SPX domain play an important role in maintaining intracellular phosphate homeostasis, and three polyphosphate synthase subunits (vacuolar transporter chaperone 2 (Vtc2), Vtc 3, and Vtc 4) all harbor the SPX domain in yeast (Secco et al., 2012). In our study, the transcript of the SPX domain-containing protein involved in vacuolar polyphosphate accumulation was significantly up-regulated in P-deficient cells (Supplementary Table S4), suggesting that S. costatum might have the ability to store P. Our results on intracellular P contents and other physiological responses in P-replete and P-resupplied cells supported this speculation (Figure 1). S. costatum might increase P allocation to polyphosphate under P deficient conditions, which was consistent with the findings in P-limited T. pseudonana (Dyhrman et al., 2012). Our result also supported the view that not all diatom polyphosphate allocation is driven by luxury uptake of phosphate (Dyhrman et al., 2012).

Phytoplankton are able to utilize organic P, and hydrolysis of cell surface P esters by AP is considered to be the most common DOP utilization mechanism (Beszteri et al., 2012; Dyhrman et al., 2012; Lin et al., 2012). In our study, both the activity and the transcript of AP was up-regulated in P-deficient cells (Supplementary Table S4), which was consistent with the findings in T. pseudonana and Karenia brevis (Dyhrman et al., 2012; Lin et al., 2012). Phosphomonoester can be hydrolyzed directly by a phosphomonoesterase, such as AP, whereas the hydrolysis of high molecular weight phosphodiesters requires phosphodiesterase and polyphosphatase (Cembella et al., 1984; Ou et al., 2015). Genes for several phosphodiesterases were significantly up-regulated in P-deficient cells, especially the PLD (Supplementary Table S4). PLD is a key regulator of cytoskeletal organization and can hydrolyze structural phospholipids (such as membrane lipids), regulating a diverse range of cellular processes, such as membrane transport and cell migration (Pleskot et al., 2013; Frohman, 2015). Its transcription and activity increase upon exposure to various stresses, such as cold, drought, and salinity (Pleskot et al., 2013). Our results indicated that S. costatum could utilize both simple organic P (i.e., G-6-P) (Figure 1; Supplementary Table S4) and complex organic P (i.e., membrane lipids) as a P source under P-deficient conditions.

Eukaryotic phytoplankton and cyanobacteria are able to replace phospholipid with non-P containing sulfolipids (sulfur containing) and betaine lipids (nitrogen containing) in a P scarce environment to decrease the requirement of cells for P (Yu et al., 2002; Van Mooy et al., 2009; Dyhrman et al., 2012). In our study, transcripts of sulfoquinovosyltransferase, UDP-sulfoquinovose synthase, and desulfoglucosinolate sulfotransferase which are involved in sulfolipid biosynthesis were significantly up-regulated in P-deficient cells (Supplementary Table S4) (Sanda et al., 2001; Yu et al., 2002). Moreover, genes encoding the betaine lipid synthase involved in betaine biosynthesis were also significantly up-regulated in P-deficient cells (Supplementary Table S4). These results indicated that S. costatum could utilize non-P containing lipids to reduce the demand for P, which might be an adaptive response of cells to ambient P deficiency. This finding, coupled with a similar response in the diatoms, T. pseudonana and Chaetoceros affinis, suggests diatoms can utilize non-P containing lipids to reduce their cellular P demand (Van Mooy et al., 2009; Dyhrman et al., 2012).

A circadian rhythm occurs ubiquitously in both prokaryotes and eukaryotes driven by a circadian clock, which is entrained by light, temperature, iron, and nitrogen signals in plants (Gutiérrez et al., 2008; Bednářová et al., 2013; Hong et al., 2013; Salomé et al., 2013; Larrondo et al., 2015). Light is regarded as the most important environmental factor involved in resetting the circadian clock (McClung, 2006; Roy and Morse, 2013; Hurley et al., 2014). Studies show that photoreceptors play a significant role in light energy capture and the circadian oscillator mechanism (Bognár et al., 1999). Phytochromes and cryptochromes are the two important photoreceptor families transducting light signal input to the circadian clock (Millar, 2004; McClung, 2006). Furthermore, many light-dependent processes controlled by phytochrome and cryptochrome are also regulated by a circadian rhythm (Millar et al., 1995; Bognár et al., 1999; Millar, 1999). In our study, transcripts of phy B and cry 1 were significantly up-regulated in P-deficient cells and declined rapidly after P resupply (Supplementary Table S5). In Arabidopsis, the light period of the phy B-deficient mutant is 1.5–2 h longer than that of the wild type; overexpression (15-fold) of phy B shortens the light period length; cry 1 also plays a similar role as phy B in regulating clock length (Somers et al., 1998). Phy B is the primary high-intensity red light photoreceptor for circadian control, while cry 1 mediates high-intensity blue light signals for the control of period length (Somers et al., 1998). Hence, the light period length might be shortened and the input pathway of the circadian clock influenced in P-deficient S. costatum (Figure 4).

The core oscillator of the circadian clock is composed of three interlocked feedback loops, and CCA1 is one of the two domain transcription factors (the other is Late Elongated Hypocotyl, LHY) and participates in each loop (Figure 4) (McClung, 2006). Furthermore, CCA1 plays a central role in the circadian control of gene expression, and is involved directly in light regulation of the gene expression in plants, providing a molecular link between phytochrome and the circadian oscillator in plant cells (Green and Tobin, 1999). Lack of the CCA1 can cause a shortened circadian period, but its overexpression results in severe disruption of the normal circadian function in plants (Wang and Tobin, 1998; Green and Tobin, 1999; Fujimori et al., 2005; Scaglioni et al., 2008). In our study, the CCA1 gene was up-regulated approximately 1.66-fold in P-deficient cells and down-regulated 1.9-fold after P resupply (Supplementary Table S5). Moreover, genes encoding CK2 subunit alpha and beta and PRR5 were also significantly up-regulated in P-deficient cells, and down-regulated in P-resupplied cells (Supplementary Table S5). CK2 can phosphorylate two components of the central oscillator in Arabidopsis, CCA1, and LHY (Sugano et al., 1998; Green and Tobin, 1999; Mulekar et al., 2012). CCA1 phosphorylation by CK2 is necessary for maintaining the normal function of the central oscillator in Arabidopsis (Daniel et al., 2004). PRR5 belongs to the PRR family, which is a negative regulator mediating the expressions of CCA1 and LHY genes (Fujimori et al., 2005; Nakamichi et al., 2010; Salomé et al., 2010). It is interesting that overexpression of the PRR5 gene is also related to circadian-associated phenotypes, such as regulation of flowering time (Fujimori et al., 2005). Our results indicated that the phosphorylation of CCA1 was affected by P deficiency, and the core oscillator of the circadian clock was re-programed in response to ambient P change (Figure 4).

In our study, transcripts of RNase MRP1 were significantly up-regulated in P-deficient cells, and down-regulated in P-resupplied-28 h cells (Supplementary Table S5). RNase P/MRP, being localized in the nucleolus and cytoplasm, plays an important role in regulating the cell cycle of yeasts (Esakova and Krasilnikov, 2010). RNase MRP mutation occurring in Saccharomyces cerevisiae can delay the cell cycle at the end of mitosis (Cai et al., 2002; Esakova and Krasilnikov, 2010). In some unicellular organisms, such as the green alga Chlamydomonas reinhardtii, the flagellate Euglena gracilis, the cyanobacterium Synechococcus elongata, and the dinoflagellate Gonyaulax polyedra, cell division is timed by a circadian mechanism (Hunt and Sassone-Corsi, 2007). However, the role of RNase P/MRP in regulating the circadian clock is poorly understood. In general, the circadian rhythm of S. costatum was disturbed by ambient P deficiency which might subsequently initiate adaptive mechanisms to the ambient P change in accordance with the signals released by the new circadian clock (Figure 4). However, the detailed regulation mechanism still needs further study.

In our study, 33 DNA-directed RNA polymerase (RNAP) I, II, and III genes were up-regulated in P-deficient cells and down-regulated in P-resupplied cells (Supplementary Table S6). DNA-directed RNAPs are complex enzymes containing multiple subunits, i.e., RNAP I, II, III, IV, and V, and are necessary for constructing RNA chains using DNA genes as templates (Jones et al., 1987; Ishihama, 2000; Iyer et al., 2003; Kwapisz et al., 2008). RNAPs regulate the process of gene transcription which allows a cell to acclimate to a changing environment (Ishihama, 2000). However, RNAP genes varied insignificantly in nitrogen-, iron-, or silicon-deficient T. pseudonana cells (Mock et al., 2008). These results suggest that variation of RNAP gene expression might be a specific response of diatoms to ambient P-deficiency. In addition, three genes encoding DNAP, including DNAP eta subunit, DNAP alpha subunit A, and DNAP I, were also up-regulated in P-deficient cells and down-regulated in P-resupplied cells (Supplemenatry Table S6). DNAP is an enzyme responsible for DNA replication by assembling nucleotides during cell division. DNAP eta can correct common defects in DNA, particularly important for accurate translesion synthesis of DNA damage resulting from ultraviolet radiation (Goodsell, 2004). Xanthine oxidase and urate oxidase are two key enzymes involved in the degradation of purine bases (Andersen et al., 2006; Abooali et al., 2014). Moreover, transcripts of enzymes involved in the synthesis of nucleotides, such as xanthine oxidase, urate oxidase, GTP synthase, UMP-CMP kinase, phosphoribosylamine-glycine ligase, IMP dehydrogenase, and adenine phosphoribosyltransferase (Van Rompay et al., 1999; Ingerson-Mahar et al., 2010; Sampei et al., 2010), were all up-regulated significantly in P-deficient cells (Supplementary Table S6). These results indicated that P deficiency increased expressions of genes related to DNA damage and RNA biosynthesis, thus the mechanism protecting the nucleotide from damage and maintaining its normal functioning might be initiated.

Phosphorus participates in the co-ordinated regulation of photosynthesis in cyanobacteria (Marcus and Gurevitz, 2000). In our study, expressions of ribulose-bisphosphate carboxylase (RuBisco) large subunit genes increased in P-deficient cells and decreased in P-resupplied-4 h cells (Supplementary Table S6). RuBisco catalyzes the first step of the Calvin cycle of photosynthesis and the oxidation of ribulose bisphosphate in the first step of photorespiration (Miziorko and Lorimer, 1983; Zhang Y.J. et al., 2015). RuBisco usually consists of two types of protein subunit, the large and the small subunits in plants, algae, and cyanobacteria, and the substrate binding sites located in the large subunit. In cyanobacteria, P binds to the RuBisco active site and to another site on the large subunit where it can influence transitions between active and less active conformations of the enzyme (Marcus and Gurevitz, 2000). Transcripts of LHC I, LHC II, P680, and P700 were up-regulated significantly in P-deficient cells, and down-regulated in P-resupplied cells (Supplementary Table S6). The LHC plays an important role in absorbing light and transferring energy to the center of the photosystem (Hiller et al., 1993, 1995; Boldt et al., 2012; Zhang Y.J. et al., 2015). P700 and P680 are the central pigment proteins of photosynthetic system I and II of eukaryotic cells. Furthermore, the expression of genes encoding key photosynthetic proteins, such as photosystem II CP43, photosystem II cytochrome c550 and cytochrome b559, cytochrome b6-f complex iron-sulfur subunit and F-type H⁺-transporting ATPase, were all up-regulated significantly in P-deficient cells and down-regulated in both P-resupplied-4 h and P-resupplied-28 h cells (Supplementary Table S6). These results indicated that P-deficiency enhanced light harvesting and photosynthesis which might be an adaptive mechanism of S. costatum to ambient P deficiency. However, diatom photosynthesis is suppressed by ambient nitrate (Bender et al., 2014) or iron limitation (Allen et al., 2008), indicating that photosynthetic responses of diatoms to different nutrients are complicated.

In the glycolysis pathway, triosephosphate isomerase is essential for efficient energy production, PGK is a major enzyme in the first ATP-generating step, and phosphofructokinase-1 is also one of the most important regulatory enzymes in glycolysis (van der Kamp, 2013). In our study, genes encoding these enzymes were up-regulated significantly in P-deficient cells. Moreover, genes encoding other important enzymes, such as fructose-1,6-bisphosphatase I, fructose-bisphosphate aldolase, pyruvate kinase, pyruvate dehydrogenase, and aldose 1-epimerase were also up-regulated in P-deficient cells (Supplementary Table S6). However, several genes related to glycolysis are down-regulated in nitrate-limited T. pseudonana, Fragilariopsis Cylindrus, and Pseudo-nitzschia multiseries (Mock et al., 2008; Bender et al., 2014). These results indicate that responses of diatoms to ambient N or P deficiency are different, which might be caused by their different storage capacity for N and P. In general, glycolysis was significantly enhanced in P-deficient cells in order to produce more energy for cells to acclimate to ambient P deficiency.

Many genes related to cell cycle regulation were identified in the transcriptome of S. costatum (Supplementary Table S6). Among them, genes encoding cyclin B and CDK were significantly up-regulated in P-deficient cells and down-regulated in P-resupplied cells (Supplementary Table S6). Cyclin B is a regulatory protein involved in mitosis and contributes to the switch-like all or none behavior of the cell in deciding to commit to mitosis (Zhuang et al., 2013). CDK-cyclin complexes are present in all eukaryotic lineages and play important roles in regulating the cell cycle and ensuring cell division (Robbens et al., 2005; Wang et al., 2013; Zhuang et al., 2013). Overexpression of both genes in P-deficient S. costatum suggested that ambient P deficiency disturbed the normal functions of cyclin B and CDK, and subsequently cell division ceased (Figure 1A).