A Reinvestigation of Multiple Independent Evolution and Triassic–Jurassic Origins of Multicellular Volvocine Algae

Keywords: volvocine algae, phylogeny, incomplete lineage sorting, divergence time, multicellularity

Phylogenetic Relationships of the Volvocales and Incertae Sedis Taxa in SV Clade

Our molecular data set included 1,362 nuclear genes from 133 Chlorophyceae species, with 5 OCC clade species as outgroups. We largely expanded taxon sampling within the Volvocales (110 species, including 81 Reinhardtinia clade species; supplementary table S1, Supplementary Material online). To reduce the bias generated by genes with fast evolutionary rates, we used 3 subsets with slow evolutionary rates of genes (1,267 genes, 1,152 genes, and 764 genes) and reconstructed the phylogenetic relationships of the Volvocales by coalescent and concatenation approaches. To better resolve the deep relationships in the Volvocales, we used different models in the concatenation-based analyses (i.e., partition strategy, C20 model, and GHOST model). These complex models have been regarded as approximating the evolutionary process better by accounting for substitution rate variation across genes, sites, and lineages. Our phylogenomic analyses showed that all the species trees were largely congruent and most phylogenetic relationships among the major clades of the Volvocales and Sphaeropleales were well-supported. The Volvocales and its main clades were recovered as monophyletic groups with full support (fig. 1 and supplementary figs. S1–S15, Supplementary Material online), which were consistent with One Thousand Plant Transcriptomes Initiative (2019) and chloroplast phylogenomic analyses (e.g., Lemieux et al. 2015; Fučíková et al. 2019; Li et al. 2023).

The Crucicarteria, Hafniomonas, and Radicarteria were recovered as the earliest diverging lineages of the Volvocales. The remaining clades of the Volvocales formed a monophyletic group, of which the Chloromonadinia was the earliest divergent clade, and Caudivolvoxa + Xenovolvoxa were recovered as the sister groups to the Reinhardtinia + Oogamochlamydinia (fig. 1 and supplementary figs. S1–S15, Supplementary Material online). However, the relative position of these earliest diverging lineages differed in various analyses. In the coalescent-based analysis, Crucicarteria is the earliest diverging lineage, followed by Hafniomonas and Radicarteria (fig. 1 and supplementary figs. S1–S3, Supplementary Material online). When using site-heterogeneous C20 model and partition strategy based on the 764 genes, Hafniomonas and Radicarteria were inferred as sister groups (supplementary figs. S4–S8, Supplementary Material online). Based on the other 3 data set of partition strategies (1,152 genes, 1,267 genes, and 1,362 genes) and GHOST model, these 3 clades clustered together as the earliest divergent clade, with Hafniomonas being the sister of Radicarteria + Crucicarteria (supplementary figs. S9–S15, Supplementary Material online). The phylogenetic conflicts in the concatenation-based analyses may imply that our data set was sensitive to the rate variation across sites at deep divergences in the Volvocales. Evaluating model-to-data fitness is an important but challenging task in phylogenetics. Methods for assessing the fitness of C20 or GHOST models are not available so far, and model fitness evaluation still needs further investigation. In addition, restricted taxon sampling could affect accuracy of phylogenetic analysis and reduce branch support (Massoni et al. 2014; Zhong and Betancur-R 2017). These uncertain relationships were also likely caused by insufficient taxon sampling. Appropriate and extensive taxon sampling is one of the most practical ways to improve the accuracy of phylogenetic inference (Heath et al. 2008; Massoni et al. 2014). Thus, intensive samplings of deep lineage and better-fitting evolutionary models are needed to improve the resolution of the phylogenetic relationship of Volvocales.

The phylogenetic distinctness of the Sphaeropleales and Volvocales appears to fade when including deeply diverging incertae sedis taxa (Lemieux et al. 2015; Fučíková et al. 2019; Li et al. 2023). Therefore, it is important to clarify the positional relationship of the incertae sedis taxa. Four incertae sedis taxa (Golenkinia, Microsporaceae, Treubarinia, and Spermatozopsis) defined in Fučíková et al. (2019) were encompassed in our analyses. Our results indicated that the incertae sedis taxa and Sphaeropleales grouped into a monophyletic lineage in the concatenation-based analysis with C20 model using 764 genes and all coalescent-based analyses, named Sphaeropleales sensu lato (Sphaeropleales s.l.), which was fully supported as sister to the Volvocales (fig. 1 and supplementary figs. S1–S4, Supplementary Material online). Conversely, in the other concatenation-based analyses, the Spermatozopsis was sister to the Sphaeropleales, and the remaining incertae sedis taxa (Golenkinia, Microsporaceae, and Treubarinia) were recovered as sister groups to the SV clade (supplementary figs. S5–S15, Supplementary Material online). Together, our results suggested that the affinities of the incertae sedis taxa were more closely related to the Sphaeropleales or SV clade than the Volvocales.

The Well-Supported Nonmonophyletic Colonial Volvocine Algae

All phylogenomic analyses recovered full support for the nonmonophyletic colonial volvocine algae. The Vitreochlamys ordinata (Chlamydomonadaceae) was sister to the Tetrabaenaceae, and they together diverged earliest in the colonial volvocine algae. The unicellular Chlamydomonas (e.g., Chlamydomonas reinhardtii, Chlamydomonas incerta, and Chlamydomonas globosa) was resolved as the sister clade of the Goniaceae + Volvocaceae (fig. 1 and supplementary figs. S1–S15, Supplementary Material online). Interestingly, the phylogenetic relationships reconstructed using slow-evolving genes could improve support for key nodes in the coalescent-based approach. In particular, the use of 764 genes with slow evolutionary rates improved support for focus node of colonizing TGV clades from 1/70 to 1/91 (fig. 1 and supplementary figs. S1–S3, Supplementary Material online). These results imply that fast evolutionary genes may introduce bias into phylogenetic reconstruction of ancient groups.

