The strain D.salina FACHB-435 was acquired from the Freshwater Algae Culture Collection at the Institute of Hydrobiology (FACHB) in Wuhan, Hubei. D. salina strain was cultured in liquid DM medium with a trace element solution diluted to 10% with Milli-Q water, and the pH value was adjusted to 7.5 after medium preparation. The cultures were incubated at a temperature of 25 ± 2°C with a photoperiod of 12/12 h light/dark cycle, and the culture bottles were shaken at least 6 times a day to ensure proper mixing. To isolate a sterile D. salina cell, techniques such as capillary separation, plate marking, and antibiotic sterilization were employed. The identification of D. salina FACHB-435, through both morphological and molecular biology methods, verified its classification within the Dunaliella genus (Figure 1A). The algal liquid cultured to the later stage of logarithmic growth was centrifuged and inoculated into the newly prepared medium with the same amount. The normal light (NL) 2,500 Lux as control groups and the high light (HL) 10,000 Lux as treatment groups were set, and three replicates were set for each experimental group. The inoculated algal liquid was placed on the light culture rack at the culture temperature of (25 ± 1) °C, the light and dark period L/D was 12 h/12 h, and the bottle was shaken 6 times a day. After inoculation, the growth state of algal cells was observed. At 0, 5, 13 and 20 days of culture, algal fluid was taken for dry weight, β-carotene content detection and RNA sequencing, respectively.

The determination of β-carotene in D. salina by high-performance liquid chromatography according to the following process: a 10 mL aliquot of algal suspension was centrifuged at 5,000 rpm for 5 min, the supernatant was discarded, and the algal biomass was transferred to a 2 mL Eppendorf tube. The sample was flash-frozen in liquid nitrogen for a few seconds before being stored in a − 80°C freezer for several hours and then freeze-dried using a lyophilizer. The freeze-dried algal powder was stored in a − 80°C freezer. To the algal powder, 2 mL of acetone was added, and the mixture was ground using a grinder at a frequency of 50 Hz for 3 min, until the algal powder turned white. After allowing the mixture to settle for 20 min, the supernatant was collected, centrifuged, and filtered through a 0.22 μm organic membrane into an HPLC sample vial. Since D. salina cells have a high lipid content that can interfere with β-carotene extraction, causing incomplete extraction or emulsification, a saponification reaction was performed prior to pigment extraction. The lipid extract was dissolved in 1 mL of toluene, and 100 μL of 0.1 g/mL KOH-methanol saponifying agent was added. The mixture was allowed to stand at room temperature for 1 h. The supernatant was then collected, centrifuged, and filtered through a 0.22 μm organic membrane into an HPLC sample vial. The sample was stored at −80°C until further analysis.

The samples were analyzed using an Agilent 7890B gas chromatograph equipped with a capillary column (HP88, 100 m, 0.25 mm, 0.2 μm film thickness) and a hydrogen/air flame ionization detector (FID). Helium served as the carrier gas at a constant flow rate of 20 mL/min. The injector temperature was maintained at 250°C, and a 1 μL sample volume was injected with a 10: 1 split ratio. The temperature program started at 140°C with a 5-min hold, then increased to 240°C at a rate of 4°C/min, followed by a 12.5-min hold. The mass spectrometer (MS) operated in electron impact ionization mode with an electron energy of 70 eV. The other parameters included: full scan mode with an m/z range of 50–550; ion source temperature at 230°C; and quadrupole temperature at 150°C. The National Institute of Standards and Technology (NIST) library database was used to identify the components. The standard curve was prepared using C8-22 fatty acid methyl esters and quantified by the internal standard method. Each experiment was independently replicated three times.

