Comparative Transcriptome Profiling Reveals Insights Into the Mechanisms Related to Explosive Growth of Alexandrium pacificum

Alexandrium pacificum is an organism that has an important impact on the aquaculture industry and human health. In this study, the digital gene expression approach was used to conduct a comparative analysis of differentially expressed genes (DEGs) that influence the explosive growth of A. pacificum following five treatment conditions: normal culture (C), high phosphorus and manganese (M), high irradiance (G), low phosphorus (P), and low nitrogen (N). Compared with the C conditions, a total of 265, 320, 185, and 150 DEGs were detected in the M, G, P, and N treatment groups, respectively. Clustering analysis suggested that A. pacificum acclimated to explosive growth using similar mechanisms in the M and G conditions. Analysis of DEGs showed that upregulation of genes associated with the pentose phosphate pathway and photosynthesis may contribute to explosive growth. Unigenes involved in the cell cycle were also found to be upregulated to promote cell division. The DEGs identified in this study may allow for the elucidation of molecular mechanisms responsible for the explosive growth of A. pacificum.


INTRODUCTION
Dinoflagellates, a large group of unicellular algae, are very important marine primary producers and grazers. Importantly, most harmful algae blooms (HAB) are caused by dinoflagellates, such as Karenia spp. and Alexandrium spp. (Zhang et al., 2013). In recent decades, the intensity and impact of HABs have increased drastically resulting in a negative impact on marine resources, water quality and human health (Genovesi et al., 2009). Dinoflagellates are also known for their unique features among eukaryotes. They contain vast amounts of DNA (from 3-250 pg · cell −1 ), and have an extremely large genome (up to 80 times of the human genome). Further, their chromosomes are permanently condensed during the cell cycle (Lin, 2006). To date, it remains unclear whether these fascinating biological characteristics are associated with HABs and explosive growth.
Alexandrium pacificum, one of the most well studied and widely distributed dinoflagellates, is known to form toxin blooms and cause paralytic shellfish poisoning (PSP). Studies have shown that not only shellfish, but also humans might be affected by PSP through the consumption of shellfish that have been contaminated by toxins (Hallegraeff et al., 2004). Thus, the development of a treatment strategy for HABs is critical.
Blooms of harmful dinoflagellates are affected by multiple environmental factors. Anthropogenic enrichment is considered to be one of the most pervasive changes altering the environment of offshore areas, especially the input of macronutrients and micronutrients, such as nitrogen and phosphorus (Li et al., 2014). Nitrogenous nutrients are thought to be important in triggering flagellate blooms (Glibert and Terlizzi, 1999). Studies have shown that the main cause of the increasing frequency of HABs is eutrophication. From the 1980s to 2000s, the annual frequency of red tide events in China had increased rapidly and the concentration of dissolved inorganic nitrogen (DIN) was increased from 2 to 5 µM (Wang et al., 2018). With respect to phosphates, the growth rate of Alexandrium tamarense increased with an increase in starting phosphate concentration (ranging from 0 to 10 µM) (Wang and Hsieh, 2005). Here, under laboratory conditions, we cultured algae using F/2 medium containing nitrogen and phosphate at concentrations of 882 and 36 µM, respectively, which are considerably higher than field conditions. Interestingly, some studies have shown that the biomass of A. tamarense observed under high concentrations of nitrogen (2,646 µM) was higher than that observed under low nitrogen concentrations when F/2 medium was used to culture the algae (Zhang et al., 2006). Likewise, when different concentrations of phosphate (3.6, 36, and 108 µM) were used to culture A. tamarense, the maximum biomass and growth rate occurred at the highest phosphate concentration (108 µM), while the lowest phosphate concentration (3.6 µM) limited growth (Yanjun et al., 2003).
Previous studies have shown that manganese could affect the photosynthetic activity and growth of dinoflagellates (Wang et al., 2019). For example, the growth rate of Alexandrium sp. was found to increase with an increasing manganese concentration (ranging from 0.9 to 2.7 µM), while in Coolia tropicalis, manganese was required to maintain the photosynthetic activity and growth rate (Paul and Hauck, 2006). Irradiance is another environmental factor that plays a vital role in stimulating growth. Laabir (2011) reported that the growth and cell yield of Alexandrium catenella increased significantly as the range of irradiance increased from 10 to 90 µmol photons m −2 s −1 . Thus, multiple environmental factors impact dinoflagellate growth and may influence HABs.
Many studies have focused on the influence of environmental factors on changes in population structure and size (Toulza et al., 2010). However, the molecular mechanisms underlying these events have been poorly studied and remain unclear. Uribe et al. (2008) developed an expressed sequence tag (EST) library and identified 6,496 unigenes from the EST database. Interestingly, genes associated with photosynthesis were highly expressed during the exponential growth phase. Similarly, an EST database developed by Toulza et al. (2010) found that the largest group of genes identified by this database were also related to photosynthesis. The large size of the dinoflagellate genome makes it difficult to sequence the whole genome. Thus, transcriptome analysis could be an efficient way to identify the underlying mechanisms of HAB development. It is required to understand the functional elements of the genome and identify the molecular constituents of cells and tissues (Laabir, 2011).
In this study, digital gene expression profiling (DGE) was used to comparatively analyze the gene expression and regulation of A. pacificum under the stage of explosive growth through high-throughput transcriptomic sequencing. Based on the environmental conditions under which HAB occur, five treatments were simulated including normal culture (C), high concentrations of phosphorus and manganese (M), high irradiance (G), low phosphorus (P), and low nitrogen (N). The aim of the study was to analyze the gene expression and regulation patterns that affect physiological and biochemical processes when A. pacificum are exposed to different environmental conditions. This database may serve as a valuable source to identify the key genes and pathways, which are activated under the stage of explosive growth of A. pacificum, as well as examine the mechanisms of HAB caused by A. pacificum.

