We obtained seven Sphaeropleales strains from National Aquatic Biological Resource Center, Institute of Hydrobiology, Chinese Academy of Sciences. The lists were Bracteacoccus aerius (FACHB-895), Bracteacoccus engadinensis (FACHB-1300), four Tetradesmus obliquus strains (FACHB-14, FACHB-276, FACHB-417 and FACHB-1235), Monoraphidium contortum (FACHB-3677). Bracteacoccus aerius and Bracteacoccus engadinensis were terrestrial, and others were collected from freshwater. All these strains were grown at 25 °C in liquid BG11 medium under a 12/12-h light/dark cycle.

DNA was extracted using a Universal DNA Isolation Kit (Axygen, Suzhou, China). A NEB Next Ultra DNA Library Prep Kit for Illumina (New England Biolabs, Ipswich, Massachusetts, USA) was used for preparing sequencing libraries which were sequenced on an DNBSEQ platform. The quality of the raw reads was initially assessed using FastQC v0.11.6 (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/ accessed on 24 June 2021). Data were trimmed using SOAPnuke software (Chen et al., 2017) and were then assembled by megahit v 1.0.3 (Li et al., 2015) with defaulted parameters. Anvio 7.0 (Eren et al., 2015) was used to remove prokaryotic sequences from the assembled genomes with defaulted parameters. Firstly, Bowtie2 (Langmead and Salzberg, 2012) was used to align the assembled genomes to the fq files, generating BAM files. Then, the ‘anvi-gen-contigs-database’ was employed to generate a contig database, followed by the use of ‘anvi-run-hmms’ to identify single-copy genes in the contig database. The ‘anvi-profile’ was used to import sample information into the database. Next, ‘anvi-interactive’ was used to visualize the results. Finally, ‘anvi-export-collection’ was used to export the eukaryotic bins, and ‘anvi-export-contigs’ was used to extract the sequences of the aimed bins, resulting in the final assembled genomes. The genomes of the seven Sphaeropleales trains were deposited in CNCB with accession numbers SAMC4339295 - SAMC4339301, respectively.

Ab initio gene prediction was performed by AUGUSTUS version 3.2.1 (Stanke et al., 2006) with four three-trained organism models: volvox, chlamy2011, and chlorella. GeMoMa-1.9 (Keilwagen et al., 2018) was employed for homolog-based gene prediction and integration of gene prediction results, with genomes and protein of Volvox carteri f. nagariensis (GCA_000143455.1), Gonium pectorale (GCA_001584585.1), Haematococcus lacustris (GCA_030144725.1), and Tetradesmus obliquus (GCA_030272155.1) from GenBank as the references. Functional annotation of predicted proteins was performed by InterProScan (Jones et al., 2014) and EggNOG-mapper (Huerta-Cepas et al., 2017) for GO term mapping and KEGG pathway analyses respectively. To assess the completeness of the genome assembly and annotation, BUSCO (Simao et al., 2015) (Benchmarking Universal Single-Copy Orthologs) was used, utilizing the chlorophyta_odb10 database for a quantitative evaluation.

Orthologous gene groups for the Sphaeropleales samples (including those newly added in this study and those available in the public database with annotation completeness greater than 80%) were estimated using Orthofinder (Emms and Kelly, 2015) with default parameters to identify single-copy or low-copy orthologues. Additionally, V. carteri f. nagariensis (GCA_000143455.1) and D. salina (GCA_002284615.2) were included as outgroups to estimate shared orthologues for subsequent divergence time estimation and gene family evolution analysis. Transporters were identified using BLASTp, with the Transporter Classification Database (TCDB) as a query, with an e-value cut-off of <1E⁻⁵. Further, repeat content for each genome was determined by RepeatModeler and RepeatMasker with default parameters. A library of repeats was first created for each assembly using RepeatModeler (Version 2.0.1) and all repeats were masked using RepeatMasker (Version 4.1.0).