Genomic DNA was isolated from 100 mL of fresh culture. The cell suspension was centrifuged at 8000 rpm for 5 min. The cell pellet was then suspended in approximately 1 mL of medium and pipetted into a 2-mL tube and centrifuged again at 8000 rpm for 5 min. Briefly, 800 μL of 2% mercaptoethanol solution was pipetted into each sample, followed by the addition of 800 μL of 10% w/v CTAB (in 0.7 M NaCl solution) and incubation at 56°C for 10 min. After extraction with an equal volume of phenol: chloroform: isoamyl alcohol (25:24:1), the mixture was centrifuged at 12,000 rpm for 10 min twice. The supernatant was dissolved in 0.4 mL of 100 μg/mL RNase and incubated at 37°C for 30 min. An equal volume of chloroform: isoamyl alcohol (24:1) was then centrifuged at 12,000 rpm for 10 min. Genomic DNA was precipitated by adding 2.5 volumes of 100% ethanol and collected by spinning at 12,000 rpm and 4°C for 10 min. After the supernatant was discarded, the resulting genomic DNA pellet was stored in 5.0 mL of 70% cold ethanol at 4°C overnight to allow the impurity to dissolve. Finally, DNA was eluted in 100 μL of 10 mM Tris HCl by centrifugation at 12,000 rpm for 1 min. The purity and concentration of DNA were analyzed using a NanoDrop 2000 Spectrophotometer (Thermo Scientific, United States). More than 5 μg of sheared and concentrated DNA was applied to size-selection by the BluePippin system. Approximately 20-kb SMRTbellTM libraries were prepared according to the released protocol from PacBio company. A total of 6.7 Gb subreads were sequenced on the PacBio sequel system, i.e., 106 × coverage of the estimated genome size. After removing the low-quality (containing 10 or more Ns and low-quality bases with quality scores ≤7) and redundant reads.

In this study, Jellyfish software (Marçais and Kingsford, 2011) was used to analyze the distribution of k-mer. FACHB-435 genome was assembled based on the third-generation PacBio data, and the FACHB-435 initial genome was corrected using the second-generation Illumina data. Based on the third-generation sequencing data, we first compared two sets of third-generation assembly software:Canu (Koren et al., 2017) and Mecat (Xiao et al., 2017).

BUSCO (Viridiplantae_odb10), EST, and Illumina alignment were used to evaluate the completeness and quality of the assembled genome (Camacho et al., 2009; Langdon, 2015; Simão et al., 2015). BUSCO software is used for the accuracy and integrity of genome assembly, gene sets, and transcripts. BUSCO software. The OrthoDB database constructs single-copy gene sets for several large evolutionary branches. The gene set is compared with the assembled genome, and the accuracy and completeness are assessed according to the ratio and completeness of the comparison. Using the complete transcript, mapping to the sequence of the spliced sample is used to assess the integrity of the gene coding region.

The repetitive sequence prediction of D. salina first uses RepeatModeler to obtain a denovo repeat sequence library, which can automatically execute the two denovo search software RECON and RepeatScout, after obtaining the consistent TE, import the result into RepeatMasker to identify and cluster the repeat sequences. Repeat sequences are classified using a process for identifying them (Flynn et al., 2020). In order to obtain high-quality genes, this project used de novo prediction (Augustus, SNAP, etc.) and homologous prediction software (GeneWise) combined with transcriptome data to predict gene structure.

First, RNAseq samples were selected and imported into Trinity de novo assembly and Trinity genome-guided assembly (Grabherr et al., 2011). Then RSEM was used to calculate transcript abundance (filter out FPKM< 1 and iso-percentage < 3%; Filtered transcripts are imported into PASA software to build more comprehensive transcripts). PASA is able to detect new transcripts using high sensitivity based on reference assembly, combined with the ability to assemble from scratch, and then train the identified transcripts as parameters. By comparison with the Uniprot plant protein database, transcripts with protein coverage greater than 95% were retained as candidate transcripts. Then, the gene prediction software SNAP, GENEMARK and AUGUSTUS based on ab initio theory were used to extract and retrain the gene model with AED score of 0 (Holt and Yandell, 2011), and the final gene prediction result was obtained by integrating MAKER software (Stanke et al., 2006). Through gene structural annotation, a total of 30,752 protein-coding genes were predicted. The predicted gene structure was evaluated by BUSCO.