Sample Preparation and Collection
Alexandrium pacificum was obtained from the Key Laboratory of Marine Genetics and Breeding, Ministry of Education of China, College of Marine Life Sciences. Stock cultures of A. pacificum were grown at 20 ± 1 • C at a photon flux density of 30 µmol m −2 s −l on a light:dark photocycle of 12:12 h and cultures were maintained in F/2 medium (without silicate).
For constructing of a reference transcriptome dataset of A. pacificum, four kinds of cultivating conditions were used including lag phase (Natural seawater with F/2 medium and the cells were collected at lag phase), log phase (Natural seawater with F/2 medium and the cells were collected at log phase) and induced log phase which contains two cultivating conditions (NaH 2 PO 4 · 2H 2 O and MnCl 2 · 4H 2 O were supplied to the final concentration of 0.144 mmol/L and 2.730 µmol/L, respectively, in the f/2 medium. The cell was collected at 7 and 12 days, then mixed for transcriptome analysis). Before treating, the algae which has already grown into logarithmic phase was placed in the dark for 48h to synchronize the cell cycle. The synchronized cells then were inoculated into three cultivating conditions at the initial concentration of 2 × 10 6 cells/L .
For preparing the growth curve, we cultured the algae for 21 days and count the algae number every 2 days using Sedgwick Counting Chamber to observe the growth trend (Figure 1) (Liu et al., 2019). For DGE analysis, five treatments were set up: M, G, C, P, and N ( Table 1). In sample M, the concentration of phosphate and manganese was 3 times than normal. In sample P and N, the concentration of nitrogen and phosphate was 10 times lower than normal. Three biological replicates were used in each treatment. Before treatment, the algae which has already grown into logarithmic phase were placed in the dark for 48 h to be synchronized of the growth. The synchronized cells then were inoculated into five treatment conditions at the initial concentration of 2 × 10 6 cells/L. All the 15 samples (5 treatments and each treatment including 3 biological replicates) were cultured for 12 days when they grow into logarithmic phase and were collected by centrifugation for 6 min at 4,000 g.  (Liu et al., 2019). Normal culture (C), high concentrations of phosphorus and manganese (M), high irradiance (G), low phosphorus (P), and low nitrogen (N) were the conditions used in the study.
To reduce bacterial contamination during collection of cells, every 100 mL of the cultures was filtered through a 10 µm pore size bolting-silk, rinsed with 300 mL sterile F/2 medium and subjected to the following treatments: the washed cells were suspended in 50 mL sterile F/2 medium containing 0.05% Tween-80 and 0.01 M EDTA (at 20 • C for 30 min), followed by ultrasonication for 1 min (50 W; ultrasonic treatment time of 5 s with an interval of 5 s). Then lysozyme (0.5 mg/mL, 20 • C for 10 min) and SDS (0.25%, 20 • C for 10 min) were added. The cells were washed three times with sterile F/2 medium to remove these reagents. The final cultures were stained with acridine orange (0.01%) for 1 or 2 min, then observed using epifluorescence microscopy (Nikon ECLIPSE 80i, Japan), only the axenic samples were harvested and used for the following experiments. Each sample contains 5 × 10 6 cells which were collected by centrifugation for 6 min at 4 000 g and then the cell pellets in 1.5 mL eppendorf tubes were frozen immediately and stored in liquid nitrogen for RNA extraction.

