Transcriptomic Analysis Reveals the Regulation Mechanism of Sporogenesis in Ulva prolifera

Ulva prolifera, the main causative species of green tide, has drawn much attention. Spore reproduction was one of the main reproduction strategies that could be induced by changing environmental factors, but the mechanism of spore formation remains obscure. Here, we cultured U. prolifera (segments) under the optimal sporulation condition, and four states in total from vegetative to reproductive were defined. Additionally, the chlorophyll fluorescence and transcriptome analysis were determined at these four states. The results showed that: (1) Compared with state I, the photosynthetic capacity (Fv/Fm, Fv′/Fm′, rETRmax) and chlorophyll content (Chl a, Chl b, carotenoids) were enhanced in state II, whereas it decreased in state III and IV (the spore formation period); (2) a total of 41,058 unigenes were expressed during the spore formation process; (3) compared with state I, the genes related with photosynthesis, terpenoid backbone biosynthesis, and carotenoid biosynthesis were significantly upregulated in states II, III, and IV whereas glycolysis was downregulated in state I; (4) some genes of the transcription factors families, such as the C3H family, may be one of the key factors that regulate genes in the spore formation; (5) 574 of the differentially expressed genes (DEGs) associated with flagella biosynthesis were annotated according to Gene Ontology and Kyoto Encyclopedia of Genes and Genomes Orthology, including 20 unigenes encoding intraflagellar transport proteins (IFTs) which had not been reported in previous transcriptome analysis in U. prolifera. This study provides a new perspective of spore formation at the gene transcriptional level, although the detailed transcription regulatory spore formation network remains to be unveiled.


