The mitochondrial phylogeny of land plants shows support for Setaphyta under composition-heterogeneous substitution models

Sampling of mitochondrial genomes

We sampled 21 taxa representing the major lineages of land plants, plus 5 green algae species as outgroup taxa. The taxa sampled for this study were: Coleochaetales (Chaetosphaeridium globosum, Coleochaete scutata), Zygnematales (Closterium baillyanum, Gonatozygon brebissonii, Roya anglica), liverworts (Aneura pinguis. Marchantia polymorpha, Pleurozia purpurea, Treubia lacunosa), mosses (Atrichum angustatum, Bartramia pomiformis, Physcomitrella patens, Sphagnum palustre, Tetraphis pellucida), hornworts (Megaceros aenigmaticus, Phaeoceros laevis), lycophytes (Huperzia squarrosa, Isoetes engelmannii), pteridophytes (Ophioglossum californicum, Psilotum nudum), and spermatophytes (Brassica napus, Cycas taitungensis, Ginkgo biloba, Liriodendron tulipifera, Oryza sativa, Welwitschia mirabilis).

Algal cultures for Coleochaete scutata (SAG 3.91) and Gonatozygon brebissonii (SAG 292) were obtained from the Culture Collection of Algae (SAG) (Georg-August-Universität Göttingen, Germany), and the algal culture of Roya anglica (ACOI 799) was obtained from the Coimbra Collection of Algae (ACOI) (Universidade do Coimbra, Coimbra, Portugal). The mitochondrial genomes of the three taxa were sequenced and assembled de-novo using standard methods as described in Civán et al. (2014), and annotated with the aid of Mitofy (Alverson et al., 2010) and NCBI BLAST analyses (Altschul et al., 1990). The remainder of the genomes were retrieved from NCBI GenBank (Bethesda, USA). The list of species samples, their classification, and the source and accession numbers of the sequences used are shown in Table 1.

Alignment and model testing

