The efficacy of single mitochondrial genes at reconciling the complete mitogenome phylogeny—a case study on dwarf chameleons

Mitogenome generation

Our sampling protocol followed the ARROW guidelines (Field et al., 2019). Specifically, we acquired institutional animal care approval from South African National Biodiversity Institute: 003/2011, 001/2013, 001/2014, 0001/2015, University of Johannesburg: 2019-10-10, and University of the Witwatersrand: 2019/10/56/B for DNA sampling. We complied with relevant national, international and institutional regulations regarding animal care during DNA sampling, ensuring animals were in our care for minimal duration and under minimal stress. We took all measures possible to follow the 3R tenets, and taxon-specific guidelines for the genetic sampling of dwarf chameleons (Herrel et al., 2012). Total genomic DNA was extracted from tissue samples (tail tips) from dwarf chameleon individuals (n = 44) using the Qiagen DNeasy blood and tissue kit (Qiagen, Hilden, Germany). Samples were stored in either 99% ethanol or RNAlater and were collected during multiple field trips in South Africa between 2010 and 2022 under permits from the provinces of Eastern Cape, KwaZulu-Natal, Limpopo, Mpumalanga, Northern Cape, Western Cape (018-CPM403-00001, 0092-CPM401-00006, AAA004-00107-0035, CPM-333-00002, CRO 3/19CR, CRO 35/15CR, CRO 36/15CR, CRO 32/20CR, CRO 33/20CR, FAUNA 110/2011, OP 4758/2010, OP 4596/2010, OP 2635, OP 1259/2014, KZN 1647/2009, MPB.5299, MPB.5371, MPB.5544, MPB.5604, CN44-59-5795, RSH 24/2021, WRO 41/03WR; WRO 15/03WR). Subsamples were exported to the USA under CITES regulations (permits #142667, 206199, 257511, 260568 and 260570 issued by the South African Department of Fisheries, Forestry and Environment), with remaining tissue samples accessioned at the South African National Wildlife Biobank, Pretoria. Sampling included 14 of the 20 described species and three additional populations of taxonomic rank not yet established (Table 1). Genomic libraries were prepared using 10 ng of DNA per sample and an Illumina DNA Library Prep Kit which made use of Transposase-mediated Tagmentation. Samples were uniquely barcoded using IDT Nextera for Illumina DNA Unique Dual 10 bp Indexes adaptors. After library preparation, Illumina re-sequencing was carried out on an Illumina NovaSeq S4. Raw sequencing reads were trimmed, and adaptors were removed using trimmomatic v0.39 (Bolger, Lohse & Usadel, 2014) according to the following settings: LEADING:20 TRAILING:20 SLIDINGWINDOW:13:20 MINLEN:23 with the adaptor set to NexteraPE-PE.fa:2:30:10:4 and a minimum phred score of 33. Trimmed reads were then passed through NOVOPlasty v4.3.1 (Dierckxsens, Mardulyn & Smits, 2017) using a publicly available Chamaeleo chamaeleon complete mitogenome (accession number: NC012427; Macey et al., 2008) as starting seed and reference genome. Assembled mitogenomes were then annotated using the MITOS v1 web server (Bernt et al., 2013).

Phylogenetic analyses

Once annotated, mitogenomes were imported into Geneious Prime 2022 (https://www.geneious.com) and published complete mitogenomes for seven other chameleon genera were downloaded from GenBank as outgroup taxa (Table 2). Using these 51 individuals (44 ingroup, seven outgroup), we created a total of 22 different alignments for subsequent phylogenetic analyses: a full mitogenome alignment (excluding the highly variable tRNAs which are difficult to align), 15 single locus alignments (13 protein-coding genes, two non-coding ribosomal genes), the short fragment of 16S (Alrefaei, 2022; Bates et al., 2013; Hay et al., 1995; Hertwig, de Sá & Haas, 2004; Lamb & Bauer, 2002; Main, Van Vuuren & Tolley, 2018; Main et al., 2022; Vences et al., 2004; Vences et al., 2005; Welton et al., 2010), a commonly used short fragment of COI (Hebert et al., 2003), a concatenation of short 16S fragment with ND2 (Tolley, Chase & Forest, 2008; Tolley, Tilbury & Burger, 2022; as these two genes are regularly concatenated in the herpetological phylogenetic literature), a concatenation of ND2 and ND5 (shown by our analyses to be among the most reliable genes for phylogenetic reconstruction–see Results), a concatenation of both ribosomal subunits, and a concatenation of all protein-coding genes (PCGs). All new sequences from this study were deposited in GenBank (Table S1).

Maximum likelihood phylogenetic analyses were run for the complete mitogenome in IQ-TREE v2.0.3 (Minh et al., 2020b) with the optimal model of molecular evolution estimated for each gene separately using the built-in ModelFinder implemented in IQ-TREE (Chernomor, von Haeseler & Minh, 2016; Kalyaanamoorthy et al., 2017). To assess whether node support would be improved by merging models of evolution for similarly evolving genes into larger partitions, a merged partitioned analysis was also run. This analysis makes use of the PartitionFinder (Lanfear et al., 2017) algorithm to generate a partitioning scheme within IQ-TREE and merge genes for model estimation to assess whether this reduces over-parameterization.

For the Bayesian analyses, we estimated the best evolutionary models and partitions for the complete mitogenome in PartitionFinder v2.1.1 (Lanfear et al., 2017). The output suggested that the dataset be partitioned into 10 partitions with GTR + I + G the optimal model of evolution for each partition (Table 3).

Bayesian phylogenetic analyses were carried out on all datasets separately and the complete mitogenome in MrBayes v3.2.7 (Ronquist et al., 2012) with partitions and models included as per the PartitionFinder results (Table 3). For each of the 22 datasets (13 separate coding genes, two ribosomal genes, 16S fragment, COI fragment concatenated 16S and ND2 alignment, concatenated ND2 and ND5 alignment, concatenated protein-coding genes, complete mitogenome), two separate runs were made, each with four independent MCMC chains of 20 million generations, sampling every 1,000 generations. For every run we assessed convergence using the R package RWTY v1.0.2 (Warren et al., 2017) with sampling considered adequate when all parameters have ESS values >200. We also made use of ASTRAL version 5.7.8 (Zhang et al., 2020) to reconcile the mitogenome tree from each of the previously estimated single mito gene trees. ASTRAL makes use of the multispecies coalescent and quartet scores to estimate the best-supported species tree assuming a model of gene tree discordance. This assumption is not necessarily relevant when estimating the mitogenome tree because the mitogenome is a single non-recombining locus and therefore each mitochondrial gene is expected to share a single coalescent history. Nonetheless, we performed this analysis as an additional measure of gene tree congruence.

Finally, gene and site concordance factors were calculated by computing individual gene trees and comparing them to the best-supported merged-partitioned mitogenome tree in IQ-TREE v2.0.3 (Minh et al., 2020b). Concordance factors are typically used to measure topological discordance among independent loci but here all loci share a single history due to a lack of recombination (Minh et al., 2020a). Therefore, we used concordance factors to measure topological differences in resolution rather than independent history.

Topological congruency and phylogenetic resolution

To assess topological congruency between different gene trees and the complete mitogenome tree, Bayesian consensus trees for each of the 22 datasets were imported into R using the package Ape v5.6-2 (Paradis, Claude & Strimmer, 2004). Using the function tree.dist.matrix in the package RWTY, Robinson-Fould (RF) topological distances between trees inferred from different gene regions were calculated pairwise to produce a distance matrix of topological distances. We used the nj function in Ape to generate a neighbour-joining dendrogram from this topological distance matrix. To visualize if topological reliability relates to the length of the gene, we constructed a broken x-axis scatterplot in R, using the package ggbreak v0.1.1 (Xu et al., 2021) and the gap.plot function, of RF distance from the complete mitogenome against gene length (in 100 bp) for each gene as well as our different concatenated partitions. We also carried out Shimodaira-Hasegawa (SH) topology tests in IQ-TREE as an additional measure of topological congruency among gene trees (Shimodaira & Hasegawa, 1999).

Saturation

To ascertain whether topological discordance across different genes might be due to saturation, we estimated saturation for each gene separately and the complete mitogenome by plotting uncorrected genetic p-distances against model corrected p-distances for each gene (Everson, 2015; Philippe et al., 1994) using the package Ape and the dist and dist.correct functions. The difference between these estimates is an approximation of the unobserved mutations hidden by multiple mutations occurring at a site along a branch in our tree. Due to an abundance of gaps in the ribosomal RNAs, which would be ignored when calculating model-corrected p-distances and thereby inflate the slope of the linear regression, we did not estimate saturation plots for these genes. In these saturation plots, the slope of the linear regression is indicative of the degree of saturation, the lower the slope of the linear regression, the more saturated the gene is assumed to be. For this reason, we tabulated the slope of the saturation linear regression for each gene along with other metrics relevant to topology: number of polytomies, mean node posterior probability (calculated based only on the supported nodes within the phylogeny), number of fully bifurcating nodes, RF distance from mitogenome, number of parsimony informative sites, gene length, and standardized information content (SIC, see Macey et al., 2004; Table 4).

Explaining RF distance

In addition to comparing the topology of each gene tree to that of the complete mitogenome, we summed the number of polytomies, mean posterior probability support across the tree, and number of fully bifurcating nodes for each gene tree and partition in our dataset and tabulated these metrics together with gene/alignment length, number of parsimony informative sites, saturation and RF-distance. We sought to assess which of these metrics best explains RF distance. We therefore ran an analysis of variance (ANOVA) on the data tabulated above using the aov function in the stats package in R, omitting number of polytomies, number of fully bifurcating nodes, and mean posterior probability because these are all intrinsically linked to the topology of the tree and therefore to RF distance.

Comparison with other published phylogenies

Our complete mitogenome topology and our ND2-16S are largely congruent with published dwarf chameleon phylogenies constructed using concatenated 16S and ND2 gene fragments (Tolley et al., 2004; Tolley, Chase & Forest, 2008; Tolley, Tilbury & Burger, 2022). However, we find that our phylogeny based on ND2 alone is more robust and has a topology more similar to our complete mitogenome phylogeny. Despite the similarities between the complete mitogenome phylogeny and the previously published phylogenies, a few differences are apparent. Our full mitogenome topology places B. setaroi and B. caffrum as sister taxa although the relationship lacks node support (PP = 0.71; BS = 47) and therefore cannot be considered a compelling difference. More notably is that the full mitogenome phylogeny supports B. ngomeense as being nested within the B. transvaalense clade (PP = 0.99; BS = 92) instead of sister to that species. This could be due to our better geographic coverage for B. transvaalense than in previous studies (Tolley et al., 2004; Tolley, Chase & Forest, 2008) allowing for a more detailed assessment of these relationships. However, none of our other (single gene or concatenated) analyses showed B. ngomeense to be nested within B. transvaalense, all of which are similar to previously published topologies, and also show node support for the monophyly of the two species. This poses the question as to whether the full mitochondrial dataset has allowed for better resolution within this particular lineage, or whether certain mito-genes have experienced directional selection, providing differing topological outcomes even within the mitochondrion.

Phylogeny

To assess if the mitogenome phylogeny can be reconciled from single mitochondrial genes, we estimated Bayesian and maximum likelihood phylogenies for the entire mitogenome as well for each single gene and various combinations of genes (all single gene trees have been archived on the Harvard Dataverse: 10.7910/DVN/RY1RQI). We found some variation in the mean node support and number of fully bifurcating nodes between the whole mitogenome tree and the single gene trees (Table 4). For example, the whole mitogenome inferred tree yielded 50 fully bifurcating nodes with a mean posterior probability support of 0.98 while a short fragment of 16S only yielded 13 fully bifurcating nodes with a mean posterior probability support of 0.86 (Table 4).

The complete mitogenome consensus phylogenies (Bayesian and maximum likelihood) recovered Bradypodion as a well-supported monophyletic clade with most internal nodes well-supported both in terms of bootstrap and posterior probability support (Fig. 1). The coalescent phylogeny estimated using ASTRAL (Fig. S1) was very similar to the Bayesian and maximum likelihood phylogenies except that the coalescent phylogeny does not recover B. caffrum and B. setaroi as monophyletic sister taxa, but this is not particularly surprising given that the sister-taxon relationship of B. caffrum and B. setaroi received very low support in the Bayesian and maximum likelihood phylogenies (PP = 0.71; BS = 47). Gene concordance factors reveal that while most nodes in the complete mitogenome consensus phylogeny are well-supported (PP ≥ 0.95; BS ≥ 75), some of these nodes do not show concordance across all gene trees (Fig. 1). For instance, the node separating the B. caffrum-B. setaroi clade from the B. caeruleogula-B. transvaalense clade is well supported by both Bayesian and maximum likelihood methods, but it has a gene concordance factor of only 33.3, meaning that only 33.3% of the gene trees in the mitogenome support that node (Fig. 1). The placement of B. caeruleogula, while well-supported in both bootstrap and posterior probability support, is only supported by 53.3% of gene trees. In both of these instances, the site concordance factors are also low (34.09% and 38.58% respectively (Fig. S2). These low concordance factors might reflect a lack of phylogenetically informative information in some of the genes, hindering phylogenetic estimation for these genes. Typically, when site and gene concordance factors agree (which they tend to in our analyses), this suggests genuine topological conflict between genes. However, given that the mitochondrion is a single non-recombining locus, this conflict is more likely an artefact of certain genes lacking the phylogenetic information to converge on the true topology.

Analysis of topological congruence and phylogenetic resolution

We found varying degrees of topological inconsistency across phylogenies inferred by different genes within the mitogenome (Fig. 2). The dendrogram of RF topological distances shows that ND2, ND5, COIII, the concatenated PCG alignment, the concatenated ND2-16S and the concatenated ND2-ND5 all cluster with the tree inferred from the complete mitogenome, suggesting that the phylogenies inferred from these single genes or concatenations most closely resemble that of the complete mitogenome. However, the concatenated ND2-16S alignment is more distant from the complete mitogenome than the ND2 gene on its own, despite the concatenated alignment consisting of more characters. Sister to this cluster is a grouping of ND4 and COI, and just outside this cluster is ND6. More distant from the complete mitogenome are the ribosomal RNAs and ND3, Cyt-b, COII, ND4l, ATP6, ATP8, the shortened 16S fragment, and ND1.

It is apparent that RF distance from the complete mitogenome tree shows an inverse correspondence with gene length; however, there is considerable variation which indicates that some genes hold more information than others (Fig. 3). For instance, our longest alignment was for the combined protein-coding genes which also had the lowest RF distance from the mitogenome tree, while the shortest alignment was the ATP8 gene which had one of the highest RF distances from the mitogenome. ND2 had one of the lowest RF distances from the mitogenome, even lower than the concatenated ND2-16S alignment. Furthermore, the concatenated rRNA alignment was longer than all individual genes and yet still had a relatively high RF distance from the mitogenome, so phylogenetic informativeness does not correlate highly with gene fragment length. Nonetheless, concatenating two of the most reliable single genes, ND2 and ND5 produced a topology even closer to the complete mitogenome than either of these genes on their own.

While the SH-tests also showed ND5 to be the best single gene proxy for the complete mitogenome, there were a few discrepancies between the RF-distances from the complete mitogenome and the results of the SH-tests (Table 4). The SH-tests supported COIII as the second-best proxy for the mitogenome, with no other single genes being congruent to the complete mitogenome. Furthermore, the SH-tests supported the concatenated alignment of ND2 and ND5 as a better proxy for the complete mitogenome than the PCG alignment.

As might be expected, the complete Bradypodion mitogenome phylogeny scores best with the fewest polytomies, the most fully bifurcating nodes, and the highest mean posterior probability support and an SH-test p-value of 1. Based on these same criteria, the short 16S fragment scored lowest with the most polytomies, lowest mean posterior probability support, and the fewest fully bifurcating nodes. It should be noted that mean posterior probability node support is calculated based only on the supported nodes within the phylogeny. If there is an unsupported polytomy, this is left out of the calculation, so these values should be consulted in conjunction with the number of polytomies and number of fully bifurcating nodes. The protein coding gene alignment, ND2, ND5, and the concatenated ND2 and ND5 alignment performed best across most of these metrics.

It is apparent from our ANOVA that saturation is not driving topological discordance within the mitogenome (p = 0.59; F = 0.30). This is probably due to strong structural constraints on the evolution of protein-coding genes within the mitochondrion, constricting the degree to which third position bases mutate. We did, however, find RF distance from the mitogenome to correlate significantly with gene/alignment length (p = 0.006; F = 11.36), and we found a significant interaction between the number of parsimony-informative sites and total gene alignment length (p = 0.02; F = 7.36).