The predicted protein sequences (PEP) sequences of each single-copy or low-copy orthologue were aligned using MAFFT v7.394 (Katoh and Standley, 2013) with the parameters -maxiterate 1000 and -globalpair. Regions with poor alignment were trimmed using TrimAl v1.2 (Capella-Gutierrez et al., 2009) with the -automated1 option. The resulting trimmed alignments of orthologous groups were then used for subsequent phylogenomic analysis. Coalescent-based analyses were employed to construct the phylogenetic tree. For these analyses, RAxML (Stamatakis, 2014) was used to perform maximum likelihood (ML) analysis of each single-copy orthologue, applying the PROTGAMMA GTR model. Additionally, ASTRAL (Zhang et al., 2018) was used to infer the coalescent-based species tree (ST) phylogeny.

The 18S rDNA sequences were aligned using MAFFT v7.0 (Katoh and Standley, 2013), and ambiguous regions were manually adjusted and refined using MEGA7 (Kumar et al., 2016). To locate the positions of 18S rDNA sequences within the sequenced genomes, nhmmer (Wheeler and Eddy, 2013) was used to align the final assembled genomes to the rDNA database, based on rDNA sequences downloaded from NCBI. Subsequently, SSU sequences were extracted using Perl scripts. To understand the phylogenetic relationships of the strains in this study, additional Sphaeropleales 18S rDNA sequences were downloaded from the public database. The 18S rDNA sequences were analyzed using jModelTest2 (Darriba et al., 2012) to select the best-fit model, which was found to be GTR + I + G. Bayesian inference (BI) method was applied to infer the phylogeny.

PhyloSuite (Zhang et al., 2020) was used to concatenate all the shared single-copy orthologous groups from the Sphaeropleales and the outgroup, resulting in a concatenated sequence of 35,684 amino acids. This concatenated alignment, along with the species tree constructed by ASTRAL (Zhang et al., 2018), was used for Bayesian divergence time estimation using the approximate likelihood method described by dos Reis & Yang, implemented in MCMCtree v. 4.9 (Yang, 2007).

The number of gene copies per family, as determined by Orthofinder (Emms and Kelly, 2015), and the timetree estimated earlier by MCMCtree (Yang, 2007). were used to analyze gene family expansion and contraction using CAFE v. 5.1 (Mendes et al., 2020). Expanded and contracted genes were then extracted for Gene Ontology (GO) functional and KEGG pathway enrichment analyses. For GO enrichment, all PEP sequences were imported into InterProScan (Jones et al., 2014) for GO term mapping. The analysis was carried out using the clusterProfiler (Yu et al., 2012) package, with a significance cutoff of p < 0.05, and the false discovery rate (FDR) method was applied to adjust for multiple testing (Benjamini and Hochberg, 1995).

The CODEML program of PAML v4.9 (Yang, 2007) was used to estimate positive selection and rapidly evolving genes based on common orthologues were described as Xiong et al (Xiong et al., 2022, 2021). The branch model was employed in the calculation of dN/dS ratios for terrestrial Sphaeropleales species and aquatic ones with the two Bracteacoccus species labeled as foreground branches. A null model (model = 0), where one dN/dS ratio was fixed across all strains, was compared with an alternative model (model = 2), where Bracteacoccus species were allowed to have a different dN/dS ratio. Likelihood ratio tests were performed to examine model fit, a chi-squared test was used to analyze p values, and multiple testing was corrected using false discovery rate (FDR) in R. The genes were considered putative rapidly evolving genes if they had an FDR-adjusted p value < 0.05 and a higher dN/dS ratio in the foreground branch than in the background branches.

Branch-site model was utilized to find genes that possibly underwent positive selection. The improved branch-site model (model = 2, Nsites = 0) was used to detect signatures of positive selection on individual amino acids in a specific branch. The two Bracteacoccus species were set as the foreground branch. The null model assumed that no positive selection occurred on the foreground branch (fix_omega = 1, omega = 1), and the alternative model assumed that sites on the foreground branch were under positive selection (fix_omega = 0, omega = 1.5). LRT were used to test model fit and Chi-square test was applied for examining the P values in R. A correction was performed for multiple testing using an FDR criterion, and BEB method was employed to statistically identify sites under positive selection. Genes with an FDR-adjusted P < 0.05 were considered as putatively selected. For the genes belonging to both under positive selection and putative rapidly evolving were performed for gene ontology (GO) functional enrichment analyses as above.