RNA Extraction and Library Construction for Sequencing
Total RNA was extracted using Trizol reagent (TaKaRa, Japan) following the manufacturer's instructions, followed by DNA digestion with DNase I (TaKaRa, Japan). The quantity and quality of the RNA samples were measured by gel electrophoresis and NanoDrop 2000 (Thermo, United States). High-quality RNA of the same quantities was used in subsequent experiments.
Fifteen libraries were constructed following the manufacturer's instructions of the Illumina Gene Expression Sample Prep Kit. The mRNA-seq libraries were sequenced 100 bp paired-end reads on the Illumina Genome Analyzer for DGE analysis. Base calling was performed by the Illumina instrument software.

Quality Control
The original data from HiSeqTM 2000 high-throughput sequencing is converted into the original sequenced sequence (Raw Data) by base calling analysis. FastQC was used for quality control. Raw data (raw reads) were cleaned by removing reads containing adapter, reads in which more than 10% of the bases were unknown, and low quality reads (where more than 50% of bases in a read had a quality value Q ≤ 5). The Q20, Q30, and GC content and sequence duplication level of the clean data were calculated. All the downstream analyses were based on clean, high-quality data.

Mapping of Reads
Transcriptome data of A. pacificum (SRX368254) were selected as the reference sequence to map of the date. Bowtie was used to map single-end clean reads to the unigene sequences. Clean reads were obtained by removing reads containing adapter, reads containing ploy-N and low quality reads from raw data. At the same time, Q20, Q30, and GC content of the clean data were calculated. All the downstream analyses were based on the clean data with high quality.

Quantification of Gene Expression Level and Differential Expression Analysis
Transcript abundance and differential gene expression were calculated with the program Cufflinks 1 . Gene expression levels were normalized with fragments per kilobase of exon per million mapped reads (FPKM) values to obtain the relative levels of expression. The false discovery rate (FDR) was used to determine the threshold of P-value in multiple test and analysis. Genes were considered to be differentially expressed only when their absolute value of log2 fold change was >1 and FDR was <0.05 between two treatments. The correlation of the detected number of counts 1 http://cufflinks.cbcb.umd.edu/ Variance-stabilized data obtained using DESeq was used to generate the heatmaps of differentially expressed genes. Clustering analysis was performed using the command hclust in R language which used K-means and SOM methods to calculate distance between clusters. The heatmap was drown using the R packages of ggplot2 and pheatmap. The heatmap and cluster analysis were performed based on the expression level of differential gene RPKM under different experimental conditions include all the 15 samples among the five treatments. Based on researches about RNA-seq analysis (Robinson and Oshlack, 2010) and the Pearson's correlation tests, we use the average value of RPKM of each treatment to complete the heatmap. Gene ontology (GO) and KEGG enrichment of differentially expressed genes were analyzed after the quantification of gene expression level. GO and KEGG enrichment was based on GOseq and KEGG database. The enrichment analysis was performed using the OmicShare tools, a free online platform for data analysis.

Quantitative Real-Time PCR Validation
Total RNA was extracted as described for DGE library preparation and sequencing. For the first-strand cDNA synthesis experiment, a Transcriptor First Stand cDNA Synthesis Kit (Roche) was used following the manufacturer's instructions. Quantitative real-time PCR was conducted using SYBR Green dye (LightCycle R 480 SYBR Green I Master). The selected genes including transketolase, phosphoglycerate kinase (PGK), phosphoribulokinase (PRK) and serine-threonine kinase (STK) were verified using the LightCycle R 480 Real-Time PCR System with the following cycling conditions: 95 • C for 5 min, followed by 45 cycles of 95 • C for 10 s, 60 • C for 10 s and 72 • C for 20 s. Specificity of the qPCR product was analyzed by melting curve analysis. The sequences of the primers used are given in Supplementary Materials. Actin and GAPDH served as internal controls. The 2 − Ct method was used to calculate relative gene expression values. The mean of fold changes was calculated for each experiment and mean >2 or <0.5 was considered as significant.

Sequencing of 15 Digital Gene Expression Libraries
All fifteen RNA samples of A. pacificum under the five treatments (each treatment contains three biological replicates) were subjected to RNA-seq. Each sample yields 11-12.6 million raw reads. After quality assessment and filtering low quality data, approximately 10.7-12.1 million clean reads were obtained in each library ( Table 2). Approximately 69.14 ∼ 80.57% reads were mapped to the transcriptome of A. pacificum (SRX368254), which also proved that the transcriptome data was a reliable reference. The mean Q20 and Q30% of five treatments were ranging from 93.92 to 95.1% and from 83.98 to 86.37% average, respectively, suggesting that the high quality sequencing data made the results more reliable. The mean GC content (%) is 63.5% average. The reproducibility of the DGE was tested by Pearson correlation analysis for every three replicates (Figure 2). Sample C was representative, other results was provided in Supplementary Materials. The rooted squares of the Pearson's correlation coefficient (R 2 ) were greater than 0.92, mostly ranging from 0.92 to 0.953 ( Table 3). The rooted squares of the Pearson's correlation coefficient (R 2 ) were close to 1, indicating both the better parallelism between samples and the reliability of the experimental results.

