Order, please! Uncertainty in the ordinal-level classification of Chlorophyceae

Bayesian analyses

Sequences of 58 protein-coding chloroplast genes were aligned using the translation-aided algorithm in Geneious v.10 (Biomatters) using the bacterial code and the remaining parameters set at default. Sites and regions of uncertain homology in variable genes were masked and trimmed prior to analyses using a custom Python script. We provide the masked alignments in supplementary files; the deploy.py script is available as supplement to Fučíková, Lewis & Lewis (2016). Phylogenetic trees were constructed using MrBayes v.3.2 (Ronquist et al., 2012) using the GTR+I+ Γ model over 6,250,000 generations and two MCMC chains in each of two parallel runs. Sites were partitioned by codon position across the entire data set. The first 20% of each run was discarded as burn-in. The data were analyzed in the amino acid (aa) form as well, implementing the aa GTR model in MrBayes.

Single-gene data sets were analyzed using MrBayes v.3.2 (Ronquist et al., 2012) as described above. A subset of these analyses (the 37 genes comprising the entire set of 68 taxa –i.e., with no missing data) was used to assess the degree to which individual genes support each clade in the phylogeny. One hundred randomly sampled post-burnin trees per gene were used for this assessment (3,700 trees total). The proportion of all 3,700 trees supporting each clade (PMT, or Proportion of Merged Trees), as well as the internode certainty (IC; Salichos & Rokas, 2013; Salichos, Stamatakis & Rokas, 2014), were calculated for the clades in Fig. 2A and plotted in Fig. 3.

ML analyses