By employing OrthoFinder (version 2.3.12), orthologous genes were found between D.salina and 12 species, including Cyanidioschyzon merolae, Porphyriopsis chorda, Gracilariopsis chorda, Porphyra umbilicalis, Isochrysis galbana, Emiliania huxleyi, Phaeodactylum tricomutum, Chlamydomonas rheinhardtii, Thalassiosira pseudonana, Chlamydomonas eustigma, Gonium pectorale, Volvox carteri. The maximum likelihood tree was constructed with rhodophyta as an outgroup. Homologous base pairs between species were identified using jcvi’s compara tool (Tang et al., 2015). Then, according to the synonymous_calculation module in the tool bio-pipeline, synonymy replacement rate (Ks) distribution of homologous base pairs obtained in the previous step was analyzed, and the whole genome replication events were inferred (Zhang et al., 2006; Rodamilans et al., 2014).

21 samples were used for RNA sequencing: Dsa-0d1, Dsa-0d2, Dsa-0d3, Dsa-NL5d1, Dsa-NL5d2, Dsa-NL5d3, Dsa-NL13d1, Dsa-NL13d2, Dsa-NL13d3, Dsa-NL20d1, Dsa-NL20d2, Dsa-NL20d3 (control: samples cultured under the normal light 2,500 Lux), and Dsa-HL5d1, Dsa-HL5d2, Dsa-HL5d3, Dsa-HL13d1, Dsa-HL13d2, Dsa-HL13d3, Dsa-HL20d1, Dsa-HL20d2 and Dsa-HL20d3 (treatment: samples cultured under the high light 10,000 Lux). Total RNA was extracted by the instruction manual of RNA extraction reagent kit Trans Zol Up Plus RNA Kit (Beijing TransGen Biotech Co., Ltd). The transcriptome library were constructed and sequenced using the Illumina Hiseq4000 platform.

The expression level of the gene is represented by the FPKM value, which represents the number of sequenced fragments included in every thousand sequencing bp per million sequencing bases. The Cufflinks v2.2.1 software was used to calculate the expression level of the genes (Roberts et al., 2011). After calculating the expression levels of all genes, we identified and statistically analyzed the differentially expressed genes (DGEs) between the control group and high-light stressed D. salina at two time points (0dVs5dNL, 5dNLVs13dNL, 13dNLVs20dNL, 0dVs5d-HL, 5d-HL Vs13d-HL, 13d-HL Vs20d-HL), searched for common DEGs, and performed clustering analysis on all DEGs. The standard for defining DEGs is: p-Value<0.05, Fold change>2 (Nikolayeva and Robinson, 2014). DEGs can be divided into two categories based on their relative expression levels: Up-regulated Genes and Down-regulated Genes.

Using the RNA from the control and experimental groups on the 5th and 13th days of culture as templates, reverse transcription was performed according to the instructions of the TransScript One-Step gDNA Removal and cDNA Synthesis SuperMix kit(Beijing TransGen Biotech Co., Ltd.) The obtained cDNA was used as a template, and fluorescence quantification detection was performed using the ABI Step One fluorescent quantitative PCR instrument, according to the instructions of the Quanshijin TransStart®TOP Green qPCR SuperMix kit. We chose three genes, DsLCYB, DsPDS and DsZDS, to validate.

We investigated the growth and β-carotene accumulation patterns of D. salina under four salinity gradients: 0.5 M, 1.5 M, 2.5 M, and 3.5 M (Figure 1B). Figure 1C shows the dry weight of algae at different salinity gradients, with the dry weight in each group increasing as culture time progressed. The dry weight of the low salinity 0.5 M group and the control group were similar, while the high salinity experimental group consistently had a lower dry weight than the control group. Under high salinity conditions, cell numbers increased slowly and cell density remained lower than in the control group. Figure 1D presents the β-carotene content of D. salina FACHB-435 at different salinities. The β-carotene content per cell initially decreased and then gradually increased during normal culture, stabilizing as culture time increased. Upon resuming culture, the β-carotene content per cell increased with salinity. On day 15, the high salt stress group’s β-carotene content per cell was significantly higher than that of the control group, with an increase of 72.15% (p < 0.01). The figure also displays the biomass of D. salina under different light intensity conditions, revealing that the biomass increases with light intensity. It can be observed that, on the 10th day of incubation, the β-carotene accumulation in individual algal cells under high light intensity was significantly higher than that in the control group. As incubation time increased, β-carotene accumulation in single cells decreased, with higher light intensities correlating to lower β-carotene content (Figure 1D).