Number of Differential Expression Unigenes
The expression level of mapped unigenes was normalized with a value of FPKM (fragments per kilobase of exon per million fragments mapped). An absolute value of log 2 fold change > 1 and false discovery rate < 0.05 was set to declare the differential expressed genes (DEG).
The fold change of the expression levels in the different treatments were compared in pairs ( Table 4). Eutrophication conditions (M) and high irradiance conditions (G) compared with C and nutrient limitation conditions (N and P) to find out the DEGs functioned in promoting the growth of the algae. The treatment N compared with P to find out the different mechanisms when the algae cope with the nutrient limitation conditions. The treatment M compared with G to figure out the different mechanisms in promoting the growth of the algae.
A hierarchical cluster was used to determine the profiles of the differentially expressed unigenes among the five treatments (Figure 3). Unigenes with the same or similar expression patterns were gathered in the same subcluster. For example, the unigenes from sub-cluster 1 show similar expression trend. The pattern under the M and G conditions were similar which indicates the similar mechanisms under the stage of explosive growth. In subcluster 1, the unigenes were found to relate to the carbohydrate metabolism and energy metabolism which enriched in the pentose phosphate pathway and carbon fixation in photosynthetic organisms. In subcluster 2, the unigenes were associated with the translation of genetic information including genes in the pathway of ribosome biogenesis. In subcluster 3, most of the unigenes were related to the amino acid metabolism, except that there were less unigenes related to the synthesis of carbohydrates (Figure 4).

Gene Ontology Functional Analysis of Differentially Expressed Genes
Gene ontology (GO) was carried out to unify the representation of gene and gene product attributed in A. pacificum under different treatment conditions. For the biological process category, the mostly highly represented terms were cellular process and metabolic process. In metabolic process, the number of DEGs under M and G were much more than under N and P indicating the metabolism of the cells were increased to deal with the rapid growth of the algae. For the cellular component category, the mostly highly represented terms were cell and cell parts. For the molecular function category, the two mostly highly represented terms were binding and catalytic activity. And the number of DEGs under M and G were much more compared with those under N and P in these two terms suggesting the cells need more carriers and enzymes to complete the biological process under M and G and the inhibition of parts of the biological process when exposed to nutrient deficiency conditions ( Figure 5). Of interest, 75 of the 172 upregulated unigenes were related to the structure of the ribosome. These findings were consistent with the KEGG analysis data. In higher plants, the PPP is associated with plant development and environmental stress responses. Here, DGE analysis revealed that 12 unigenes involved in the PPP were upregulated under M conditions including 6-phosphogluconate dehydrogenase (6PGDH), glucose-6-phosphate dehydrogenase (G6PDH), glyceraldehyde-3-phosphate dehydrogenase (GAPDH), transketolase and transaldolase.
Upregulation of unigenes involved in carbon assimilation such as ribulose-3-phosphate isomerase, phosphoglycerate kinase and fructose-1, 6-bisphosphate aldolase was also observed under M conditions. Upregulation of such genes would be associated with an increased rate of photosynthesis.
Finally, M treatment conditions led to an upregulation of genes that can affect the development of organisms, including PCNA, calmodulin (CaM), glycine-rich protein and ubiquitin. Nine transcription factors were also found to be differentially regulated included basic leucine zipper (bZIP), basic transcription factor 3 (BTF3), and myeloblastosis oncogene homolog (MYB).
Taken together, these findings indicate that eutrophic conditions promote expression of genes associated with nucleotide synthesis, cell growth and photosynthesis.

