Fossil-calibrated molecular clock data enable reconstruction of steps leading to differentiated multicellularity and anisogamy in the Volvocine algae

Abstract

Background

Throughout its nearly four-billion-year history, life has undergone evolutionary transitions in which simpler subunits have become integrated to form a more complex whole. Many of these transitions opened the door to innovations that resulted in increased biodiversity and/or organismal efficiency. The evolution of multicellularity from unicellular forms represents one such transition, one that paved the way for cellular differentiation, including differentiation of male and female gametes. A useful model for studying the evolution of multicellularity and cellular differentiation is the volvocine algae, a clade of freshwater green algae whose members range from unicellular to colonial, from undifferentiated to completely differentiated, and whose gamete types can be isogamous, anisogamous, or oogamous. To better understand how multicellularity, differentiation, and gametes evolved in this group, we used comparative genomics and fossil data to establish a geologically calibrated roadmap of when these innovations occurred.

Results

Our ancestral-state reconstructions, show that multicellularity arose independently twice in the volvocine algae. Our chronograms indicate multicellularity evolved during the Carboniferous-Triassic periods in Goniaceae + Volvocaceae, and possibly as early as the Cretaceous in Tetrabaenaceae. Using divergence time estimates we inferred when, and in what order, specific developmental changes occurred that led to differentiated multicellularity and oogamy. We find that in the volvocine algae the temporal sequence of developmental changes leading to differentiated multicellularity is much as proposed by David Kirk, and that multicellularity is correlated with the acquisition of anisogamy and oogamy. Lastly, morphological, molecular, and divergence time data suggest the possibility of cryptic species in Tetrabaenaceae.

Conclusions

Large molecular datasets and robust phylogenetic methods are bringing the evolutionary history of the volvocine algae more sharply into focus. Mounting evidence suggests that extant species in this group are the result of two independent origins of multicellularity and multiple independent origins of cell differentiation. Also, the origin of the Tetrabaenaceae-Goniaceae-Volvocaceae clade may be much older than previously thought. Finally, the possibility of cryptic species in the Tetrabaenaceae provides an exciting opportunity to study the recent divergence of lineages adapted to live in very different thermal environments.

Background

Major evolutionary transitions occur when multiple autonomous units (e.g., genes) combine to form an interdependent autonomous unit (e.g., chromosomes) capable of storing and transmitting information in a novel way [1]. Over the past four billion years a relatively small number of such transitions have resulted in a myriad of innovations that have contributed to the diversification of life on Earth [1] [2]. Among the most conspicuous of these is the transition from organisms whose individuals consist of one cell to individuals that consist of many.

Multicellularity has independently evolved from a unicellular ancestor at least 45 times [3], and has been documented in all three domains of life: Bacteri a[4], Archaea [5], and Eukarya [6]. In most cases, this transition opened the door to evolution of division of labor, which in turn paved the way for cellular differentiation, where cells in a multicellular body take on specific tasks. Task specialization has the potential to boost a multicellular organism’s fitness, provided that the programmed partitioning of functions provides metabolic, structural, and/or genetic advantages over retaining all functions in all cells. Cellular differentiation occurs in most but not all multicellular eukaryotes.

Though multicellularity and differentiation have repeatedly evolved, we have limited knowledge of what selects for these features and the genetic steps that enable their initial evolution [7]. To fill this knowledge gap, we can use comparative genomics to study evolutionary lineages, or clades, in which some species are unicellular and others multicellular, and in which the multicellular species exhibit varying degrees of cellular differentiation. Several eukaryotic clades fulfill these requirements [8, 9], notably the volvocine algae [10], a small extant group of freshwater green algae nested within Viridiplantae in the Chlorophycean order Volvocales (Fig. 1).

Fig. 1.

Cladogram and photographs of volvocine genera. Depicted evolutionary relationships derived from Lindsey et al. [11] with Goniaceae modified to be monophyletic based upon more recent evidence from Ma et al. [12]. The three multicellular volvocine families (Tetrabaenaceae, Goniaceae, and Volvocaceae) are represented by orange, purple, and green font and borders. Unicellularity is denoted by black font and borders. Blue lines correspond to anisogamous lineages, and blue to pink gradient lines represent lineages in which oogamy has evolved. A single asterisk (*) denotes somatic cell differentiation, and double asterisks (**) denote complete germ-soma differentiation. Superscript letters correspond to their respective photographs. Photographs are arranged in the following order: (A) Tetrabaena socialis, (B) Chlamydomonas reinhardtii, (C) Gonium pectorale, (D) Astrephomene gubernaculifera, (E) Platydorina caudata, (F) Colemanosphaera charkowiensis, (G) Pandorina morum, (H) Volvulina compacta, (I) Yamagishiella unicocca, (J) Eudorina elegans, (K) Pleodorina starrii, (L) Volvox carteri. Photos are not to scale. Figure Credit for (B) and (C): Deborah Shelton

The volvocine algae are facultatively sexual [13]; thus, in addition to being a model system for the study of multicellularity and differentiation, they can also provide insight into another major transition: the evolution of male and female sexes via gametic differentiation. The definition and elaboration of the sexes, as well as the evolution of sexually selected traits, is predicated on the advent of anisogamy, gamete types that differ in size. Oogamy constitutes a specialized form of anisogamy wherein small male gametes are motile and large female gametes are not (e.g., vertebrates and land plants are oogamous). Anisogamy is hypothesized to have evolved from isogamy [14], in which the genetic determinants of mating type segregate during meiosis but the resulting gametes are identical in size and shape [15]. Because multicellular lineages such as animals and plants lack extant isogamous forms, the origins of anisogamy cannot be deduced using these groups. By contrast, the volvocine algae range from unicellular, isogamous species like Chlamydomonas to multicellular, oogamous species like Volvox, which exhibits complete germ-soma differentiation. Following in the vein of previous volvocine authors, terms such as “multicellular”, “colonies”, and “colonial” are used interchangeably. The use of this terminology (e.g., “colonies” to describe multicellular volvocine taxa) for this group can be traced back to the early 20^th century [16]. We refer to any volvocine species with more than one cell as “colonial” or “multicellular”. For additional context on terminology used to describe volvocine algae, see footnote in Kirk, 1998 pp. 115-116) [17].

Multicellular volvocine species are grouped within three families: Tetrabaenaceae, Goniaceae, and Volvocaceae. According to recent phylogenetic evidence [11], these families are paraphyletic. The unicellular Vitreochlamys ordinata is sister to the Tetrabaenaceae. This clade of V. ordinata + Tetrabaenaceae is sister to several Chlamydomonas species plus the multicellular Goniaceae + Volvocaceae. Both Vitreochlamys and Chlamydomonas are isogamous, undifferentiated unicells. The genera composing the Tetrabaenaceae, Tetrabaena and Basichlamys, are 4-celled colonies that are isogamous and undifferentiated. The Goniaceae comprises two genera, Gonium and Astrephomene. Gonium consists of flattened 8-32 celled colonies that are isogamous and undifferentiated. Astrephomene, also isogamous, includes 32-128 celled spheroidal colonies that exhibit somatic cell differentiation. Somatic cells engage in specific cellular functions, principally motility, are mortal, and do not pass on genetic information to subsequent generations. Among the 8 genera that comprise the Volvocaceae, Pandorina, Volvulina, and Yamagishiella are isogamous and undifferentiated; Colemanosphaera, Eudorina, and Platydorina are anisogamous and undifferentiated; Pleodorina is anisogamous with somatic differentiation, and Volvox is oogamous with specialized germ and somatic cells. Specialized germ cells pass on genetic information and are distinct from undifferentiated cells in that they do not significantly contribute to colony motility.