Kirk (2005) suggested that the volvocine algae followed an evolutionary process of a simple, progressive increase in size and complexity. From an ancestral unicellular to multicellular, Kirk (2005) identified out of 12 developmental changes, none of which have obvious analogs in the unicellular C. reinhardtii (Umen 2020). Thus, in agreement with Lindsey et al. (2021), our phylogenetic results overturned previous knowledge of the monophyletic origin of the colonial TGV clade (supplementary table S2, Supplementary Material online) and suggested that multicellularity independently evolved at least twice in the volvocine algae: one originated in the Tetrabaenaceae and another in the ancestor of Goniaceae and Volvocaceae (fig. 1). The unicellular taxa (e.g., C. incerta, Chlamydomonas debaryana, Chlamydomonas schloesseri, and V. ordinata) were found nested within the colonial clade. The morphological study has shown the close relationship between Vitreochlamys and the colonial volvocine algae, such as the expanded middle layer of the cell wall in Vitreochlamys (Nakazawa et al. 2001). The transformation of the cell wall to the extracellular matrix (ECM) is an important characteristic in the colonial algae (Kirk 2005). Herron and Michod (2008) indicated that the expanded gelatinous middle layer of the cell wall in Vitreochlamys might be an early step in the conversion of the cell wall into ECM. These results implied that unicellular relatives are critical for understanding the early steps of multicellular evolution, especially the unicellular Chlamydomonas and Vitreochlamys, which are the closest relatives of colonial volvocine algae.

The monophyly of the Goniaceae has been widely discussed (Nozaki et al. 2000; Herron and Michod 2008; Herron et al. 2009; Nakada et al. 2010; Nozaki et al. 2014; Hu et al. 2019, 2020; Lindsey et al. 2021). In our study, Goniaceae was recovered as monophyletic in the coalescent-based analysis (fig. 1 and supplementary figs. S1–S3, Supplementary Material online). Astrephomene of Goniaceae was sister to the Volvocaceae in the concatenation-based analyses with C20 model using 1,362 genes (supplementary fig. S7, Supplementary Material online), whereas Astrephomene was sister to Volvulina or Volvox that embedded in the Volvocaceae using the partition strategy, GHOST model, or other 3 subsets of C20 model (supplementary figs. S4–S6 and S8–S15, Supplementary Material online). Previous studies supported that the Goniaceae was not monophyletic using coalescent and concatenation approaches, but these studies had limited sampling of genes and species (Yamashita et al. 2021; Lindsey et al. 2021, supplementary table S2, Supplementary Material online). The use of multiple genes has indisputably increased the resolution of phylogenetic inference (Delsuc et al. 2006; Roure et al. 2013). Furthermore, Nute et al. (2018) showed that a large number of genes are necessary to achieve a high degree of accuracy in reconstructing phylogenetic relationships based on the coalescent approach. Lindsey et al. (2021) used 40 single genes to reconstruct the phylogeny, and a small number of genes are unlikely to provide enough phylogenetic information.

Gene Tree Conflicts and Coalescent Simulations

Although phylogenetic relationships of the colonial volvocine algae and main clades of the Volvocales were largely resolved with high support by dense gene/taxon sampling and by employing site-heterogeneous and heterotachy models, several inconsistent relationships of interest were observed: 1) the phylogenetic relationship among the three earliest diverging clades within the Volvocales (Radicarteria, Hafniomonas, and Crucicarteria), 2) the phylogenetic position of incertae sedis taxa, and 3) the monophyly of Goniaceae. These inconsistent relationships were all accompanied by short internal branches in coalescent analyses, and topological conflicts among gene trees were widely observed by using quartet supports among three topologies (q1–q3) of the focal internal branch. Of the 1,362 gene trees, 42% (572 genes trees) supported Radicarteria as sister to the remaining Volvocales clades (excluding the earliest 3 clades), 29% (395 gene trees) supported the remaining Volvocales clades as sister to Radicarteria + Hafniomonas, and the other 29% (395 gene trees) supported Hafniomonas as sister to the remaining Volvocales clades (fig. 2 A , node I). Similarly, 36% (491 gene trees) supported the main topology ([Microsporaceae, Sphaeropleales], Golenkinia), 33% (449 gene trees) supported a first alternative topology ([Microsporaceae, Golenkinia], Sphaeropleales), and remaining 31% (422 gene trees) supported a second alternative topology ([Sphaeropleales, Golenkinia], Microsporaceae) (fig. 2 A , node II). As well, 42% (572 gene trees) supported Gonium as sister to Astrephomene, with 29% (395 gene trees) each supporting the remaining 2 alternative topologies (fig. 2 A , node III). In these analyses, the frequencies of gene trees favoring the two alternative topologies (q2–q3) were almost equal.

To further explore whether incomplete lineage sorting (ILS) could better explain the gene tree conflicts, we assembled three data sets (five-taxon (I), five-taxon (II), and seven-taxon) around three conflicting nodes (node I, node II, and node III) and performed coalescent simulations. We estimated the theta parameter for internal branches, which could reflect the level of ILS (high theta value means large ancestor population size and high ILS level) (fig. 2 B ). Then, we simulated 100,000 gene trees from species trees for each data set using the estimated theta values under the coalescent model. The considerable agreement between the observed gene trees and simulated ones with ILS was detected in five-taxon (I) and five-taxon (II) (Pearson's correlation coefficient = 0.998, P < 0.01 for five-taxon (I) data set; Pearson's correlation coefficient = 0.988, P < 0.01 for five-taxon (II) data set; fig. 2 C ). There was a moderate correlation between the simulated and observed topologies in seven-taxon (Pearson's correlation coefficient = 0.90, P < 0.01; fig. 2 C ). We further simulated the gene trees without ILS for each data set, and these gene trees showed a relative low agreement with our observed gene trees (Pearson's correlation coefficient = 0.58/0.65/0.77, for five-taxon (I)/five-taxon (II)/seven-taxon; fig. 2 C ). These results increased the potential for ILS and indicated that ILS may account for the topological discordance of gene trees in the earliest diverging lineages of the Volvocales and in the ancestral branch between incertae sedis taxa and Sphaeropleales. In contrast, for the colonial volvocine algae Goniaceae, the moderate correlation implied that ILS may partly contribute to the observed discordance in gene trees. In this situation, the multispecies coalescent model is likely to perform better than the concatenation approach, as the coalescent methods can better accommodate gene tree discordance due to ILS (Zhong et al. 2013; Liu et al. 2015).