The seven newly added Sphaeropleales strains were collected prior to 2000 and were deposited in National Aquatic Biological Resource Center, Institute of Hydrobiology, Chinese Academy of Sciences. Phylogenetic analysis based on 18S rDNA sequences extracted from the genomes was performed, and the results are shown in Supplementary Figure S1 . The phylogenetic tree with high support values at the basal nodes confirmed species identification. We also performed a phylogenetic analysis using low-copy orthologues across the entire Chlorophyta, including core Chlorophyta and prasinophytes by Astral ( Figure 1 ), which enables highly accurate phylogenomic estimation, even in the presence of high levels of gene tree conflict because of incomplete lineage sorting (Mirarab et al., 2014) or horizontal gene transfer (Davidson et al., 2015). The results were almost identical with the previous studies, supporting the phylogenetic relationships within core Chlorophyta and between prasinophytes. Notably, our analysis confirmed that P. coloniale (CCMP 1413) was distinct from both core Chlorophyta and prasinophytes, suggesting it belongs to a newly identified phylum, Prasinodermophyta (Li et al., 2020).

Using Illumina platforms, we sequenced and assembled draft genomes for B. aerius, B. engadinensis, M. contortum, S. obliquus, and three T. obliquus strains. These genomes, along with seven previously sequenced Sphaeropleales genomes, were summarized in Table 1 . The contig N50 values of the seven newly sequenced genomes were nearly identical, ranging from 20.5 to 22.34 kb. The BUSCO completeness for all genomes exceeded 90%, indicating that the genome assemblies were sufficiently complete for comparative analyses of genomic and gene constituents. The genome sizes of the seven sequenced species ranged from 45.62 Mb to 109.8 Mb, with the GC content varied from 55.72% to 57.63%. Gene counts varied significantly among species, ranging from 6,874 to 13,568. The repeat composition covered 3.06% to 18,88% of the genome.

Among all the Sphaeropleales species, the genomes size of Scenedesmus species exhibited significant variation, ranging from 39.8 Mb to 151.90 Mb, representing both the smallest and largest genomes within this group. In contrast, Monoraphidium, Raphidocelis, and Bracteacoccus displayed moderate genome sizes, ranging from 46.6 Mb to 69.5 Mb. And Tetradesmus species possessed relatively larger genomes that exceed 100 Mb. Most species have a GC content around 56%, exception for Monoraphidium and Scenedesmus sp. PABB004 greater than 70%. The repeat content varied considerably among all genomes examined, Scenedesmus sp. PABB004 and R. subcapitata showed the lowest repeat content about 1.3%, and T. obliquus (SAMC4339298) displayed the highest at 18.8%. Bracteacoccus and Monoraphidium have moderate repeat content, ranging from 3% to 4%, whereas Tetradesmus and most Scenedesmus species exhibit higher ranging 6.78% to 18.88%. minutum,

It is generally accepted that there is some degree of correlation between the genome size, the proportion of repetitive sequence, and the number of genes (Feng et al., 2017; Hou and Lin, 2009; Nishitsuji et al., 2020; Wang et al., 2021). These observations support the idea that larger genomes generally having more gene count, fewer repeats and lower GC content. While, it seemed larger genomes exhibited more gene count at the genus level of Sphaeropleales, the relationship the genomes size among the GC and repeat content showed no obvious pattern, more high-quality genomes will contribute to explore the relationships.

We conducted orthologous gene family analysis based on the newly sequenced seven Sphaeropleales strains along with seven high-quality previously sequenced genomes at different taxonomic levels ( Figure 2 ). The orthologous analysis revealed that 3,761 gene families were shared or conserved across all fifteen genomes ( Figure 2A ). Additionally, each strain contained unique gene families, comprising 2-3.5% of the total gene count: 2,009 families were unique to T. obliquus, 554 to S. obliquus, 411 to B. aerius, 257 to Scenedesmus sp. NREL 46B-D3, 247 to B. engadinensis, 151 to M. minutum, 144 to R. subcapitata, and 138 to M. contortum. At the genus level ( Figure 2B ), 6,733 gene families were shared or conserved across the five genera. Furthermore, 5,587 gene families were unique to Tetradesmus, 2,391 to Bracteacoccus, 1,432 to Monoraphidium, and 720 to Scenedesmus. At the family level ( Figure 2C ), 7,304 gene families were shared or conserved among the three families. Additionally, 10,656 gene families were unique to Scenedesmaceae, 3,149 to Selenastraceae, and 2,667 to Bracteacoccaceae. In whole, the Tetradesmus species, in particular, had the most unique gene families (5,587), likely due to their larger genomes.