Differentially Expressed Genes in Response to High Irradiance Conditions
Under G conditions, the growth rate of A. pacificum was significantly increased and the number of algae in the log phase was significantly higher compared with normal conditions (Table 4). A total of 320 DEGs were identified under G conditions relative to C conditions ( Table 4) including 228 upregulated unigenes and 92 downregulated unigenes. Based on GO enrichment analysis, there were 197 and 57 genes associated with the enriched GO terms in the up-and downregulated unigenes, respectively. KEGG analysis revealed that the upregulated unigenes were highly enriched in ribosome biogenesis and carbon fixation in photosynthetic organisms.
Digital gene expression analysis revealed that unigenes associated with electron transfer reactions in photosynthesis were upregulated including PsbO, PsbU, cytochrome b6f complex, ferredoxin and thioredoxin. Among the 320 DEGs, unigenes related to carbon fixation including ribulose-3phosphate isomerase, phosphoglycerate kinase, ribulose-1,5biphosphate carboxylase and fructose-bisphosphate aldolase were upregulated under G conditions. Similar to M conditions, unigenes encoding PCNA, CaM and ubiquitin, which are related to the development of organisms, were upregulated. In addition, three unigenes encoding transcription factors were identified including MYB. Thus, high irradiance promotes expression of genes associated with growth and photosynthesis.
A total of 150 DEGs were identified under N conditions compared compared with C conditions ( Table 4). After treatment with low nitrogen, 46 unigenes were upregulated and 104 unigenes were downregulated. KEGG analysis revealed that the DEGs were highly enriched in nitrogen metabolism. Based on GO enrichment analysis, 21 upregulated unigenes and 84 downregulated unigenes were enriched in GO terms. The downregulation of six unigenes involved in nitrogen metabolism including nitrate reductase, which plays a vital role in nitrogen assimilation, was observed. In addition, that amino acid metabolism is disrupted. Thus, in contrast to other conditions, nutrient limitations promote expression of genes associated with metabolism and stress.

Validation of Gene Expression
To validate the expression profiles obtained by RNA-Seq, RT-PCR was performed on four genes. Normalized gene expression according to 2 − CT method was shown in Figure 6. Expression of the four genes under M and G were all higher than those under N and P (p < 0.05). The trend in RT-PCR expression was in agreement with the DGE data (Figure 7).

The Role of the Pentose Phosphate Pathway Under the Stage of Explosive Growth
Cell growth requires biosynthesis of critical intermediates such as nucleotides and amino acids. During cellular proliferation, a key feature of the metabolic transformation is the enhancement of biosynthetic capacity (Levine and Puzio-Kuter, 2010). The PPP, a fundamental component of cellular metabolism, can be subdivided into oxidative and non-oxidative PPPs. In this study, 35 unigenes associated with the PPP were found in the transcriptome data (SRX368254). Furthermore, a total of 14 DEGs involved in the PPP were identified under different treatments including G6PDH, 6PGDH, transketolase, transaldolase and ribose 5-phosphate isomerase ( Table 5). The upregulation of transketolase and transaldolase under M conditions was consistent with the comparative transcriptome data of A. pacificum (SRX368254). The oxidative branch converts glucose 6-phosphate into ribulose 5-phosphate through the consecutive reactions of G6PDH and 6GPDH (Miclet et al., 2001). In the non-oxidative branch, transketolase and transaldolase are the key enzymes of the non-oxidative PPP, which connect glycolysis and the PPP. They catalyze the conversion of ribose 5-phosphate and sedoheptulose 7-phosphate to complete the PPP (Kochetov and Sevostyanova, 2005;Samland and Sprenger, 2009).
Several studies have demonstrated that the expression levels of G6PDH, transketolase and transaldolase were correlated with the progression of cancer, suggesting that enhanced PPP activity is associated with cell proliferation (Deberardinis et al., 2008). In The "x-" axis represents the treatments. The "y-" axis represents the value of the relative expression level [log2 (ratio)]. In the subcluster, the gray line represents the value of the relative expression levels, the black line represents the mean value of the relative expression level. The subclusters in figure represent the subclusters in Figure 3. Normal culture (C), high concentrations of phosphorus and manganese (M), high irradiance (G), low phosphorus (P), and low nitrogen (N) were the conditions used in the study.
our study, under M conditions, the unigenes encoding G6PDH, 6PGDH, transketolase, transaldolase, ribose 5-phosphate isomerase and pyruvate kinase were upregulated ( Table 5), suggesting increased PPP activity, which may enhance the cell biosynthetic capacity to adapt to rapid proliferation. In addition, transketolase and transaldolase were also found to be upregulated in the induced log phase in the A. pacificum comparative transcriptome dataset (SRX368254).
Under nutrient limiting conditions (P, N), transketolase, transaldolase and ribose 5-phosphate isomerase were downregulated, which reduced the conversion of ribose 5phosphate and resulted in the inhibition of nucleic acid synthesis and potentially growth arrest. Interestingly, normal expression of the unigenes encoding G6PDH and 6GPDH was observed under N and P conditions compared with C conditions ( Table 5). G6PDH and 6GPDH, key enzymes in the oxidative PPP, yield NADPH, which maintains the redox balance under stress conditions (Slekar et al., 1996). Hence, under the nutrient limiting-induced stress, algae still need to express G6PDH and 6GPDH to yield NADPH in order to maintain the redox balance for survival.