Besides ILS, ancestral hybridization can similarly result in conflicting gene genealogies. To explore the contribution of ILS and hybridization to phylogenetic inconsistency, we performed the χ² test of the quartet frequencies of two minor topologies between empirical gene trees and simulated gene trees with ILS. Our analyses suggested that ILS rather than ancestral hybridization was a major factor causing gene tree heterogeneity surrounding nodes I, II, and III (supplementary fig. S16, Supplementary Material online). In addition, the topological conflicts of most nodes were mainly due to ILS (supplementary fig. S16, Supplementary Material online).

Previous studies have shown that the multispecies coalescent model can accommodate gene tree heterogeneity and produce more accurate species trees in the presence of ILS (Liu et al. 2015; Rannala et al. 2020); thus, the multispecies coalescent model is likely to be appropriate for our data set. In addition, we found that using slow-evolving genes could improve support for key nodes in our coalescent analyses. The phylogenetic tree inferred using 764 genes and the coalescent model supported that Crucicarteria is the earliest diverging lineage, followed by Hafniomonas and Radicarteria within the Volvocales; the four incertae sedis taxa and Sphaeropleales grouped into a monophyletic lineage; and the family Goniaceae is monophyletic. Notably, the monophyly of the Goniaceae was also supported by morphological studies (Kirk et al. 1986; Nozaki 1990). Each cell in that Astrephomene and Gonium is surrounded by a tripartite boundary, which is the important feature that distinguishes their mode of colony formation from all other colonial algae in the Volvocaceae.

Cellular Differentiation Independently Arose at Least Five Times in the Volvocine Algae

The transition from unicellular to multicellular organisms symbolizes the increase in the level of complexity and provides opportunities for cellular differentiation. Our analyses supported that cellular differentiation evolved independently and repeatedly in the colonial volvocine algae. The cellular differentiation independently evolved at least five times: once in Astrephomene, once in section Volvox, and at least three times in the EVP (Eudorina + Volvox + Pleodorina) clade (fig. 1), but we cannot rule out the possibility of an independent origin of Pleodorina thompsonii in the EVP clade (fig. 1). Previous studies suggested at least three separate origins of cellular differentiation: in Astrephomene, section Volvox, and at least once in the EVP clade (Herron and Michod 2008; Grochau-Wright et al. 2017). Conversely, our results implied that cellular differentiation arose independently up to six times. In addition, Volvocine algae has a relatively simple multicellular organization with only two cell types (germ cells and somatic cells). Species with somatic cells and germ cells are considered completely differentiated, such as Volvox within the EVP clade, and their germ cells are used for reproduction, and somatic cells are used for motility; those with soma but no specialized germ are partially differentiated, such as Astrephomene, section Volvox, and Pleodorina, and their somatic cells perform the function of motility. Multicellularity and cellular differentiation appear to be an evolutionarily dynamic trait with multiple origins in the volvocine algae.

Impacts of Different Factors on the Divergence Time of the Volvocales

In attempting to establish a robust timeline of the Volvocales evolution, we explored the impacts of tree topology, data set size, and node prior distributions on time estimations. Firstly, we evaluated the effects of different phylogenetic hypotheses on the time estimations of the major Volvocales clades. These trees accounted for alternative phylogenetic relationships in the Goniaceae and uncertain position of the incertae sedis taxa (the concatenation-based tree using C20 model and the coalescent-based tree based on the 764 genes). Our dating results showed that differences in tree topology had small effects on the time estimates of most nodes (fig. 3 E ), such as the time of the crown Volvocales, TGV clade, and ancestor of Goniaceae + Volvocaceae (fig. 3 E and supplementary table S3, Supplementary Material online). In contrast, nodes with uncertain phylogenetic relationships exhibited some variations in age estimates, like the time estimation of Goniaceae (supplementary table S3, Supplementary Material online). Clock-like genes could minimize errors associated with molecular clock model misspecification. We then explored the effects of various numbers of clock-like genes (400 and 100) for the time estimations. Our dating results showed little difference between the estimates based on 764 genes, 400 genes, and 100 genes. Most nodes had highly congruent age estimates and 95% credibility intervals (CI) (fig. 3 C and D and supplementary table S3, Supplementary Material online). Thus, it is a useful strategy to select a subset of genes under the assumption of the molecular clock, especially in the face of computational burden.

In the absence of reliable fossil records, secondary calibration is an efficient approach for time estimation. Recently, Hou et al. (2022) investigated the green algae evolutionary timeline with the one-billion-year-old Proterocladus fossil (Tang et al. 2020) and proposed several hundred million years of earlier origin of the green algae. Thus, we employed secondary calibrations (table 1) from Hou et al. (2022) to estimate the divergence time of the Volvocales. Accounting for uncertainty in the calibration by using a prior distribution is important. The 95% CI of node age estimates from the primary study was usually used to build a uniform prior distribution. Recent studies suggested that normal distribution is preferred over the uniform prior (Morris et al. 2018; Swenson et al. 2019), with most of the probability density around the mean value. The normal distribution allows for bidirectional uncertainty, and it does not impose unjustified boundaries on the potential node age. Thus, we applied multiple secondary calibrations on internal nodes based on posterior times from Hou et al. (2022) and explored the effect of using uniform or normal prior distributions for secondary calibrations on the time estimates. Molecular dating analyses of the two prior distributions produced almost identical age estimates of the Volvocales and Tetrabaenaceae, as well as the ancestor of Goniaceae and Volvocaceae (fig. 3 A and B and supplementary table S3, Supplementary Material online). These results implied that using different distributions at secondary calibrations had little impact on age estimates for the Volvocales.