To further investigate gene functions within the Sphaeropleales, we classified all annotated proteins from 14 genomes into functional categories using the Gene Ontology (GO) database ( Supplementary Figure S2 ). The predicted proteins were categorized into three primary GO domains: molecular function (MF), cellular component (CC), and biological process (BP). We compared the contents of each category and their corresponding percentages. Among the 11 strains, unique GO domains were identified, with T. obliquus FACHB-276 exhibiting the highest diversity, containing 278 distinct domains across the three GO categories. For the top ten functional categories at the CC and MF levels, most species shared nearly identical terms. We also calculated the PFAM domain categories and their percentages ( Supplementary Figure S3 ) and KEGG pathway categories ( Supplementary Figure S4 ) for all species, which showed that all the species shared similar top-ten categories and associated percentages.

Based on the Transporter Classification Database, we analyzed the transporters in the Sphaeropleales genomes. A total of 320 transporters were identified, as shown in Supplementary Table S1 . The transporters across the Sphaeropleales species exhibited similar patterns in terms of the top ten transporter categories and their respective numbers. As additional species were included in the analysis, the results were consistent with previous studies (Suzuki et al., 2018), which showed that Sphaeropleales species possess notably higher numbers of genes encoding for H+/hexose transporters (2.A.1.1, TCAD ID), amino acid permeases (2.A.18.2), peptide transporters (2.A.17.3), aquaporins (1.A.8.8), and metal-nicotianamine transporters (2.A.67.2) compared to Chlamydomonas reinhardtii ( Table 2 ). The Sphaeropleales are a dominant group of freshwater algae, well-adapted to diverse environmental conditions (Baldisserotto et al., 2012; Fawley et al., 2004). These species also exhibit high sensitivity to exogenous substances (Nagai et al., 2016). Previous studies have suggested that Sphaeropleales have the ability to adapt to various environmental conditions for possessing a significantly greater number of genes related to H+/hexose transporters, amino acid permeases, peptide transporters, aquaporins, and metal-nicotianamine transporters comparing to C. reinhardtii (Suzuki et al., 2018). Aquaporins help algae adapt to high salt stress by facilitating the transport of small polar molecules, such as water, across cell membranes, thereby regulating intracellular osmotic pressure (Anderberg et al., 2011). Additionally, the presence of multiple metal-nicotianamine transporters, ABC transporters, and genes involved in heavy metal ion and xenobiotic transport ( Supplementary Table S2 ) suggested that Sphaeropleales may have a high sensitivity to metals, positioning them as potential phytoremediation organisms for removing heavy metal pollution from aquatic environments (Benderliev and Ivanova, 1994, 1996; Murata et al., 2006; Schaaf et al., 2004). Genes related to H+/hexose cotransport, amino acid/peptide transporters, and nitrate/nitrite transporters are likely key to their rapid growth under varying nutrient conditions (Cho et al., 1981; Sauer et al., 1983).