Photosynthesis Under the Stage of Explosive Growth
Photosystems are functional and structural units of protein complexes involved in photosynthesis that carry out the primary photochemistry of photosynthesis: the absorption of light and the transfer of energy and electrons. Of the two photosystems [photosystem I (PS I) and PS II], the cytochrome b 6 f complex mediates electron transfer by oxidizing lipophilic plastoquinone and reducing plastocyanin (Kurisu et al., 2003). In this study, unigenes encoding PsbO, PsbU and the cytochrome b6f complex were found to be upregulated after G treatment ( Table 5). The PsbO protein is an extrinsic subunit of PS II and is thought to play a central role in stabilization of the catalytic manganese cluster   in which water is split into four protons and O 2 via a series of four redox steps (Murakami et al., 2002). PsbU, a 12 kDa protein of PS II, is thought to regulate PS II function in cyanobacteria (Shen et al., 1998).
The ferredoxin/thioredoxin system (FTS), the major redox system responsible for photosynthesis in plant chloroplasts, is composed of three proteins, thioredoxin (Trx), ferredoxin thioredoxin reductase (FTR) and ferredoxin (Fd). In the FTS, Fd transfers electrons received from PS I to FTR, which then reduces Trx, providing a biochemical link between light reactions and the regulation of metabolism (Holmgren, 1995;Schürmann, 2003). In chloroplasts, most thioredoxin-regulated enzymes are associated with carbon metabolism including phosphoribulokinase (PRK), fructose-1,6-bisphosphatase and G6PDH (Ruelland and Miginiac-Maslow, 1999).
Among the DEGs, the unigenes encoding Trx and Fd were upregulated under G conditions ( Table 6), suggesting that the FTS was active. Furthermore, together with PsbO, PsbU and the cytochrome b6f complex, Trx and Fd promote the transport of photosynthetic electrons and provide ATP and NADPH for carbon assimilation and the PPP.
For algae, many studies have suggested that the Calvin-Benson cycle (C3 pathway) was predominant compared with the C4 and crassulacean acid metabolism (CAM) pathways (Tsuji et al., 2009). However, a recent study suggested that C4 pathways may also be used in algae. Consistent with this, 72 and 52 unigenes associated with the C3 and C4 pathways, respectively, were found in the A. pacificum transcriptome dataset (SRX368254). Among these unigenes, 16 and 9 unigenes associated with the C3 and C4 pathways, respectively, were differentially expressed under different treatment conditions ( Table 6). In the C3 pathway, three DEGs encoding components of the ribulose 1, 5-bisphosphate carboxylase (Rubisco) complex, which is the central carboxylation enzyme for CO 2 fixation, Negative numbers indicated down-regulated and positive number mean up-regulated. The dashes (-) mean that the unigenes was not significantly differentially expressed under these conditions. Normal culture (C), high concentrations of phosphorus and manganese (M), high irradiance (G), low phosphorus (P), and low nitrogen (N) were the conditions used in the study.
were upregulated under G and P treatment conditions. The C3 pathway of carbon fixation is catalyzed by Rubisco, producing two three-carbon molecules of 3-phosphoglyceric acid, through the carboxylation of the five-carbon ribulose-1,5-biphosphate (Wang et al., 2011). In addition, G treatment resulted in the upregulation of DEGs associated with the C3 pathway including fructose-bisphosphate aldolase, phosphoglycerate kinase, ribulose-3-phosphate isomerase and PRK. All of the five afore-mentioned genes were involved in the rate-limiting step of the C3 pathway, suggesting that the C3 pathway may be involved in the carbon assimilation process under these conditions. Unigenes invon the A. pacificum transcriptome dataset (SRX368254). Among these unigenes, the differential expression of enzymes that catalyze the C4lved in the C4 pathway were also found i pathway was observed including phosphoenolpyruvate carboxykinase (PEPCK), pyruvate orthophosphate kinase (PPDK) and pyruvate kinase (PK) ( Table 6). In the C4 pathway, PEPCK is a key enzyme that catalyzes the release of CO 2 from oxaloacetate to produce phosphoenolpyruvate (PEP) to complete the C4 pathway (Hibberd and Covshoff, 2010). Since PPDK catalyzes the conversion of pyruvate to PEP (Parsley and Hibberd, 2006), upregulated PPDK expression under G treatment conditions could promote the C4 pathway. The function of PK is similar to that of PPDK. However, PK catalyzes an irreversible reaction, and does not consume ATP (Baud et al., 2007). Under high temperature, high light, and the current CO 2 concentration in the atmosphere, the C4 pathway is more efficient than C3 photosynthesis because it increases the CO 2 concentration around the major CO 2 fixating enzyme Rubisco (Bräutigam and Gowik, 2016). When HAB occurs, a rapid increase in the size of the population requires efficient photosynthesis. Hence, the C4 pathway may exist in A. pacificum to promote carbon assimilation and support increased growth. Negative numbers indicated down-regulated and positive number mean up-regulated. The dashes (-) mean that the unigenes was not significantly differentially expressed under these conditions. Normal culture (C), high concentrations of phosphorus and manganese (M), high irradiance (G), low phosphorus (P), and low nitrogen (N) were the conditions used in the study.
There is a significant amount of transcriptomic information about C4-related enzyme variations among various algae (Zhao and Su, 2014). In Bacillariophyta, Phaeodactylum tricornutum and Thalassiosira pseudonana have developed a C4-like photosynthesis pathway in addition to the carbon-concentrating mechanism (Armbrust et al., 2004;Kroth et al., 2008), while in Chlorophyta, Ulva linza and Ulva prolifera have been shown to express the C4 pathway . Taken together, these findings suggest that the C4 pathway may be more complex in A. pacificum. Furthermore, the transcriptome data provides the foundation to elucidate the mechanisms underlying the complete carbon fixation system.