Based on the 17-K-mer results, it showed that the haploid genome size of of D. salina, we FACHB-435 was 386.49 Mb. Genome homozygosity (100%- genome heterozygosity) was 99.713%, homozygosity > 99%, indicating high homozygosity of genome (Supplementary Table S1). Only a single peak indicates low genome complexity, and there is no obvious peak behind the main peak, which further indicates that the repeat sequence is low, indicating good sequencing results (Supplementary Figure S1).

In this study, genomes were assembled with three-generation data and corrected with second-generation data. After two analysis of three-generation assembly software, it is obvious that the assembly of MECAT software (Genome Size: 472 Mb N50:458 kb) is better than the assembly results of CANU software (Genome Size: 511 Mb N50:185 kb) (Table 1). Our assessment of genomic assembly complete showed that the overall complete BUSCOs was 91.0%, where complete and single-copy BUSCOs 42.2%, complete and duplicated BUSCOs 48.8%. The proportion of complete and duplicated BUSCOs is relatively high. Fragmented BUSCOs was 6.3%, Missing BUSCOs 2.7% (Supplementary Table S2).

The annotation results were evaluated by BUSCO, as shown Supplementary Table S3. After the whole genome annotation is completed, 303 lineal homologous gene sets in eukaryota_odb9 are evaluated. The results showed that completeness was 74.2%, complete and single-copy BUSCOs was 27.7%, and complete and duplicated BUSCOs 46.5%, which was slightly higher. The fragments were 3.3% and the missing BUSCOs 22.5%. The integrity of gene annotation is high, which can be used for metabolic studies related to D. salina (Supplementary Table S3).

The repetitive sequence prediction of D. salina first uses RepeatModeler to obtain a de novo repeat sequence library, which can automatically execute the two de novo search software RECON and RepeatScout, after obtaining the consistent TE, import the result into RepeatMasker to identify and cluster the repeat sequences. The process of identifying repeat sequences was used to classify them, and the resulting repeat sequence classification was shown in Supplementary Table S4.

RNAseq samples were selected and imported into Trinity de novo assembly and Trinity genome-guided assembly. Then RSEM was used to calculate transcript abundance (filter out FPKM<1 and iso-percentage < 3%). Filtered transcripts are imported into PASA software to build more comprehensive transcripts. PASA is able to detect new transcripts using high sensitivity based on reference assembly, combined with the ability to assemble from scratch, and then train the identified transcripts as parameters. By comparison with the Uniprot plant protein database, transcripts with protein coverage greater than 95% were retained as candidate transcripts. Then, the gene prediction software SNAP, GENEMARK and AUGUSTUS based on ab initio theory were used to extract and retrain the gene model with AED score of 0, and the final gene prediction result was obtained by integrating MAKER software. Through gene structural annotation, a total of 30,752 protein-coding genes were predicted. The predicted gene structure was evaluated by BUSCO, and the results were shown in Supplementary Table S5. By BUSCO analysis (Supplementary Table S5), the total completeness was 75.9%, with 44.9% gene replication accounting for a slightly higher proportion. Fragmented is only 4% and Missing 20.1%. On the whole, if the deletion of conserved genes within the species is not taken into account, the quality of the genome assembled by us is medium.

According to the evaluation results, the replication integrity of D.salina takes up 44.9% of the total, which is analyzed by RepeatMasker software, as shown in Supplementary Table S5. The Total repeat fraction is 43.75% of the genome, accounting for about half of the genome, among which the Class-I of TEs accounts for a relatively high proportion, accounting for 30.25% of the genome. The specific proportion distribution of TEs in repeated sequences is shown in Supplementary Table S4, where class-I accounts for 69%, Class II 25% and Unknown 2%. The above data all indicate that D. Salina genome is a highly repetitive genome (Supplementary Table S4). 77,287 ESTs were constructed by means of transcriptome analysis without parameter assembly, and 96.21% of the ESTs were covered by the EST database. The results of comparison were shown in Supplementary Table S5, indicating that the integrity of the genome coding region was high.