Furthermore, Sphaeropleales species possessed a gene encoding pyrroline-5-carboxylate reductase, which synthesizes proline to alleviate osmotic stress under cold conditions (Liu et al., 2020). And pyrroline-5-carboxylate reductase is also associated with halotolerance (Arora et al., 2019; Kishor et al., 1995; Pérez-Arellano et al., 2010). Previous study indicated that copies of the pyrroline-5-carboxylate reductase gene of the halotolerant microalga Scenedesmus sp. NREL 46B-D3 were upregulated (Calhoun et al., 2021) in the cold stress, and pyrroline 5-carboxylate reductase was upregulated under the salt stress in the cold tolerant M. minutum 26B-AM and S. obliquus UTEX393 (Calhoun et al., 2022). They also exhibited genes encoding antioxidant enzymes (e.g., catalase, ascorbate peroxidase) and antioxidant biosynthesis pathways (e.g., glutamate, β-carotene) ( Supplementary Table S3 ), which helped mitigate oxidative stress from excess reactive oxygen species (O₂ and H₂O₂) under cold stress (Fryer et al., 1998; Liu et al., 2020; Prasad, 1997; Van Breusegem et al., 1999). The findings from this study, based on a larger set of sequenced Sphaeropleales genomes, supported with adaptation for different environmental condition.

All Sphaeropleales species showed the presence of genes involved in vitamin B biosynthetic pathways, including the vitamin B6 biosynthetic process, thiamine (vitamin B1) biosynthesis and metabolism, biotin (vitamin B7) synthase activity, and biotin biosynthesis ( Supplementary Table S4 ), indicating the ability to synthesize these vitamins de novo, similar to C. reinhardtii, Cyanidioschyzon merolae, and Thalassiosira pseudonana (Croft et al., 2006). A previous study on the marine diatom Skeletonema costatum has shown that a mixture of vitamin B compounds plays a crucial role in mitigating the harmful effects of hypersalinity (Gede et al., 2017). Furthermore, the possession of a complete pathway for thiamine biosynthesis contributes to enhanced biotic and stress resistance (Almutairi et al., 2021). The biosynthesis of the essential amino acid methionine can occur via both B12-dependent and B12-independent isoforms of methionine synthase (MetH and MetE, respectively) (Croft et al., 2005, 2006; Helliwell et al., 2011). In this study, all the Sphaeropleales species showed no genes about cobalamin (vitamin B12) biosynthesis, while exhibited genes about cobalamin (vitamin B12) metabolic process, and cobalamin binding ( Supplementary Table S4 ), indicating that methionine synthesis occurred solely via the VB12-independent pathway, namely MetE isoform in Sphaeropleales. A previous survey of 306 species aimed to determine whether these algae require vitamin B12 (Croft et al., 2006). The results indicated that none of the Sphaeropleales species showed a requirement for cobalamin, which acts as a cofactor for enzymes involved in rearrangement-reduction reactions and methyl transfer reactions. Vitamin B12 (VB12) can only be produced by bacteria (both eubacteria and archaea) in nature, and its concentration in the natural environment is typically lower than required in culture (Croft et al., 2006). Therefore, these species acquire VB12 or its precursors through a symbiotic relationship with bacteria. Such symbiotic interactions between bacteria and algae are widespread, as many algae species are capable of acquiring vitamin B12 from their bacterial partners (Croft et al., 2006; Daisley, 1969).

All species contained genes involved in the assembly, movement, and organization of cilia ( Supplementary Table S5 ), as well as genes associated with meiosis ( Supplementary Table S6 ). In the Sphaeropleales cell cycle, the stage with motility flagellates/) or meiosis are either not dominant or not well understood (Trainor and Burg, 1965; Yamagishi et al., 2017). Suzuki et al. proposed that immobility may force cells to adapt to different environmental conditions aided by their numerous transporters (Suzuki et al., 2018). Another possibility that a cryptic sexual cycle or previously unobserved motile life cycle stage with flagella may exist in these organisms, or these genes are nonfunctional in Sphaeropleales. This phenomenon has been reported in C. zofingiensi (Blanc et al., 2010; Roth et al., 2017).