Because the volvocine algae are closely related, but vary in cellularity, differentiation, and sexuality, comparative genomics offers the prospect of understanding how and when each of these traits evolved. Nearly 15 years ago Herron et al [18]. estimated divergence times for this group while analyzing the temporal sequence of certain developmental traits. Since then, a wealth of molecular, paleontological, and phylogenetic data has come to light, motivating us to reevaluate their conclusions. For example, multiple volvocine genomes [19,20,21,22,23] have been published, as well as the transcriptomes of all known extant volvocine genera [11, 24], including a new genus [25]; furthermore, genomes and transcriptomes have been sampled among different strains within many volvocine species [11, 24]. Also, the fossil record has been enriched by the discovery of a billion-year-old multicellular chlorophyte [26], firm dates of more than a billion years have been established for simple, multicellular red algae [27], and red algal fossils have been observed that may be 1.6 billion years old [28]. Using a dataset of 40 nuclear-protein coding genes, Lindsey et al [11]. recently reported that (i) multicellularity evolved at least twice in this group, (ii) cell differentiation evolved at least four times, (iii) anisogamy evolved from isogamous ancestors at least three times, and (iv) oogamy evolved at least three times. The first two of these findings have since been bolstered by Ma et al. [12].

When divergence times were estimated by Herron et al. [18], their molecular dataset consisted of 5 chloroplast protein-coding genes plus the 18S nuclear, small ribosomal subunit across 35 volvocine taxa. Subsequent studies undertaking ancestral-state reconstruction [29,30,31,32] of this group have relied solely on the 5 chloroplast gene dataset for their analyses. Unsurprisingly, the resulting trees exhibited essentially the same branching order, none of which reflect the new insights gained by Lindsey et al. [11] As noted, Ma et al. [12] also inferred divergence times of the volvocine algae using large nuclear datasets, employing a single relaxed clock model. However, these researchers did not use their phylogenetic inferences to perform ancestral-state reconstruction of sexual and developmental traits.

Here, we provide a geologically-calibrated roadmap of when multicellularity, differentiation, and anisogamy arose in the volvocine algae. Since no reliable fossils exist for this group, we sampled fossils across the Archaeplastida, or kingdom Plantae (sensu lato), where reliable fossils are plentiful. Specifically, we sampled 14 fossil taxa across the three major Archaeplastida clades (Rhodophyta, Streptophyta, and Chlorophyta) selecting them so as to calibrate our time-tree over an interval of at least one billion years. Our molecular data consist of amino acid sequences for 263 single-copy nuclear genes drawn from 164 taxa across the Archaeplastida. Sequences were obtained both from publicly available genomic and transcriptomic datasets as well as from our own published RNA-Seq dataset [11]. Our goal was to reanalyze the divergence times and the gain and/or loss of traits related to multicellularity, differentiation, and sexuality for the volvocine algae. Our data, presented as ancestral-state reconstructions, lead us to conclude that there were two independent origins of multicellularity in the volvocine algae, as predicted by Lindsey et al. [11] and Ma et al. [12] Moreover, under one of our four relaxed clock models reassessment of divergence times indicates that multicellularity may have originated in the Goniaceae + Volvocaceae sometime between the Carboniferous and Triassic, much older than previous inferrences [12, 18]. Among the Tetrabaenaceae we find that multicellularity may have arisen as early as the Carboniferous or as recently as the Cretaceous. Using our chronogram, we also reassessed the developmental milestones proposed in David L. Kirk’s 12-Step Program [33] to differentiated multicellularity in the volvocine algae. Through ancestral state reconstructions, we show that the temporal sequence of events he hypothesized almost 20 years ago is essentially correct. To investigate multicellularity’s role in the evolution of sex in the volvocine algae, we reconstructed 7 sexual traits. When the gain and loss of these traits are mapped onto our time tree, multicellularity’s role as a driver of sexual evolution becomes evident, consistent with the results of Hanschen et al [32]. Finally, our analyses of morphological, molecular, and divergence time data suggest the possibility of a cryptic species in the multicellular Tetrabaenaceae.

Results and Discussion

Selection of orthologues in taxa across the Archaeplastida provides the basis for inferring how multicellularity and cellular differentiation evolved in the volvocine algae. We sampled a total of 164 taxa representing all major clades of Archaeplastida: Rhodophyta, Streptophyta, and Chlorophyta. Specifically, we sampled 12 rhodophytes across 3 red algal subclades, 45 streptophytes representing 13 major green algal and land plant subclades, and 107 chlorophytes, representing prasinophytes, Ulvophyceae, Trebouxiophyceae, and Chlorophyceae, including 68 unicellular and multicellular volvocine algae (Fig. 2, Additional File 1: Figure S1, Figure S2, and Figure S3). We chose the red algae for our outgroup, as multiple studies indicate that they are sister to the green algae + land plants (Embryophyta), and that the red and green algae share a common plastid ancestor [34, 35]. Although our chief aim was to discern, on a geological timescale, the evolution of traits leading to multicellularity and differentiation in the volvocine algae, we included other green algal genera such as Coleochaete and Tetraselmis to further illuminate the history of Viridiplantae. The former is generally considered to be relatively closely related to land plants [36, 37], while the latter is believed to be sister to the three major Chlorophyta clades [38]. Because all multicellular lineages necessarily evolved from unicellular ancestors, clarifying which lineages are sister to which multicellular clades will enable a deeper understanding of how this major transition occurred.

Fig. 2.

Time-calibrated phylogeny of the Archaeplastida. Branching order of the tree was inferred under maximum-likelihood analysis from an aligned amino acid, concatenated dataset of 263 nuclear genes. Numbers on branches represent bootstrap and posterior probability values, respectively. Branch lengths, corresponding to time, were inferred under the CIR relaxed clock model using 16 most clock-like genes as determined by Sortadate. Blue bars correspond to the inferred 95% HPD interval for each node. Red bubbles correspond to calibrated nodes (Table 1), and blue bubbles correspond to key divergences among the volvocine algae (Figure S4B). Members of the multicellular volvocine algae (Tetrabaenaceae, Goniaceae, and Volvocaceae) are denoted in orange, purple, and green. Taxa in black font are unicellular volvocine algae

A single dataset consisting of 263 single-copy, protein coding genes was compiled and analyzed under Maximum-Likelihood (ML), Bayesian inference (BI), and coalescent-based (CB) phylogenetic methods. These 263 genes were conserved across the three main clades of Archaeplastida. Our concatenated alignment for ML and BI analyses represents an aggregate of 79,844 amino acids, equivalent to 239,532 nucleotide positions, with a total of 62,106 parsimony-informative sites. All raw reads used to complete our single-gene and concatenated alignments encompassing 164 total taxa were mined from previously published data located in public repositories (Additional File 2: Table S1) [11, 20,21,22, 24, 39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83].

Phylogenetic analyses of the volvocine algae indicate multiple independent origins of multicellularity and differentiation. Our ML, BI, and CB analyses all indicate two independent origins of multicellularity among the volvocine algae: one in the lineage leading to the Tetrabaenaceae and another in the lineage leading to the Goniaceae + Volvocaceae (Fig. 2, Additional file 1: Figure S1, Figure S2, and Figure S3). These results bolster the findings of Lindsey et al. [11] and Ma et al. [12] As before, our ML and BI analyses have identical branching orders. The Tetrabaenaceae + Vitreochlamys ordinata are shown to be sister to several Chlamydomonas species + the Goniaceae and Volvocaceae in our ML and BI trees (Maximum-likelihood bootstrap support (MLBS) = 100, Bayesian posterior probablility (BPP) = 1.0) (Fig. 2, Additional file 1: Figure S1, and Figure S2), and this topology both corroborates and expands on with those presented in several earlier studies [11, 12]. Our CB tree, however, indicates a slightly different branching order for the Tetrabaenaceae. The resulting coalescent-based tree shows the Tetrabaenaceae + V. ordinata forming a clade with Chlamydomonas reinhardtii and its relatives, and this clade is shown to be sister to the Goniaceae + Volvocaceae (Coalescent posterior probability (CPP) = 1.0) (Additional file 1: Figure S3). Of note, the sister relationship between the Tetrabaenaceae + V. ordinata and the Chlamydomonas clade is poorly supported in our CB analysis (CPP = 0.34).