C. eustigmais, an acidophilic extremophile, is also an extreme organism. Both C. eustigmais and D. salina are representative of extremophiles within the Volvocales order (Hirooka et al., 2017). G.pectoral and V. carteri are representative multicellular taxa in Volvocales (Prochnik et al., 2010; Hanschen et al., 2016). Using C. rheinhardtii as a neutrophilic relative reference for polar organisms, and based on current DNA sequencing analysis, it is the most likely ancestral species that evolved from single-celled to multicellular groups. So we selected these most representative species for our next comparative genomics analysis (Merchant et al., 2007), G. pectorale and V. carteri, along with the current genome, enabled us to analyze the evolutionary relationships of species within the Volvocales family using OrthoFinder software (Figure 2A) (Emms and Kelly, 2015). This provides a reference for inferring the genome content and structure of the most recent common ancestor of all extant Volvocales, thereby improving our understanding of their origins and evolution. The results showed that 10,934 gene families were clustered, including 1,020 single-copy orthologous gene families (Figure 2A). After extracting the gene families of each species, the total number of gene families for the five species was 5,180, with 41 (0.4%) specific to C. eustigma, 21 (0.2%) specific to C. reinhardtii, 35 (0.3%) specific to V. carteri, 18 (0.2%) specific to G. pectorale, and 147 (1.3%) specific to D. salina (Figure 2B).

From the gene family function annotation results and considering the data in Figures 3A,B, all gene family functions in D. salina are primarily concentrated in cellular processes, metabolic processes, multicellular organismal processes, catalytic reactions in molecular functions, membrane in cellular components, and cell parts in cellular components within the biological process when compared to the other four species. A comparison of functional annotation of specific gene families in D. salina, V. carteri, and G. pectoral was conducted. Since both D. salina and C. eustigma are extreme organisms, with D. salina being a halophilic microorganism and C. eustigma, an acidotolerant microorganism, we aimed to explore the function of gene families shared by extreme organisms in Volvocales. To achieve this, we performed GO functional enrichment. During biochemical processing, the highest number of gene families were enriched by cellular processes, single-organism processes, and metabolic processes. This was followed by response to stimulus, biological regulation, and detoxification. In comparison to molecular functional enrichment, the number of gene families enriched in molecular functions was significantly smaller, mainly including catalytic activity, binding, and antioxidant activity. Regarding cellular components, the primary focus was on cells, cell components, organelles, organelle parts, and membranes. These functions facilitate the normal growth of extreme organisms in their corresponding environments (Figures 3A–C).

We constructed a high-confidence phylogenetic tree and estimated the divergence times of the 13 species using genes extracted from a total of 1,020 single-copy families. The results of species evolutionary tree clustering showed one branch of the multicellular group, where the two extremophiles clustered together, further demonstrating that the two representative multicellular algae, G. pectorale and V. carteri, evolved from the ancestral species of C. reinhardtii, and differentiated later than D. salina and C. eustigma. The expansion and contraction of gene families allow us to observe the changes in gene families during the evolution of species. We discovered 79 expanded and 5 contracted gene families in D. salina (Figure 4A).

According to the syntenic blocks detected in C. eustigma, C. reinhardtii, D. salina, G. pectorale and V. carteri. The D. salina whole-genome duplication (WGD) took place 13 MYA (Ks ≈ 0.0809) and C. reinhardtii whole-genome duplication (WGD) that took place 19 MYA (Ks ≈ 0.1169) (Figure 4B). We construct a phylogenetic tree. In combination with R8s software, the reported fossil evidence of these three species is C. reinhardtii (2.5 million years), G. pectorale (1.5 million years) and V. carteri (65–99 million years). Volvocales C. eustigma, D. salina, and other species diverged about 250 million years ago, C. eustigma and multicellular algae diverged about 232.5 million years ago, and G. pectorale and Volvox diverged about 150 million years ago (Figure 4A). Chlorophyta and haphtophyta diverged about 1,278 million years ago. And rhodophyta diverged with the chlorophyta and haphtophyta about 1,406 million years ago.

We used synonymous substitution rates (Ks) between collinear paralogous genes to identify potential whole-genome duplication (WGD) events, based on the assumption that the number of silent substitutions per site between two homologous sequences increases in a relatively linear manner with time. A density plot of Ks values for the collinear gene pairs of D. salina and C. reinhardtii suggested a more recent genome-wide replication event between D. salina and C. reinhardtii occurred about 46 million years ago, and an earlier genome-wide replication event 122 million years ago. The D. salina genome had a large-scale gene replication event about 1.3 million years ago, while the C. reinhardtii genome had a large-scale gene replication event about 1.9 million years ago (Figure 4B).