To gain a comprehensive understanding of the orthogroups differences between species, particularly those with assignable functions, we performed gene family expansion and contraction analyses using CAFÉ, based on 29,073 gene families. As shown in Figure 3 , all species and most ancestral nodes exhibited a substantial number of gene family expansions, with contractions being more prevalent. Specifically, B. engadinensis and B. aerius exhibited the highest number of contracted gene families, with 17,786 and 17,499 respectively, surpassing all other nodes, while T. obliquus (FACHB-276) displayed the most expanded gene families, with 2,142. In terms of ancestral nodes, the common ancestor of the Selenastraceae family had the largest number of expansions and contractions, with 10,832 contracted gene families and 471 expanded gene families. In contrast, the common ancestor of the Bracteacoccaceae family showed the smallest number, with only 11 expanded gene families and no contractions. The number of expanded or contracted gene families can be influenced by various factors, including gene duplication, de novo gene creation, gene loss, the functional roles of gene families, and environmental changes (Albalat and Canestro, 2016; Guo, 2013; Prachumwat and Li, 2008). Consequently, the observed differences in gene family expansions and contractions in Sphaeropleales are likely due to multiple causes. A previous study suggested that gene family contraction could contribute to genome reduction (Qiu et al., 2015). This may help explain why Scenedesmus sp. PABB004exhibited the largest number of gene family contractions, while simultaneously having the smallest genome size within the Scenedesmaceae, and the Bracteacoccus species showed relatively small genome size and almost the fewest gene count.

To investigate whether gene family evolution correlates with habitat adaptation, we performed Gene Ontology (GO) enrichment analysis based on the expansions and contractions in the common ancestor of Bracteacoccaceae. The GO enrichment analysis revealed that the expanded families were enriched in six GO terms ( Figure 4A ), primarily related to methionine biosynthetic, cobalamin binding, tRNA Modification, pentosyltransferase activity, and metal ion binding.

The branch model of PAML was used to compare dN/dS ratios between terrestrial and aquatic Sphaeropleales species based on 3757 common orthologues. Among these, 3032 orthogroups showed significantly higher dN/dS ratios in the two terrestrial Bracteacoccaceae species, indicating the occurrence of rapid evolution ( Supplementary Table S7 ). Positive selection analysis was performed based on branch-site model, and the null and alternative models were compared. The null model considered that the foreground branch only has dN/dS = 1, and the alternative model assumed that sites on the foreground branch have dN/dS > 1 (positive selection). When the two terrestrial Bracteacoccaceae species were labelled as the foreground branch, 13 orthogroups showed the FDR-adjusted P value less than 0.05, which indicated that these 13 orthogroups may possibly under positive selection ( Supplementary Table S8 ), and they also under rapid evolution based on the result of Branch model. GO enrichment analysis showed these orthogroups were enriched in 87 GO terms ( Figure 4B ), among which 55 GO terms belonged to biological process category, 22 belonged to cellular component category and 10 belonged to molecular function category. Among all the three categories, there were some common functions such as oxidation-reduction process, the function of mitochondrial, the biosynthetic and metabolism of starch, polysaccharide and other organics ( Figure 4B ).

Bracteacoccaceae included only one genus Bracteacoccus, which are coccoid green algae that occurs in a wide range of soil types worldwide, spanning climates from polar to tropical and not avoiding even heavily polluted localities (Patova and Dorokhova, 2008; Ronquist et al., 2012; Stibal et al., 2006). Gene Ontology (GO) enrichment analysis, based on the expansions and contractions in the ancestral lineage of Bracteacoccaceae, revealed that the enriched GO terms were primarily associated with methionine biosynthesis, cobalamin binding, tRNA modification, pentosyltransferase activity, and metal ion binding. A previous study about C. reinhardtii indicated that methionine biosynthesis is an essential cellular mechanism for adaptation to thermal stress (Xie et al., 2013). Additionally, Zeng et al. identified key pathways in regulating natural variations in phenylpropanoid content, including flavone C-pentosyltransferase proteins, which are involved in UVB protection in Qingke (Tibetan hulless barley) (Zeng et al., 2020). It has also been reported that tRNA modifications is correlated with cold temperatures, drought, increased soil salinity, and developmental stages in vascular plants (Mmbando, 2024). Furthermore, the evolutionary analysis based on PAML showed the two species of Bracteacoccus exhibited genes experienced both rapidly evolution and positively selected. Among these, common functional categories included oxidation-reduction processes, mitochondrial functions, and the biosynthesis and metabolism of starch, polysaccharides, and other organic compounds. We considered these expanded gene families, rapidly evolution and positively selected were likely related to the adaption to the terrestrial habitat for Bracteacoccaceae.