In accordance with the findings of Lindsey et al. [11], all three of our phylogenetic analyses indicate a minimum of 4 independent origins of somatic cellular differentiation and a minimum of 3 origins of anisogamy among the volvocine algae. Origins of somatic cell differentiation occur in the following lineages: (i) Astrephomene, (ii) section Volvox, (iii) Pleodorina thompsonii, and (iv) the Pleodorina japonica + Volvox carteri clade. Anisogamy evolved from isogamous ancestors at least three times in the following lineages: (i) Astrephomene, (ii) section Volvox, and (iii) in the Eudorina + Volvox + Pleodorina (EVP) clade.

In contrast to several recent volvocine studies [11, 23], all three of our phylogenetic analyses support the major conclusion that the Goniaceae are monophyletic, albeit with varying support values (MLBS = 95, BPP = 1.0, CPP = 0.67). This conclusion was reached by multiple earlier investigations [12, 18, 31, 32, 84,85,86,87,88,89,90,91]. Consistent with results of Lindsey et al. [11] and Ma et al. [12], we find that Volvox section Volvox is sister to the Pandorina + Volvulina + Colemanosphaera (PVC) and the (EVP) clades within the Volvocaceae. Our ML + BI and CB trees contain minor branching order differences in the EVP clade, but these have no bearing on our major findings.

Variation in divergence times inferred under different relaxed molecular clock models necessitate validation tests. All divergence time estimates were inferred by Phylobayes 4.1b [92] under a Bayesian approach for our inferred ML and CB species trees. For each topology, a total of 14 nodes were calibrated across Rhodophyta, Streptophyta, and Chlorophyta (Fig. 2 and Table 1), and four relaxed clock models: autocorrelated lognormal (LN) [93] and Cox-Ingersoll-Ross (CIR) [94] models, and the uncorrelated gamma (UGAM) and white -noise (WN) models [95]. Autocorrelated models (i.e., LN and CIR) allow the rate of evolution to vary across branches, and the rate of evolution is more similar along branches for closely related species compared to distantly related taxa [93, 94]. Uncorrelated models (i.e., UGAM and WN) assume each branch of will have its own unique rate of evolution irrespective of rates across branches for closely and distantly related species [94, 95]. Inferred divergence times established under all clock models are largely consistent across the three major red and green algal clades (Figure S4). However, there are nodes such as the root age, earliest rhodophyte divergence, and major divergences within the volvocine algae where one model infers a date markedly younger or older than the others (Fig. 3).

Table 1 Fossils used to calibrate nodes in Archaeplastida tree

Fig. 3.

Estimated divergence times of the volvocine algae. Branching order of the tree was inferred under maximum-likelihood analysis from an aligned amino acid, concatenated dataset of 263 nuclear genes. Numbers on branches represent bootstrap and posterior probability values, respectively. Branch lengths, corresponding to time, were inferred under the CIR relaxed clock model using 16 most clock-like genes as determined by Sortadate. Blue bars correspond to the inferred 95% HPD interval for each node. Green bubbles correspond to a developmental trait gain, and red bubbles corresponds to loss of a trait. The figure table lists the 12 developmental traits identified by Kirk in their original order. Blue bubbles indicate key divergences in the volvocine algae. Members of the multicellular volvocine algae (Tetrabaenaceae, Goniaceae, and Volvocaceae) are denoted in orange, purple, and green. Taxa in black font are unicellular volvocine algae

The estimated mean root age of all tested clock models varies by ~630 million years (MY) with the WN model estimating the earliest red algal divergence as ~2016 million years ago (MYA), and the CIR model estimating it as late as ~1385 MYA (Additional file 1: Figure S4B). Averaging the mean root ages of the LN, CIR, and UGAM models produces an average root age of ~1503 MYA, reducing the variance in estimated mean root age between WN and the others to ~500 MY. Similarly, for the earliest rhodophyte divergence, the WN model inferred a mean date significantly older than all other models by at least 300 MY, whereas the LN, CIR, and UGAM models estimated mean ages within <100 MY of each other. Large date discrepancies such as these may be solely due to differences in clock algorithms [104, 105].

For all volvocine algae divergences, the UGAM model estimated ages markedly younger than the three other models we tested (Additional file 1: Figure S4A, B (divergences 9-11)). It is noteworthy that the inferred ages of the UGAM model for the volvocine algae are very similar to the dates inferred by Ma et al. [12], who used a single relaxed clock model. When taking the average of the mean dates inferred by the LN, CIR, and WN models, there is a consistent ~130 MY difference between the average and a mean node date inferred by the UGAM model for this group. Significantly, the UGAM 95% highest posterior density (HPD) intervals for the volvocine algal divergence do not overlap with any of the 95% HPD intervals inferred by other models, conversely the other three models see a comfortable overlap in their 95% HPD intervals for volvocine divergence estimates.

Given that the WN and UGAM models inferred outlier dates for certain key divergence events, we decided to exclude both as clock models. Ages inferred by the CIR model in this study were never observed to be outliers for key divergence events, and the ages inferred by the CIR model for the volvocine algae are nearly identical to the WN model’s estimates for this group. Additionally, divergence estimates inferred under the CIR model for the volvocine algae were the most conservative among the LN, WN, and CIR results. Altogether, the foregoing considerations prompted us to report dates for major divergence events using data produced by the CIR relaxed clock model (Additional file 1: Figure S4A, B). Lepage et al. [94] performed model comparisons of the CIR, LN, UGAM, and WN relaxed clock models using real, as opposed to simulated, molecular datasets that varied in type (nuclear vs mitochondrial and nucleotide vs amino acid) and number of taxa. Their tests concluded that autocorrelated clock models performed best against the various datasets used in their study. Furthermore, LN and CIR models outperformed UGAM and WN models when taxonomic sampling was high, and the CIR model was recommended by the authors for dense taxonomic datasets. Fossil cross validation tests were performed as described earlier for each relaxed clock model, with the objective of identifying a best-performing analysis. However, no relaxed clock model performed markedly better than any of the others (Additional file 2: Table S2).

Fossil selections and Archaeplastida phylogenetic results and divergence times. Since no reliable fossils exist for the volvocine algae, we selected 14 fossil taxa across the Archaeplastida where reliable fossils are abundant (Fig. 2 and Table 1). For each primary fossil calibration used in this study, each calibration point was constrained to a range rather than a fixed-point estimate, acknowledging the inherent uncertainty in fossil ages. With each date range, we specified a soft bound where 2.5% of the total probability mass is positioned outside of the specified lower and/or upper bounds. Detailed information regarding each fossil taxon and calibration point may be found in Additional file 3: Supplementary Methods [18, 26, 27, 75, 96, 97, 99,100,101,102,103,104,105,106,107,108,109,110,111,112,113,114,115,116,117,118,119,120,121,122,123,124,125,126,127,128,129,130,131,132,133,134,135,136,137,138,139,140,141,142]. Moreover, our results indicate that divergence times among the Archaeplastida inferred under the CIR model are largely congruent with those of prior studies. Maximum-Likelihood (ML), Bayesian inference (BI), and coalescent-based (CB) phylogenies for Archaeplastida are largely congruent, well-supported, and, overall, confirm relationships inferred during past studies (Additional file 1: Figure S1, Figure S2, and Figure S3). The ML and BI phylogenies have identical branching orders (Additional file 1: Figure S1 and Figure S2), and the only differences between the ML + BI and CB trees occur in Chlorophyta, the major green algal and seaweed clade. With the volvocine algae being our principal clade of interest, we elected to report general Archaeplastida phylogenetic inferences and divergence times in Additional file 3: Supplementary Methods.