The Evolutionary Timescale of the Volvocales

We used 100 clock-like genes with normal distribution under the coalescent-based tree to estimate the timeline of the Volvocales. Our dating analyses indicated that the Volvocales originated in the Neoproterozoic (696.6–551.1 Ma; fig. 4 and supplementary table S3, Supplementary Material online), earlier than the time estimated by Del Cortona et al. (2020) (586–426 Ma) and narrower than the interval estimated by Hou et al. (2022) (695–417 Ma). The common ancestor of colonial TGV clade was diverged in the Early Triassic to Early Jurassic (251.4–180.7 Ma) (fig. 4 and supplementary table S3, Supplementary Material online), similar to the estimated time by Herron et al. (2009) (260–209 Ma). Two independent origins of multicellularity occurred in the volvocine algae, once in the Tetrabaenaceae, originating at 225–140.1 Ma, and another in the ancestor of Goniaceae and Volvocaceae, emerging at 213.1–153 Ma (fig. 4 and supplementary table S3, Supplementary Material online). The Goniaceae and Volvocaceae originated almost simultaneously, with the Goniaceae diverging at 199.1–135.5 Ma and the Volvocaceae occurring at 196–139.3 Ma (fig. 4 and supplementary table S3, Supplementary Material online). As demonstrated by Herron et al. (2009), most morphological and developmental changes occurred in the early and rapid radiation soon after multicellular divergence from unicellular. Similarly, our time estimation indicated that after a time span of about 33 Ma, the Volvox algae diverged from its closest unicellular relatives and completed the initial innovations: transformation of cell wall into ECM and genetic control of cell number.

Taken together, our time-calibrated phylogeny provides a framework for understanding the Volvocales evolution. Molecular dating analyses indicated that the volvocine algae occurred in the Cryogenian–Ediacaran period, which was characterized by two widespread glaciations, resulting in extremely cold global climate and dramatic fluctuations in biogeochemical cycling (Hoffman et al. 2017). These conditions were unfavorable for planktonic eukaryotes, resulting in a significant reduction in photosynthesis and decreasing overall diversity of eukaryotes (Cohen and Macdonald 2015). Since the Cambrian, the decline of glaciation led to increase in habitat space and increased nutrients and oxygen (Lyons et al. 2014), and the break-up of Pangaea during the Triassic to Cretaceous led to the sharply increased continental rainfall (Chaboureau et al. 2014). Combined with the global climate change, we speculate that the decrease in glaciation and subsequent increase in nutrients and oxygen provide external conditions for the origin of the multicellularity in the volvocine algae. Subsequently, the rapid radiation of volvocine may provide a source of food and oxygen for the survival of the zooplankton and promote the diversity of higher plants and animals.

Taxon Sampling and Data Collection

We collected and analyzed nuclear protein-coding genes from 133 species mined from 18 genomes and 115 transcriptomes. The public genome data were downloaded from Phytozome (https://phytozome.jgi.doe.gov/pz/portal.html) and GenBank (http://www.ncbi.nlm.nih.gov/genbank/). The 49 transcriptomes were obtained from the "1,000 plants" project database (http://www.onekp.com/public_data.html) and 60 previously generated transcriptomes (Hu et al. 2019; Lindsey et al. 2021). In addition, we also assembled six transcriptomes of the Volvocales using Trinity with default settings (downloaded from NCBI SRA database). Our data set includes representatives from the major lineages of Sphaeropleales and Volvocales. Five species from the OCC clade were designated as outgroups. We sampled 110 species representing the major clades of the Volvocales, 81 from the Reinhardtinia clade that contains colonial volvocine algae and unicellular species. Our data set includes 60 colonial volvocine algae, covering all genera of the 3 families (Tetrabaenaceae–Goniaceae–Volvocaceae). The complete list of taxa and information of data sources are provided in supplementary table S1, Supplementary Material online.

Ortholog Selection, Alignment, and Trimming

Data from 7 published whole genomes of the Volvocales (C. reinhardtii, Volvox carteri, D. salina, and Gonium pectorale) and Sphaeropleales (Chromochloris zofingiensis, Raphidocelis subcapitata, and Tetradesmus obliquus) were used to identify single-copy orthologs using OrthoFinder v2.4.0 with default parameters. The resulting 2,742 of single-copy orthogroups were used as references in HMMER v3.1b2 (Eddy 2011) to select orthologous genes (OGs) from other 126 genomic or transcriptomic data.

Amino acid sequences of each OG were aligned by MAFFT v7.310 (Katoh and Standley 2013) using the L-INS-I algorithm. We excluded poorly aligned regions using Gblocks 0.91b (Castresana 2000) with half gaps allowed and other default parameters. The short sequences of each OG were removed using trimAl v1.4 (Capella-Gutiérrez et al. 2009) with the parameters -resoverlap 0.5 and -seqoverlap 50. In addition, the sequences shorter than 100 amino acids and taxon occupancy below 80% were excluded, resulting in 1,362 OGs for phylogenetic analysis.

Phylogenetic Inferences

The data set of 1,362 OGs was analyzed using coalescent and concatenation approaches. For the coalescent approach, individual gene tree was reconstructed by using RAxML v8.2.12 (Stamatakis 2014) with the best fitting model (-m PROTGAMMAAUTO -auto-prot = bic) and 200 rapid bootstrap replicates. We collapsed low support branches (<10% bootstrap support) in gene trees with Newick Utilities v1.6 to minimize potential impacts of gene tree error. The species tree was inferred using ASTRAL v5.7.3 (Zhang et al. 2018) with support values estimated by local posterior probabilities and multilocus bootstrapping (gene and site resampling). For the concatenation approach, sequences of 1,362 genes were concatenated into a single alignment representing an aggregate of 341,331 amino acids of 133 taxa by Geneious v10.2.3 (Kearse et al. 2012). The maximum likelihood (ML) tree was reconstructed using IQ-TREE v2.0.5 (Minh et al. 2020) with three strategies: 1) the gene-wise partitioned analysis with optimal partitioning schemes and appropriate substitution models, 2) using site-heterogeneous C20 model, and 3) using heterotachy GHOST model. Nodal support values were estimated using SH-aLRT test (Guindon et al. 2010) and ultrafast bootstrap (Minh et al. 2013) with 1,000 replicates. To reduce the bias generated by genes with fast evolutionary rates for constructing phylogenetic relationships, we calculated the average evolutionary rates of 1,362 genes using PhyKIT (Steenwyk et al. 2021) and ranked according to the rate of evolution. Finally, we divided 1,362 genes into 3 sub-data sets based on the relatively low rates. The 3 sub-data sets were composed of 317,806 (1,267 genes), 289,550 (1,152 genes), and 205,840 (764 genes) amino acids, respectively. The phylogenetic analyses described above were carried out for these three data sets.