INTRODUCTION
Ulva spp., the most ubiquitous and important primary producers around the world's oceans (Raffaelli et al., 1998), can be used as animal feed in aquaculture or an important protein food source for humans because they are highly rich in proteins (He et al., 2021). They were also the domain species of the green tide, especially for the Ulva prolifera, which is the main causative species for their rapid growth rate and complex reproductive abilities (Cui et al., 2018;Wang et al., 2018a). The reproductive strategies of U. prolifera include parthenogenesis, sexual and asexual, and the asexual spores, sexual male and female gametes (Lin et al., 2008). Recently, many studies about the effects of ecological factors, such as temperature, light intensity on the growth, physiological and biochemical characteristics of U. prolifera have been reported (Cui et al., 2015;Wu et al., 2018a;Wang et al., 2020a). However, the life history of U. prolifera is still very controversial, and whether gametophytes or sporophytes dominate the life cycle has not been determined (Wang et al., 2020b).
For macroalgae, the spore reproduction could be affected by biotic and abiotic factors, such as the thallus fragmentation, temperature, and light intensity (Han et al., 2008;Gao et al., 2010Gao et al., , 2017. Previous studies have shown the following: about 0.5 mm diameter of U. prolifera fragments were almost entirely converted to sporangia (Gao et al., 2010), abrupt temperature changes could stimulate the reproduction of U. rigida (Gao et al., 2017), and the minimum light intensity required for the spore release was 5-10 µmol m −2 s −1 in U. pertusa (Han et al., 2008). The spore release of U. pseudocurvata was also spontaneously triggered as soon as the light cycle began (Smith, 1947). Moreover, to understand better the changes in the reproductive process of Ulva, some studies divided the reproductive progress, from the transformation of vegetative cells into reproductive cells, into two to four stages (Wichard and Oertel, 2010;Carl et al., 2014;Wang et al., 2016;Chavez-Sanchez et al., 2018). Although the previous study showed that spore formation was regulated by the redox status of the plastoquinone pool that could be affected by photosynthetic electron transport chain (PET) (Wang et al., 2016), the regulatory mechanism of spore formation of U. prolifera is not yet clear.
RNA-sequencing, as a revolutionary tool for transcriptomics, has been wildly used for biological inquiries (Wang et al., 2009;Du et al., 2021). The previous study of transcriptome analysis showed some of the potential molecular mechanisms of growth have been revealed, such as the specific photoprotective mechanism associated with LhcSR and PsbS . Expressed sequence tags analysis found that some sequences code for proteins involved in the formation of flagella (axonemaltype dynein light chain), the length of the flagella (LF4), and the movement of the flagella (CAM kinase and PF16) (Stanley et al., 2005). Flagella was the key for spores to swimm in the water to survive and proliferate under different environmental conditions (Silflow and Lefebvre, 2001). In U. prolifera and U. linza, the existence of C 4 -type metabolism pathway also has been confirmed by transcriptome analysis (Xu et al., , 2013. Compared with the other three Ulva species, U. prolifera demonstrated noticeably larger numbers of total transcription factors (TFs) which were prominent in many TF types, and several growth-related genes were enriched (Wang et al., 2019). However, few studies were related to spore formation.
So, in this study, the global mRNA expression patterns of U. prolifera at multiple sampling time points from vegetative cells to matured spores were examined by using RNA-seq technology to understand better the morphological conversion from vegetative cells to matured spores of U. prolifera.

Species and Culture Conditions
Ulva prolifera seeding (3 mm) was obtained from Xiangshan Xuwen Seaweed Development Co. Ltd., Xiangshan County, Zhejiang Province, China (29 • 05 065N, 121 • 56 153E), and it was identified through morphology and genotyped by sequencing the 18S and ITS rRNA genes. They were cultured at 20 • C, 100 µmol m −2 s −1 condition until ca. ≥18 cm length and were used in the following experiments. No specific permits were required for the using this species.

Morphology Observation
The process of spore formation and release was observed by using the optical microscope with digital camera (ToupCam TM , Touptek, China). According to Carl et al. (2014) and Chavez-Sanchez et al. (2018), four states of thalli segments, from vegetative to reproductive, were defined. Moreover, the synchronization of different thalli was defined as the number of more than 50% of different thalli that changed significantly. The thalli with a high synchronization than 90%, were cultured in the laboratory after being cut into small segments. So, when the state was confirmed, the samples were harvested for transcriptome analysis.

Chlorophyll Fluorescence Parameters Measurement
After defining the state, a PSI fluorometer (AquaPen-C, Photo System Instruments, Drasov, Czechia) was used to determine the effective quantum yield of PSII (F v /F m ), the non-photochemical quenching (NPQ), and the rapid light curves (RLC). After adapting to dark conditions for 15 min, NPQ and F v /F m were measured with the cultivated light intensity as the actinic light. RLC was fitted as rETR = I/(aI 2 + bI + c) (Eilers and Peeters, 1988), where I is the photon flux density of activity light (µmol m −2 s −1 ), and a, b, c are the adjustment parameters. The initial slop (α), the maximum electron transport rate (rETR max ), and the light saturation parameter (E k ) were expressed as a function of the parameters a, b, and c as follows: α = 1/c; rETR max = 1/(b + 2(a × c) 1/2 ); E k = c/(b + 2(a × c) 1/2 ).

Chlorophyll Content Determination
Chlorophyll a (Chl a), Chlorophyll b (Chl b), and carotenoid (Car) are important photosynthetic chlorophyll in green macroalgae. Approximately 0.01 g freshly weighed thalli were extracted in 3 mL methanol at 4 • C in darkness overnight, and then the absorption spectrum of the supernatant was obtained by scanning the sample from 250 to 750 nm with a scanning spectrophotometer (Yuanxi Instrument Co., Ltd, Shanghai, China). The chlorophyll concentration, including Chl a, Chl b. and Car, was calculated according to Wellburn (1994).

RNA Isolation and cDNA Library Construction
Total RNA of the samples was extracted with an E.Z.N.A. R Plant RNA Kit (Omega Bio-tek, Norcross, GA, United States) according to the manufacturer's instructions. The concentration and the quality of RNA were determined using a NanoDrop ND1000 spectrophotometer (Thermo Scientific, Waltham, MA, United States) and Agilent Bioanalyzer 2100 (Agilent, Santa Clara, CA, United States).
RNA samples and a pooled RNA from these 12 samples (three biological replicates, four states) were submitted for differential expression analysis and de novo transcriptome analysis, respectively. The 13 cDNA libraries were constructed using mRNA library preparation kit (MGIEasy TM , MGI, Shenzhen, China). The mRNA with ploy(A) structure was enriched with oligo magnetic adsorption. The enriched mRNAs were fragmented and then reversed into a double-strand cDNA (dscDNA) with an N6 random primer. Sequencing adaptors were linked to the end of the purified dscDNA, and the ligation product was amplified by PCR with specific primers. The PCR product was thermally denatured into a single strand, and then the single-stranded DNA was circularized with a bridge primer to obtain a single-stranded circular DNA library.

RNA-Sequence, de novo Assembly of Transcriptome, and Gene Identification
The library for de novo transcriptome was sequenced by 150 bp double ends sequencing method while the other 12 libraries for differential expression analysis were sequenced by 50bp single end sequencing method on a DNBSEQ platform, performed by Wuhan Genomic Institution (BGI, Shenzhen, China). Additionally, although the genome of U. prolifera was reported (the accession number in NCBI was SDUY00000000.1; He et al., 2021), the annotation of this genome has not been found in the NCBI database. So, the de novo transcriptome analyzed was chosen in this study. The details were similar to Hu et al. (2019).

Differentially Expressed Genes and Their Gene Ontology and Kyoto Encyclopedia of Genes and Genomes Pathway Enrichment Analyses
Differentially expressed genes (DEGs) were identified using the DESeq method based on the negative binomial distribution . There are two standards to select DEGs, namely, using the Q value < 0.001 and the absolute value of fold change ≥ 2 as two threshold values to screen out DEGs in different comparison groups. The Q value was a corrected p-value calculated by using KOBAS 2.0. method (Mao et al., 2005). The Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway and Gene Ontology (GO) enrichment analysis of DEGs were performed using the KEGG database and GO database (He et al., 2019a). A Q value cutoff of 0.05 was set for the identification of remarkably enriched GO and KEGG pathways for DEGs.

Expression Validation by Quantitative Real-Time PCR Analysis
Ten DEGs were randomly selected for quantitative realtime PCR (qRT-PCR) to validate the RNAseq results. The single-strand cDNAs were synthesized from 1 µg of the total RNA using HiScript II 1st strand cDNA synthesis Kit (Vazyme, Nanjing, China) according to the manufacturer's instructions. qRT-PCR was performed on a Mastercycler Realplex2 realtime PCR instrument (Eppendorf, Germany) using a ChamQTM Universal SYBR qPCR Master Mix Kit (Vazyme, China). The qRT-PCR program was as follows: 95 • C for 30 s, 40 cycles of 95 • C for 15 s, 60 • C for 15 s, and 72 • C for 20 s. The primers were designed using Primer Premier 5 software (Premier Biosofe International, San Francisco, CA, United States) and are listed in Supplementary  Table 1. Four technical replicates were performed for each sample. The relative expression levels of the genes were normalized by comparison with 18S rDNA expression and calculated using the 2 − Ct relative quantitative method (He et al., 2019a).

Statistical Analysis
Three replicates for each treatment were used for all the treatments. The data were shown as means ± SD. The data were analyzed using SPSS software (version 22.0 USA) and assessed by one-way analysis of variance (ANOVA). A turkey post hoc test (Turkey HSD) was performed to show the difference between different states, such as chlorophyll content and the gene expression levels obtained from qRT-PCR among treatments. The significant levels were set as p < 0.05.

Morphology Change of Ulva prolifera in the Spore Formation
According to previous studies (Pettett, 2009;Carl et al., 2014;Park, 2020) and based on the observations, the reproductive process of U. prolifera was divided into four states (Figure 1). In detail, first, initially, the chloroplasts were visible to the interior of the green vegetative cell (state I; Figure 1A). Second, the chloroplasts were prepared to displace to one side of the cell (state II; Figure 1B). The spores filled the interior of the cell (state III; Figure 1C) and were matured with a clear shape and an eyespot (state IV; Figure 1D). Previous studies showed that sporulation can be identified in Ulva thalli through a change in coloration (Pettett, 2009;Carl et al., 2014;Park, 2020). During this study, the color of the thalli segments converted to yellow from green, which was significantly observed in states III and IV, indicating the spore formatted in these two states. Moreover, an eyespot was also observed in the states III and IV, and especially in state IV, the eyespot has a clear shape. Therefore, we speculated that states III and IV were the crucial periods of spore formation.

Chlorophyll Fluorescence Parameters and Chlorophyll Content
With the change of algae cell morphology, photosynthesis changed significantly at the same time. Compared with state I, the relative electron transport rate (rETR) and effective quantum yield (F v '/F m ') were higher in the states II and III, and state III has higher NPQ than I and II (Figures 2D-F). NPQ is one of the adaptive mechanisms employed by plants and algae for coping with complex environmental change that works by two independent mechanisms, with the faster-activated quenching catalyzed by the monomeric light-harvesting complex (LHCII) proteins and the slowly activated quenching catalyzed by LHCII trimers, both processes depending on zeaxanthin but to a different extent . Previous study showed that the induction of NPQ was regulated by the generation of zeaxanthin . In this study, state I was observed at the end of the light period whereas state II was observed at the end of the dark period, the genes encoding zeaxanthin synthesis were upregulated from state I to state II, which was consistent with the previous study, showing that the content of zeaxanthin could be enhanced by high light (Eismann et al., 2020). But compared to state II, the expression of genes encoding for violaxanthin de-epoxidase (VDE; EC: 1.23.5.1) decreased by 2-4 folds, whereas the expression of genes encoding for betaring hydroxylase (LUT) and zeaxanthin epoxidase (ZEP; EC: 1.14.15.21) was upregulated (Supplementary Table 2) in state III. The NPQ in state III was higher than that in state II, but because we did not measure the content of zeaxanthin, the relationship between NPQ and zeaxanthin, as well as the role of zeaxanthin during the sporulation should be further studied. The other chlorophyll fluorescence parameters (F v /F m , F v '/F m ' , NPQ, rETR max , and α) were significantly decreased in state IV (Figure 2). The rETR max and α provide an approximation of the rate of electron flow via PS II into the photosynthetic electron transport chain and the efficiency of energy transfer from the light-harvesting antenna to PSII reaction centers, respectively, and these values also reflect the overall photosynthetic capacity of algae (Serodio et al., 2006).
Moreover, the highest chlorophyll content (Chl a, Chl b, and Car) was observed in the state II cells (Figure 3A), but gradually decreased in the state III and IV. Although there were no significant differences in the ratio of Chl a to Car or Chl b to Car between state III and IV (Figure 3B), the state IV cells showed the lowest chlorophyll content but highest Car percentage (23%) in the total chlorophyll content (Chl a + Chl b + Car) ( Figure 3A). This indicated that the proportion of carotenoids was synchronized with the formation of spores, which was consistent with the previous study (Park, 2020). This may be the reason why the algae color was yellow in states III and IV. Additionally, carotenoid as a type of photoprotection pigment, which indicated higher antistress ability of algae (Siefermann-Harms, 1987), can protect the membranes from oxidants, such as reactive oxygen species (ROS), generated in the spore formation process of U. prolifera.

De novo Assembly of Transcriptome and Gene Functional Annotation
RNA-sequencing using the DNBSEQ platform generated a final 10.56 Gb of clean data. De novo assembly obtained 41,779 unigenes with an N 50 length of 2,619 bp by the Trinity method. The total length was 69,331,130 bp, and the average length was 1,659 bp. After BLASTx searches (E-value ≤ 10 −5 ) against seven public databases, including the NCBI non-redundant protein sequences (Nr), NCBI Nucleotide Sequence Database (Nt), EuKaryotic Orthologous Groups database (KOG), Swiss Prot database (Swiss Prot), Pfam database (Pfam), KEGG, and GO (Supplementary Table 3), the transcriptome expression profile analysis of U. prolifera was performed using 12 cDNA samples, produced 21.68M clean reads each sample in average, and the average mapping rate to the transcriptome described above was 92.54%, and a total of 41,058 unigenes were identified.
To validate the reliability and reproducibility of the RNAseq analysis, ten unigenes were randomly selected from different metabolism pathways for qRT-PCR validation. The expression patterns of these selected unigenes in our transcriptome data were highly consistent with the qRT-PCR results ( Supplementary  Figure 1), indicating that the RNA-seq data were credible and sufficient for subsequent analysis.

Principal Component Analysis and Correlation Analysis
Principal component analysis (PCA) was performed to get an overview of the transcriptome difference between different states (Supplementary Figure 2A). The first principal component (PC1), accounting for 34.5% of the total variance, clearly separated samples of II, III, and IV states from those of state I, indicating that spore formation could significantly change the transcriptional component of algal cells. Likewise, correlation analysis results was consistent with that of PCA. Correlation coefficients among samples of the three biological replicates were all greater than 0.9 (Supplementary Figure 2B).

Transcriptomic Changes of Thalli in Different Spore Formation States
The transcriptomic profiles of state II with I, III with II, and IV with III were compared. In detail, compared with state I, a total of 21,149 unigenes, including 11,645 upregulated unigenes (URGs) and 9,504 downregulated unigenes (DRGs) were founded in state II (Supplementary Figure 2C). The enrich result of URGs and DRGs are shown in Supplementary Figures 3, 4. The ribosome would encode amino acids to form polypeptide chain by utilizing mRNA as the model and produce protein for microtubule protein, transporter, cell skeleton protein, and so on. The upregulation of genes in the photosynthesis pathway explained the increased photosynthetic activity in state II with I (Figure 2). Photosynthesis also could provide several stress-resistant mechanisms to affect photosynthesis activity and provide power for other metabolism pathways. A previous study found that ubiquinone and other terpenoid-quinone biosynthesis pathways could provide precursors for carotenoid biosynthesis (Lohr et al., 2012), and the carotenoid could affect photosynthesis (Eismann et al., 2020).
Compared with state II, a total of 20,856 DEGs, including 9,539 URGs and 11,317 DRGs were found in state III (Supplementary Figure 2C). Purine metabolism may be associated with the repair of PS II because defective purine biosynthesis impairs chloroplast biogenesis and alters the induction of high light stress pathways in Arabidopsis (Woo et al., 2011).
Compared with state III, a total of 15,185 DEGs, including 6,473 URGs and 8,712 DRGs were found in state IV (Supplementary Figure 2C). Hormones play essential roles in the germination, growth, and development of plants by signaling pathways (Li et al., 2016), such as abscisic acid (ABA), jasmonate (JA), ethylene, and cytokinin (CK) (Nakashima and Yamaguchishinozaki, 2013).

Differentially Expressed Genes During the Sporulation Process
Photosynthesis Photosynthesis products are the main energy sources for algae that include light and dark reactions, each of which consists of a series of redox reactions related to energy capture and conversion (Liu, 2016). At the initial step of photosynthesis, light energy is captured by the peripheral light-harvesting complexes, and then the synthesis of ATP and NADPH was affected by regulated F v /F m (Yang et al., 2020). Compared with state I, the genes related to photosynthesis participating in the constitution of photosystem I (PS I) (PsaA, PsaB, PsaD, etc.) ( Figure 4A, PS I), photosystem II (PS II) (PsbD, PsbC, PsbB, etc.) (Figure 4A, PSII), cytochrome b6/f complex (cytob6/f) (PetA, PetC, PetN, etc.), ATP synthase (F-type ATPase), and photosynthethic electron transport (PET) were upregulated in the state II but were downregulated in state III ( Figure 4A and Supplementary Table 4). A previous study showed PsaB encodes P700 chlorophyll A2 apoproteins, one of PS I reaction, and the upregulation of this gene could significantly increase photosynthetic electron transport (Qian et al., 2010), resulting in an increase in photosynthesis (Figure 2). However, the gene expression of photosynthesis was also upregulated in the state D while the chlorophyll fluorescence parameters were decreased (Figure 2). We speculated that the polar growth of algae segment maybe one reason, but it needed to be further studied. Additionally, compared with state I, the gene encoding for ferredoxin NADP reductase (FNR) of NADPH synthesis was upregulated 5. 3-, 1. 6-, and 2-folds in the state II, III, and IV, respectively, indicating higher NADPH and photosynthesis electron transport of PS I in state II, but lower in III and IV. NADPH could enhance CO 2 assimilation by providing power (Sui et al., 2015).
The DEGs among four states in carbon fixation pathway are shown in Figure 4B. The rbcl gene, the key gene in carbon fixation, encoding ribulose-bisphosphate carboxylase (Rubisco), was upregulated (Figure 4B), indicating the increase in photosynthesis in U. prolifera (Qian et al., 2010). The rbcl gene could utilize CO 2 to catalyze ribulose-1,5-bisphosphate (RuBP) to 3-phosphoglycerate (3-PGA), and during this process, the role of carbonic anhydrases (CA), C 4 enzymes [phosphoenolpyruvate carboxylase (PEPC; EC: 4.1.1.31), and phosphoenolpyruvate carboxykinase (PEPCK; EC: 4.1.1.49)] cannot be ignored . In detail, atmospheric CO 2 could be transferred into bundle sheath cells by PEPC and PEPCK, ensuring efficient transfer of CO 2 to Rubisco to participate in 3-PGA synthesis (Hatch, 1987), and finally generate 3-phosphoglyceraldehyde (3-PGAld) ( Figure 4B). However, the gene expression of PEPCK insignificantly changed from state II to III, whereas the rbcl gene expression was downregulated 2-4 folds, indicating that redundant CO 2 existed in the cell, which could decrease the pH of the thylakoid membrane, and which could produce NPQ . Moreover, oxaloacetate also could be catalyzed by malate dehydrogenase (MDH; EC: 1.1.1.37) to form malic acid which could be catalyzed by NADP-malate dehydrogenase (NADP-ME; EC: 1.1.1.40) to provide CO 2 to Rubisco (Yang et al., 2020). The expression of the gene encoding for MDH was upregulated 2.6 folds while NADP-ME was downregulated 2.1 folds from state II to state III (Supplementary Table 5), indicating the malic acid synthesis increased. Malic acid is an organic acid that could help cells to transduce signal in the stress condition by upregulating gene expression (Kumar et al., 2016) and stimulate spore formation (He et al., 2019b). Additionally, compared with state II, rbcl gene was downregulated in state III, which might decrease 3-PGA synthetic and inhibit the electrons transport of Cytochrome b6/f complex and PET, but the F v '/F m ' has no change, which was measured from PS II, and the PsaB expression was upregulated 2.1 folds. The PsaB is the key gene of PS II (Qian et al., 2010), that might enhance the oxidation of the plastoquinone (PQ), which was caused by plastid terminal oxidase (PTOX). We also found the expression of gene encoding PTOX was upregulated from state II to state III. The role of PTOX was transferring electrons from PQH 2 to oxygen during the oxidation of PQ pool to enhance the oxidation state of the plastoquinone pool (Nawrocki et al., 2015). The enhanced oxidation of the plastoquinone pool could promote the formation of spores (Wang et al., 2016) and catalyze NADP to NADPH to provide power for MDH catalyze reaction (Chen et al., 2020). Then, we speculated that PTOX was one of the key genes for regulating spore formation.
The 3-PGAld was a photosynthetic product that could be converted to starch in the chloroplast stroma or enter the cytoplasm via triose phosphate translocator and synthesize sucrose (Yang et al., 2020). It was also the raw material for glycolysis and terpenoid backbone biosynthesis. A previous study showed that glucose-1-phosphate adenylyltransferase (APG; EC: 2.7.7.27), starch synthase (EC: 2.4.1.21), and 1,4-alpha-glucan branching enzyme (GBE; EC: 2.4.1.18) were involved in starch synthesis (Stark et al., 1992), and in this study, the genes responsible for these enzymes were upregulated in state II ( Figure 5A). Additionally, the synthesis of 3-PGA was increased in carbon fixation pathway in state II (Figure 4B), which was consistent with the previous study, that is the APG could be regulated by 3-PGA and P i as positive and negative effectors, respectively (Preiss, 1991). Glucose 1-phosphate (G1P) as the substrate of ADP-glucose (ADPG) was also increased, catalyzed by fructose-bisphosphate aldolase (EC: 4.2.1.13), fructose-1,6bisphosphatase I (EC: 3.1.3.11), glucose-6-phosphate isomerase (EC: 5.3.1.9), and phosphoglucomutase (EC: 5.4.2.2). The expression of genes of these enzymes were upregulated by 7-9.2, 2.6-3.5, 2.8, and 3.5 folds, respectively. The starch catabolism genes, alpha-amylase (EC: 3.2.1.1) and glycogen phosphorylase (EC: 2.4.1.1), were downregulated 3.2-4.6, 2.1-36.8 folds, implying that the starch content might increase, and the enzymes related to starch synthesis also could be used to synthesize sucrose ( Figure 5B). Starch, as the major carbohydrate, could provide power for spore formation (Wang et al., 2016). Notedly, the chloroplast disappeared in state III, but the spore would produce chromatophore at the anterior end (Chen et al., 2011). Chromatophore was the product of primary endosymbiotic integration of a cyanobacterium into a eukaryotic (Reyes-Prieto et al., 2007). Therefore, we suggested that chromatophores have the photosynthesis power as chloroplast and that polar growth occurred in state IV, which may be the reason why the gene expression of starch synthase was downregulated in state II to III and III to IV, but APG was upregulated in the state III to IV, which might prepare the precursor for starch synthesis for the growth of later seedlings.

Glycolysis
Glycolysis is a fundamental metabolic pathway that breaks down glucose into two three-carbon metabolites and generates energy in the cytoplasm. The DEGs in glycolysis are shown in Figure 5B. Compared with state I, the genes encoding for hexokinase (EC: 2.7.1.1) and 6-phosphofructokinase 1 (EC: 2.7.1.11), which are the rate-limiting enzymes in glycolysis pathway (El-Gharbawy and Koeberl, 2014), were downregulated in state II while the expressions of genes encoding fructose-bisphosphate aldolase, glyceraldehyde 3-phosphate dehydrogenase (EC: 1.2.1.13), and phosphoglycerate kinase (EC: 2.7.2.3) were upregulated. So, we speculated that the 3-PGA increased in the glycolysis pathway because the algae can transfer 3-PGA from carbon fixation pathway to glycolysis via triose phosphate translocator. The increase of 3-PGA could replace glucose as the raw material downstream for glycolysis to produce pyruvate and produce a negative effect for upstream, causing downregulation of the gene expression of hexokinase and 6-phosphofructokinase and upregulation of pyruvate kinase (EC: 2.7.1.40), and all this is shown in the Figure 5B. These results also showed that fructose-6 -phosphate would increase, which was the substrate of sucrose, and the gene expression of sucrose-phosphate synthase (SPS; EC: 2.4.1.14) and sucrose-6-phosphatase (SPP; EC: 3.1.3.24) was upregulated in state II, indicating the glycolysis pathway was downregulated to prepare energy for spore formation in states III and IV, which also was observed in previous study (Wang et al., 2016).

Carotenoid Biosynthesis Pathway
Moreover, the expressions of genes encoding for Car biosynthesis pathway were significantly upregulated in three states (II, III, and IV) (Figure 6B), as compared with those in state I, such as phytoene synthase (CrtB), phytoene desaturase (PDS), zetacarotene desaturase (ZDS), zeta-carotene isomerase (Z-ISO), prolycopene isomerase (crtISO), lycopene beta-cyclase (CrtLb), and beta-carotene 3-hydroxylase (CrtZ) were upregulated 8. 6, 2-22.6, 6.5-16, 2.5, 3.7-34.2, 3.5-64, and 3.2-10.6 folds, respectively ( Figure 6B). CrtB is the limiting enzyme for carotenoid biosynthesis which could catalyze two GGPPs to form phytoene (Eismann et al., 2020), and compared with state I, the gene expression was upregulated 8.6, 22.6, 6.5 folds in states II, III, and IV, respectively. These results verified the changes in Car content, which was observed in Figure 3, and Car, in addition to photosynthesis, protects membranes from excess of photon energy by producing NPQ by xanthophyll cycle  and affects spore formation (Park, 2020). Moreover, some studies showed that carotenoid biosynthesis is related to the accumulation of astaxanthin in Haematococcus pluvialis and Chlorella zofingiensis (Zhang et al., 2016;Du et al., 2021), whereas the expression of CrtZ related to the astaxanthin synthesis pathway was upregulated during the sporulation process ( Figure 6B). However, it is still uncertain whether astaxanthin exists in U. prolifera, and the further research is needed.

Transcriptional Factors
Transcription factors are the main regulators that control the expression of genes involved in plant growth, development, and responses to external stimuli (Ahammed et al., 2020). They activate or inhibit genes by binding to promoter functional elements or interacting with other proteins (Ng et al., 2018). In our study, a total of 219 of unigenes involved in 29 TF families were detected during the process of sporulation. Among them, many genes were involved in the AP2-EREBP family, the MYB family, the G2-like family, and the WRKY family whereas the ARF family, Alfin-like family, bHLH family, and bZIP family have only one gene ( Figure 7A). The identified 42 of the DEGs, covering 19 families, were constantly significantly upregulated during the process of the spore formation ( Figure 7B, label 1), including four G2-like genes, two C2C2-CO-like genes, two PLATZ genes, one bHLH gene, six WRKY genes, seven AP2-EREBP genes, etc. (Supplementary Table 6), whereas 23 of the DEGs were significantly downregulated ( Figure 7B, label 2). Additionally, most of the orthologs for the 29 URGs in other species were hypothetical or predicted proteins. Five TFs unigenes which involved four families seemed to be unique in U. prolifera because no orthologs could be found in other species (Supplementary Table 7), and the expressions of CL2970.Contig1_mix12 and Unigene3886_mix12 were significantly downregulated in the process of spore formation. The expressions of Unigene4015_mix12 and Unigene17383_mix12 were significantly upregulated in states III and IV, respectively, as compared with A (Supplementary Table 7). Notedly, states III and IV were the crucial periods of spore formation according to Figure 1. Therefore, we speculated that those two genes would affect spore formation at different times, but need to be further studied.
The previous study showed that the bZIP TF is involved in regulating adventitious shoot development in Paulownia kawakamii (Low et al., 2001). In our study, the expression of Unigene14540_mix12, which belonged to bZIP family, was highly upregulated in the state IV as compared with other states (including state I, II, and III), indicating that polar growth may have occurred in state IV, and this also was consistent with the results of the light saturation parameter (E k ), which increased in state IV ( Figure 2C). The TFs of RWP-RK family have profound and broad influence on the early transcriptome response to nitrate signaling and assimilation (Marchive et al., 2013). Moreover, some TFs of RWP-RK family regulate the storage of starch by controlling carbon partition in Chlamydomonas (Chardin et al., 2014). MYB and WRKY TFs also played important roles in regulatory signal transduction and responses to biotic and abiotic stresses (Smita et al., 2015;Wu et al., 2018b).
The functions of the remaining TFs can be divided into three primary GO categories according to GO annotation, including biological process, cellular component, and molecular function (Supplementary Figure 5). Notedly, the functions of CL1966.Contig2_mix12 of C3H family in the biological process was reproduction and reproduction process, in which the gene expression was upregulated in states II and III whereas it was downregulated in stage IV (Supplementary Figure 6). Moreover, the expression of this gene was the highest in stage III, which is also the main stage of spore formation. All of these results indicated that the main process of reproduction occurred in state III and ended in state IV, which also verified the result of Figure 1 and our speculation that the polar growth might occur in state IV. These results implied that the CL1966.Contig2_mix12 may be one of the key genes in the spore formation. However, the involvement of specific functions of these TFs in the spore formation of U. prolifera merits further in-depth investigations. Additionally, Unigene5846_mix12 has an identity above 40% with diguanylate cyclase, and its overexpression could decrease the sporulation frequency and early sporulation gene transcription of Clostridioides difficile (Edwards et al., 2021), but it did not significantly change in this study, indicating that the spore formation was not inhibited.

Flagella Biosynthesis
The spore of U. prolifera has four flagella , and flagella are one of the most ancient cellular organelles, constituted by skeleton protein, providing motility for living in water which is vital for algal cells to survive and proliferate under different environmental conditions (Silflow and Lefebvre, 2001). In this study, 574 of the DEGs associated with flagella biosynthesis in U. prolifera were annotated according to GO annotation (Supplementary Table 8). To understand the expression of genes associated with flagella synthesis in different states for better, we arranged and divided these genes into different groups, noting labels as shown in Figure 8A. Among these genes, 37 of the DEGs were constantly upregulated in the process of spore formation (Figure 8A, label 1), but 12 of the DEGs were constantly downregulated ( Figure 8A, label 6); 20 of the DEGs were upregulated in II state but downregulated in states III and IV ( Figure 8A, label 3). Although the flagella did not synthesize in state II (Figure 1), the ribosome pathway was significantly upregulated in state II (Supplementary Figure 4A). Therefore, these genes may prepare material for flagella synthesis in the later period, such as protein and cell skeleton, and assemble only the 146 DEGs that were upregulated in state III to produce flagella ( Figure 8A, label 4) because the spore was observed in state III according to morphology (Figure 1). Only 44 of the total DEGs that were upregulated in state IV may be associated with spore movement (Figure 8A, label 5), the previous study showed that spore could keep rotating depending on the flagella with much higher speed and collided with the cytocyst wall, and a hole was formed at the top of the cytocyst where germ cells were released one by one in a kind of amoeboid movement . Moreover, 243 of the DEGs were upregulated in states II and III state while they were downregulated in IV may have the preparation and assembly effects to flagella synthesis ( Figure 8A, label 2). We also found that including 20 of the DEGs encoding intraflagellar transport proteins (IFTs) has the highest expression level in state III, but decreased in IV ( Figure 8B). The assembly and maintenance of flagella require the involvement of IFTs, which are also closely related to flagella movement, signal reception, and signaling (Croft et al., 2018;Wang et al., 2018b). However, those genes were annotated in Chlamydomonas reinhardtii and were not reported in previous transcriptome analysis of U. prolifera (He et al., 2019a;Wang et al., 2019; Table 1). Combining the phenomenon of state III, we suggest that state III is a critical period for sporulation and provides an observable expression gene for subsequent exploration of how environmental factors affect sporulation.

CONCLUSION
The spore formation of U. prolifera was divided into 4 states, cells in the state II accumulated energy by enhancing photosynthesis and starch synthesis and inhibiting catabolism, and provide power for the spore production in the state III and IV. The synthesis of carotenoid had a high degree correlation with spore formation. According to the expression of gene associated with flagella and intraflagellar transport proteins biosynthesis in U. prolifera, state III is a critical period for sporulation. Moreover, some genes of TF families, such as C3H family, may be one of the key regulating genes in spore formation. These results provide a new perspective of spore formation at the gene transcriptional level although the detailed transcription regulatory spore formation network remains to be unveiled.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are publicly available. This data can be found here: https://www.ncbi.nlm.nih. gov/sra/?term=PRJNA794305.

AUTHOR CONTRIBUTIONS
YL and NX: conceptualization, funding acquisition, supervision, and writing-review and editing. JJ: data curation and writingoriginal draft. JJ, YL, CH, WZ, and NX: formal analysis. JJ and YL: visualization. All authors: contributed to the article and approved the submitted version.