Maximum Likelihood (ML) analyses were carried out on the concatenated nt and aa data sets to complement the Bayesian analysis results. RAxML (Stamatakis, 2014) was chosen as analysis tool because of the size of our data set, but because RAxML does not allow the GTR amino acid model (which we consider the most realistic and therefore chose to implement it in MrBayes), we instead used the LGF model for the amino acid analysis, allowing for a gamma distribution of rates among sites (PROTGAMMALGF). The LGF model was selected using the script ProteinModelSelection.pl available from the RAxML webpage (https://cme.h-its.org/exelixis/web/software/raxml/hands_on.html). Two hundred bootstrap pseudoreplicates were carried out to assess branch support.

Supplementary phylogenetic analyses

Because the initial concatenated analyses yielded conflicting topologies (nt vs. aa, Fig. 2), we used Shimodaira’s Approximately Unbiased (AU) test (Shimodaira, 2002) within PAUP* (Swofford, 2002) to see whether the two topologies were significantly different from each other. The test and its results are described in full in the Supplements.

On the concatenated data, a suite of analyses was conducted to probe the effects of analysis type and model selection on the final topology. To allow site-specific evolutionary processes in tree inference, we analyzed the nt data using SVDquartets (Chifman & Kubatko, 2014; Chifman & Kubatko, 2015) with 200 bootstrap pseudoreplicates implemented in PAUP* (Swofford, 2002).

We also carried out a supplementary set of Bayesian analyses with ambiguously recoded Serine, Isoleucine and Arginine codons (as described in Fučíková, Lewis & Lewis, 2016) to probe the possible effects of convergent codon usage. Further, we analyzed the 1st and 2nd positions only, and in a separate analysis the 3rd positions, to determine whether there was conflict between these data partitions. The data sets with analysis specifications and the resulting consensus trees are in the Supplements. The full concatenated data set was also analyzed using PhyloBayes v.4.1 and implementing the CAT-GTR model, allowing (in addition to gamma-distributed rates) site-specific substitution processes and thereby attempting to further mitigate branch attraction and other systematic bias issues in the data (Lartillot & Philippe, 2004; Lartillot, Lepage & Blanquart, 2009). The commonly used 18S nuclear ribosomal gene was analyzed to provide a chloroplast-independent estimate of the phylogeny. Because of the unique alignment issues associated with this gene, we selected the program BAli-Phy v.3, which estimates the optimal alignment as well as the phylogeny (Suchard & Redelings, 2006). The GTR+I+ Γ model was implemented and 12,000 iterations were run after pre-burnin. Parameter stability was checked using Tracer v1.7 (Rambaut et al., 2018) and 20% of the run were discarded as burnin. Lastly, an analysis combining the 18S and plastid nucleotide data was conducted, and is described in full in the supplementary methods.

We used TREESPACE (Jombart et al., 2017) to visualize the variability of posterior trees among genes in 2-to 3-dimensional Euclidean space. In this approach, the pairwise distances between posterior trees were computed (as Robinson-Foulds unweighted metric, Robinson & Foulds, 1981) from the package phangorn, and decomposed into a low-dimensional Euclidean space using metric multidimensional scaling (MDS). Before applying MDS, TREESPACE transformed the unweighted Robinson-Foulds tree distances into Euclidean distances using Cailliez’s transformation (Cailliez, 1983). The analysis was performed on posterior trees from 37 genes comprising an entire set of Chlorophycean taxa (68 spp.) sampled in the study, corresponding to the data used for Fig. 3. TREESPACE was performed on the collections of 3,700 posterior trees (representing 100 randomly sampled trees/gene) obtained from MrBayes MCMC analysis.

Chloroplast genome evolution in the Chlorophyceae

The enhanced taxon sampling in our data set was useful in showing several genomic features with peculiar evolutionary histories, despite the fact that most of our new genome data were not complete. In Fig. 4 (see Table S1 for detailed information), the psaC gene appears to have acquired a group II intron multiple times. This intron is inserted at position 25 in the chlorophyceans Borodinellopsis texensis, Jenufa perforata and Carteria sp. and is quite large, spanning over 6.6 kb in Borodinellopsis, 2.4 kb in Carteria, and 1.5 kb in Jenufa despite containing no detectable open reading frames (cutoff 500 bp) or protein domains. A group II intron in the same psaC position was also detected in Dictyochloris and Elakatothrix, but due to the partial nature of the data it is unclear how large the intron is, or if it is cis-spliced. PsaC also has a trans-spliced group II intron at position 25 in all studied OCC taxa –which is one of the genomic-structural synapomorphies for the clade. In all other included Chlorophyceae the psaC gene is intact, suggesting that these intron insertions have occurred independently several times. Such independent acquisition of an intron is consistent with the findings of others, e.g., McManus et al. (2012). Additional examples of dynamic intron movement are likely to emerge from chloroplast data, although it is in some cases difficult to determine whether the data represent multiple acquisitions or multiple losses of introns.

Perhaps most curious is the case of psaA, which is trans-spliced in the SV clade. The one salient exception is Carteria cerasiformis, which has a completely intact (or orthodox) psaA. An analysis of the psaA gene data did not reveal any indication of horizontal gene transfer unless from a close relative, and examination of the published genome did not reveal introns or other extra-genic DNA that may indicate joining of formerly separated portions. Thus, this apparent re-joining of exons remains unexplained, though perhaps not unique in Chlorophyceae (e.g., the mitochondrial rns gene in Neochloris in Fučíková et al., 2014). While most SV taxa have three trans-spliced psaA exons, the first two exons appear re-joined in Tetradesmus obliquus and Golenkinia longispicula. A retroprocessing mechanism behind genomic reconnection was recently proposed by Grewe, Zhu & Mower (2016), who discussed an analogous example in a mitochondrial gene of a vascular plant (Pelargonium, Geraniaceae). In general, trans-spliced arrangement is considered prohibitive of intron loss, but green algae might provide further exceptions to this rule.

The loss of infA is common to Volvocales and to Chaetophorales, providing a further example of convergent gene loss from the plastome. The most striking example, however, is the loss of the light-independent protochlorophyllide oxidoreductase (LIPOR) gene cassette (chlB, chlL and chlN) from Hafniomonas laevis and Neochloris aquatica, and also possibly from Golenkinia longispicula. The incomplete nature of the data for the latter taxon leave some room for doubt, but all other expected genes were recovered and reported from Golenkinia by Lemieux et al. (2015). The loss of the LIPOR genes from plastids of diverse algal groups, including red algae and other classes of green algae (Hunsperger, Randhawa & Cattolico, 2015) appears correlated with multiple copies of the nuclear-encoded por gene (light-dependent protochlorophyllide oxidoreductase). Nuclear genomes of the pertinent Chlorophyceae are unfortunately not available at this time, so we can only speculate whether the LIPOR genes were lost from Hafniomonas, Neochloris, and Golenkinia completely, and if por genes occur as single or multiple copies in these algae.

Increased character and taxon sampling as a phylogenetic remedy

The obvious advantage of whole plastome data is the great number of nucleotide characters available for phylogenetic inference. The obvious drawback is that there are still relatively few green algal taxa with whole cp genome sequences available. The case is reverse for the 18S marker: many taxa have an 18S accession in GenBank, but the entire gene is only ca. 1,800 base pairs long. In addition, not all of the gene’s nucleotides are phylogenetically informative, and some cannot be confidently aligned across distantly related taxa. Further, many 18S GenBank accessions contain partial data, only comprising a few hundred nucleotides. We have attempted a direct comparison of 18S and plastid data using the same taxon set (Fig. 2 and Fig. S3) and found the 18S phylogeny to be quite different from the plastid trees (e.g., Table 2), though often consistent with previous 18S-based studies. The character-rich plastid data set provided strong (often absolute) support for most relationships in our trees—but different analyses yielded different, and well-supported, topologies. Thus, simply adding more characters does not necessarily increase confidence in a phylogeny—much of the informative variation appears restricted to the shallow splits of the tree (Fig. 3), and the high support shown in concatenated trees may raise false confidence in the deep splits.

The problematic relative arrangement of Pediastrum, Neochloris and Chlorotetraedron (representatives of the families Hydrodictyaceae and Neochloridaceae) was explored in Fučíková, Lewis & Lewis (2016) without a satisfying explanation for the conflicting topologies in different analyses. Similarly, here we recovered all three possible groupings of the three taxa (Fig. 2 and Figs. S1–S5). Most recently this issue received extensive attention in McManus et al. (2018), who expanded the taxon sampling in Hydrodictyaceae to 14 species in five genera, but still recovered conflicting relationships between Hydrodictyaceae and Neochloridaceae in different single-gene and concatenated analyses. Expanding the sampling of Neochloridaceae, especially to include the genus Tetraedron, is the next logical step in the effort to resolve this systematic problem. However, given the varying signal from individual chloroplast and mitochondrial genes reported in McManus et al. (2018), it is possible that another confounding factor needs to be considered, such as horizontal gene transfer or insufficient model complexity in analyses. For the scope of our study, positioning of the two families is a minor issue, but the results of McManus et al. (2018) exemplify that increased taxon sampling may not always be a sufficient solution to all phylogenetic conflict; adding species helped within Hydrodictyaceae but did nothing to resolve the relationships among closely related families. Indeed, our main goal was to add resolution to the chlorophycean tree by increasing taxon sampling especially of deeply diverging incertae sedis lineages, but instead of converging on a single topology, we present several viable phylogenetic hypotheses. To our knowledge, we have sampled the known incertae sedis thoroughly, but perhaps species that would ‘break up’ additional long branches in the tree are yet to be discovered, or they may be extinct.

Unexpected and problematic lineage placements

The filament-formers Parallela transversalis and Microspora sp. grouped strongly together (Fig. 2, Supplementary Data). Our analyses either placed Microsporaceae along with their apparent coccoid relative Dictyochloris in the proximity of Scenedesminia (all aa cp analyses and the 18S analysis), or within the incertae sedis clade (all nt cp analyses). Previously, Parallela was placed in the proximity of Volvocales by Novis et al. (2010) based on data from three plastid genes, but because our study offers a better sampling of deeply diverging chlorophyceans combined with genome-scale sequence data, our results are likely more realistic, even if ambiguous. The family Microsporaceae was listed as member of Sphaeropleales by Tsarenko (2005), but as discussed in Fučíková, Lewis & Lewis (2014a), verifying this placement will be difficult due to the lack of type material for Microspora. Indeed the strain UTEX LB472 used in our study is not a type, and does not even have a species-level identification, so it should be viewed with caution as representative of Microsporaceae.

The respective affiliations of Spermatozopsis and Dictyochloris are fairly robust based on our analyses, but seem curiously at odds with the current understanding of these taxa based on morphology and ultrastructure. Spermatozopsis is traditionally viewed as member of Volvocales because of the clockwise orientation (CW) of its flagellar apparatus (Melkonian & Preisig, 1984). Among our analyses, only the BAli-Phy and PhyloBayes nt analyses placed Spermatozopsis at the base of Volvocales, separate from Sphaeropleaceae (Figs. S3 and S5), consistently with 18S analyses of Lewis (1997), and Němcová et al. (2011). Curiously, Marin’s (2012) nu rDNA analysis, using a closely related species Spermatozopsis exsultans, recovered this genus sister to Sphaeropleaceae, like our cp analyses did, but his cp rDNA analysis and the combined nu and cp rDNA analysis recovered the genus at the base of Volvocales with representatives of Carteria (also see simplified trees in Fig. 1) The presence of the gene infA in the cp genome of Spermatozopsis contrasts to its absence in all other volvocalean taxa examined, including the deeply diverging Hafniomonas and Carteria. Dictyochloris, in turn, invites placement in Sphaeropleales based on its flagellar ultrastructure nearly identical to that of Bracteacoccus, as reported by Watanabe & Floyd (1992). However, the unusual parallel configuration of the basal bodies in both taxa, similar to the volvocalean Heterochlamydomonas, led to further investigations by Shoup & Lewis (2003) who, again using 18S, placed Dictyochloris squarely in Sphaeropleales but unrelated to Bracteacoccus, and interestingly also showed Treubarinia as sister to Sphaeropleaceae. The phylogenetic tree of Němcová et al. (2011), however, places Dictyochloris sister to Microspora, consistently with our Fig. 2.

Another “straggler” taxon subtended by a long, weakly placed phylogenetic branch has previously been the coccoid genus Mychonastes (e.g., Krienitz et al., 2011). Whereas Fučíková, Lewis & Lewis (2014a) showed the genus at the base of Sphaeropleales, in none of our scenarios was Mychonastes the deepest diverging branch of the order. The lineage appears as a long-branched member of Scenedesminia without a particularly close alliance to other families.

Sphaeropleales—greatly broadened or greatly narrowed?

Though the above described phylogenetic results are intriguing, our primary focus is the inconsistently recovered monophyly of Sphaeropleales, which has far-reaching taxonomic implications (Fig. 1). As shown in Table 2, the clade designated here as Scenedesminia is phylogenetically quite robust. In fact, aside from Treubarinia it is our only focal group that appears unaffected by analysis type, and is supported by 18S data as well. Sphaeropleaceae, while morphologically quite different from Scenedesminia, share similarities in flagellar apparatus ultrastructure with them (Hoffman, 1984; Watanabe, Floyd & Wilcox, 1988; Wilcox & Floyd, 1988). Yet, molecular support for a monophyletic Sphaeropleales incl. Scenedesminia has been weak in previous studies, especially where 18S was the sole marker used (e.g., Shoup & Lewis, 2003). The support improved with the use of whole plastome protein coding data (Fučíková, Lewis & Lewis, 2016), even though not all genes agreed on the order’s monophyly. Despite building on the study of Lemieux et al. (2015) by including additional incertae sedis Chlorophyceae, our study does not bring the desired order to Sphaeropleales. Instead, Sphaeropleaceae, Treubarinia, and various other deeply diverging chlorophycean taxa vary in their placement across analyses (Table 2). Most of our cp nt phylogenies agree on Sphaeropleales s. l., which would include Microsporaceae, Dictyochloridaceae and Sphaeropleaceae as previously noted (Tsarenko, 2005) along with Scenedesminia. The hypothetical broadened order further includes Treubarinia (previously similarly placed by Shoup & Lewis (2003) and corroborated by Lemieux et al., 2015), Jenufa and Golenkinia, and surprisingly also Spermatozopsis similis. The cp nt trees disagree on the branching order within Sphaeropleales s. l., but the Bayesian amino acid phylogeny (Fig. 2B), both PhyloBayes trees (Figs. S4, S5), and the BALi-Phy 18S tree (Fig. S3) contradict this arrangement entirely, rendering Scenedesminia order-less and instead grouping Sphaeropleaceae at or near the base of Volvocales.

The PhyloBayes results hint at the possibility of site-specific evolutionary patterns in the nucleotide data, which the CAT model specifically tries to address (Lartillot & Philippe, 2004). SVDquartets, in turn, aims to account for incomplete lineage sorting and the resulting conflict among signals from different genes (Chifman & Kubatko, 2014; Chifman & Kubatko, 2015). Even though plastid genes are often assumed to evolve as a single locus, among-gene conflict has been demonstrated in previous studies (e.g., Fučíková, Lewis & Lewis, 2016; McManus et al., 2018), suggesting the possibility of more complex inheritance patterns. Lateral gene transfer is a possible explanation for among-gene conflict, but our data are at present insufficient to support or reject this possibility. Our SVDquartets tree (Fig. S2) supports the most “traditional” view of distinct Volvocales and Sphaeropleales, with Sphaeropleaceae and Scenedesminia grouping together, though also including Spermatozopsis as sister to Scenedesminia. However, neither of the two orders is well supported—both receive less than 50% BS—and therefore this result cannot be viewed as definitive. The SVDquartets analysis also strongly groups all incertae sedis taxa: Treubarinia, Microsporaceae, Dictyochloris, Golenkinia and Jenufa receive 97% BS as a clade. This grouping appears in several other analyses but mostly embedded within Sphaeropleales s.l. (Table 2, Fig. 2A, Figs. S1–S6).

It is possible that other factors contribute to the observed conflicts in phylogenetic signal. Saturation of 3rd positions is often blamed for phylogenetic issues, for example. However, in our case, this problem is unlikely to apply: saturation implies that the rate of evolution is so high that no historical signal is present. Saturated sites represent noise and thus are expected to reduce support overall. They should not strongly support any particular clade, nor should they conflict strongly with others sites that have a strong signal. We show this in our supplementary analyses: 3rd positions do not have much signal at the critical deep nodes (manifested by low BPP and polytomies), so they are unlikely to cause conflict. On the other hand, 3rd positions are very informative at shallow divergences, and thus likely contribute to good resolution in many parts of the concatenated tree.

One concern is that saturated 3rd position sites are not as affected by stabilizing selection and thus might have a different base composition than 1st or 2nd position sites. This is effectively handled by partitioning (as implemented in our analyses), which allows 3rd position sites to have their own nucleotide frequencies. Moreover, assessing saturation is problematic in itself, as “saturation plots” comparing model-corrected vs. uncorrected distances rely on pairwise comparisons. These comparisons are misleading because the largest distances span the entire tree using data from only two taxa. Bayesian and maximum-likelihood methods for estimating phylogeny instead estimate much shorter individual edges of the tree using all the data. Genes or data subsets (e.g., 3rd codon positions) that appear to be saturated using pairwise saturation plots may be evolving at a rate ideal for maximum-likelihood or Bayesian approaches. Evaluation of saturation is an area much in need of better methods. In summary, given the disparity of results among different analyses, the class Chlorophyceae clearly has a history of complex evolutionary processes that currently cannot all be accommodated in a single analysis.

It is entirely possible that with the addition of even more data (both in terms of taxa and genes), Sphaeropleales will be shown as a robust order. At present, however, it is at least equally possible that Sphaeropleales contains only one family, the Sphaeropleaceae. This would mean that Scenedesminia, and perhaps other lineages, be considered distinct orders as well. In our view, given the phylogenetic uncertainty, such drastic taxonomic changes do not seem justified. The problems with solidly anchoring the order-defining family Sphaeropleaceae in the chlorophyte phylogeny could lead to years of taxonomic changes and redefinitions in the future. Therefore, we adopt a practical non-Linnaean solution with well-supported clades as taxonomic currency instead (Nakada, Misawa & Nozaki, 2008). We consider Scenedesminia and Treubarinia useful names for the two well-supported groups. In addition, within Scenedesminia the clade of mostly colony-forming families (Hydrodictyaceae, Neochloridaceae, Scenedesmaceae and Selenastraceae, plus the coccoid Rotundellaceae) was well supported in nearly all our analyses. In the future, this ‘colonial clade’ might warrant a formal recognition, but here we merely note it as an interesting hint at a pattern in morphological evolution. In general, rather than make formal taxonomic changes in the Linnaean framework, we prefer the approach to name clades, as has been the practice in phycology for a long time, e.g., in Trebouxiophyceae, where order- and family-level relationships have long been unclear and a clade system has been in use (e.g., Neustupa et al., 2013; Fučíková, Lewis & Lewis, 2014b).