Molecular clock analysis reveals that the crown Tetrabaenaceae-Goniaceae-Volvocaceae clade arose sometime between the Carboniferous and Triassic, and ancestral state reconstructions confirm our prior inference that multicellularity arose twice. According to our ancestral state reconstructions (ASR) (Additional file 1: Figures S5 and 6), multicellular lineages share a common unicellular ancestor ~298 MYA (CIR relaxed clock model, 95% HPD interval 349-237 MYA), marking the emergence of the crown Tetrabaenanceae-Goniaceae-Volvocaceae (TGV) clade sometime in the Carboniferous to the Triassic periods (Figs. 2 and 3, node A). Our 95% HPD interval for the TGV clade partly overlaps with ranges reported in previous studies [12, 18]. For example, the range reported by Herron et al. [18] is mostly confined to the Triassic (260-209 MYA), whereas that reported by Ma et al. [12] extends across the Early Triassic and Early Jurassic (251-180 MYA).

From the TGV clade’s common ancestor (Fig. 3, node A), our ASR results show that two unicellular lineages diverged (Additional file 1: Figures S4 and 5), each of which led to independent origins of multicellularity in the volvocine algae (Figs. 2 and 3). From one branch of this split, a multicellular ancestor of the Tetrabaenaceae diverged from a unicellular ancestor of Vitreochlamys ~267 MYA (316-206 MYA), and from the other branch a multicellular ancestor of the Goniaceae + Volvocaceae diverged from a unicellular ancestor of Chlamydomonas ~274 MYA (323-214 MYA (Fig. 3B and Fig. 3). These results confirm two independent origins of multicellularity among the volvocine algae as indicated by previous phylogenies [11, 12, 86]. According to our data, multicellularity arose in Tetrabaenaceae sometime between the Jurassic and Cretaceous (126-49 MYA), while multicellularity arose in the Goniaceae + Volvocaceae during the Permian-Jurassic (297-190 MYA).

Ancestral state reconstructions suggest that the steps leading to differentiated multicellularity in the volvocine algae occur in the order envisioned by Kirk. In 2005, David L. Kirk, a developmental and cell biologist, hypothesized how differentiated multicellularity, as manifest in Volvox carteri, might have evolved from a unicellular, Chlamydomonas-like ancestor [33]. In his essay, attention was brought to 12 ontogenetic traits (Fig. 3) that occur in the developmental cycle of V. carteri as it matures from a single-cell to a >500-celled organism. These 12 traits were then shown to be present in various extant Goniaceae and Volvocaceae, allowing them to be positioned in a stepwise manner along a simplified volvocine phylogeny, beginning with the most recent common ancestor of C. reinhardtii and V. carteri. Through the ordered mapping of these traits, Kirk suggested a developmental program that could explain how multicellularity with a division of labor might evolve. This hypothesis, known as Kirk’s 12-Step Program to differentiated multicellularity, has been formally tested once by Herron et al. [18] using ancestral-state reconstructions from a previous study [84]. Herron et al. [18] suggested that the 12-Step Program was not as orderly as hypothesized due to the phylogenetic position of the Tetrabaenaceae, which are absent from Kirk’s hypothesis. Below, we show that, using ancestral-state reconstructions based on our more robust phylogeny, Kirk’s 12-Step Program occurs temporally much as originally hypothesized.

For continuity, Kirk’s original order of 12 developmental traits has been retained in our analysis (Fig. 3). Approximately 298 MYA (349-237 MYA), two lineages emerged from the common unicellular ancestor of the colonial volvocine algae (Fig. 3, node A). One of these lineages eventually gave rise to the Tetrabaenaceae (Fig. 3). In this lineage, Steps 1, 3, 5 and 6 evolved in the multicellular ancestor of Basichlamys and Tetrabaena (Fig. 3 and Additional file 1: Figures S7-14). According to our ASR, this multicellular ancestor was established ~89 MYA (126-49 MYA); no other traits identified by Kirk are present either in this ancestor or in any of the known extant Tetrabaenaceae. Importantly, the number of Kirk’s traits gained in the lineage leading to the Tetrabaenaceae is twice that reported by Herron et al. [18] This conclusion is supported by recent studies confirming incomplete cytokinesis (Step 1) and basal body rotation (Step 3) in Tetrabaena [143, 144].

The second lineage diverged from a Chlamydomonas-like ancestor ~274 MYA (323-214 MYA), and eventually gave rise to the multicellular Goniaceae and Volvocaceae (Fig. 3, node B). ASR indicates that the common ancestor to the Goniaceae and Volvocaceae was multicellular (Additional file 1: Figures S5 and 6), that it had evolved 5 of the 12-Step Program traits over a period of ~27 MY (Fig. 3, node C) (Additional file 1: Figures S7-16), and was most likely non-spheroidal (Additional file 1: Figures S17 and 18). By ~247 MYA (297-190 MYA), the non-spheroidal, multicellular ancestors to the Goniaceae and Volvocaceae had evolved: (Step 1) an incomplete cytokinesis where daughter cells were left connected to each other via cytoplasmic bridges, (Step 3) peripheral colony cells with rotated flagellar basal bodies beating in parallel allowing the colony to efficiently mobilize through water, (Step 4) central-to-peripheral polarity similar to what is seen in Gonium, (Step 5) daughter cells embedded into an extracellular matrix with two distinct boundaries, and (Step 6) colony cell number controlled by its genetic code rather than the environment (Fig. 3).

For some volvocine algae, a developmental process known as “inversion” must occur. During development, certain bowl-shaped and spheroidal colonies have their flagella pointing inside the bowl or spheroid, preventing the colony from being able to swim. The process of inversion flips the curvature for bowl-shaped volvocines and turns spheroidal colonies inside out; thus, allowing the flagella to provide motility to the colony. In the volvocine algae, inversion occurs in the following three states: i) does not occur, ii) partial inversion, and iii) complete inversion. For unicellular volvocines, the Tetrabaenaceae, and Astrephomene, inversion does not occur in their developmental cycle. Due to the way volvocine cells divide, embryos resulting in Gonium and Volvocaceae colonies must undergo “partial” and “complete” inversion, respectively, so that the eventual colonies have the ability swim. For Gonium, embryos are in the shape of a bowl with their flagellar ends positioned in the concave region, and volvocacean embryos, in similar fashion, have their flagellar ends located inside the developing colony. By ~146 MYA (206-106 MYA), the ancestor to Gonium had developed (Step 2) “partial” inversion whereby a curvature reversal in the embryo occurs, allowing the flagellar ends to be situated on the exterior of the colony. By ~228 MYA (279-175 MYA), volvocacean embryos developed "complete" inversion (Step 7), where cells, in a coordinated fashion, move and curl outward along the anteroposterior axis, resulting in an embryo whose flagellar ends point outside the colony (Additional file 1: Figure S19 and 20) (Step7). In both examples, cells in colonies are held in position by cytoplasmic bridges formed via incomplete cytokinesis (Step 1). Along with Step 7, by ~228 MYA (279-175 MYA) spheroidal volvocacean ancestors had evolved body plans with increased volume. Diminished colony volume occurred secondarily in lineages leading to Platydorina and the paraphyletic Pandorina (Fig. 3).

The spheroidal body plan in the Volvocaceae and in the lineage leading to Astrephomene likely consisted of cells positioned along a colony’s periphery, resulting in a transparent organism with an interior filled with extracellular matrix. Despite the morphological similarity among spheroidal colonies, earlier studies [145, 146] and our ASR results indicate that such a body plan arose independently twice in the volvocine algae (Additional file 1: Figures S17 and 18). Yamashita et al. [145] elucidated the developmental mechanisms that underlie the two independently evolved forms. During Astrephomene embryogenesis daughter protoplasts gradually rotate their apical ends so that flagella extend towards the posterior of the growing colony after each successive cell division. Development proceeds differently in the Volvocaceae. In volvocacean early embryos, chloroplasts are positioned at one end of each protoplast, facing outward. As the embryo matures, these cells change shape, with the chloroplast ends forming acute edges [146]. This cell shape change results in bending the epithelium at an opening of the embryo known as the phialopore [147]. Epithelial bending causes the colony to turn itself inside-out allowing flagella to be located on the exterior of the developing colony. This process in the Volvocaceae is known as “complete” inversion [145]. Yamashita and Nozaki [146] showed that neither of these developmental mechanisms occur in Tetrabaena or Gonium embryogenesis and thus concluded that spheroidal body plans in the volvocine algae are independently derived via separate mechanisms.