To further investigate the mechanism of β-carotene accumulation in D. salina under high light stress, this study analyzed the differences in gene expression in the carotenoid metabolic pathway of D. salina under high light stress and normal light incubation using high-throughput transcriptome sequencing. The evaluation of the sequencing output data for the seven samples during high light stress of D. salina is shown in the Supplementary Table S6.

The percentage of clean reads was high, indicating good sequencing quality. Additionally, the GC content of the bases in all seven libraries was between 50 and 60%, which is within the normal range. The Figure 5A shows the accumulation of β-carotene in single algal cells of D. salina under high light and control conditions. The accumulation of β-carotene in the high light and control groups increased and then decreased with the increase of incubation time. The β-carotene content in the high light experimental group was significantly higher than that in the control group at day 13 (p < 0.05). Afterward, it gradually became lower than that of the control group. The results showed that high light could promote the accumulation of β-carotene in D. salina at a certain stage of culture (Figure 5A).

The annotation success rates for the four major databases are shown in Supplementary Table S7. The percentage of unigenes successfully annotated in eggNOG is 61.87%, in GO database is 34.31%, in KOG database is 61.87%, and in KEGG database is 21.30%. These percentages indicate that the gene function annotation of the transcriptome is comprehensive (Supplementary Table S7).

As shown in Supplementary Figure S2 and Supplementary Tables S7, S8, GO categories include: A total of 35,230 unigene notes were classified into cell component classification, among which the largest proportion of unigene related to cells and cell components accounted for 22.7%, and the unigene related to organs and organ components accounted for 29.4%. These categories make up the majority of cell classification notes. A total of 32,412 unigene annotations were classified into the classification of biological processes, of which 12.6% were related to metabolic processes, 14.2% were related to cellular biological processes, and 13.5% were related to single biological processes. These categories make up the majority of cell metabolism annotations. A total of 628 unigene annotations were classified into the classification of molecular processes, among which binding related unigene accounted for 29.6% and catalytic activity related unigene accounted for 43.9%, which accounted for the majority of molecular process annotations. In summary, most of the genes in the classification of cell components were related to organelles, cell membranes and macromolecular complexes. In the process of biological metabolism, most genes are related to stress response, signal transduction pathway and metabolic response. In terms of molecular function, most genes are related to transport, catalysis and binding activities (Supplementary Figure S3 and Supplementary Tables S7, S8).

KOG was divided into 26 annotation groups, and the genes successfully annotated by KOG were classified according to the KOG group. The results were shown in Supplementary Table S9, as shown in Supplementary Figure S3. A total of 7,955 unigene were annotated in the KOG classification. Among them, the highest gene annotation included general function prediction (13%), signal transduction mechanism (10%), and post-translational modification (11%). 274 genes related to lipid transport and metabolism were found in KOG annotation. The above data indicated that in the early stage of stress, a large number of algal cells carried out biological processes such as transcriptional activation and protein translation in response to stress, accompanied by lipid anabolism (Supplementary Figure S4 and Supplementary Table S9).

DEGs between two adjacent time points (AVsB) of stress treatment group and control group were detected by DEGseq software, and DEGs were divided into up-regulated and down-regulated genes. The statistical results of DEGs in the high-light treatment group and control group are shown in Supplementary Table S10 and Figures 5C,D. The results showed that the gene expression of D. salina under high light stress was different from that of normal culture. The most obvious difference was that in the control group, 6,142 up-regulated genes were lower than 0dVs5d (6250), while 3,949 down-regulated genes were higher than 0dVs5d (2068) in 5dVs13d differential genes (Figures 5C,D; Supplementary Table S10).