The sequences of each of 43 mitochondrial protein-coding genes were aligned using the program MAFFT (vers. 6.8; Katoh & Toh, 2008). Nucleotide alignments were manually edited in Geneious (vers. 9; http://www.geneious.com) to remove out-of-frame indels, misaligned portions, premature stop codons, positions with less than 50% taxon occupancy, and ambiguities. Genes that were missing within the algal outgroup taxa (Coleochaetales and Zygnematales) or were under 200 bp in length were discarded. The final data set consisted of 36 genes for 26 taxa, with 9.5% of gene sequences missing from the alignment. Amino acid alignments were generated by translation from each of the 36 nucleotide matrices using SeaView (vers. 4.5.4; Gouy, Guindon & Gascuel, 2009). The best-fitting substitution models for the 36 amino acid alignments were inferred in PartitionFinder (Lanfear et al., 2012) using the BIC selection criterion. The stmtREV (Liu et al., 2014) with a gamma-distribution of among-site rates (+G) model was the best-fitting model for all genes but one, for which the chosen model was JTT+G (Jones, Taylor & Thornton, 1992).

Gene tree estimation and monophyly test

Gene trees were estimated from individual nucleotide matrices, using the general time-reversible model (Rodriguez et al., 1990) with a gamma rate distribution and estimated base frequencies (GTR+G+Fest), and from amino acid matrices, using the best-fitting models inferred inPartitionFinder (stmtREV+G+F_est, JTT+G+F_est). Bayesian MCMC analyses were performed with P4 (Foster, 2004) using tree-homogeneous (henceforth referred to as CV1, i.e., one composition vector) and tree-heterogeneous (NDCH) composition models. Each analysis had two independent runs which were assessed for convergence by calculation of the marginal likelihood (chains were considered to have converged if they differed by <10 units) and of the average standard-deviation of split support between trees sampled from the posterior distributions (chains were considered to have converged if <0.01). Posterior predictive simulations of the X² test statistic of composition homogeneity was used to assess composition fit (Foster, 2004). For each gene, 50% majority-rule consensus trees were generated from the best-fitting analyses of the nucleotide and amino acid data.

Gene tree topologies inferred from each of the 36 amino acid alignments were tested using “gene genealogy interrogation” (GGI; Arcila et al., 2017) to ascertain whether the non-monophyly of the five major lineages under scrutiny (hornworts, liverworts, mosses, tracheophytes, and the outgroups) was statistically supported by the data. We were to consider any gene that supported the non-monophyly of one of these clades as aberrant and not suitable for inclusion in the combined analyses, as the monophyly of hornworts, liverworts, mosses, tracheophytes and embryophytes has been consistently recovered in molecular phylogenies (e.g., Qiu et al., 2006; Liu et al., 2014). Optimal trees for each gene were compared to each of fifteen constraint trees representing all possible topologically resolved combinations of the five monophyletic groups, using a nonparametric bootstrapping test. The results were assessed for statistical significance with the Approximately Unbiased (AU) test (Shimodaira, 2002) in CONSEL (vers. 0.1k; Shimodaira & Hasegawa, 2001). Optimal trees were estimated in RAxML (MPI-compiled vers. 8.2.8; Stamatakis, 2014) using the model stmtREV+G+F_est and starting from 20 random trees. Constraint trees were written in Newick format with internal branches within each of the five major clades collapsed to a polytomy. Each constraint tree was optimized in RAxML under the stmtREV+G+F_est model. Constraint topologies were considered statistically supported by the data if the p-value of the AU test was equal or greater than 0.05, meaning that the monophyly of each of the five lineages could not be rejected. In every gene, at least one of the constraint topologies was supported by the data, meaning that the monophyly of each clade could not be rejected, thus all 36 genes were included in downstream analyses.

Analyses of concatenated nucleotide data

A nucleotide alignment with 24,864 characters was obtained from the concatenation of the 36 individual genes. A second concatenated nucleotide matrix was constructed by codon-degenerate recoding of the data where ambiguity codes are used to negate synonymous substitutions (Criscuolo & Gribaldo, 2010; Cox et al., 2014). Bayesian MCMC analyses were performed on the concatenated and codon-degenerate nucleotide matrices using a tree-homogeneous composition model (CV1; 2 replicates) and the tree-heterogeneous composition (NDCH2; 4 replicates) model, as implemented on P4. In contrast to the original NDCH model (Foster, 2004) that allows an a priori defined number of compositions to evolve on the tree, the NDCH2 model estimates a separate composition for each node of the tree, constrained by a sampled Dirichlet prior on how much the composition vectors may differ from the empirical composition. Model fit to composition heterogeneity was inferred during the Bayesian MCMC with posterior predictive simulations of the X² statistic of composition homogeneity, where p-values equal or greater than 0.05 indicate acceptance of the model. The GTR+G+F_est model of substitution was used for for all MCMC of nucleotide and degenerate nucleotide data. Marginal likelihoods were estimated in P4 according to the eq16 method of Newton & Raftery (1994).

Analyses of concatenated amino acid data

An amino acid alignment with 8,288 characters was obtained by concatenation of the amino acid translations of the 36 genes. Bayesian MCMC analysis was performed on the concatenated amino acid data using both tree-homogeneous composition (stmtREV+G+F_{CV 1}; 2 replicates) and tree-heterogeneous composition (stmtREV+G+F_NDCH2; 4 replicates) models, with the fit of the composition evaluated by posterior predictive simulations as described for the nucleotide data. In addition, a Bayesian MCMC analysis was also performed using PhyloBayes (MPI-compiled vers. 1.6; Lartillot, Lepage & Blanquart, 2009) under the model stmtREV+G+F_CAT.

Alignments of individual genes, the concatenated data, and the resulting tree files of each analysis are available on Zenodo (doi: 10.5281/zenodo.3554149). All ML analyses of the concatenated nucleotide and amino acid data sets were consistent with the homogeneous Bayesian MCMC analyses and are not reported here, but the resulting tree files are also available on Zenodo.

Nucleotide and codon-degenerate data

All individual genes were best-fit by a tree-heterogeneous composition model with two composition vectors on the tree (CV2; Table S1). Majority-rule consensus trees resulting from the best-fitting Bayesian MCMC analyses of individual genes had low resolution in general, but liverworts were supported (>95% posterior probability (PP)) as the earliest-branching lineage in two genes (nad 3 and nad 5), whereas the mosses were supported as the earliest-branching lineage in one gene (ccm C). All other resolutions of the bryophyte lineages relative to the tracheophyte clade were not statistically supported, and the Setaphyta clade was not resolved in any gene tree.

Bayesian MCMC analysis of the concatenated nucleotide data set, assuming a homogeneous composition (CV1), resulted in a tree with mosses as the sister-group to the remaining land plants (PP = 1.0; Fig. 1A; Fig. S1), and hornworts as sister-group to the tracheophytes (PP = 1.0). In contrast, the analysis of the degenerate data set under a homogeneous base composition (Fig. 1C; Fig. S2) returned a tree with liverworts the earliest-branching lineage of embryophytes (PP = 1.0) and hornworts the sister-group to tracheophytes (PP = 0.93). However, the homogeneous model was rejected for both data sets by the posterior predictive simulation of the X² statistics of homogeneity, with a tail-area probability of 0.0, thereby indicating that the data were not composition homogeneous. The tree-heterogeneous composition analysis (NDCH2) of the concatenated nucleotide data also resulted in mosses supported as the earliest-diverging land plant lineage (PP = 0.98; Fig. 1B; Fig. S3), but placed liverworts as the sister-group to tracheophytes (PP = 0.94). The NDCH2 model was a good fit to the data according to the posterior predictive simulation (X² p = 0.99). When analysing the codon-degenerate data, the NDCH2 model recovered hornworts as the sister-group to the remaining land plants with full branch support (PP = 1.0; Fig. 1D; Fig. S4), and mosses fully supported (PP = 1.0) as the sister-group to liverworts (i.e., the clade Setaphyta). The NDCH2 model was not a good statistical fit to the data for the best scoring MCMC run, according to the X² posterior predictive test (X² p = 0.038; Fig. S4).

Amino acid data

Individual mitochondrial protein alignments were best-fit by both homogeneous and heterogeneous composition models, with some being best-fit by a model with up to four compositions (CV4), indicating that they are highly heterogeneous in composition among lineages (Table S1). Majority-rule consensus trees of best-fitting Bayesian MCMC analyses of individual proteins were poorly supported regarding relationships among bryophyte lineages. Indeed, only one tree (ccm C) showed any statistically supported resolution and placed the mosses as the earliest-branching lineage of embyrophytes. Although not statistically supported, one amino-acid tree showed bryophytes as a monophyletic group (rps 7), and the clade Setaphyta was present in three others (atp 4, rpl 2, and sdh 4).

When the concatenated amino acid data were analysed with a tree-homogeneous composition model (CV1), the resulting tree recovered liverworts as the sister-group to the remaining embryophytes without statistically significant support (PP = 0.82), and hornworts as the sister-group to tracheophytes, also without support (PP = 0.82; Fig. 1E; Fig. S5). The homogeneous composition model did not fit the data (X ² p = 0.0). The Phylobayes stmtREV+G+F_CAT analysis (stationary, data-heterogeneous composition model) of the concatenated amino acid data resulted in an unsupported placement of liverworts as the sister-group to the remaining embryophytes (PP = 0.89; Fig. S6). A posterior predictive composition homogeneity test using Phylobayes (readpb_mpi parameter -comp) showed that the data rejected the model and that the data were therefore composition tree-heterogeneous (p = 0.0). When the tree-heterogeneous composition NDCH2 model was used to model the concatenated protein data, the separate analyses did not converge on the same tree topology, although the NDCH2 model was a good fit (X² p = 0.1022). The tree obtained from the analyses with the best marginal likelihood (-lnLh 142829.1129) supported hornworts as the earliest-branching lineage (PP = 1.0), with the liverworts as the sister-group to mosses (PP = 1.0; Figs. 1F, 2, Fig. S7).

Gene tree discordance in mitochondrial data is likely due to mis-modeling and insufficient phylogenetic signal

Alternative hypothesis testing using nonparametric and parametric bootstrapping has been applied before to the mitochondrial land plant phylogeny to test the fit of the data (Liu et al., 2014). Here we chose a different approach, and used the optimized likelihood of constraint trees to identify genes that did not support the monophyly of the major embryophyte lineages (hornworts, liverworts, mosses, and tracheophytes) and of the outgroup. The optimal trees of the 36 mitochondrial genes, inferred under maximum-likelihood, show varied topologies, among which none is predominant. Because the four major land plant lineages are known to be monophyletic (as shown by many studies, e.g., Wickett et al., 2014) our concern was to identify gene trees that showed non-monophyly of one of these groups. The strategy we adopted allowed us to discern whether such topologies truly reflect underlying data or whether they are one among different topologies supported by the data. The AU test indicated, in all genes, that the monophyly of each land plant lineage was not statistically contradicted. This result suggests that the observed phylogenetic conflict among mitochondrial gene trees is unlikely to be explained by biological processes, such as horizontal gene transfer or duplication-loss, affecting specific lineages within each of the four major groups, and that any observed paraphyly of major groups on gene trees is probably the result of inadequate data modeling or paucity of phylogenetic signal.

Synonymous substitutions are responsible for the placement of mosses as the earliest-branching lineage of embryophytes

The tree inferred from the concatenated nucleotide data set of 36 mitochondrial genes shows mosses as the sister-group to the remaining land plants, as previous analyses of mitochondrial nucleotide data have shown (Liu et al., 2014). However, the mosses are replaced by the liverworts in the same position when analysing codon-degenerate recoded data. Codon degenerate recoding is used to eliminate synonymous substitutions, which are unconstrained by selection at the protein level and therefore can be subject to high rates of substitution, and ultimately saturation and loss of phylogenetic signal. Indeed, as the exclusion of synonymous substitutions is sufficient to eliminate phylogenetic signal that supports mosses as the sister-group to the remaining land plants, these results illustrate that despite being the slowest evolving genomic compartment in plants, phylogenetic inferences from highly divergent mitochondrial genomes are also affected by substitutional saturation due to the effect of high substitution rates. Moreover, this observation implies that caution should be taken when invoking biological explanations (e.g., hybridisation, incomplete lineage sorting) to explain incongruence among phylogenies when inadequate modeling of the substitution process may represent a simpler explanation.

Codon degenerate nucleotide data and protein data support Setaphyta under tree-heterogeneous models of substitution

With the nucleotide data, both the tree-homogeneous and NDCH2 tree-heterogeneous models support mosses as the earliest-diverging group in the embryophytes. The likely incorrect placement of mosses as the earliest-diverging group that is recovered with the best-fitting tree-heterogeneous NDCH2 model suggests that homoplasy driven by high nucleotide substitution rates (saturation) may overwhelm the ability of the model to correct for composition bias, hence the need to use codon-degenerate recoded data in combination with tree-heterogeneous models. Indeed, when codon-degenerate recoded data are used, contrasting supported relationships are obtained under tree-homogeneous and tree-heterogeneous composition models. Whereas using a homogeneous model for the analysis of the codon-degenerate data shows liverworts well supported as the sister-group to other embryophytes, the tree-heterogeneous analysis (NDCH2) model places liverworts as the sister-group to the mosses (clade Setaphyta), with maximum support, and hornworts as the sister-group to all other embryophytes, also with maximum branch support. These results demonstrate that the phylogenetic signal contained in non-synonymous sites is also subject to composition biases and that tree-heterogeneous composition models are required to model the data effectively. Contrasting results were also obtained when the amino acid data were analysed with tree-homogeneous and tree-heterogeneous models. The tree inferred under the poorly-fitting (P = 0.0) homogeneous model (CV1) resolves liverworts as the sister-group to the remaining land plants (PP = 0.82). By contrast, the Bayesian MCMC run with the highest marginal likelihood under the NDCH2 (fitting) model shows strong support for the Setaphyta clade (PP = 1.0 –mosses plus liverworts) with the hornworts as the earliest-branching lineage of embryophytes (PP = 1.0). Liu et al. (2014) observed topological congruence between mitochondrial nucleotide and protein data that placed liverworts as the sister-group to all other embryophytes, but this placement of liverworts received low branch support in different analyses, and thus no firm conclusions regarding these cladogenic events were put foward. Importantly, in that study no tree-heterogeneous analyses of the codon-degenerate data were performed, nor was the protein data analysed with more than two composition vectors on the tree. In contrast, our analyses of the codon-degenerate nucleotide data and amino acid data using a better-fitting tree-heterogeneous model resulted in well supported, congruent topologies, strengthening the argument in favour of the analyses presented here, that show support for the clade Setaphyta.

The land plant mitochondrial phylogeny is partially congruent with nuclear and chloroplast phylogenies

In contrast to previous analyses of the land plant mitochondrial phylogeny, we show that both nucleotide and amino acid data carry signal that joins mosses and liverworts as sister lineages (clade Setaphyta). This phylogenetic signal is typically obscured due to substitution saturation in the nucleotide data and among-lineage composition bias in both the nucleotide and amino acid data. In the nucleotide data, phylogenetic signal supporting mosses as the sister-group to the remaining land plants is eliminated when codon-degenerate recoded data is analysed, and instead the liverworts are found as the sister-group to all the remaining embryophytes under tree-homogeneous composition models. However, it is only when a combination of codon-degenerate recoding and a better-fitting tree-heterogeneous composition model is used that the mosses and liverworts appear resolved as sister taxa, therefore suggesting that both substitutional saturation and among-lineage composition heterogeneity are important evolutionary processes to be modeled in the nucleotide data. Similarly, the unsupported placement of liverworts as the earliest-branching lineage is obtained using tree-homogeneous composition models with the amino acid data, but better-fitting tree-heterogeneous composition models again support the mosses plus liverwort clade.

Support for the moss-liverwort sister-group relationship has been found in trees previously inferred from nuclear and chloroplast protein-coding data (e.g., Nishiyama et al., 2004; Cox et al., 2014; Puttick et al., 2018; Sousa et al., 2019). The clade can be resolved by mitochondrial data with our analyses, and therefore avoids the necessity of calling upon biological explanations, such as hybridisation, to account for incongruence among the phylogenies of the three plant genomes regarding the placement of mosses and liverworts. However, if the placement of the hornworts as the earliest-branching lineage of embryophytes does indeed reflect the true mitochondrial topology, then it is in conflict with the nuclear and chloroplast data which suggest the bryophytes are likely monophyletic. A biological process involving a rapid divergence of the hornworts from other bryophytes, after the tracheophyte-bryophyte split, and the retention of a copy of the mitochondrion that was lost in all other embryophyte lineages, could be invoked to explain the observed phylogenetic conflict. However, the incongruence of the mitochondrial data could, of course, still be a result of mis-modeling or lack of phylogenetic signal. It is likely that further sampling of mitochondrial genomes from hornworts and other bryophyte lineages may aid resolution of the phylogeny, but such analyses would only be informative if they were in combination with the heterogeneous composition models that have been shown here to be necessary to correctly model the underlying processes of mitochondrial evolution.