Cell Division
Among the DEGs, 92 unigenes were enriched in the ribosome biogenesis pathway. Ribosome biogenesis underlies a cell's capacity to grow as cell growth requires a large numbers of ribosomes, the molecular factories that carry out protein synthesis (Lempiäinen and Shore, 2009;Donati et al., 2012). Based on the A. pacificum transcriptome data (SRX368254), 16 unigenes related to DNA replication including PCNA were upregulated, especially in the induced log phase. PCNA was also upregulated under M and G treatment conditions (Table 6), reflecting the proliferative state of A. pacificum.
Proline-rich and glycine-rich proteins are required for the synthesis of cell walls, which are required before division and are involved in many important regulatory processes in plants (Josè and Puigdomènech, 1993;Bocca et al., 2005). Among the DEGs, unigenes encoding proline-rich and glycine-rich proteins were found to be upregulated under G and M treatment conditions. However, under the P and N treatment conditions, these unigenes were not differentially expressed due to their involvement in the environmental stress response (Li et al., 2012;Kim et al., 2005).
Calmodulin is the most important Ca 2+ receptor in plants and participates in transcriptional regulation by acting on transcription factors, especially during responses to environmental changes. In this study, CaM and CaM-related proteins such as serine-threonine kinase were upregulated under M and G conditions ( Table 7). CaM and the CaM-dependent signaling system is essential for the proliferative capacity of cells (Chafouleas et al., 1982). In mammalian cells, overexpression of CaM revealed a linear relationship between CaM concentration and rate of G1 progression (Joseph and Means, 2000). In contrast, reduction of CaM levels led to the majority of cells being arrested in G2, suggesting that CaM has a critical role in mediating entry into and progression through mitosis (Keith et al., 1988). In addition, three microRNAs in A. pacificum including miR169a, miR169c, and miR169f, which target CaM-dependent kinase II (CaMKII), were found to be upregulated in the log phase (Geng et al., 2015). In A. pacificum, CaM was shown to activate CaMKII. In animals, CaMKII can promote cell proliferation and growth by facilitating replication of the centriole. At the same time, the CaMKII inhibitor can arrest cells in the G2 phase of the cell cycle. Hence, CaMKII may play a major role in regulating the cell cycle and rapid growth during the log phase in A. pacificum.
The high mobility group box (HMGB) proteins are an important component of chromatin and share some functions with histone H1 (Johns, 1982). In A. pacificum, histone-like proteins such as H2A.X are rarely expressed (Hackett et al., 2004). Based on the transcriptome dataset (SRX368254), histones including all four core nucleosomal histones, histone H1, linker histone H1 and H5 family proteins, were found to be expressed at low levels . However, HMGB binding affinity is much lower than H1, suggesting that HMGB may interact with the linker region only when histone H1 is absent or weakly expressed (Pallier et al., 2003). In this study, the DEGs encoding HMGB2 were upregulated under M and G conditions, suggesting that HMGB2 may play a critical role in DNA packaging. HMGB2 has been shown to be a component of mitotic chromosomes and is associated with condensed chromosomes in HeLa cells during mitosis (Scaffidi et al., 2002). Thus, HMGB2 may be associated with the condensed chromosomes of A. pacificum during the cell cycle (except DNA replication).