In the high-light treatment group, the number of up-regulated genes in 5d-HLVs13d-HL was 8,466, higher than that of 0dVs5d-HL (6123) in the previous stage, while the number of down-regulated genes was 2,969, lower than that of 0dVs5d-HL (4886). At the same time, it was also observed that compared with the previous stage of 13dVs20d and 13d-HLVs20d-HL, the number of up-regulated genes decreased and the number of down-regulated genes increased, and the number of down-regulated genes was more than the number of up-regulated genes. According to the data analysis, the transcriptome of D. salina underwent drastic changes in the early stage, and the transcriptome was in an activated state, and then showed a downward trend, and the number of down-regulated genes in each differentially expressed gene set was lower than that of up-expressed genes, except for 13-20d. The number of DEG in the treatment group was higher than that in the control group, and the number of genes up-regulated by 5d-HLVs13d-HL was relatively higher, and the number of genes up-regulated by 20d-HL was at least 2,101. It can be seen from the above that the number of differential genes before and after 5 to 13 days of culture was the highest in both high-light stress and normal culture of D. salina, and the number of up-regulated genes was higher than that of down-regulated genes (Supplementary Table S10).

Through the enrichment of differentially expressed genes in the KEGG pathway, 52 genes in the metabolic pathway were annotated in the biosynthesis of secondary metabolites, 325 genes were annotated in lipid metabolism, and 87 genes were annotated in the metabolism of terpenes and polykeones. KEGG enrichment of differential genes is shown in Figure 5E. The predicted up-regulated and down-regulated gene enrichment numbers of KEGG pathway are shown in Supplementary Tables S11, S12. It notes that the predicted up-regulated genes of carotenoid metabolic pathway include DsCrtB, DsPDS, DsZ-ISO, DsZDS, DsCRTISO, DsLUT5, DsCrtL-B and DsCCD8, while the predicted down-regulated genes include DsCrtF and DsLUT1. The four genes that were both up-regulated and down-regulated were DsZEP, DsCrtR-b, DsCruA/P and DsCrtZ 4 (Supplementary Table S13).

Through the enrichment of 5d Vs13d differentially expressed genes in KEGG pathway, 40 genes in metabolic pathway were annotated in secondary metabolite biosynthesis, 220 genes were annotated in lipid metabolism, and 70 genes were annotated in terpenoid and polyketo metabolism. The down-regulation of carotenoid metabolic pathway prediction genes is shown in Supplementary Table S14. There are 12 down-regulated genes, namely DsCrtB, DsPDS, DsZ-ISO, DsZDS, DsCrtISO, DsLUT5, DsCrtF, DsLUT1, DsZEP, DsCrtR-b, DsCruA/P and DsZEP. DsCrtL-b and DsCCD8 are both upward and downward forecast.

KEGG enrichment of differential genes is shown in Figure 5F. The predicted up-regulated and down-regulated gene enrichment numbers of KEGG pathway are shown in Supplementary Tables S15, S16. It can be seen from the results that compared with the results of 5d-HLVs13d-HL, the genes enriched by 5D-VS13D KEGG differential genes in the carotenoid metabolic pathway were mostly down-regulated genes, which may be closely related to the culture under high light stress condition. Moreover, on the 13th day of culture, the β-carotene content in the high light group was significantly higher than that in the control group (Supplementary Tables S15, S16). It was further indicated that the high light stress could promote the synthesis of carotene in the salt algae.

When comparing data from 5 days versus 13 days, the genes expressed differently in the carotenoid metabolism pathway during the 5-day high light versus 13-day high light comparison were predominantly up-regulated, leading to increased β-carotene levels in individual cells compared to the control. This further supports the conclusion that high light stress enhances β-carotene synthesis. As shown in the Figure 5B, the expression of DsPDS and DsZDS genes was higher on day 13 under high light treatment than on day 5. The expression of both DsPDS and DsZDS genes increased on day 13 compared to day 5, which was consistent with the predicted results. The β-carotene accumulation of individual cells in the high light experimental group and the control group was compared on day 5 and day 13. The β-carotene content of single cells was higher in the high light group than in the control group on day 13, which was consistent with the statistical results of differentially expressed genes. Since the DsLCYB gene is one of the key enzymes for β-carotene synthesis, and there was no DsLCYB gene among the predicted up-and down-regulated genes, the expression of the DsLCYB gene was verified. The results showed that there was no significant change in the expression of the DsLCYB gene in the samples from the two time periods, which was consistent with the predicted results (Figure 5B).