Large Variation in the Ratio of Mitochondrial to Nuclear Mutation Rate across Animals: Implications for Genetic Diversity and the Use of Mitochondrial DNA as a Molecular Marker

Abstract

It is commonly assumed that mitochondrial DNA (mtDNA) evolves at a faster rate than nuclear DNA (nuDNA) in animals. This has contributed to the popularity of mtDNA as a molecular marker in evolutionary studies. Analyzing 121 multilocus data sets and four phylogenomic data sets encompassing 4,676 species of animals, we demonstrate that the ratio of mitochondrial over nuclear mutation rate is highly variable among animal taxa. In nonvertebrates, such as insects and arachnids, the ratio of mtDNA over nuDNA mutation rate varies between 2 and 6, whereas it is above 20, on average, in vertebrates such as scaled reptiles and birds. Interestingly, this variation is sufficient to explain the previous report of a similar level of mitochondrial polymorphism, on average, between vertebrates and nonvertebrates, which was originally interpreted as reflecting the effect of pervasive positive selection. Our analysis rather indicates that the among-phyla homogeneity in within-species mtDNA diversity is due to a negative correlation between mtDNA per-generation mutation rate and effective population size, irrespective of the action of natural selection. Finally, we explore the variation in the absolute per-year mutation rate of both mtDNA and nuDNA using a reduced data set for which fossil calibration is available, and discuss the potential determinants of mutation rate variation across genomes and taxa. This study has important implications regarding DNA-based identification methods in predicting that mtDNA barcoding should be less reliable in nonvertebrates than in vertebrates.

mitochondria, mutation rate, substitution rate, polymorphism, natural selection

Introduction

In animals, mitochondrial DNA (mtDNA) is assumed to experience a higher mutation rate than nuclear DNA (nuDNA) (Ballard and Whitlock 2004). Surprisingly, despite the popularity of mtDNA as a marker in evolutionary studies, this assumption only relies on a handful of comparisons involving mostly vertebrates species (Brown etal. 1979; Miyata etal. 1982; Nabholz etal. 2009). The situation is much less clear in nonvertebrates. Depending on species, a mtDNA mutation rate much higher (Thomas and Wilson 1991; Blouin etal. 1998; Oliveira etal. 2008; Willett 2012), similar (Powell etal. 1986; Vawter and Brown 1986; Brower and DeSalle 1998), or lower (Shearer etal. 2002; Hellberg 2006; Huang etal. 2008) than the nuDNA mutation rate has been reported. Taking a comparative approach, Lynch etal. (2006) reported contrasted values for the ratio of mtDNA to nuDNA mutation rate (μ_mit/μ_nuc) among animals. In this analysis, μ_mit/μ_nuc varied from 0.57 and 1.43 in corals and mosquitoes, respectively, to around 18 in crickets and aphids and up to 25 in vertebrates. More recently, Havird and Sloan (2016) also reported a large range of μ_mit/μ_nuc in animals and plants with flies (Diptera) showing a clearly lower μ_mit/μ_nuc than mammals. These analyses, however, were based on a limited number of data sets (less than ten nonvertebrate taxa), so that the extent of the variation in μ_mit/μ_nuc and its distribution among taxonomic groups is currently unclear.