ILS Simulation and Hybridization

Normalized quartet supports of the internal branches in the species tree were performed in ASTRAL v5.7.3 with the parameter -t 8. To investigate if gene tree discordance around three conflicting nodes can be explained by ILS, we carried out coalescent simulations following the approach of Cloutier et al. (2019) and Morales-Briones et al. (2021).

First, our data set of 133 taxa was filtered to 3 small data sets: 5-taxon (I), 5-taxon (II), and 7-taxon. The five-taxon (I) data set was selected around the internal branch between the Radicarteria and remaining Volvocales clades (excluding the earliest three clades), the five-taxon (II) data set was selected around the internal branch between the Microsporaceae and Sphaeropleales, and the seven-taxon data set was selected around the internal branch between the Gonium and Astrephomene. The OGs of each data set were aligned and trimmed according to the above "Ortholog selection, alignment and trimming" method. Taxon occupancy below 100% and length shorter than 100 amino acids of OGs were removed, resulting in 1,512 genes for the 5-taxon (I) data set, 1,638 genes for the 5-taxon (II) data set, and 1,391 genes for the 7-taxon data set. Then, gene trees for each data set were reconstructed by using RAxML v8.2.12 (Stamatakis 2014) with the best fitting model; species trees were estimated using ASTRAL v5.7.3, with internal branch lengths assigned in coalescent units. An ultrametric species tree with branch lengths in mutational units (uT) was estimated using LG + GAMMA model and strict molecular clock in PAUP v4.0a (Swofford 2002). The branch lengths in mutational units and coalescent units (t = T/4Ne) were used to estimate the population size parameter theta (theta = uT/t; Degnan and Rosenberg 2009) for internal branches. Theta for terminal branches was set to the constant 1. We simulated 100,000 gene trees from species trees for each data set under the multispecies coalescent model using the function sim.coaltree.sp in Phybase v.1.5 and the estimated theta values (Liu and Yu 2010). In addition, we also simulated 100,000 gene trees from species trees by setting theta to 1/200 of the minimum theta value to represent the absence of ILS or extremely low ILS. Finally, we calculated the topological frequencies of the observed gene trees and the simulated gene trees for each data set. And the correlation between the observed frequencies and the simulated frequencies of gene trees was performed using the cor.test function in R package.

Furthermore, we explored whether ancestral hybridization event during the diversification of volvocine algae was a main cause of gene tree discordance. We estimated the population size parameter theta for internal branches based on the data set of 764 OGs. The 20,000 gene trees were simulated from species trees under the multispecies coalescent model using the function sim.coaltree.sp in Phybase v.1.5 and the estimated theta values. We then calculated quartet frequencies of 2 minor topologies between the empirical gene trees and the simulated gene trees with ILS based on the data set of 764 OGs. Finally, we implemented a two-sided χ² test on those quartet frequencies. If P value is >0.1, phylogenetic discordance surrounding that internal node is due to a significant level of ILS. If P value is lower than 0.1, it indicates phylogenetic discordance surrounding that internal node results from hybridization and other factors (Ma et al. 2021).

Divergence Time Estimation

Divergence times were estimated based on 764 genes using MCMCTree in PAML v4.9 (Yang 2007). The approximate likelihood (Reis and Yang 2011) calculations in MCMCTree were implemented using CODEML under the LG + Γ4 + F model. The MCMC process was run for 15 million generations sampled every 1,000 generations after a burn-in of 1.5 million iterations. We ran 2 independent chains to evaluate convergence and confirmed that the effective sample sizes (ESS) of all parameters were above 200 using Tracer v1.7.1 (Rambaut et al. 2018).