In our updated phylogeny, the most recent common ancestor of the Volvocaceae exhibited 7 of the 12 traits identified by Kirk (Fig. 3). Moreover, there is only a single conflict in the order in which they evolved according to his hypothesis. This conflict arises from ambiguity as to when inversion first evolved in a volvocine ancestor; indeed, ASR yields alternative inferences depending upon which analytical assumptions are used. When we consider this conflict as an unordered 3-state problem (no inversion, “partial” inversion, or “complete” inversion) (Additional file 1: Figure S19-20), which assumes that no state is a prerequisite for another to evolve, the ancestor to the Goniaceae and Volvocaceae is shown not to have evolved inversion (Fig. 3, node C). However, when the analysis is performed as an “ordered” 3-state problem, where “partial” inversion is required for “complete” inversion to evolve, the common ancestor to the Goniaceae and Volvocaceae is shown to have partial inversion (Additional file 1: Figure S21). Lastly, when inversion is treated as a 2-state problem (e.g., no inversion, inversion), the common ancestor of the two clades is shown to have evolved inversion (Additional file 1: Figure S22-23). This last result conflicts with Yamashita and Nozaki’s [146] explanation for how inversion in Gonium and the Volvocaceae are mechanistically distinct. Inversion should be treated as different states (e.g., “partial” and “complete”), as in our first example of an “unordered” 3-state problem, because it is unknown whether “partial” inversion is necessary for “complete” inversion to evolve. Thus, we report that only the ancestor to Gonium had evolved “partial” inversion (Step 2) by ~146 MYA (206-106 MYA).

The evolution of somatic cellular differentiation (Step 9) is the first step towards a complete reproductive division of labor where a specific cell line (the soma) evolved to house cells (germ-line) that transmit genetic information to subsequent generations. Somatic cell lines observed in Astrephomene, Pleodorina, and Volvox are distinct in that they assume specific cellular functions, notably motility, exhibit a finite lifespan, and do not pass on genetic information. Two other cell types are present in the volvocine algae: undifferentiated cells and fully-differentiated germ cells. Undifferentiated cells, the ancestral cell state in volvocine algae, participate in motility and can undergo mitotic division to produce the next generation. By contrast, fully-differentiated germ cells, observed only in Volvox, also undergo mitotic division to generate progeny, but these cells never significantly participate in motility. By comparing our ASR results to our phylogenetic results we can infer specific number of independent origins of sterile somatic cells. Thus, the gain of sterile somatic cells (Step 9) occurred five times in three different genera of the volvocine algae (Fig. 3). According to our results, true somatic differentiation arose in (i) Pleodorina thompsonii sometime after it diverged from the common ancestor it shared with V. carteri (169-99 MYA) and in the common ancestors of: (ii) Astrephomene (170-75 MYA), (iii) section Volvox (136-57 MYA), (iv) the Volvox + Pleodorina japonica clade (166-98 MYA), and (v) P. starrii + V. gigas (140-82 MYA) (Fig. 3 and Additional file 1: Figures S24 and 25). In addition to this trait being gained, it was also lost once in the Eudorina clade sister to V. gigas + V. powersii.

The advent of true germ-soma differentiation in the volvocine algae is characterized by evolution of specialized germ cells (Step 10) that do not significantly contribute to colony motility. In genera like Astrephomene and Pleodorina that have not made this step, colonies consist of somatic cells as well as biflagellate gonidia, cells that contribute to motility and can undergo mitotic division. Indeed, in these taxa gonidia greatly outnumber somatic cells. In A. gubernaculifera a 32 or 64-celled colony may have 30 or 60 gonidia [148], respectively, and 128-celled P. californica colonies can have >70 gonidia [149]. Gonidia in these instances contribute significantly to motility until cell division is initiated, which ultimately leads to hatched daughter colonies and subsequent death of the mother. In contrast, Volvox colonies are all composed of >500 somatic cells and far fewer flagellated or aflagellate gonidia. For example, at the end of V. rousseletii’s embryogenesis, a colony may consist of ~2000-14,000 biflagellate, nearly indistinguishable cells with only 1-16 gonidia [17, 150]. In V. rousseletii, gonidia contribute to motility for <20% of a colony’s asexual life cycle before resorbing their flagella and undergoing mitotic division [17]. Volvox carteri colonies are typically composed of ~2000-6000 cells with ~8 gonidia [150]. Somatic and gonidial cells are clearly defined during V. carteri’s embryogenesis. Gonidial cells are aflagellate; thus, in this species gonidia never contribute to colony motility during its life cycle [17]. Therefore, Volvox spp. are the only extant volvocine algae to have evolved specialized germ cells (Step 10). In this group, specialized germ cells evolved by ~96 MYA (136-57 MYA) in section Volvox, ~46 MYA (68-28 MYA) in the ancestor to V. gigas and V. powersii, and ~129 MYA (166-98 MYA) in V. aureus and V. carteri’s most recent common ancestor (Fig. 3 and Additional file 1: Figures S24 and 25). This trait was lost once in P. japonica (Fig. 3 and Additional file 1: Figures S24 and 25).

The final two steps of the 12-Step Program are coupled, consisting of asymmetric cell division (Step 11) and a bifurcated cell division program (Step 12) (Additional file 1: Figures S26-29). Both traits evolved twice in the genus Volvox: (i) V. africanus sometime after ~106 MYA and (ii) ~98 MYA for the ancestor of V. obversus and V. carteri. During V. carteri embryogenesis, anterior embryo cells undergo asymmetric divisions producing both small and large cells. Size controls cell fate during embryogenesis. Cells <8μm may divide up to 6 times before they arrest, whereas larger cells undergo 1-2 asymmetric cell divisions. Smaller cells, and the remaining embryonic cells, eventually become soma, whereas the larger cells are destined to become specialized germ cells. Through these two coupled processes (Steps 11 and 12), complete division of labor is achieved between somatic and germ cell lines.

Ancestral state reconstruction shows that evolution did not progress linearly within the volvocine algae from a unicellular organism like Chlamydomonas to one that is multicellular and differentiated like V. carteri. Numerous trait gains and losses occurred in the history of this clade, exemplifying the non-linear nature of evolution. However, Kirk’s 12-Step Program does not attempt to linearize the evolutionary history of this group, but rather to provide a roadmap of steps needed for true differentiated multicellularity to evolve in the volvocine algae. Our results indicate that the temporal sequence in which Kirk organized his 12 traits aligns well with the inferred evolutionary history of this group.

Ancestral-state reconstructions support previous findings that multicellularity drives evolution of sexual traits in the volvocine algae. The volvocine algae in their various forms have given rise to multicellular lineages that express a variety of sexual traits. For this group, asexual reproduction via mitosis is the typical means by which the next generation of unicells or colonies is produced. Sexual reproduction is facultative, and is precipitated by reduction in environmental nitrogen and, in some cases, a sex-inducing hormone [151]. Regardless, a diploid zygote is formed after fertilization, resulting in a spore. Meiosis only occurs in the volvocine algae when spores germinate to produce haploid progeny [151].

Unicellular forms, such as Chlamydomonas (sensu Pröschold et al. [152]), and in some multicellular forms (Astrephomene, Basichlamys, Gonium, Pandorina, Tetrabaena, Volvulina, and Yamagishiella), plus and minus type gametes are morphologically identical (isogamy). However, in many multicellular forms (Colemanosphera, Eudorina, Platydorina, and Pleodorina) the two gamete types differ in size, with the smaller of the two recognized as being male (anisogamy) [153]. In certain instances (Volvox), gametic differentiation may be further exaggerated, with the small male gamete being motile, and the large female gamete being sessile (oogamy). One longstanding question is: why does anisogamy evolve in the first place? Another is: what role does multicellularity play in gametic differentiation? Because the volvocine algae encompass both unicellular and multicellular forms, and because they exhibit the full range of known gamete types, they are especially well-suited to addressing these two fundamental questions.