Understanding this variation has several important implications. mtDNA is by far the most frequently used marker for DNA-based species identification (the so-called “DNA barcoding” (Hebert etal. 2003; Hebert and Gregory 2005) and biodiversity monitoring using environmental DNA (Thomsen etal. 2012). The success rate of species identification using barcoding techniques is known to vary among taxonomic groups (Meier etal. 2006; Elias etal. 2007); whether this pertains to the variation in mtDNA mutation rate is worth investigating. Full mitochondrial genome data, typically obtained from next-generation sequencing, are also widely used in phylogenetic analyses (Botero-Castro etal. 2013; Barker 2014; Fabre etal. 2014; Tilak etal. 2014). Characterizing the taxonomic variation in mtDNA vs. nuDNA evolutionary rate should help assessing the performance of the mtDNA marker, compared to nuclear genes, for phylogenetic purposes.

Importantly, the μ_mit/μ_nuc ratio is also relevant to the interpretation of patterns of within-species polymorphism across genomes and species (Ellegren and Galtier 2016). In a meta-analysis of 1,683 species of animals, Bazin etal. (2006) found that the level of mtDNA polymorphism was surprisingly constant among taxonomic groups and uncorrelated to nuDNA polymorphism. This was interpreted as reflecting the effect of positive selection on mtDNA evolution. Under neutrality the genetic polymorphism is governed by the mutation/drift balance, but the fixation of a newly arising, positively selected allele (that is, a selective sweep) leads to complete loss of genetic diversity at linked loci—a process known as genetic hitchhiking Maynard-Smith and Haigh (1974). If pervasive, genetic hitchhiking is expected to homogenize the amount of polymorphism among species and erase the relationship between polymorphism and N_e (Gillespie 2001). This hypothesis appeared plausible because the animal mitochondrial genome, which is nonrecombining and gene dense, is expected to be more strongly affected by linked selection than nuclear loci (Hurst and Jiggins 2005; Bazin etal. 2006; Berlin etal. 2007). This interpretation, however, was debated because a correlation between mtDNA and nuDNA polymorphism was reported within specific groups such as mammals (Mulligan etal. 2006; Nabholz etal. 2008b; Piganeau and Eyre-Walker 2009).

Much of this discussion has focused on the effects of N_e and selection, in absence of large scale data on the relative mtDNA vs. nuDNA mutation rate. This is, we believe, in large part due to the fact that to explain the lack of relationship between mtDNA and nuDNA diversity, the μ_mit/μ_nuc ratio would have to exactly compensate for the difference in average N_e between taxa, a hypothesis deemed implausible by Bazin etal. (2006). However, more recently, Piganeau and Eyre-Walker (2009) have suggested the existence of a negative relationship between mitochondrial mutation rate per generation and N_e in mammals. Theoretically, the mutation rate is expected to evolve towards a lower limit determined by the power of random genetic drift (1/N_e), and therefore could be negatively correlated to N_e (Lynch 2010). There is no obvious reason, however, why this relationship should apply specifically to the mitochondrial genome, and not the nuclear genome.

In this study, we take advantage of the availability of phylogenetic data sets in a large number of taxa to revisit the mutation rate hypothesis. Analyzing 121 multilocus and four phylogenomic data sets, we estimate the μ_mit/μ_nuc ratio in a variety of animal phyla and show that the variation in μ_mit/μ_nuc is indeed sufficient to explain the absence of correlation between mtDNA and nuDNA polymorphism. Taxa with a high average nuclear genetic diversity (such as insects or mollusks) experience a mtDNA mutation rate that is close to the nuDNA mutation rate, leading to similar level of polymorphism in the two genomic compartments. In contrast, taxa with a low average nuclear genetic diversity (such as birds, mammals, or scaled reptiles) have a mtDNA mutation rate markedly higher that the nuDNA mutation rate, leading to a mtDNA polymorphism much more elevated that the nuDNA polymorphism. Controlling for the variation in μ_mit/μ_nuc, we were able to recover a positive correlation between nuDNA and mtDNA polymorphism. We conclude that there is no need to invoke positive selection to explain the results of Bazin etal. (2006), and discuss the causes and consequences of μ_mit/μ_nuc variation across animal phyla.

Results

Literature Review for Multilocus Data sets

We performed a systematic review of volumes 62–97 of the journal Molecular Phylogenetics and Evolution (January 2012–April 2016) and retrieved 122 data sets comprising at least one mitochondrial and one nuclear alignment covering a common set of at least five closely related species (supplementary table S1, Supplementary material online). We sorted these data sets according to 12 taxonomic groups, each of them represented by a variable number of data sets: Squamata N = 20; Teleostei N = 20; Insecta N = 19; Amphibia N = 13; Aves N = 11; Mollusca N = 10; Mammalia N = 8; Arachnida N = 8; Crustacea N = 6; Testudines N = 4; Chondrichthyes N = 3; Porifera N = 1. These groups are similar to those used by Bazin etal. (2006) except that we split “Pisces” in Teleostei (bony fishes) and Chondrichthyes (cartilaginous fishes) and Sauropsida in Aves (birds), Squamata (scaled reptiles) and Testudines (turtles). We also retrieved data in Arachnida, but did not find any suitable Echinoderm data set. Our single sponge (Porifera) data set yielded a ratio of mitochondrial to nuclear mutation rate lower than 1 (μ_mit/μ_nuc = 0.52), which is consistent with Huang etal. (2008). Porifera were excluded from further analyses due to insufficient sample size (N = 1). The final data sets comprised an average 38.6 species, 1.4 mitochondrial genes and 2.1 nuclear genes.

Extensive Variation in the Ratio of Mitochondrial to Nuclear Mutation Rate

The expected molecular divergence at neutrally evolving sites is equal to the product of mutation rate by divergence time (Kimura 1968; Birky and Walsh 1988). We selected the same set of species for mtDNA and nuDNA, such that the ratio of molecular divergences estimates the ratio of mutation rates, assuming that the times of divergences are similar for mtDNA and nuDNA (i.e., neglecting the variation in ancestral coalescence time among loci). To limit the impact of natural selection, we only analyzed third codon positions and estimated branch lengths using classical phylogenetic maximum likelihood methods. Molecular divergence was defined as the sum of branch lengths, and the μ_mit/μ_nuc ratio was estimated by dividing the mitochondrial molecular divergence by the nuclear one. The average number of substitutions per phylogeny was 3,377 and 443 for mitochondrial and nuclear loci, respectively.

The estimated ratio of mutation rates (μ_mit/μ_nuc) was highly variable among taxonomic groups (fig. 1, ANOVA testing the effect of taxonomical groups on μ_mit/μ_nuc variation, R² = 0.63, P < 2.2e−16). In Arachnida, the average mtDNA mutation rate (μ_mit) was only 3.1 times higher than the average nuDNA mutation rate (μ_nuc), whereas in Squamata (scaled reptiles) mitochondrial genes had a mutation rate 26.4 times higher than nuclear genes, on average. Among the 22 pairwise comparisons revealing a significant difference, 21 involved a comparison between a vertebrate and a nonvertebrate taxon (Tukey Test, adjusted p < 0.01). No statistical difference could be detected among nonvertebrate groups, whereas only one comparison was significant among vertebrate groups, namely, the Teleostei–Squamata comparison (Tukey Test, P = 0.003, Teleostei have a ratio of 11.6 on average). Between vertebrates and nonvertebrates, 21 comparisons out of 28 showed statistically significant differences (Tukey Test, adjusted P < 0.01).

Fig. 1.

Across-taxa distribution of the ratios of mitochondrial to nuclear mutation rate in animals. Boxes give the quartiles; whiskers extend to 1.5 times the interquartile range. Dots indicate ratio of individual data sets. Colors vary according to the value of the median of the taxon.

A number of surprising outliers were uncovered in several groups. The data set with the highest μ_mit/μ_nuc ratio in nonvertebrates was a mollusk (μ_mit/μ_nuc = 39.9), whereas all the other mollusks showed a ratio close to 4 (fig. 1). This extreme data set is composed of 95 species of land snails (Pulmonata: Bradybaenidae, Hirano etal. 2014) and we have no strong reason to question the reliability of this estimate. Similarly, in insects, two data sets yielded a μ_mit/μ_nuc ratio above 20, with the highest ratio obtained from a data set composed of 127 species of large Australian cicada (μ_mit/μ_nuc = 24.3, Hemiptera: Cicadidae; Owen etal. 2015). The data set with the largest ratio of all is composed of only ten species of Eurasian toads (Anura: Bufonidae; μ_mit/μ_nuc = 79.2; Recuero etal. 2012) and therefore might be affected by sampling error.

Extensive Variation of μ_Mit/μ_Nuc Is Supported by Genome-Wide Data sets

Each of the data sets analyzed above included just a handful of genes. To assess the genome-wide variation in μ_mit/μ_nuc, we used four phylogenomic data sets—two in vertebrates (primates and birds, 5,566 and 3,146 genes, respectively) and two in nonvertebrates (Drosophila and Bivalvia, 8,359 and 398 genes, respectively). In this analysis, we computed a μ_mit/μ_nuc per nuclear gene by dividing the mitochondrial divergence (averaged across all mitochondrial genes) by the divergence of each individual nuclear gene. Some among-nuclear gene variation in μ_mit/μ_nuc was detected in all four taxonomic groups, but the distribution of μ_mit/μ_nuc was clearly influenced in the first place by taxonomy (fig. 2, ANOVA, R² = 0.96, P < 2.2e−16). The mean μ_mit/μ_nuc per gene varied from 0.8 in Drosophila and 1.8 in bivalves to 7.4 in birds and 32.5 in primates. This genome-wide analysis therefore confirms the existence of a conspicuous difference in μ_mit/μ_nuc between vertebrate and nonvertebrate taxa, and demonstrates that the effect reported in figure 1 is not specific to the genes typically used in phylogenetic analyses.

Fig. 2.

Across-genes distribution of the ratio of mitochondrial to nuclear mutation rate in four phylogenomic data sets. Boxes give the quartiles; whiskers extend to 1.5 times the interquartile range. Dots indicate ratio of individual genes above or below the whiskers limit. Colors vary according to the value of the median of the taxon. Y-axis is log-transformed.

Correcting for Mutation Rates Recovers the Correlation between mtDNA and nuDNA Coalescence Time

Our data clearly demonstrate that mtDNA and nuDNA mutation rates do not vary in a similar way across metazoan taxa. This result calls for a reappraisal of the relationship between mtDNA and nuDNA polymorphism (Bazin etal. 2006). To this aim, we computed the corrected mtDNA polymorphism, π_mit*, defined as the mtDNA polymorphism π_mit divided by the ratio of mutation rates μ_mit/μ_nuc. Assuming selective neutrality and a balanced sex ratio, the corrected mtDNA polymorphism is expected to equal 1/4 of the nuDNA polymorphism (π_nuc). We obtained the π_mit and π_nuc values used by Bazin etal. (2006) and use the μ_mit/μ_nuc ratio estimated above.

As reported by Bazin etal. (2006), the average π_mit varied little among taxa and was not correlated to the average π_nuc (Spearman's rank correlation ρ = 0.52, P = 0.20, fig. 3). Once corrected for μ_mit/μ_nuc, the correlation coefficient increased and became statistically significant (Spearman's rank correlation ρ = 0.90, P = 0.004, fig. 3), revealing a higher average mtDNA coalescence time in species with an elevated π_nuc. Our results therefore suggest that the absence of correlation between π_mit and π_nuc can be explained by the variation in μ_mit/μ_nuc among animal taxa.

Fig. 3.

Relationship between mitochondrial and nuclear average synonymous diversity across matazoan taxa. Y-axis: uncorrected (red) and mutation rate-corrected (blue) average level of synonymous polymorphism in mtDNA X-axis: average level of synonymous polymorphism in nuclear DNA. Y-axis and X-axis are log-transformed. Data from Bazin etal. (2006). Am: Amphibia, Ar: Arachnida, C: Crustacea, I: insecta, M: Mammalia, Mol: Molluska, P: Pisces, S: Sauropsida.

We also reproduced the main figure of Bazin etal. (2006, their fig. 1) by plotting the average π_mit* and π_nuc against the average allozyme heterozygosity (fig. 4). π_mit* was not significantly correlated to the allozyme heterozygosity (Spearman's rank correlation ρ = 0.69, P = 0.07) but the correlation coefficient increased markedly compared to the uncorrected π_mit (π_mit vs. allozymic heterozygosity, Spearman's rank correlation ρ = 0.02, P = 0.98) and became close to the correlation between allozymic and π_nuc.

Fig. 4.

Relationship between mitochondrial synonymous diversity and allozyme heterozygosity across matazoan taxa. Y-axis: uncorrected (blue/circles) and mutation rate-corrected (red/trianlges) average level of synonymous polymorphism in mtDNA and nuclear DNA (purple/squares) X-axis: average level of allozyme heterozygosity Y-axis and X-axis are log-transformed. Abbreviations: see legend to figure 3.

In our data sets, π_mit* is on average 1.5 times lower than π_nuc. This is substantially higher than the neutral expectation (π_nuc = 4 π_mit*). This could be explained by (i) an underestimation of the μ_mit/μ_nuc ratio, (ii) and unbalanced sex ratio with an excess of females.

Absolute Mutation Rates

The large variation in μ_mit/μ_nuc among metazoan taxa could result from a variable μ_mit, a variable μ_nuc or variation in both μ_mit and μ_nuc. To investigate these possibilities, we estimated the per year mutation rate for 28 data sets in which at least one fossil calibration point was available (Sauropsida N = 8; Insecta N = 7; Mammalia N = 5; Amphibia N = 2; Testudines N = 2, Aves N = 1; Chondrichthyes N = 1; Mollusca N = 1; Teleostei N = 1).

Including all data sets, μ_mit was not correlated to μ_nuc (Spearman's rank correlation ρ = 0.31, P = 0.11; fig. 5), but a significant correlation was detected within vertebrates, that is, excluding insects and mollusks (Spearman's rank correlation ρ = 0.71, P < 0.001; fig. 5). The absolute μ_mit was higher than μ_nuc in all taxa but the difference was much more pronounced in vertebrates (fig. 6A). The per year nuclear mutation rate was significantly lower in vertebrates than in insects and mollusks (median μ_nuc = 0.0012 per site per 10⁶ years in vertebrates, 0.0044 in nonvertebrates, Wilcoxon' test P = 0.002), whereas no significant difference was detected between vertebrates and nonvertebrates as far as the mitochondrial mutation rate was concerned (median μ_mit = 0.0199 per site per 10⁶ years in vertebrates, 0.009 in nonvertebrates, Wilcoxon' test P = 0.07).

Fig. 5.

Relationship between the per year nuclear and mitochondrial mutation rates. Vertebrate data sets are in blue (squares). Nonvertebrate are in orange (Insects, circles) and green (Mollusk N = 1, triangle). The blue line is the regression line computed using only vertebrate data.

Fig. 6.

Distribution of the per year and per generation mutation rates across taxa and genomes. (A) Per year mitochondrial (circles with crosses) and nuclear (diamonds) mutation rate in vertebrates and nonvertebrates. (B) Per generation mitochondrial (circles with crosses) and nuclear (diamonds) mutation rate in vertebrates and nonvertebrates. Vertical grey lines joint mitochondrial and nuclear estimates obtained from the same data set. Y-axis is log-transformed.

Polymorphism levels are determined by mutation rate per generation, not per year. We obtained rough estimates of the average per generation mutation rates using age at sexual maturity or number of generation per year as a proxy for generation time. This considerably altered the distribution of mutation rate across vertebrates and nonvertebrates (fig. 6B). The average estimated per generation μ_mit was higher in vertebrates than in nonvertebrates (median μ_mit = 26.75 × 10⁻⁹ per site per generation in vertebrates, 3.90 × 10⁻⁹ in nonvertebrates, Wilcoxon' test P < 0.001) whereas the average per generation μ_nuc was not significantly different between nonvertebrates and vertebrates (median μ_nuc= 1.60 × 10⁻⁹ per site per generation in vertebrates, 1.31 × 10⁻⁹ in nonvertebrates, Wilcoxon' test P = 0.58, fig. 6B). Therefore, the shorter generation time of nonvertebrate species seems to result in a lower mitochondrial mutation rate per generation. The lower per generation mutation rate of nonvertebrate species apparently compensates for their higher N_e leading to a similar π_mit between vertebrates and nonvertebrates and explains the homogeneity of mitochondrial polymorphism levels across animals.

Estimation of Neutral Divergence and the Effect of GC Content

In order to estimate mutation rate variation from divergence data it is important to limit the impact of natural selection. The divergence of neutral sites is strictly proportional to mutation rate (Kimura 1968). This is true even for neutral sites genetically linked to selected sites (e.g., third codon position) (Birky and Walsh 1988). In this respect, a common practice is to use synonymous sites as proxy of neutral sites (but see Chamary etal. 2006). However, the estimation of the number of synonymous sites could be an issue. Specifically, the method of Goldman and Yang (1994) (GY94) appears to be influenced by GC content and non-stationary base composition (Bierne and Eyre-Walker 2003; Guéguen and Duret 2017). In our data set, we detected a strong relationship between the proportion of synonymous sites (i.e., number of synonymous sites/total number of sites) and GC content (supplementary fig. S5, Supplementary Material online). In the GY94 model, synonymous and nonsynonymous sites are counted as mutational opportunities, that is, the relative number of synonymous and nonsynonymous sites per codon depends on the estimated rates of substitution to its neighbor codons. We use the F3x4 model, where the equilibrium codon frequencies are empirically estimated from observed nucleotides frequencies. Equilibrium codon frequencies influence substitution rates between codons under GY94 and, consequently, the relative number of synonymous and nonsynonymous positions. In GC-rich alignments, the substitution rate between G- and C-ending codons, which are often nonsynonymous, are assumed to be higher than in alignments with balanced GC content—and similarly for substitutions between A- and T-ending codons in AT-rich alignments. This leads to low estimated numbers of synonymous sites in alignments with extreme GC content. This is illustrated by the fact that the proportion of synonymous sites greatly increases if we assume a common equilibrium frequency for all codons (f = 1/61) instead of F3 × 4. For example, the most AT-rich mitochondrial alignment (COI in Melitta, Hymenoptera; GC = 7%) has an estimated 9% of synonymous sites with the F3x4 parameter but this proportion increases to 30% when equal equilibrium codon frequencies are assumed. Similarly, in the most GC-rich nuclear alignment (shaw in Aleyrodidae, Hemiptera; GC = 63%), the proportion of synonymous sites increase from 9% to 38%.

It is unclear whether the F3 × 4 parameter reflects a true biological feature, that is, whether in AT-rich alignments, transversion between AT-rich codons are truly more frequent that in alignments with a more balanced GC content. Nevertheless, we replicated all our analyses using either the synonymous divergence (estimated under GY94) or 4-fold degenerated sites, instead of all third-codon positions in our main analysis. The 4-fold degenerated sites have the advantage to be both synonymous and straightforward to count. However, the number of 4-fold degenerated sites can be limiting especially when the divergence is high—we only included sites that were 4-fold degenerated in all the analyzed species. With both type of sites, the main results remained unchanged. All the results obtained using the synonymous divergences and 4-fold degenerated sites are presented as supplementary material text S1 and supplementary figures S1–S4, Supplementary Material online.

Finally, the distribution of GC content greatly varies between genomes and taxa. Particularly, the most AT-rich loci are nonvertebrate mitochondrial loci and the most GC-rich loci are nuclear loci. It is possible that GC content is biologically linked to substitution rate via GC dependent mutation rate and/or GC-biased gene conversion (Duret and Galtier 2009). However, the average mitochondrial GC content and nuclear GC content have not effect on the variation of μ_mit/μ_nuc as tested in a multiple-factor type II ANOVA accounting for the effect of taxonomy (P > 0.10 for GC content effect and P = 7.7e−08 for taxonomy).

Discussion

High μ_Mit/μ_Nuc in Vertebrates

Our analysis of 121 data sets encompassing 4,676 species of animals revealed extensive variation in μ_mit/μ_nuc among phyla of animals. Although somewhat continuous, this variation was mostly partitioned between vertebrates (with μ_mit/μ_nuc typically above 10) and nonvertebrates (with μ_mit/μ_nuc typically below 5). Our results confirm early reports, based on very small samples, of a similar order of magnitude for μ_mit and μ_nuc in nonvertebrate taxa (Powell etal. 1986; Vawter and Brown 1986; Brower and DeSalle 1998). We also detected a few non-vertebrate taxa with a high μ_mit/μ_nuc, similar to that of vertebrates, but these are exceptions.

Analyzing 28 data sets for which fossil calibrations were available, we estimated the per year mutation rate separately for mtDNA and nuDNA. We found that the per year μ_nuc is higher in nonvertebrates than in vertebrates, perhaps reflecting an effect of generation time. The per year μ_mit, in contrast, did not differ significantly between vertebrates and nonvertebrates. This is surprising, knowing the extensive variation in per year μ_mit that was previously reported, and here confirmed, within mammals (Nabholz etal. 2008a; Welch etal. 2008) and within birds (Nabholz etal. 2009, 2016; Nguyen and Ho 2016). In particular, the strong effect of life-history traits, such as body mass, longevity, and generation time, on the per year μ_mit that was documented in mammals and birds (Nabholz etal. 2008a; Welch etal. 2008; Nabholz etal. 2009, 2016) does not seem to hold at the Metazoa scale. Using a proxy of generation time, we observed that the average per generation μ_mit, which matters as far as polymorphism is concerned, is much higher in vertebrates that in nonvertebrates.

Robustness to Gene Sampling

Many of the μ_mit/μ_nuc ratios were estimated based on data sets comprising just one locus per genome (supplementary table S1,Supplementary Material online), which might not be representative of the genome-wide ratio. To test whether the variation in μ_mit/μ_nuc is robust to gene sampling, we gathered four genome-wide data sets of vertebrates and non-vertebrates species. In all cases, we identified substantial variation in μ_mit/μ_nuc among genes. This could be due to a variable nuclear mutation rate across the genome (Ellegren etal. 2003; Hodgkinson and Eyre-Walker 2011), but part of the variation could also come from stochastic effects (although we excluded alignments shorter than 200 bp) and alignment errors. Despite this large variation, a strong taxonomic effect was detected, confirming that the variation in μ_mit/μ_nuc among taxonomic groups is not an artifact of gene sampling. Many of our Arachnida and Mollusca data sets all comprise a single nuclear locus, namely Histone H3. This might have influenced our result if third codon positions were particularly fast-evolving in this gene. However, in the bivalves phylogenomic data set, the μ_mit/μ_nuc ratio for histone H3 is 3.6 which is higher than the mean μ_mit/μ_nuc (1.8). This demonstrates that histone H3 is rather a slow evolving gene (leading an high μ_mit/μ_nuc) compared to other nuclear protein-coding genes and that our results, if anything, underestimate the true effect.

Reconsidering the Genetic Draft Hypothesis

Our analysis demonstrates that the variation in μ_mit/μ_nuc among animal clades explains the absence of correlation between π_mit and π_nuc reported by Bazin etal. (2006). The μ_mit/μ_nuc-corrected π_mit* was well correlated to π_nuc, indicating that the variation in effective population size across animals affects the nuclear and mitochondrial genomes in a comparable way, as previously reported within some vertebrate clades (Mulligan etal. 2006; Popadin etal. 2007; Nabholz etal. 2008a; Piganeau and Eyre-Walker 2009). Bazin etal. (2006) invoked a more pronounced effect of genetic draft (Gillespie 2001) on mtDNA in large-N_e taxa to interpret the lack of correlation between π_mit and π_nuc. Our analysis provides an alternative explanation to the results of Bazin etal. (2006). It shows that the constancy of mean π_mit among animals groups could be explained by a negative correlation between μ_mit/μ_nuc and effective population size.

Our results imply that the action of natural selection is not pervasive enough to erase the effect of genetic drift on mitochondrial diversity. They do not, however, exclude that positive selection sometimes occurs in mtDNA evolution (da Fonseca etal. 2008). Recently, James etal. (2016) have estimated that ∼26% (5.7–45%) of amino-acid substitutions affecting mitochondrial proteins are adaptive in animals, and some of the reported instances of mitochondrial/nuclear phylogenetic incongruence have been interpreted as reflecting adaptive mitochondrial introgression (Toews and Brelsford 2012). In our analysis, we report a significant correlation between π_mit* and π_nuc but the correlation between π_mit* and allozyme heterozygosity was not statistically significant. It should be noted that here the values of μ_mit/μ_nuc (this study) and polymorphism (Bazin etal. 2006) were estimated from distinct data sets, including different loci and species. As a consequence, the correlation coefficients we report are likely to be underestimated.

Why Is μ_Mit/μ_Nuc Correlated to N_e?

There is no obvious theoretical reason to expect a relationship between μ_mit/μ_nuc and effective population size. This partly explains why Bazin etal. (2006) did not consider mutation rate as a plausible confounding factor, in absence of large-scale data on μ_mit/μ_nuc. To gain insight into the causes of μ_mit/μ_nuc variation, we estimated absolute mutation rates per year. Our results indicate that both μ_mit and μ_nuc are variable among species and higher level taxa. The variation between vertebrates and nonvertebrates seems more pronounced for the nuclear genome than for the mitochondrial genome. Finally, μ_mit appears correlated to μ_nuc in vertebrates, suggesting that the per year mitochondrial and nuclear mutation rates are influenced by common determinants. This is not surprising as several life-histories traits such as body-mass and sexual maturity are known to influence both mitochondrial and nuclear mutation rates (Bromham 2009). The absence of correlation between the two rates in insects could simply result from insufficient power due to the limited number of data sets available (N = 7). We were able to estimate absolute mutation rates in only 28 data sets in total, which limits our ability to test alternative hypotheses that could explain the variation in μ_mit/μ_nuc. Below, we discuss four nonexclusive hypotheses that could potentially be tested with an extended data set.

First, the metabolic rate could have an impact on the mtDNA mutation rate (Martin etal. 1992; Martin and Palumbi 1993; Rand 1994). Because of its location, the mitochondrial genome is expected to be impacted by mutagenic byproducts of aerobic respiration such as reactive oxygen species (ROS) (Turrens 2003). Insects have a lower mass-specific metabolic rate than endothermic vertebrates (Makarieva etal. 2008) and therefore, insects could have a reduced production rate of ROS that could limit their mtDNA mutation rate. However, many ectothermic vertebrates have a mass-specific metabolic rate similar or lower than that of many insects (Makarieva etal. 2008). Moreover, scaled reptiles have a similar μ_mit/μ_nuc than birds and mammals despite their lower mass-specific metabolic rate (White etal. 2006). Finally, the relationship between mass-specific metabolic rate and ROS is currently unclear and may not be linear (Galtier etal. 2009).

Second, the number of genome replications per generation could differ between the mitochondrial and nuclear genomes, and among taxa of animals. A high μ_mit/μ_nuc could be explained by a higher ratio of mtDNA/nuDNA replication cycles per generation in vertebrates than in nonvertebrates. In mammals, it has been shown that oocyte maturation involves a dramatic increase in the number of mitochondrial genomes (Mishra and Chan 2014). A large increase in mtDNA amounts, however, was also reported in the in oocytes of some nematodes (Tsang and Lemire 2002). Whether the number of mtDNA replications per generation generally differs between vertebrates and nonvertebrates is therefore currently unclear.

Third, the mitochondrial genome is replicated by a specific DNA polymerase, namely DNA polymerase γ. A single event involving a reduction in the fidelity of DNA polymerase γ in the ancestor of vertebrates would be sufficient to explain most of the μ_mit/μ_nuc variation we report.

Fourth, Lynch (2010) provided some theoretical justification and empirical evidence for a negative correlation between per generation mutation rate and effective population size. Natural selection is expected to push down the mutation rates to a value at which further decrease would imply a negligible fitness increase, compared to the strength of genetic drift (1/N_e) (Lynch 2010). The relationship between N_e and the fitness effect of the germline mutation rate might be different for the mitochondrial and the nuclear genomes therefore leading to a variable μ_mit/μ_nuc scaling with N_e. Direct estimation of mutation rates through sequencing of genomes in mutation accumulation or pedigree experiments are still rare (Lynch 2010; Smeds etal. 2016) and a formal test of this hypothesis may not be possible in the short term. However, indirect evidence for or against the above hypotheses could be obtained by confronting phylogenetic estimation of μ_mit/μ_nuc with species' life-history traits taken as proxies of metabolic rate (Makarieva etal. 2008), generation time (Bromham etal. 1996), and effective population size (Romiguier etal. 2014).

Implications for DNA Barcoding

Finally, our results have implications for DNA-based identification methods. An ideal barcoding marker should form robust monophyletic clades reflecting species boundaries. This ability depends on the mutation rate of the considered marker, which determines the expected number of mutations diagnostic of a given species, and the strength of genetic drift (1/N_e), which determines the age of the most recent common ancestor of individuals in a species. If species boundaries were predictable from nuDNA divergence (Roux etal. 2016), our results would imply that mtDNA barcoding is more powerful in vertebrates than in non-vertebrates species (fig. 7). For a given level of nuDNA divergence, one should sequence 3–10 times more mtDNA sites in a nonvertebrate than in a vertebrate to obtain the same level of phylogenetic resolution. To some extent, this is also true of bony fishes and amphibians compared to birds, mammals and scaled reptiles. This is in line with reports of a relatively low success of the barcoding technique in some non-vertebrate taxa (Meier etal. 2006; Elias etal. 2007). Finally, direct comparison of species delimitation methods and biodiversity monitoring using the same mtDNA markers across distinct phyla of animals may be biased knowing the large variation in μ_mit/μ_nuc we report.

Fig. 7.

Schematic representation of mitochondrial DNA genealogies in two fictive pairs of vertebrate and non-vertebrate species. T_split is the age of population split (without gene flow) between the two species in number of generations. The blue and red parts of the genealogies correspond to within-species polymorphism. 2N_e is the average time to the most recent common ancestor in a demographically stable population. Stars represent mutations. N_e is assumed to be ten times higher, and μ ten times lower, in nonvertebrates than in vertebrates. In vertebrates mtDNA, a small N_e and a high mutation rate facilitate the identification of monophyletic clades that match species boundaries. In nonvertebrates, a high N_e and a low mutation rate decrease the expected number of diagnostic mutations.

Materials and Methods

Multilocus Data sets

We extracted all data sets from « Molecular Phylogenetics and Evolution » by systematically scanning the table of content between volume 62 (January 2012) and 97 (April 2016). All selected data sets have at least one mitochondrial gene and one nuclear gene. Any species represented by only nuclear or mitochondrial sequence were excluded. We also selected only one sequence per species, following the taxonomy adopted in the publications. This step aims at excluding polymorphism data from the data sets. The data sets including less than five species were excluded. The data sets were split in 12 taxonomical groups: Amphibia; Arachnida; Aves; Chondrichthyes; Crustacea; Insecta; Mammalia; Molluska; Porifera; Squamata; Teleostei; Testudines. Because the literature review was very time consuming, we limited the research to 20 data sets per taxonomical group.

Sequences were extracted from NCBI/GenBank (http://www.ncbi.nlm.nih.gov/nucleotide/; last accessed April 2016) using the accession numbers available from the publications. Details of the data sets are presented in supplementary table S1, Supplementary Material online.

Molecular Divergence Estimation

Coding DNA Sequences (CDS) were extracted using custom made BioPython scripts (Cock etal. 2009). Sequences were aligned gene by gene using MACSE (Ranwez etal. 2011). Phylogenetic trees were inferred using RAxML with the substitution model GTRCAT (Stamatakis 2006). One phylogeny per locus was inferred. This will avoid potential discordance between gene trees that could lead to a spurious excess of inferred substitutions produced by incomplete lineage sorting (Mendes and Hahn 2016).

Next, we extracted third codon positions in order to limit the effect of natural selection. Sites containing more than 20% of missing data were excluded. Molecular divergence was estimated using the BASEML program (Yang 2007) and the F84 substitution model without rate variation among sites. Total third codon positions Divergence (TD) was estimated gene by gene by computing the sum of branch lengths.

Assuming that all mutations are neutral at third codon positions, we have: TD = μ × t; where μ is mutation rate and t the divergence time. For every data set, the total divergence time is the same for the nuclear and mitochondrial phylogenies because the same set of species was analyzed for both markers. As a result, we are able to estimate the ratio of mitochondrial and nuclear mutation rate μ_mit/μ_nuc by computing the ratio of TD. When several mtDNA and/or nuDNA markers were available, we took the mean ratio across all mitochondrial versus nuclear pairs of markers.

Besides third codon position analysis, we also calculated branch lengths in unit of per synonymous site synonymous substitutions using CODEML (Yang 2007) and using 4-fold degenerated sites. Four-fold degenerated sites were extracted from the alignment using custom C ++ program based on the Bio ++ libraries (Guéguen etal. 2013). These analyses were conducted the same way as explained above.

Alignments are available from Figshare repository (URL: https://doi.org/10.6084/m9.figshare.5117467.v1).

Polymorphism Data

Mitochondrial and nuclear synonymous polymorphism (π_mit and π_nuc, respectively) were obtained from Bazin etal. (2006) (Eric Bazin, pers. commun.). We defined the mutation-corrected polymorphism, π_mit*, as:

$${\pi }_{\text{mit}}^{\text{*}}{\text{= π}}_{\text{mit}}\hspace{0.17em}{\mu }_{\text{nuc}}{/\text{μ}}_{\text{mit}}$$

μ_nuc/μ_mit being defined as the median μ_nuc/μ_mit of each taxonomical group of Bazin etal. (2006). For ‘Pisces’ group, we used the median μ_nuc/μ_mit of the Teleostei and Chondrichthyes.

Assuming neutrality and a balanced sex-ratio, π_mit* is expected to be proportional to π_nuc.

$${\pi }_{\text{mit}}^{\text{*}}=N_{\mathrm{e}}\hspace{0.17em}{\mu }_{\text{mit}}\hspace{0.17em}{\mu }_{\text{nuc}}{/\text{μ}}_{\text{mit}}\hspace{0.17em}=\hspace{0.17em}{N}_{\mathrm{e}}\hspace{0.17em}{\mu }_{\text{nuc}}\hspace{0.17em}=\hspace{0.17em}{\pi }_{\text{nuc}}/4$$

Phylogenomic Data sets

Four phylogenomic data sets were analyzed: two vertebrates (Birds and Primates) and two nonvertebrates (Drosophila and Bivalvia). In birds, we used the alignment of Jarvis etal. (2014). We focused on six species of passerines + parrots: ((Melopsittacus, Nestor), Acanthisitta (Corvus, Geospiza, Taeniopygia)). In primates, alignments were obtained from the ORTHOMAM database (version 9, http://www.orthomam.univ-montp2.fr/; last accessed April 2016; Douzery etal. 2014). The data set included 12 species: (((((((Homo, Pan), Gorilla), Pongo), Nomascus), ((Macaca, Papio), Chlorocebus)), Callithrix), Tarsius, (Otolemur, Microcebus)). The Bivalvia data sets were obtained from González etal. (2015). We focused on the 11 species of Neoheterodontei: ((Cycladicama, Diplodonta), ((Mercenaria, (Glossus, Arctica)), (Corbicula, Cyrenoida))))))). Finally, the Drosophila data sets were obtained from Romiguier etal. (2014). The phylogenetic relationship between the 12 species is ((((((D. melanogaster, (D. Simulans, D. sechellia)), (D. erecta, D. yakuba)), D. ananassae), (D. pseudoobscura, D. persimilis)), D. willistoni), ((D. mojavensis, D. virilis), D. grimshawi)).

Nucleotide alignments were available for all the data sets except for the Bivalvia in which only amino-acid alignments were produced. For these data sets, we extracted the nucleotide sequences from assembled transcriptomes (Vanessa L. Gonzalez, pers. commun.). Sequences were aligned gene by gene using MACSE (Ranwez etal. 2011).

Full mitochondrial genomes were available for birds, primates and Drosophila. For Bivalvia, we used only the Cytochrome Oxidase 1 gene. Mitochondrial sequences were aligned gene by gene using MACSE (Ranwez etal. 2011).

For all the data sets, alignments were cleaning using Gblocks option: −t = c −b4 = 5 −b5 = h (Talavera and Castresana 2007) and sites containing more than 20% of missing data were excluded. Alignments including less than 200 sites were excluded. Molecular divergence was estimated using BASEML program (Yang 2007) and the F84 substitution model without rate variation among sites. Total third codon positions Divergence (TD) was estimated gene by gene by computing the sum of branch lengths. μ_mit/μ_nuc was estimated by dividing the TD of a mitochondrial gene by the mean TD of the nuclear genes.

Absolute Mutation Rates and Generation Times

When they were available, fossil calibration points were extracted from the publications. We re-estimated a single phylogeny based on mitochondrial and nuclear sequences using RAxML with the substitution model GTRCAT partitioned by gene (Stamatakis 2006). Molecular dating was performed gene by gene using MCMCTREE (Yang 2007) with the HKY85 + Gamma model. The absolute mutation rate was computed as the average rate across genes.

Generation times were approximated using age at sexual maturity or number of generations per year. For vertebrates, we mostly relied on the AnAge database (version 13, Tacutu etal. 2012). For nonvertebrates, generation times were obtained from the literature and, in rare cases, from web-sites (see supplementary table S2, Supplementary Material online for details). In most cases, no information was available for the very species of our data sets and we used the median value of a taxonomical group including our species (mostly at the family level). For example, for a data set composed of Trioceros Chameleo species (Ceccarelli etal. 2014), we used the median value of all Chameleonidae in AnAge.

Statistical Analysis

All the statistical analyses were performed using R (R Core Team 2013). The effect of taxonomy on μ_mit/μ_nuc was tested with an Analysis of variance (ANOVA). The effect of the average GC content of mitochondrial (GCmito) and nuclear (GCnucl) GC content was tested with a multiple-factor type II ANOVA in a model including taxonomy (log (μ_mit/μ_nuc) ∼ taxonomy + GCmito + Gcnucl). Variables were log transformed. Correlations were assessed using Spearman’s non-parametric method.

Acknowledgments

This work is dedicated to the memory of Eric Bazin. We are grateful of Vanessa L. Gonzalez for provided the transcriptomes of Bivalvia species. We thank Sylvain Glémin and the anonymous reviewers for helpful discussion. This work was supported by Agence Nationale de la Recherche grants ANR-15-CE12-0010 “DarkSideOfRecombination” to B.N. and N.G, ANR-14-CE02-0002-01 “BirdIslandGenomic” to B. N and ANR-10-BINF-01-01 “Ancestrome” to N.G. This is ISEM publication number ISEM 2017-141.

Supplementary Material

Supplementary data are available at Molecular Biology and Evolution online.

Author notes

Associate editor: Nadia D. Singh