1. Buchheim MA, Kirkwood AE, Buchheim JA, Verghese B, Henley WJ. 2010. Hypersaline soil supports a diverse community of Dunaliella (Chlorophyceae). J Phycol. 46:1038–1047. [Google Scholar]
2. Capella-Gutiérrez S, Silla-Martínez JM, Gabaldón T. 2009. Trimal: a tool for automated alignment trimming in large-scale phylogenetic analyses. Bioinformatics 25(15):1972–1973. [DOI] [PMC free article] [PubMed] [Google Scholar]
3. Castresana J. 2000. Selection of conserved blocks from multiple alignments for their use in phylogenetic analysis. Mol Biol Evol. 17(4):540–552. [DOI] [PubMed] [Google Scholar]
4. Chaboureau AC, Sepulchre P, Donnadieu Y, Franc A. 2014. Tectonic-driven climate change and the diversification of angiosperms. Proc Natl Acad Sci USA. 111:14066–14070. [DOI] [PMC free article] [PubMed] [Google Scholar]
5. Cloutier A, et al. 2019. Whole-genome analyses resolve the phylogeny of flightless birds (Palaeognathae) in the presence of an empirical anomaly zone. Syst Biol. 68(6):937–955. [DOI] [PMC free article] [PubMed] [Google Scholar]
6. Cohen PA, Macdonald FA. 2015. The proterozoic record of eukaryotes. Paleobiology 41(4):610–632. [Google Scholar]
7. Degnan JH, Rosenberg NA. 2009. Gene tree discordance, phylogenetic inference and the multispecies coalescent. Trends Ecol Evol. 24(6):332–340. [DOI] [PubMed] [Google Scholar]
8. Del Cortona A, et al. 2020. Neoproterozoic origin and multiple transitions to macroscopic growth in green seaweeds. Proc Natl Acad Sci U S A. 117(5):2551–2559. [DOI] [PMC free article] [PubMed] [Google Scholar]
9. Delsuc F, Brinkmann H, Chourrout D, Philippe H. 2006. Tunicates and not cephalochordates are the closest living relatives of vertebrates. Nature 439:965–968. [DOI] [PubMed] [Google Scholar]
10. Eddy S. 2011. Accelerated profile HMM searches. Plos Comput Biol. 7:e1002195. [DOI] [PMC free article] [PubMed] [Google Scholar]
11. Fučíková K, Lewis PO, Neupane S, Karol KG, Lewis LA. 2019. Order, please! Uncertainty in the ordinal-level classification of Chlorophyceae. PeerJ 7:6899. [DOI] [PMC free article] [PubMed] [Google Scholar]
12. Graur D, Martin W. 2004. Reading the entrails of chickens: molecular timescales of evolution and the illusion of precision. Trends Genet. 20(2):80–86. [DOI] [PubMed] [Google Scholar]
13. Grochau-Wright ZI, et al. 2017. Genetic basis for soma is present in undifferentiated volvocine green algae. J Evolution Biol. 30(6):1205–1218. [DOI] [PMC free article] [PubMed] [Google Scholar]
14. Guindon S, et al. 2010. New algorithms and methods to estimate maximum-likelihood phylogenies: assessing the performance of PhyML 3.0. Syst Biol. 59(3):307–321.
15. Guiry MD, Guiry GM. 2023. AlgaeBase. World-wide electronic publication, National University of Ireland, Galway. https://www.algaebase.org; searched on 2023年02月06日
16. Hanschen E, Herron M, Nozaki H, Michod R. 2018. Multicellularity drives the evolution of sexual traits. Am Nat. 192(3):E93–E105. [DOI] [PMC free article] [PubMed] [Google Scholar]
17. Heath TA, Hedtke SM, Hillis DM. 2008. Taxon sampling and the accuracy of phylogenetic analyses. J Syst Evol. 46(3):239–257. [Google Scholar]
18. Hedges SB, Kumar S. 2004. Precision of molecular time estimates. Trends Genet. 20(5):242–247. [DOI] [PubMed] [Google Scholar]
19. Herron MD. 2016. Origins of multicellular complexity: Volvox and the volvocine algae. Mol Ecol. 25(6):1213–1223. [DOI] [PMC free article] [PubMed] [Google Scholar]
20. Herron MD, Hackett JD, Aylward FO, Michod RE. 2009. Triassic origin and early radiation of multicellular volvocine algae. Proc Natl Acad Sci U S A. 106(9):3254–3258. [DOI] [PMC free article] [PubMed] [Google Scholar]
21. Herron MD, Michod RE. 2008. Evolution of complexity in the volvocine algae: transitions in individuality through Darwin's Eye. Evolution 62(2):436–451. [DOI] [PubMed] [Google Scholar]
22. Hoffman PF, et al. 2017. Snowball Earth climate dynamics and Cryogenian geology-geobiology. Sci Adv. 3(11):e1600983. [DOI] [PMC free article] [PubMed] [Google Scholar]
23. Hou Z, et al. 2022. Phylotranscriptomic insights into a Mesoproterozoic–Neoproterozoic origin and early radiation of green seaweeds (Ulvophyceae). Nat Commun. 13(1):1610. [DOI] [PMC free article] [PubMed] [Google Scholar]
24. Hu Y, et al. 2019. Evolutionary analysis of unicellular species in Chlamydomonadales through chloroplast genome comparison with the colonial volvocine algae. Front Microbiol. 10:1351. [DOI] [PMC free article] [PubMed] [Google Scholar]
25. Hu Y, Xing W, Song H, Hu Z, Liu G. 2020. Comparison of colonial volvocine algae based on phylotranscriptomic analysis of gene family evolution and natural selection. Eur J Phycol. 55(1):100–112. [Google Scholar]
26. Jackson C, Knoll AH, Chan CX, Verbruggen H. 2018. Plastid phylogenomics with broad taxon sampling further elucidates the distinct evolutionary origins and timing of secondary green plastids. Sci Rep. 8(1):1523. [DOI] [PMC free article] [PubMed] [Google Scholar]
27. Jarvis ED, et al. 2014. Whole-genome analyses resolve early branches in the tree of life of modern birds. Science 346(6215):1320–1331. [DOI] [PMC free article] [PubMed] [Google Scholar]
28. Katoh K, Standley DM. 2013. MAFFT multiple sequence alignment software version 7: improvements in performance and usability. Mol Biol Evol. 30(4):772–780. [DOI] [PMC free article] [PubMed] [Google Scholar]
29. Kearse M, et al. 2012. Geneious Basic: an integrated and extendable desktop software platform for the organization and analysis of sequence data. Bioinformatics 28(12):1647–1649. [DOI] [PMC free article] [PubMed] [Google Scholar]
30. Kirk DL. 2005. A twelve-step program for evolving multicellularity and a division of labor. Bioessays 27(3):299–310. [DOI] [PubMed] [Google Scholar]
31. Kirk DL, Birchem R, King N. 1986. The extracellular matrix of Volvox: a comparative study and proposed system of nomenclature. J Cell Sci. 80(1):207–231. [DOI] [PubMed] [Google Scholar]
32. Lemieux C, Vincent A, Labarre A, Otis C, Turmel M. 2015. Chloroplast phylogenomic analysis of chlorophyte green algae identifies a novel lineage sister to the Sphaeropleales (Chlorophyceae). BMC Evol Biol. 15:264. [DOI] [PMC free article] [PubMed] [Google Scholar]
33. Li X, et al. 2023. Chloroplast phylogenomics of unicellular and colonial Volvocales provides perspectives on the evolution of morphological characters. J Syst Evol. 61(1):127–142. [Google Scholar]
34. Lindsey CR, Rosenzweig F, Herron MD. 2021. Phylotranscriptomics points to multiple independent origins of multicellularity and cellular differentiation in the volvocine algae. BMC Biol. 19(1):182. [DOI] [PMC free article] [PubMed] [Google Scholar]
35. Liu C, Huang X, Wang X, Zhang X, Li G. 2006. Phylogenetic studies on two strains of antarctic ice algae based on morphological and molecular characteristics. Phycologia 45:190–198. [Google Scholar]
36. Liu L, Xi Z, Davis CC. 2015. Coalescent methods are robust to the simultaneous effects of long branches and incomplete lineage sorting. Mol Biol Evol. 32(3):791–805. [DOI] [PubMed] [Google Scholar]
37. Liu L, Yu L. 2010. Phybase: an R package for species tree analysis. Bioinformatics 26(7):962–963. [DOI] [PubMed] [Google Scholar]
38. Lyons TW, Reinhard CT, Planavsky NJ. 2014. The rise of oxygen in Earth's early ocean and atmosphere. Nature 506(7488):307–315. [DOI] [PubMed] [Google Scholar]
39. Ma J, et al. 2021. The Chloranthus sessilifolius genome provides insight into early diversification of angiosperms. Nat Commun. 12:6929. [DOI] [PMC free article] [PubMed] [Google Scholar]
40. Massoni J, Forest F, Sauquet H. 2014. Increased sampling of both genes and taxa improves resolution of phylogenetic relationships within Magnoliidae, a large and early-diverging clade of angiosperms. Mol Phylogenet Evol. 70:84–93. [DOI] [PubMed] [Google Scholar]
41. Minh BQ, et al. 2020. IQ-TREE 2: new models and efficient methods for phylogenetic inference in the genomic era. Mol Biol Evol. 37(5):1530–1534. [DOI] [PMC free article] [PubMed] [Google Scholar]
42. Minh BQ, Nguyen MA, von Haeseler A. 2013. Ultrafast approximation for phylogenetic bootstrap. Mol Biol Evol. 30(5):1188–1195. [DOI] [PMC free article] [PubMed] [Google Scholar]
43. Morales-Briones DF, et al. 2021. Disentangling sources of gene tree discordance in phylogenomic data sets: testing ancient hybridizations in Amaranthaceae s.l. Syst Biol. 70(2):219–235. [DOI] [PMC free article] [PubMed] [Google Scholar]
44. Morris JL, et al. 2018. The timescale of early land plant evolution. Proc Natl Acad Sci U S A. 115(10):E2274–E2283. [DOI] [PMC free article] [PubMed] [Google Scholar]
45. Nakada T, Misawa K, Nozaki H. 2008. Molecular systematics of Volvocales (Chlorophyceae, Chlorophyta) based on exhaustive 18S rRNA phylogenetic analyses. Mol Phylogenet Evol. 48(1):281–291. [DOI] [PubMed] [Google Scholar]
46. Nakada T, Tomita M, Nozaki H. 2010. Volvulina compacta (Volvocaceae, Chlorophyceae), new to Japan, and its phylogenetic position. J Jap Bot. 85:364–369. [Google Scholar]
47. Nakada T, Tsuchida Y, Tomita M. 2019. Improved taxon sampling and multigene phylogeny of unicellular chlamydomonads closely related to the colonial volvocalean lineage Tetrabaenaceae-Goniaceae-Volvocaceae (Volvocales, Chlorophyceae). Mol Phylogenet Evol. 130:1–8. [DOI] [PubMed] [Google Scholar]
48. Nakazawa A, Krienitz L, Nozaki H. 2001. Taxonomy of the unicellular green algal genus Vitreochlamys (Volvocales), based on comparative morphology of cultured material. Eur J Phycol. 36(2):113–128. [Google Scholar]
49. Němcová Y, Eliáš M, Škaloud P, Hodač L, Neustupa J. 2011. Jenufa gen. nov.: a new genus of coccoid green algae (Chlorophyceae, incertae sedis) previously recorded by environmental sequencing. J Phycol. 47(4):928–938. [DOI] [PubMed] [Google Scholar]
50. Nozaki H. 1990. Ultrastructure of the extracellular matrix of Gonium (Volvocales, Chlorophyta). Phycologia 29(1):1–8. [Google Scholar]
51. Nozaki H, et al. 2000. Origin and evolution of the colonial volvocales (Chlorophyceae) as inferred from multiple, chloroplast gene sequences. Mol Phylogenet Evol. 17:256–268. [DOI] [PubMed] [Google Scholar]
52. Nozaki H, Itoh M, Watanabe MM, Kuroiwa T. 1996. Ultrastructure of the vegetative colonies and systematic position of Basichlamys (Volvocales, Chlorophyta). Eur J Phycol. 31(1):67–72. [Google Scholar]
53. Nozaki H, Kuroiwa T. 1992. Ultrastructure of the extracellular matrix and taxonomy of Eudorina, Pleodorina and Yamagishiella gen. nov. (Volvocaceae, Chlorophyta). Phycologia 31(6):529–541. [Google Scholar]
54. Nozaki H, Yamada TK, Takahashi F, Matsuzaki R, Nakada T. 2014. New "missing link" genus of the colonial volvocine green algae gives insights into the evolution of oogamy. BMC Evol Biol. 14:37. [DOI] [PMC free article] [PubMed] [Google Scholar]
55. Nute M, Chou J, Molloy EK, Warnow T. 2018. The performance of coalescent-based species tree estimation methods under models of missing data. BMC Genomics. 19(Suppl 5):286. [DOI] [PMC free article] [PubMed] [Google Scholar]
56. Nye E, Feist-Burkhardt S, Horne DJ, Ross AJ, Whittaker JE. 2008. The palaeoenvironment associated with a partial Iguanodon skeleton from the Upper Weald Clay (Barremian, Early Cretaceous) at Smokejacks Brickworks (Ockley, Surrey, UK), based on palynomorphs and ostracods. Cretaceous Res. 29(3): 417–444. [Google Scholar]
57. One Thousand Plant Transcriptomes Initiative . 2019. One thousand plant transcriptomes and the phylogenomics of green plants. Nature 574(7780):679–685. [DOI] [PMC free article] [PubMed] [Google Scholar]
58. Pröschold T, Darienko T, Krienitz L, Coleman AW. 2018. Chlamydomonas schloesseri sp. nov. (Chlamydophyceae, Chlorophyta) revealed by morphology, autolysin cross experiments, and multiple gene analyses. Phytotaxa 362:21. [Google Scholar]
59. Rambaut A, Drummond AJ, Xie D, Baele G, Suchard MA. 2018. Posterior summarization in Bayesian phylogenetics using Tracer 1.7. Syst Biol. 67(5):901–904. [DOI] [PMC free article] [PubMed] [Google Scholar]
60. Rannala B, Edwards SV, Leaché A, Yang Z. 2020. The multispecies coalescent model and species tree inference. Phylogenet. Genom. Era. book section 3.3, pp. 3.3: 1-21.
61. Reis M, Yang Z. 2011. Approximate likelihood calculation on a phylogeny for Bayesian estimation of divergence times. Mol Biol Evol. 28(7):2161–2172. [DOI] [PubMed] [Google Scholar]
62. Roure B, Baurain D, Philippe H. 2013. Impact of missing data on phylogenies inferred from empirical phylogenomic data sets. Mol Biol Evol. 30:197–214. [DOI] [PubMed] [Google Scholar]
63. Sauquet H. 2013. A practical guide to molecular dating. C R Palevol. 12(6):355–367. [Google Scholar]
64. Smith S, Brown J, Walker J. 2018. So many genes, so little time: a practical approach to divergence-time estimation in the genomic era. PLoS One 13(5):e0197433. [DOI] [PMC free article] [PubMed] [Google Scholar]
65. Stamatakis A. 2014. RAxML version 8: a tool for phylogenetic analysis and post-analysis of large phylogenies. Bioinformatics 30(9):1312–1313. [DOI] [PMC free article] [PubMed] [Google Scholar]
66. Steenwyk JL, et al. 2021. PhyKIT: a broadly applicable UNIX shell toolkit for processing and analyzing phylogenomic data. Bioinformatics 37(16):2325–2331. [DOI] [PMC free article] [PubMed] [Google Scholar]
67. Swenson U, Havran JC, Munzinger J, Mcloughlin S, Nylinder S. 2019. Metapopulation vicariance, age of island taxa and dispersal: a case study using the pacific plant genus planchonella (Sapotaceae). Syst Biol. 68(6):1020–1033. [DOI] [PMC free article] [PubMed] [Google Scholar]
68. Swofford D. 2002. PAUP*. Phylogenetic analysis using parsimony (*and other methods) (version 4). Sunderland: Sinauer Associates. [Google Scholar]
69. Tang Q, Pang K, Yuan X, Xiao S. 2020. A one-billion-year-old multicellular chlorophyte. Nat Ecol Evol. 4:543–549. [DOI] [PMC free article] [PubMed] [Google Scholar]
70. Tao Q, Barba-Montoya J, Kumar S. 2021. Data-driven speciation tree prior for better species divergence times in calibration-poor molecular phylogenies. Bioinformatics 37(Suppl_1):i102–i110. [DOI] [PMC free article] [PubMed] [Google Scholar]
71. Tappan HN. 1980. The paleobiology of plant protists. San Francisco: Freeman WH. [Google Scholar]
72. Teyssèdre B. 2007. Precambrian paleontology in the light of molecular phylogeny an example: the radiation of the green algae. Biogeosci Discuss. 4:3123–3142. [Google Scholar]
73. Umen JG. 2020. Volvox and volvocine green algae. Evodevo 11(1):13. [DOI] [PMC free article] [PubMed] [Google Scholar]
74. Vigran JO, Mangerud G, Mork A, Bugge T, Weitschat W. 1998. Biostratigraphy and sequence stratigraphy of the Lower and Middle Triassic deposits from the Svalis Dome, central Barents Sea, Norway. Palynology 22(1):89–141. [Google Scholar]
75. Yamashita S, et al. 2021. Genome sequencing of the multicellular alga Astrephomene provides insights into convergent evolution of germ-soma differentiation. Sci Rep. 11:22231.
76. Yang ZH. 2007. PAML 4: phylogenetic analysis by maximum likelihood. Mol Biol Evol. 24(8):1586–1591. [DOI] [PubMed] [Google Scholar]
77. Zhang Z, et al. 2020. Adaptation to extreme antarctic environments revealed by the genome of a sea ice green alga. Curr Biol. 30(17):3330–3341.e7. [DOI] [PubMed] [Google Scholar]
78. Zhang Z, et al. 2022. Origin and evolution of green plants in the light of key evolutionary events. J Integr Plant Biol. 64:516–535. [DOI] [PubMed] [Google Scholar]
79. Zhang C, Rabiee M, Sayyari E, Mirarab S. 2018. ASTRAL-III: polynomial time species tree reconstruction from partially resolved gene trees. BMC Bioinformatics. 19:153. [DOI] [PMC free article] [PubMed] [Google Scholar]
80. Zhong B, Betancur-R R. 2017. Expanded taxonomic sampling coupled with gene genealogy interrogation provides unambiguous resolution for the evolutionary root of angiosperms. Genome Biol Evol. 9(11):3154–3161. [Google Scholar]
81. Zhong B, Liu L, Yan Z, Penny D. 2013. Origin of land plants using the multispecies coalescent model. Trends Plant Sci. 18(9):492–495. [DOI] [PubMed] [Google Scholar]

Supplementary Materials

Data Availability Statement

Supplementary figures S1–S16, Supplementary Material online, and supplementary tables S1–S3, Supplementary Material online, are available for this paper. The alignment data, phylogenetic trees, and time trees are available from https://doi.org/10.6084/m9.figshare.21652112.