Hanschen et al. [32] systematically investigated whether multicellularity is a prerequisite for the evolution of anisogamy and whether multicellular volvocines exhibit greater complexity in the form of more elaborate sexual traits. Their study, like ours, concludes that multicellularity is correlated with the acquisition of anisogamy, as all anisogamous volvocines are multicellular, and anisogamy appears to evolve from isogamous ancestors (Fig. 4) [32]. Using our updated phylogeny, with its newly inferred relationships, we conducted ASR for 7 sexual traits identified in previous studies [31, 32]. Except where noted (Fig. 4), we focused on the same strains used in the two prior studies of how sexual traits are distributed among the volvocine algae. Like those earlier studies, we excluded Vitreochlamys, Pleodorina thompsonii, and Volvox ovalis from certain ASR analyses due to a lack of information about how they reproduce sexually.

Fig. 4.

Sexual trait gain and loss in the volvocine algae. Branching order is of a collapsed tree inferred under maximum-likelihood analysis from an aligned amino acid, concatenated dataset of 263 nuclear genes. Numbers on branches represent bootstrap and posterior probability values, respectively. Branch lengths, corresponding to time, were inferred under the CIR relaxed clock model using 16 most clock-like genes as determined by Sortadate. Blue bars correspond to the inferred 95% HPD interval for each node. Green bubbles correspond to a sexual trait gain, and red bubbles corresponds to loss of a trait. The figure table lists the 7 mapped sexual traits and their acronyms. Pink branches denote lineages where anisogamy evolved. Members of the multicellular volvocine algae (Tetrabaenaceae, Goniaceae, and Volvocaceae) are denoted in orange, purple, and green. Taxa in black font are unicellular volvocine algae

Unicellular ancestors of the volvocine algae, like Chlamydomonas, express the following sexual phenotypes: (i) isogamy, (ii) full meiotic hatching of zygospores, (iii) dioecy, (iv) outcrossing, (v) external fertilization of the egg, and no evidence of sexual dimorphisms such as (vi) “extrafertile” females or (vii) “dwarf” males. Prior to the evolution of anisogamy, the number of "gone cells” released from a zygospore decreases. (Note: a gone cell constitutes the haploid meiotic product arising from a diploid zygospore [154].) Chlamydomonas reinhardtii, Tetrabaena, and Gonium produce all four haploid gone cells per diploid zygospore [88]. According to our ASR results full meiotic hatching is lost either once or twice in the volvocine algae (Additional file 1: Figures S30 and 31) Our MrBayes ASR indicates a single loss, but this could be an artifact arising from low sampling of taxa that exhibit full meiotic hatching. In contrast, and consistent with Hanschen et al. [32], our Phytools ASR indicates two losses of full meiotic hatching in the following lineages: once in Astrephomene sometime after diverging from Gonium ~220 MYA (274-170 MYA), and again in the Volvocaceae ~228 MYA (279-175 MYA) (Fig. 4 and Additional file 1: Figures S30 and 31). Each loss represents a reduction in meiotic products from four haploid gone cells to a single haploid gone cell and three polar bodies. Among our identified sexual traits, no other gain or loss occurs until the evolution of anisogamy.

The most significant difference between our findings and those of Hanschen et al. [32] is the number of times we find anisogamy to have evolved independently. Their results inferred two independent origins of anisogamy from isogamous ancestors, whereas we conclude three, consistent with Lindsey et al. [11] (Fig. 4 and Additional file 1: Figures S32 and 33) [11, 32]. Section Volvox’s position as sister to all other Volvocaceans underlies this distinction. The three origins of anisogamy in our analysis occur in the ancestors of the following lineages: (i) section Volvox (136-57 MYA), (ii) Platydorina + Colemanosphaera (113-46 MYA) and (iii) the Eudorina + Volvox + Pleodorina (EVP) clade (196-117 MYA) (Fig. 4 and Additional file 1: Figure S32 and 33). It is noteworthy that we infer three independent origins of anisogamy, regardless of whether this trait is treated as a 2-state problem (e.g., isogamy and anisogamy) or a 3-state problem (e.g., isogamy, anisogamy, or oogamy) (Additional file 1: Figures S32-35).

Consistent with our mapping of developmental traits, anisogamy evolved in lineages that had already acquired a diverse suite of characters differing from a unicellular ancestor. In each case, the eventual anisogamous ancestor had evolved a spheroidal body plan with increased colony size, polarity, and in the case of section Volvox, may also have acquired true germ-soma differentiation. It is possible that in some anisogamous ancestors, the body plan was larger with distinct cell types compared to a Pandorina-like or Volvulina-like morphology, as hypothesized by Hanschen et al. [32] In section Volvox and the EVP clade, internal fertilization was gained around the same time (sometime after ~228 MYA (279-175 MYA) for section Volvox and ~152 MYA (196-117 MYA) for the EVP clade) as anisogamy, suggesting that anisogamy might be a prerequisite for this trait (Fig. 4 and Additional file 1: Figures S36 and 37).

After gaining anisogamy and internal fertilization, oogamous lineages appear rapidly to gain several other sexual traits (Fig. 4). This suggests that oogamy, an extreme form of anisogamy, may drive further sexual evolution. By ~95 MYA (136-57 MYA), section Volvox’s ancestor had become monoecious (i.e., colonies could produce both sperm and ova) and capable of selfing. Necessarily, such an ancestor had lost dioecy (i.e., colonies able to produce only sperm or ova) as well as the capacity to outcross (Additional file 1: Figures S38-41). By the same time, this ancestor gained a sexual dimorphism in the form of “extrafertile” females (Fig. 4). During sexual reproduction, an “extrafertile” female Volvox colony exhibits a doubling in the number of eggs inside its colony, thus giving sexual reproduction the potential to numerically outperform an asexual reproductive cycle.

MrBayes and PhyTools ASR results indicate different ancestral states for when the “dwarf” male phenotype was gained. Results from MrBayes ASR (MBASR) indicate “dwarf” males were gained independently in the lineages leading to V. africanus, V. tertius, and in the common ancestor to V. carteri f. nagariensis and V. carteri f. weismannia (Additional file 1: Figure S44). Phytools ASR, however, indicate that the common ancestor to V. carteri and V. africanus had evolved the “dwarf” male phenotype with two losses occurring in the lineages leading to V. obversus and V. dissipatrix (Additional file 1: Figure S45). These differences arise due to model selection. MBASR results are replicated in Phytools if the “equal rates” model is exclusively used (Additional file 1: Figure S46). When an ANOVA test is performed in Phytools comparing the “equal rates” and “all-rates-different” models, the “all-rates-different” model has a markedly higher weight (0.70 vs 0.29) indicating that it is highly favored over the “equal rates” model. Therefore, in the common ancestor to V. carteri and V. africanus, two sexual dimorphisms, “extrafertile” females and “dwarf” males, were gained by ~118 MYA (152-90 MYA) (Fig. 4 and Additional file 1: Figures S42-43 and 45). The “dwarf” male phenotype is expressed as an exaggeratedly reduced male colony size when compared to female colonies. Within this same clade, V. africanus and V. dissipatrix independently evolved monoecious, selfing colonies, in addition to the previously mentioned sexual dimorphisms. Sexual trait gains occur in various other lineages of the volvocine algae, but none are as clear as the ones occurring in the two oogamous lineages of Volvox.

In summary, we conclude that multicellularity was a likely driver to the evolution of anisogamy. Once anisogamy is established, a variety of other sexual traits may arise, notably the evolution of oogamy, which in turn opens the door for more elaborate forms of sexual dimorphism such as extrafertile females and dwarf males. These derived sexual traits are positively associated with traits related to the degree of multicellular complexity, such as increased cell number, expanded extracellular matrix, organismal polarity, inversion, and cellular differentiation [32, 33]. Thus, our more robust phylogeny of the volvocine algae, based on not 5 but on 263 protein coding loci, confirms previous findings [31, 32] that multicellularity ultimately drives the evolution and elaboration of diverse sexual traits.