Nitrogen Metabolism Under the Stage of Explosive Growth
Nitrogen metabolism can be regulated by changes in the abundance of mRNA of several components involved in nitrogen uptake and nitrogen assimilation. Nitrate reductase (NR), which can catalyze NO 3− reduction to NO 2− is critical in nitrogen metabolism. NR activity could be induced by NO 3 − (Lillo, 2008). As for some research, when the algae experiences N stress, the cell upregulates the enzymes to promote translocation of external nitrogen for supplementing (Li et al., 2020). In this study, the DEGs encoding NR were down regulated Negative numbers indicated down-regulated and positive number mean up-regulated. The dashes (-) mean that the unigenes was not significantly differentially expressed under these conditions. Normal culture (C), high concentrations of phosphorus and manganese (M), high irradiance (G), low phosphorus (P), and low nitrogen (N) were the conditions used in the study. Negative numbers indicated down-regulated and positive number mean up-regulated. The dashes (-) mean that the unigenes was not significantly differentially expressed under these conditions. Normal culture (C), high concentrations of phosphorus and manganese (M), high irradiance (G), low phosphorus (P), and low nitrogen (N) were the conditions used in the study. under nitrogen-limiting conditions (Table 8), mainly because the external nitrogen has been used up when the cell was collected at 12th day, further resulting in the inhibition of nitrogen assimilation. Photosynthesis plays a critical role in nitrogen assimilation by providing ATP and reducing equivalents. Furthermore, large amounts of inorganic nitrogen are used in photosynthesis to produce sugars and amino acids. In C3 plants, PEPC is essential for the provision of carbon skeletons used for amino acid biosynthesis (Oaks and Hirel, 2003). The downregulation of PEPC under nitrogen-limiting conditions could lead to the inhibition of carbon assimilation and affect the efficiency of photosynthesis, which in turn would result in insufficient energy and provision of carbon skeletons for nitrogen assimilation (Champigny and Foyer, 1992).

Transcription Factors
In this study, transcription factors including bZIP, MYB, BTF-3, and EF-Tu were found. The upregulation of BTF-3 and EF-Tu was detected under M and G conditions. BTF-3 is a 27 kDa protein that forms a stable complex with RNA polymerase IIB and is required for transcriptional initiation (Kanno et al., 1992). In higher plants, gene silencing of BTF-3 caused leaf yellowing and abnormal leaf morphology without altering the overall growth of the plant (Yang et al., 2007), suggesting that BTF3 plays a vital role in tissue development. EF-Tu plays a role in protein synthesis by promoting the GTP-dependent binding of aminoacyl-tRNA to the A site of the ribosome (Riis et al., 1990). The mitochondrial EF-Tu gene was found to be differentially expressed during flower development, with the highest EF-Tu expression occurring during relatively early stages of flower development (Lee et al., 2002).
The expression of bZIP and MYB were upregulated under N and P treatment conditions. bZIP and MYB are critical for regulating defense gene expression in plant stress tolerance. In Arabidopsis, most of the transcription factor genes induced under P starvation conditions belong to the MYB superfamily (Bender and Fink, 1998). In addition, in Lotus japonicas and soybean, MYB genes were found to be upregulated under nitrogen-limiting conditions, indicating that MYB has a role in regulating flavonoid biosynthesis in response to nitrate starvation (Rubio et al., 2001;Miyake et al., 2003).

CONCLUSION
The variety and number of unigenes as well as the differences in gene expression patterns under distinct treatment conditions suggest that complicated and diverse regulatory mechanisms are activated under the stage of explosive growth of A. pacificum. Several common metabolic pathways are induced in A. pacificum, such as the pentose phosphate pathway and photosynthesis under eutrophication conditions and the nitrogen metabolism pathway under nutrition deficient conditions. Based on the DGE analysis, PPP and photosynthesis were activated which enhance the biosynthetic capacity and provide intermediates for cell proliferation. In addition, unigenes related to cell division such as PCNA, CAM, Proline-rich and glycine-rich proteins were detected which can facilitate the cell cycle and promote the growth of A. pacificum (Figure 8). This sequencing dataset and analysis may serve as a valuable resource to study the mechanisms of explosive growth of A. pacificum and to further elucidate the process of HAB.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: NCBI (accession: PRJNA715858).

AUTHOR CONTRIBUTIONS
YL, ZS, and SZ contributed to conception and design of the study. ES organized the database. ZN performed the statistical analysis.
YL wrote the manuscript. ZZ and JQ contributed to the validation experiment. All authors contributed to manuscript revision, read, and approved the submitted version.

FUNDING
This work was supported by National Natural Science Foundation of China (Grant Nos. 41176098 and 41676091).