We would be remiss if anisogamous/oogamous unicellular chlamydomonads outside the volvocine algae were not briefly discussed. Using morphological and SSU rRNA data, Pröschold and colleagues have proposed placing oogamous Chlamydomonas strains in two other genera: Oogamochlamys and Chloromonas [152]. In their taxonomic revision, Chlamydomonas reinhardtii is proposed as the type species of the genus [152]. Pröschold et al. [152] contend that because these other Chlamydomonas species are not monophyletic with C. reinhardtii (i.e., do not fall within the Reinhardtinia clade), they should be reclassified from Chlamydomonas to other genera [152]. Thus, while we recognize that anisogamy and oogamy exist in unicellular genera outside the volvocine algae, discussion of these taxa is beyond the scope of our investigation.

Molecular clock analysis also suggests cryptic volvocine species in the Tetrabaenaceae. Much like the longstanding debate on what constitutes a biological species [155, 156], consensus is lacking for how we should recognize cryptic species [157]. Most simply, cryptic species are those that have been binned into a single taxon at the species rank. Discriminating cryptic species can enlarge estimates of biodiversity and provide insight into the mechanisms that underlie adaptation and speciation. Following suggestions and criticisms offered by Struck et al. [158], we provide here divergence times and molecular evidence that point to the existence of cryptic species within the genus Tetrabaena. This insight is reinforced by previously published biogeographic [82, 159], phenotypic [159], and genetic data [160] on this enigmatic genus.

Tetrabaena is a monotypic, colonial volvocine alga and one of two genera composing the multicellular Tetrabanaeceae. Several strains of this genus have been described as far back as 1841 [161]; all described Tetrabaena spp. have been of freshwater origin except one Antarctic strain, NIES-691. Following collection of the Antarctic T. socialis from King George Island [162], Nozaki and Ohtani [159] provided morphological data and temperature growth profiles that demonstrated phenotypic idiosyncrasies of this strain in comparison to the freshwater NIES-571. Phenotypic traits distinctive to NIES-691 included: (i) vegetative colonies measuring up to 50 μm in diameter (~18 μm larger than NIES-571), (ii) immobile colonies that adhere to surfaces at 20 °C, and (iii) failure to grow at 25 °C. In contrast to NIES-691, NIES-571 was observed to produce smaller colonies and to grow normally at all tested temperatures ranging from 5-25 °C, with no phenotypic changes that produce immobile colonies or ones adhering to surface s[159]. These phenotypic differences, which may be local environmental adaptation, point towards the possibility of a cryptic species in Tetrabaena, arguing for the examination of detailed molecular data on these two strains.

In 1997, Mai and Coleman [160] produced complete Internal Transcribed Spacer 2 (ITS-2) sequences for 111 Volvocales. including Tetrabaena socialis strains NIES-571 and NIES-691. When these two T. socialis ITS-2 sequences are subjected to pairwise alignment, the percent identity between them is 92.47% (Additional file 2: Table S3). Similar to the percent identity (91.25%) between ITS-2 sequences of V. gigas UTEX 1895 and V. powersii UTEX 1864, two sister taxa recognized as separate species, when these are subject to pairwise alignment (Additional file 2: Table S3). When the same test is applied to ITS-2 sequences of Volvox carteri f. nagariensis strains (Poona and 72-52), the percent identity between them is 98.06% (Additional file 2: Table S3). Considered in relation to previously reported phenotypic data, these molecular differences lead us to hypothesize that NIES-571 and NIES-691 are cryptic species within the genus Tetrabaena.

This hypothesis is bolstered by our clock analysis and molecular data. The two strains appear to have diverged from one another ~51 MYA (82-25 MYA) (Fig. 3). This divergence time estimate stands in contrast to that estimated between strains of Volvox carteri f. nagariensis (NIES-865 and HK10) as well as to that estimated between strains of Pleodorina starrii (NIES-1362 and NIES-1363). In each of these two examples, strains are thought to have diverged from one another <1 MYA (Fig. 3). And in both cases, there is good reason to believe that each strain pair is a biological species sensu May r[163]. V. carteri NIES-865 is an isolate of V. carteri HK1 0[164, 165], and mating experiments [166] between Pleodorina starrii strains (NIES-1362 and NIES-1363) produce fertile offspring. Exploration of conserved regions from our 16 amino acid sequence alignments used in our molecular clock analyses further highlight differences in sequence percent identity. Pairwise percent identity between V. carteri and P. starrii strains was 100% and 99.65-100%, respectively, whereas the percent identity between each set of T. socialis sequences ranges between 88-99.27% indicating substantial evolutionary differences in the protein sequences (Additional file 2: Table S4).

Given the divergence times indicated by our dataset, and the phenotypic and molecular features known to distinguish these T. socialis strains, we propose recognizing Tetrabaena socialis NIES-571 as a cryptic species to clarify socialis N-691’s current status as the type species. Formal elevation of NIES-571 to species status will require additional data. For example, complete plastid sequences should be obtained to complement the pairwise alignments of nuclear genes provided here. Because chloroplast genomes are known to be highly conserved [167], marked differences between plastid sequences of these two T. socialis strains would strengthen the case for one to be recognized as a cryptic species. Also, NIES-571 and NIES-691 should be crossed. And, even if mating experiments result in viable F1 progeny, additional crosses should be undertaken to assess the fertility of F1 hybrids in relation to their parents and to one another.

Conclusions

Through the use of a more comprehensive molecular and taxonomic dataset, we have confirmed the major findings of Lindsey et al.’s phylogeny [11], and further supported them via ancestral-state reconstructions using fossil-calibrated molecular clock analyses. We find that the origin of the crown Tetrabaenaceae-Goniaceae-Volvocaceae (TGV) clade may be much older than previously thought, possibly emerging as early as the Carboniferous Period. We also find that multicellularity in the Goniaceae + Volvocaceae evolved sometime between the Carboniferous-Triassic Periods, and that multicellularity in Tetrabaenaceae may have evolved as recently as the Cretaceous Period. Given the wide range of time in our 95% HPD intervals for volvocine divergence times (e.g., second independent origin of multicellularity occurring between 297-190 MYA), we hesitate to speculate about specific geological and ecological changes that may have attended the evolution of multicellularity in this clade. However, even a conservative view of our divergence times shows that volvocine diversification occurred in the context of ecosystem reorganization in the wake of end-Permian mass extinction [168]. Armed with our time-tree and ancestral-state reconstructions, we reassessed Kirk’s 12-Step Program to differentiated multicellularity using a new, more robust phylogeny of the volvocine algae; we determined that the essential elements of his hypothesis to be correct. Lastly, our data point to the possibility of a cryptic species in the multicellular Tetrabaenaceae, a possibility that could change its monotypic status, and provide an exciting opportunity to investigate the genomics of two diverged species that have adapted to very different environments.

In ancestral character state reconstructions, it is nearly always the case that multiple, nearly equally likely reconstructions exist, and this limits the confidence we can have in the specifics of the (marginally) most parsimonious. For example, in the clade at the bottom of Fig. 4, which includes P. indica, P. starrii, V. powersii, V. gigas, and some E. elegans strains, our reconstruction shows a gain of sterile soma at the base of the clade and a loss in E. elegans. However, it would be only slightly less parsimonious to show independent gains in P. indica + P. starrii and in V. powersii + V. gigas, with no loss. There is, however, no way to reconstruct the history of this trait that doesn't involve several independent gains and/or losses.

The distinction between mortal somatic cells, which do not pass on genetic information to subsequent generations, and potentially immortal germ cells, which do, has been recognized as a fundamentally important evolutionary trait as far back as Weismann's germ-plasm theory [169]. The evolution of anisogamy is the first step in the differentiation of male and female sexes, setting the stage for sexual selection to produce all of the spectacular sexual dimorphisms we see throughout the animal kingdom [14, 170, 171] as well as the more modest dimorphisms in some species of Volvox (dwarf males and extrafertile females). Although nearly equally likely reconstructions preclude a high degree of confidence in the details of these traits' histories, this much is clear: given the long-recognized fundamental evolutionary significance of these traits, they are surprisingly evolutionarily labile, with multiple independent losses and/or gains, within the relatively recently-evolved multicellular volvocine algae.

Methods

Command-line arguments and software repository links for each program have been uploaded to a single file that may be downloaded from our Dryad repository (DOI: 10.5061/dryad.mcvdnck6b) [172].

Retrieval of genomic and RNA-Seq data for Archaeplastida phylogeny. A total of 164 taxa were sampled across Archaeplastida based upon publicly available genomic and RNA-Seq data. Specifically, longest primary transcript files were downloaded for 20 taxa from the Phytozome database, and protein files for 11 taxa were downloaded from Ensembl, NCBI, OrcAE, and Tokyo Institute of Technology. RNA-Seq data was downloaded for 133 taxa, and de novo transcriptomes were assembled for each taxon. Accession numbers and/or database information for each strain may be found in Additional file 2: Table S1.

Quality control of RNA-Seq reads. Quality of all RNA-Seqs was initially checked using FastQC v.0.11.8 with an additional FastQC assessment post-trimming. Quality trimming was conducted with Trimmomatic v.0.3 9[135]. Command line parameters and details for Trimmomatic in addition to adapter file information may be found in Additional file 3: Supplementary Methods.

De novo transcriptome assembly. Quality filtered, paired-end reads were used to assemble de novo transcriptomes with SOAPdenovo-Trans v1.0.4 [173] using k-mer size of 25, GapCloser was used to close gaps in each de novo transcriptome, and CDHIT v4.8.1 [174] was used under default parameters to cluster redundant transcripts from each assembled transcriptome.

Annotation of de novo transcripts. Each assembled transcript had its longest open reading frame predicted using TransDecoder v5.5.0. Diamond [136] BLASTP and HMMER [137] protein domain identification searches were conducted using each .pep file generated in the previous step. Using the results of Diamond and HMMER, nucleotide CDS and amino acid PEP files were then generated by Transdecoder v5.5.0. Detailed information for Diamond BLASTP and HMMER searches may be found in Additional file 3: Supplementary Methods.

Orthologous gene identification infers 284 single-copy genes shared among red and green algae. A training set of 10 Archaeplastida taxa (1 rhodophyte, 5 streptophytes, and 4 chlorophytes) were chosen from the Phytozome database to infer a putative list of single-copy genes present across Rhodophyta, Streptophyta, and Chlorophyta (Additional file 2: Table S5). The primary transcript file for each taxon in our training set was used to determine a putative list of single-copy genes using Orthofinder v2.3.8 [175]. Orthofinder returned 284 orthologous groups where all species were present and all genes single-copy. These 284 orthologous genes across all species were aligned using MAFFT v7.487 [176]. HMMER was subsequently utilized to build profile Hidden Markov Model (HMM) files for each gene alignment. These profile HMMs were then used by Orthofisher [177] to mine our genomes and de novo transcriptomes for those single-copy genes. Putative organellar sequences from the inferred single-copy dataset were ruled out, and detailing about this process can be found in Additional file 3: Supplementary Methods.

Gene sequence alignments. Once the amino acid sequences were confirmed to be of nuclear origin and single-copy, they were aligned using MUSCLE v3.8.31 [178], resulting in 284 untrimmed gene alignments. Furthermore, each alignment was manually inspected and aligned in Aliview v1.26 [179]. Further details are found in Additional file 3: Supplementary Methods. Following manual quality control of alignments, ambiguously aligned regions from each alignment were trimmed using Trimal v1.4.1 [180] allowing for only conserved and reliably aligned regions to later undergo phylogenetic analysis. A total of 21 gene alignments from the original 284 alignments were discarded and not used in any downstream analyses. These 21 gene alignments were discarded due to poor sequence alignment and/or many missing sequences.

Parameters for Phylogenetic analyses. Each single-gene alignment had its best evolutionary model predicted by ProtTest v3.4.2 [181] under the Akaike Information Criterion (AIC). Each predicted evolutionary model and its associated alignment may be found in Additional file 2: Table S6. All single-gene alignments underwent Maximum-Likelihood (ML) analysis using IQtre e[182]. Phyutility v2.7.1 [183] was utilized to concatenate all gene alignment files for ML and BI analyses. IQtree was used to run a concatenated, gene-partitioned ML analysis with 1000 rapid bootstrap replicates [184]. The gene-partition strategy was based upon the best evolutionary model predicted by ProtTest. The Bayesian inference (BI) analysis was conducted with MrBayes 3.2.7a [142] under the same gene-partitioning strategy previously described. A total of 4 runs each with 4 Markov chains (3 heated and 1 cold). Trees were sampled every 1000 generations over 1,000,000 generations during which convergence (average standard deviation of split frequencies = <0.01) between runs was reached, and a burn-in of 25% was used (ngen=1000000 nruns=4 samplefreq=1000 nchains=4 starttree=random). Lastly, ASTRAL [185] was utilized to perform a coalescent-based analysis using the 263 single-gene trees produced by IQtree. Each single-gene tree had branches with low bootstrap support (<10) collapsed using Newick Utilities v1.6 [186].

Parameters for molecular clock analyses. For all clock analyses, Sortadate [138] was utilized to assemble concatenated datasets of 8 and 16 genes to offset high computational cost. For each of the smaller datasets, genes were ranked with an emphasis on their Sortadate bipartition score with the scores ranging from: 48-56% (8 genes) and 45-56% (16 genes) (Additional file 2: Table S7). Additional details concerning Sortadate may be found in Additional file 3: Supplementary Methods.

Phylobayes 4.1b [92] was used to perform relaxed molecular clock analyses on each dataset under the CIR (-cir), lognormal autocorrelated (-ln), uncorrelated gamma (-ugam), and white-noise (-wn) models with the LG protein substitution model applied (-lg), discrete gamma rate across sites (-dgam 4), and birth-death prior on divergence times (-bd). The red algae were specified as the outgroup taxa according to manual instructions (-r <outgroup.txt>). For our main analyses a total of 14 nodes were calibrated using fossil taxa with either a minimum bound or minimum + maximum bounds (Fig. 2) with the soft bounds flag (-sb) allowing date sampling to occur slightly outside of the specified bounds (5% on a pure minimum date, and 2.5% on lower and upper bounds). No secondary node calibrations were used in this study. The root age (red and green algal divergence) was set to a maximum bound of 2000 MY. Separate runs for our maximum-likelihood and coalescent-based species tree topologies were conducted. A total of two MCMC chains were conducted per main run for a minimum of 10,000 generations, and the convergence of the two chains was determined in Phylobayes. A summarization of the chains was done after discarding the first 2000 generations (~20%) as burn-in.

Ancestral state reconstructions of discrete volvocine characters Two programs were used to estimate the probability of ancestral character states at internal nodes in the volvocine algae: the Phytools [139] R package and MrBayes [142] via the MBASR toolkit [141]. For each of the two programs, our tree topology was constrained to the maximum-likelihood, molecular-clock tree inferred using 16 protein-coding genes under the CIR model. For each ASR method, the following discrete character states were ancestrally reconstructed in the volvocine algae: (i) cellularity, (ii) spheroidal body plan, (iii) developmental traits identified by Kirk in his 12-Step Program [33], and (iv) 7 sexual traits from two previous studies [31, 32]. To estimate the probabilities of ancestral character states, stochastic character mapping ("simmap”) was applied using Phytools. R scripts used to run each simmap analysis may be found at Dryad [172]. In MBASR, a Markov chain Monte Carlo (MCMC) simulation was conducted for 1,000,000 generations (MBASR n.samples=10,000) with a sample taken every 100 generations to estimate ancestral character state probabilities. Additionally, information for runs conducted in Phytools and MBASR may be found in Additional file 3: Supplementary Methods.