Stronger Evidence for Relaxed Selection Than Adaptive Evolution in High-elevation Animal mtDNA

Abstract

Mitochondrial (mt) genes are the subject of many adaptive hypotheses due to the key role of mitochondria in energy production and metabolism. One widespread adaptive hypothesis is that selection imposed by life at high elevation leads to the rapid fixation of beneficial alleles in mtDNA, reflected in the increased rates of mtDNA evolution documented in many high-elevation species. However, the assumption that fast mtDNA evolution is caused by positive selection, rather than relaxed purifying selection, has rarely been tested. Here, we calculated the d_N/d_S ratio, a metric of nonsynonymous substitution bias, and explicitly tested for relaxed selection in the mtDNA of over 700 species of terrestrial vertebrates, freshwater fishes, and arthropods, with information on elevation and latitudinal range limits, range sizes, and body sizes. We confirmed that mitochondrial genomes of high-elevation taxa have slightly higher d_N/d_S ratios compared to low-elevation relatives. High-elevation species tend to have smaller ranges, which predict higher d_N/d_S ratios and more relaxed selection across species and clades, while absolute elevation and latitude do not predict higher d_N/d_S. We also find a positive relationship between body mass and d_N/d_S, supporting a role for small effective population size leading to relaxed selection. We conclude that higher mt d_N/d_S among high-elevation species is more likely to reflect relaxed selection due to smaller ranges and reduced effective population size than adaptation to the environment. Our results highlight the importance of rigorously testing adaptive stories against non-adaptive alternative hypotheses, especially in mt genomes.

elevation, altitude, selection, adaptation, mitochondria

Introduction

Altitudinal gradients are natural laboratories for ecology because they encompass overlapping clines in many environmental factors. Organisms must cope with concordant global changes in temperature, oxygen availability, and UV radiation with elevation, as well as local clines in precipitation and primary productivity. Accordingly, adaptation to high altitude is apparent in the genetic, physiological, and ecological traits of both animals and plants (Projecto-Garcia et al. 2013; Read et al. 2014; Zhu et al. 2018; Storz and Scott 2019; Gutierrez-Pinto et al. 2021; Ye et al. 2020; Storz 2021). Adaptations to hypoxia and low temperature at increased elevation have been extensively studied, particularly in birds and mammals, because of the presumed selection pressure these factors exert on the high aerobic metabolism of active endotherms (Chappell et al. 2007; Cheviron and Brumfield 2012; Scott et al. 2015; Storz and Scott 2019; Gutierrez-Pinto et al. 2021; Nabi et al. 2021).

Several phenotypes strongly associated with elevational adaptation, including changes in hemoglobin activity, body insulation, or behavior, are likely to be largely or entirely linked to variation in the nuclear genome. However, the process of oxidative phosphorylation, the key component of aerobic metabolism, involves the coordinated activity of genes from both the nuclear and mitochondrial genomes (Rand et al. 2004; Sloan et al. 2018; Moran et al. 2024). The electron transport system which powers oxidative phosphorylation generates both the majority of cellular ATP and a significant portion of endotherm body heat via non-shivering thermogenesis. Thus, electron transport system proteins, and mitochondrial structures and processes more generally, have emerged as key foci in studies of adaptation to elevation (Luo et al. 2013; Murray and Horscroft 2016; Scott et al. 2018; Hill 2019).

While some studies have proposed mechanistic adaptive hypotheses for specific substitutions in mitochondrial genes at high elevations (see Scott et al. 2011; Ji et al. 2012; Kostin and Lavrenchenko 2018), studies more often report broad molecular signatures of selection and then assume their functional salience for high-altitude living without further tests. For instance, mitochondrial adaptation to elevation has repeatedly been inferred from unique or excess nonsynonymous substitutions in mammals (Xu et al. 2005; Luo et al. 2008; Hassanin et al. 2009; Ning et al. 2010; Yu et al. 2011), birds (Zhou et al. 2014), reptiles (Jin et al. 2018), amphibians (Wang et al. 2021), freshwater fishes (Li et al. 2013), and arthropods (Zhang et al. 2013; Li et al. 2018; Yuan et al. 2018). This has often been reported via elevated ratios of nonsynonymous to synonymous substitution rates—the d_N/d_S ratio. While positive selection can explain excessive nonsynonymous substitutions, so can relaxed selection, a distinction which is rarely acknowledged in even the most exhaustive studies of this genre (e.g. Wang et al. 2021; see Zwonitzer et al. 2023). The d_N/d_S ratio alone cannot distinguish positive and relaxed selection from one another and additional metrics are needed (Wertheim et al. 2015).

Studies have long highlighted the sensitivity of d_N/d_S to relaxed selection, particularly due to reduced effective population size (N_e) (Kliman et al. 2000; Woolfit and Bromham 2003). Relaxed selection provides an alternative explanation for substitution patterns at high-elevation because high-elevation species should often have lower N_e than their low-elevation relatives. Typically, less habitat is available at higher elevations (see Eakins and Sharman 2012), and, all else being equal, species with less habitat should have smaller ranges, lower population sizes, and experience stronger genetic drift and less efficient selection. In addition, high-elevation populations should be more fragmented by topography, which might further reduce N_e (e.g, Hahn et al. 2012; Polato et al. 2017). Among possible markers, animal mitochondrial genes may be especially affected by changes in N_e. In contrast to the nuclear genome, animal mtDNA lacks recombination, is uniparentally inherited, and is effectively haploid, making it more susceptible to genetic drift than nuclear loci according to classic theory (Lynch and Blanchard 1998; Neiman and Taylor 2009; but see Bazin et al. 2006; Edwards et al. 2021). Formal tests that can disentangle positive from relaxed selection are now available (Wertheim et al. 2015) and can be used to investigate molecular patterns of adaptation related to altitude (e.g. Gutiérrez et al. 2023; Fig. 1a).

Fig. 1.

a) Schematic representation of selection inference from both d_N/d_S ratios (calculated by codeML in PAML) and k-values (calculated in RELAX). b) Illustration of different models assessed with PAML.

To foster a mechanistic understanding of adaptive and relaxed mitochondrial evolution across elevation, patterns of molecular evolution should be evaluated with regard to taxonomic group, metabolic strategy (endotherm vs. ectotherm), and body size. The impact of temperature should also be parsed from that of oxygen limitation by examining latitude as a co-variate. To that end, we examined signatures of positive and relaxed selection in all mitochondrial-encoded (mtDNA) protein-coding genes across nearly 800 species that span a range of elevations. We combined our analyses with high-quality ecological data to characterize the genetic consequences of elevation, latitude, range size, and body size across terrestrial vertebrates, freshwater fishes, and arthropods. We took two approaches; first, we compared high-elevation species to their low-elevation sister taxa and outgroups in 153 phylogenetically-independent triads, estimating selection on high and low branches and comparing it to ecological data. We then examined 6 family or sub-order level clades totaling 380 species, evaluating selection on tips at progressively higher elevations and latitudes or smaller range sizes. Together, these analyses provide multiple streams of evidence for distinguishing between adaptive evolution and neutral processes in relation to elevation and other ecological factors.

Methods

Selection of mtDNA in Focal High-Elevation Species

Identification of High-Low-Outgroup Triads

We first tested for selection in the mtDNA of high-elevation focal species compared with their low-elevation relatives. We constructed phylogenetic trees of three species (“triads”) consisting of a high-elevation taxon, a low-elevation sister taxon, and a low-elevation outgroup (supplementary table S1, Supplementary Material online). To be included in our analyses, the high- and low-elevation taxa had to differ in maximum elevation by at least 500 m. The outgroup was the most closely related taxon for which a mitogenome could be recovered which also had a maximum elevation 500 m or more lower than the high-elevation focal species. High- and low-elevation taxa were typically members of the same genus but sometimes were closely related genera, subspecies, or discrete populations. The high-elevation taxon in each comparison represented what was inferred to be an independent origin of higher-elevation living compared to the other triads based on the distribution of elevation ranges among its close relatives and any available information on the history of the group. Each comparison was thus internally phylogenetically controlled and both high- and low-elevation taxa were independent of high- or low-elevation taxa in any other comparison. Care was taken to pair species that were matched for latitude, aridity, migratory behavior, and other potentially important ecological factors (supplementary methods S1.1, Supplementary Material online).

Analysis of Selection on High-elevation Species

Triads were organized into six taxonomic groups: mammals (n = 55 triads), birds (n = 34), non-avian reptiles (n = 13), amphibians (n = 23), bony fishes (n = 13), and arthropods (n = 16). Mitogenomes were downloaded from GenBank or assembled de novo from Sequence Read Archive (SRA) data with MitoFinder v 1.4 (Allio et al. 2020) using default settings and the most closely related species with a published mitogenome as a template (supplementary file S1, Supplementary Material online). The 13 mitochondrial protein-coding genes for all species of a particular taxonomic group were aligned by MUSCLE (Edgar 2004), cleaned of start and stop codons, and concatenated in alphabetical order. Then, the three species of each triad were parsed into 154 individual alignments in which selection was analyzed.

Individual genes and concatenated sets of 13 genes (i.e. the mt coding sequence) were analyzed using codeML in PAML v. 4.7 (Yang 2007) to calculate branch-specific values of the d_N/d_S ratio, a measure of nonsynonymous substitution bias. Commonly, d_N/d_S ratios less than one are thought to indicate purifying selection, ratios above one to indicate positive selection, and ratios close to one to represent neutrality, although this is conservative when considering an entire coding sequence (CDS) because all mt genes are expected to be under strong purifying selection (d_N/d_S < 1; Fig. 1a). As such, our interest was mainly in relative increases or decreases in d_N/d_S among related species and the relationship between relative d_N/d_S and ecological variables, rather than absolute d_N/d_S values. For each triad, we calculated d_N/d_S on the concatenated CDS under five different models (Fig. 1b): (1) a three-rate model giving the high and low taxa independent rates relative to the combined ancestral and outgroup branch, (2) a two-rate model grouping the low taxon with the outgroup and ancestral branch, (3) a “two-rate-reversed” model grouping the high taxon with the outgroup and ancestral branch, (4) a “two-rate-equal” model grouping the high and low taxa and their ancestral branch together, relative to the outgroup, and (5) a one-rate (null) model with the same d_N/d_S applied to all branches (Fig. 1b). We removed any triad where d_N/d_S of the whole CDS on any branch was calculated to be greater than 1, an unrealistically high value for a whole mt CDS. Each triad had a minimum of 5/13 genes available for analysis in all 3 species. Model suitability was assessed with likelihood ratio tests and Akaike's Information Criterion (AIC). We also used the three-rate model to calculate d_N/d_S for each gene individually within each triad. We removed all runs for individual genes with predicted d_N/d_S values for any species, under any model, that were either 0.001 or 999, indicating insufficient genetic variation for the model to converge on an accurate value.

The concatenated CDS was then used as the basis for running RELAX (Wertheim et al. 2015), a phylogenetic test for relaxed selection, on the Datamonkey webserver (https://datamonkey.org/relax). The high-elevation taxon was selected as the test branch and all other branches were reference branches. RELAX examines the distribution of d_N/d_S values among sites on the test branch compared to the distribution in the reference branches. The assumption is that when selection is relaxed, d_N/d_S values converge toward one, while when selection is intensified, some sites experience intensified positive selection and others experience intensified purifying selection. Thus, the distribution of d_N/d_S values moves away from one and towards both extremes on test branches when intensified selection is responsible for elevated d_N/d_S values. RELAX calculates a k-value, with k-values below one indicating relaxed selection on test branches, and values above one indicating intensified selection (Fig. 1a). The program then performs a likelihood ratio test for the preferred k-value relative to a null model where test and reference branches have the same distribution of d_N/d_S values. The test is agnostic as to the identity and physicochemical properties of amino acid substitutions, only considering whether they are synonymous or nonsynonymous.

Analysis of Ecological Data

For each species in each triad in all six taxonomic groups, elevational range limits were taken from the International Union for the Conservation of Nature (IUCN) and/or several other expert sources or determined from analysis of Global Biodiversity Information Framework (GBIF) records (www.gbif.org), range maps, and habitat information (supplementary table S1, Supplementary Material online). We used the highest maximum and lowest minimum reputable elevational limits for analyses and calculated midpoint elevations and latitudes, with midpoint latitude for species with ranges that straddle the equator taken as the point midway between the equator and the absolute maximum latitude. For mammals, birds, reptiles, and amphibians, IUCN shapefiles were downloaded from https://iucnredlist.org when available and latitudinal limits and range sizes were extracted using a custom Python script (see Data Availability; supplementary file S1, Supplementary Material online). If IUCN shapefiles were unavailable, latitudinal limits and range size were estimated from available range maps, habitat distributions, and descriptive information where possible. Latitudinal limits and range size were not calculated for fishes and arthropods because shapefiles and expert range maps were typically unavailable. Body size was recorded from one of several published databases or drawn from the literature. Body size was recorded as mass for mammals and birds and log-transformed for analysis, and as snout-vent length for reptiles and amphibians and standard length for fishes. Size was not recorded for arthropods as it was typically unavailable.

Statistical Analyses

We compared d_N/d_S between high and low branches under the three-rate model (Fig. 1b) using paired, one-sided Wilcoxon signed-rank tests. Comparisons of the full mt CDS and individual genes were made within taxonomic groups with significance evaluated both before and after Bonferroni correction by the number of independent tests. We used Phylogenetic Generalized Least Squares (PGLS; Grafen 1989) implemented with the R packages APE (v. 5.7.1) and phangorn (v. 2.11.1) to test for associations between d_N/d_S values or k-values and various ecological factors. For the correlation matrix in PGLS, we used a phylogeny of the 138 high-elevation vertebrate species from each triad inferred via maximum likelihood based on a MUSCLE-generated alignment (Edgar 2004) of amino acid sequences from all mt CDS input into RAxML v. 8 using default parameters under a PROTGAMMAWAG model (Stamatakis 2014; supplementary file S2, Supplementary Material online). We also used PGLS to test for relationships between both absolute and relative d_N/d_S or k and ecological factors. In other words, for relative values, the predictor variable was the ratio between high- and low-taxon values for the ecological factor, and the response was the ratio between high- and low-taxon values for d_N/d_S or k. Ecological factors examined were absolute maximum, minimum, and midpoint elevations, maximum, minimum, and midpoint latitude, range size, and body size of the high-elevation species, as well as these values relative to values for the low-branch species. Variables were naturally log-transformed when necessary to meet assumptions of normality.

Selection on mtDNA Across Wide-ranging Clades

Assembly of Clades

To complement our analyses on closely related triads and more generally examine selective pressures on mtDNA across elevation, latitude, and range size, we also surveyed 6 family- or sub-order-level clades. These groups were characterized by abundant mitogenomes in NCBI and Eurasia-centered distributions typically including both the Tibetan plateau and a high degree of latitudinal variability. Unlike in the triad-based analyses of focal high-elevation taxa, species were only included if they had complete mitogenomes in NCBI, IUCN shapefiles, and elevational data (supplementary table S2, Supplementary Material online). This meant that no freshwater fish or arthropod groups were analyzed. There were three endotherm groups, Cercopithecidae (Old World monkeys; n = 72 species), Bovidae + Moschidae (sheep, goats, cattle, and musk deer; n = 88), and Phasianoidea (Pheasants, New World quail, guineafowl, and curassows; n = 67), as well as three ectotherm groups, Cryptobranchoidea (primitive salamanders; n = 34), Ranoidea (true frogs and relatives; n  = 65), and Lacertoidea (wall lizards and whiptails; n = 55). We downloaded mitogenomes and aligned them by clade as above. We used RAxML v. 8 to construct phylogenies for each group from an alignment of the thirteen mt protein-coding genes as above, except using a GTRGAMMA model on a nucleotide alignment for the Bovidae + Moschidae (due to low amino acid divergence, last common ancestor: 16 to 19 Ma; Bibi 2013) with six to seven outgroup taxa included for each group (supplementary table S2; file S2, Supplementary Material online).

Analysis of Ecological Data

In addition to extracting elevational limits from the literature and extracting latitudinal limits and range size from IUCN shapefiles as in the triad analyses, we calculated a measure which we call the Climatic Index (CI), which combines latitude and elevation to express a species’ position in both climatic dimensions relative to the rest of its clade:

$$CI=100*\sqrt{((L_{mid}/L_{Max})*(E_{mid}/E_{Max}))}$$

where L_mid is the latitudinal midpoint of a given species’ range and E_mid is its elevational midpoint, while L_Max and E_Max are the maximum latitude and elevation attained by any species in the group. Thus, the dimensionless CI expresses with one number how close a species’ midpoints are to the theoretical maxima for its clade (e.g. a species with elevational and latitudinal midpoints at 50% of the theoretical maximum would have a CI of 50.)

Analysis of Selection Across Clades

We ran codeML in PAML for each group by binning species into rate categories based on midpoint elevation, midpoint latitude, CI, or range size (6 groups × 4 parameters = 24 runs). Terminal branches of the phylogeny were assigned to bins of width 500 m for midpoint elevation, 5° for midpoint latitude, 5 dimensionless units for climatic index, and base-e exponents for range size in km². These multi-rate analyses thus had between 7 and 13 rate bins with 1 and 27 species (terminal branches) per bin. All internal branches were unlabeled, making them one binned “background” rate for the analysis. We also ran codeML on a one-rate (null) model with no labeled branches, which produced a value close to that for the “background” bin of the multi-rate analysis. d_N/d_S values for each bin were divided by the tree-wide d_N/d_S values for the null, one-rate model and plotted as a regression over the parameter of interest. The significance of the binning scheme relative to the null model was assessed using likelihood ratio tests. We then performed RELAX on each bin of the aforementioned phylogenies, treating the binned taxa as “test” branches and the internal branches—corresponding to the unlabeled “background” d_N/d_S category—as the “reference” branches, to test for relaxed selection as a function of different ecological factors.

Results

High-elevation Taxa Have Elevated d_N/d_S Relative to Low-elevation Sister Taxa

In total, we parsed, aligned, and concatenated 791 mt CDS including 41 novel mt CDS derived from publicly available transcriptomes, ultra-conserved element (UCE), and whole-genome sequencing (WGS) datasets (supplementary tables S1 and S2; file S1, Supplementary Material online). We first calculated d_N/d_S and k-values on these CDS for 153 three-species triad comparisons (supplementary table S3; file S3, Supplementary Material online). In 95 out of 153 triads (62.1%), the three-rate model was preferable to the one-rate model by likelihood ratio test at P = 0.05. In 58/95 of these triads (61.1%), the high-elevation taxon had a higher d_N/d_S value. After Bonferroni correction (P = 0.0003) the three-rate model was preferred in 40 out of 153 triads (26.1%), with 21/40 (52.5%) having higher d_N/d_S in the high taxon. d_N/d_S values for concatenated alignments under the three-rate model were always less than 1.0, indicating purifying selection in line with expectations for animal mitochondrial genomes as a whole (Fig. 2c, “All”;  supplementary table S4, Supplementary Material online). Individual genes showed greater variability, including values over 1.0 (Figure 2b, supplementary table S4, Supplementary Material online). Across all triad comparisons, the median d_N/d_S value on the high branch (0.0458) was 12.0% greater than that on the low branch (0.0408; P = 0.008, V = 7211.5, d.f. = 152; paired, one-sided Wilcoxon signed-rank test; Fig. 2c; supplementary table S4, Supplementary Material online). A similar trend was found among the 95 comparisons where the three-rate model was preferred over a one-rate model at P = 0.05 (7.1% higher on the high branch; P = 0.019, V = 2837, d.f. = 152; Fig. 2c  “LRT”). Results from other models confirmed the finding of higher d_N/d_S in the high-elevation taxon and supported the use of the three-rate model in subsequent analyses (supplementary results, Supplementary Material online).

Fig. 2.

d_N/d_S ratios for paired high- and low-elevation species (light blue and brown, respectively) under the three-rate model. a) d_N/d_S ratios for the whole mitochondrial CDS by taxon. b) d_N/d_S ratios for all taxa by individual gene. c) d_N/d_S ratios for the whole mitochondrial CDS under different constraints: All: 153 comparisons, LRT: only those in which the three-rate model is preferred to the 1-rate (null) model by likelihood ratio test, Mid-1000: only those in which the high and low species differ by at least 1,000 m in midpoint elevation, Mid-2000: only those in which they differ by 2,000 m, No-overlap: only those in which high and low species have absolutely no overlap in elevational ranges, Larger hi-range: only those in which the high-elevation species has a larger range than the low, Temperate: only those where both taxa have midpoint latitudes above 30°, Tropical: only those where both taxa have midpoint latitudes below 23.5°, Ectotherm: only those among reptiles, amphibians, fishes, and arthropods, and Endotherm: only those among mammals and birds. Black lines connect medians while blue lines connect comparisons where the high-elevation species has higher d_N/d_S and brown lines connect pairs where the opposite is true. Asterisks indicate uncorrected significance by paired Wilcoxon signed-rank test (* = P < 0.05; ** = P < 0.005).

There was some variation in overall d_N/d_S ratios among taxonomic groups, with reptiles tending to display the highest values (Kruskal–Wallis test on one-rate model: P = 0.0008, X² = 20.897, d.f. = 5). The trend of moderately higher d_N/d_S ratios in high-elevation taxa under the three-rate model held across terrestrial groups, but not in fishes (Fig. 2a, supplementary table S5, Supplementary Material online). Only in amphibians was the difference in medians statistically significant (P = 0.022, V = 204, d.f. = 22), and in no group was it significant after Bonferroni correction (P < 0.008). When examining differences in relative high-taxon d_N/d_S between genes, there was a trend toward weakly elevated d_N/d_S ratios across the mitogenome, but this was not statistically significant in any particular gene when analyzed alone (Fig. 2b, supplementary table S6, Supplementary Material online). ATP6 and ATP8 tended to display higher d_N/d_S ratios than other genes overall (Kruskal–Wallis test on one-rate model: P < 0.0001, X² = 760.74, d.f. = 12; supplementary table S6, Supplementary Material online). We also examined d_N and d_S separately to help understand variation in d_N/d_S ratios. Both d_N and d_S of the concatenated, 13-gene alignments were log-normally distributed (supplementary figs. S1 and S2, Supplementary Material online), with no significant differences among taxonomic groups (supplementary fig. S3 and table S7, Supplementary Material online). Considered separately, d_N, but not d_S, was higher in the high-elevation branches, suggesting that the trend in d_N/d_S is driven by differences in selection rather than mutation rate (supplementary table S8, Supplementary Material online).

We constrained our dataset of all triads to several subsets to test whether d_N/d_S was consistently higher in high-elevation species (Fig. 2c). In all but two cases, median d_N/d_S was statistically significantly higher among the high-elevation species prior to Bonferroni correction (supplementary table S9, Supplementary Material online). This included subsets where the difference in midpoint elevations was at least 1,000 m (“Mid-1000”) or 2,000 m (“Mid-2000”) or the pairs did not overlap at all in elevation (“No-Overlap”), and subsets of only temperate taxa (“Temperate”), tropical taxa (“Tropical”), or endothermic taxa (“Endotherm”). The only subsets in which the difference in median d_N/d_S was not statistically significantly different was across pairs where the high-elevation species had a larger range size than its low-elevation sister species (“Larger hi-range”; P = 0.386, V = 316, d.f. = 33) and among the ectotherms (“Ectotherm”). After Bonferroni correction (P < 0.004), no comparison was statistically significant except for that within temperate species (P = 0.0036, V = 527, d.f. = 36; supplementary table S9, Supplementary Material online).

Relative Elevation and Range Size Predict Relative d_N/d_S

We used high- and low-taxon d_N/d_S values from the three-rate model as the response variable in phylogenetic least-squares regressions against ecological variables. Among the high-elevation species, the range-wide elevational midpoint was not a significant predictor of d_N/d_S (P = 0.7056, residual d.f. = 135, T = 0.3785; PGLS). However, when examining relative rather than absolute values, there was a weak but statistically-significant phylogenetic correlation between the log-transformed ratio of the high/low-taxon elevational midpoints and the ratio of their d_N/d_S values (PGLS; P = 0.0377, residual d.f. = 135, T = 2.098; Fig. 3a). Like elevational midpoint, there was no correlation between high-taxon latitudinal midpoint and d_N/d_S under the three-rate model (P = 0.2626, residual d.f. = 104, T = 1.126; PGLS). There was also no correlation between the ratio of latitudinal midpoints and the ratio of d_N/d_S values (P = 0.6855, residual d.f. = 104, T = −0.4061; PGLS). There was no correlation between log-transformed range size and d_N/d_S values (P = 0.472, residual d.f. = 99, T = 0.7219; PGLS). However, there was a weak but significant negative correlation between the log-transformed ratio of high/low-taxon range sizes and the ratio of d_N/d_S values (P = 0.0212, residual d.f. = 99, t = −2.342; PGLS; Fig. 3b). The square-root of range size decreased strongly with minimum elevation among high-elevation taxa, as predicted (P < 0.0005, T = −3.845, d.f. = 99; Fig. 3c). We also tested for correlations between d_S and latitude, elevation, and range-size in the same dataset, none of which were significant (supplementary table S10, Supplementary Material online).

Fig. 3.

Molecular evolution in animal mtDNA as predicted by ecological correlates. a) Relative high-species d_N/d_S increases with relative high-species midpoint elevation in a phylogenetically controlled correlation (PGLS) using whole-CDS d_N/d_S ratios. b) Relative high-species d_N/d_S decreases with relative high-species range size under a PGLS model. c) The square-root of range size decreases with minimum elevation under a linear model, illustrating how higher-living species have smaller ranges. d) k-values decrease with relative high-species d_N/d_S under a PGLS model. “rse” = residual standard error, a more appropriate measure of goodness-of-fit than r² for PGLS models. Shading around lines of best fit shows 95% confidence intervals.

Higher Relative d_N/d_S Values are Associated With Relaxed Selection

There was a nearly-significant negative correlation between the natural log of k, the parameter for relaxed selection, and the high/low ratio of d_N/d_S values (P = 0.518, residual d.f. = 130, T = −1.962; PGLS; Fig. 3d), suggesting that with relatively greater high-species d_N/d_S, selection was more likely to be relaxed rather than intensified. Log-transformed k-values were not predicted either by the log-transformed high/low ratio of midpoint elevations (P = 0.3337, d.f. = 130, T = −0.9702; PGLS) or by the log-transformed high/low ratio of midpoint latitudes (P = 0.6981, residual d.f. = 103, T = −0.389; PGLS). Log-transformed k-values were also not predicted by the log-transformed high/low ratio of range sizes (P = 0.969, residual d.f. = 98, T = 0.0389; PGLS) or the log-transformed ratio of high/low body sizes (P = 0.3049, residual d.f. = 113, T = 1.0307; PGLS). Log-transformed k-values were also not predicted by the ratio between high- and low-elevation d_S values (supplementary table S10, Supplementary Material online).

Body Size Predicts Absolute d_N/d_S Among Endotherms

Because body size was measured as mass for endotherms and length for ectotherms, we did not perform a regression of absolute body size with d_N/d_S from the three-rate model across the whole dataset. However, among the endotherms, the relationship between log-transformed body mass and d_N/d_S was positive and highly statistically significant (P < 0.0005 residual d.f. = 77, T = 4.094; PGLS; Fig. 4a). Among the ectotherms, the relationship between body length and d_N/d_S was not significant (P = 0.5348, residual d.f. = 35, T = 1.921; PGLS; Fig. 4c). Finally, we tested for correlations between body size in both ectotherms and endotherms and d_S; we found a significant negative relationship between body mass and d_S in endotherms, but no relationship in ectotherms (supplementary fig. S4 and table S10, Supplementary Material online).

Fig. 4.

d_N/d_S ratios relative to body size among high-elevation a) endotherms and c) ectotherms. Endotherms have a highly significant positive relationship between log body mass and d_N/d_S while ectotherms have a non-significant negative relationship between body length and d_N/d_S. Box plots show lack of a difference in body size between high- and low-elevation b) endotherms and d) ectotherms, suggesting trends are not a consequence of confounding elevation with body size. “rse” = residual standard error, a more appropriate measure of goodness-of-fit than r² for PGLS models.

K-Values Decrease While d_N/d_S Increases With Elevation in Most Clades

We also examined how selection varied with ecological factors across 6 families and suborders, finding that the relationship between ecological factors, d_N/d_S, and k-values varied among and within clades (supplementary table S11; file S4, Supplementary Material online). In nearly all analyses (23/24), multi-rate d_N/d_S models were preferable to the null, 1-rate model by likelihood ratio test (LRT) at P < 0.05, and after Bonferroni correction (P < 0.002; 22/24 cases). There was, however, a large degree of variability (P-values varied by over 200 orders of magnitude). Most of the variation was among the 6 taxonomic groups, with the Cercopithecidae displaying the most support for multi-rate models and Phasianoidea the least. The comparisons with the most support for the multi-rate model, such as Cercopithecidae binned by CI or Ranoidea by elevation, tended to have one or several highly variable bins, suggesting the highly significant likelihood ratio tests reflect the accommodation of this variability. These were not, however, necessarily the comparisons with the best evidence for linear trends between molecular evolution and ecological parameters. The significance of k-values is presented in terms of each bin, reflecting the results of the likelihood ratio test for relaxed selection among those taxa.

We calculated relative d_N/d_S of each equal-width bin as the d_N/d_S for that bin divided by the d_N/d_S for the whole phylogeny under the null (one-rate) model. We visualized the trend in relative d_N/d_S across bins within a clade with a linear regression weighted by the number of taxa in that bin. Because the binning scheme was not compatible with PGLS, we discuss the change in d_N/d_S in terms of trends without testing for the statistical significance of these linear regressions. For the same reason, we also discuss trends in k-values across the parameter space weighted by sample size rather than the statistical significance of these linear models.

Relative d_N/d_S ratios tended to increase or remain constant with midpoint elevation for all groups, supporting the triad analyses suggesting moderately elevated d_N/d_S ratios among high-elevation taxa (Fig. 5a). k-values tended to decrease or remain constant with increasing elevation in 5/6 groups and increased in one, the Ranoidea (Fig. 5b). Based on k-value likelihood ratio tests, selection was significantly intensified at low elevations in Cercopithecidae and significantly intensified at high elevations in Ranoidea; by contrast, selection was significantly relaxed at high elevation in 3/6 groups, the Bovidae + Moschidae, Cryptobranchoidea, and Lacertoidea.

Fig. 5.

Trends in d_N/d_S and k-values across midpoint elevation in six family- or sub-order-level clades. a) d_N/d_S ratios (y-axis) for species binned by midpoint elevation, relative to the null rate. The dashed line at y = 1 indicates equivalency between the bin rate and the null rate; points above this line are relatively higher d_N/d_S, indicative of intensified positive selection or relaxed purifying selection. The size of the points is proportional to the number of species in the bin. Solid lines are linear trend lines weighted by the number of species. b) k-values (y-axis) of each bin relative to the internal branches of the phylogeny. Points above the dashed line at y = 1 are indicative of intensified selection (intensified purifying and/or positive selection) while points below the line are indicative of relaxed selection. The size of points is proportional to the number of species in the bin. White points had non-significant k-values, while gray points were significant at P < 0.05 and black points were significant after Bonferroni correction (P < 0.002). Silhouettes by Nat Jennings.

With increasing midpoint latitude, relative d_N/d_S ratios tended to increase in Bovidae + Moschidae and Lacertoidea while decreasing or remaining constant in the other groups (supplementary fig. S5a, Supplementary Material online). k-values tended to increase with latitude in Cercopithecidae, Phasianoidea, and Ranoidea, but decrease in the other three groups (supplementary fig. S5b, Supplementary Material online). Based on k-value likelihood ratio tests, selection was significantly relaxed at high latitudes in Bovidae + Moschidae, Cryptobranchoidea, and Lacertoidea, while being significantly intensified at high latitudes in Phasianoidea and Ranoidea. Ranoidea also showed evidence for significantly relaxed selection at low latitudes, while in Cercopithecidae, both relaxed and intensified selection was significant at low latitudes.

Relative d_N/d_S ratios tended to increase with increasing climatic index (incorporating both elevation and latitude) in Bovidae + Moschidae and Lacertoidea, while remaining relatively flat in the other groups (supplementary fig. S6a, Supplementary Material online). k-values tended to increase with climatic index in Phasianoidea, Cryptobranchoidea, and Ranoidea, were relatively static in Cercopithecidae, and tended to decrease in Bovidae + Moschidae and Lacertoidea (supplementary fig. S6b, Supplementary Material online). Based on k-value likelihood ratio tests, selection was significantly intensified at high climatic indexes in Ranoidea while being significantly relaxed at high climatic indexes in Bovidae + Moschidae and Lacetoidea. In Phasianoidea, selection was significantly relaxed at low climatic indexes, while in the other two groups, trends were more ambiguous.

With increasing range size, relative d_N/d_S ratios tended to increase in Cercopithecidae, Cryptobranchoidea, and Ranoidea while decreasing in the other 3 groups (supplementary fig. S7a, Supplementary Material online). k-Values tended to increase with range size in Cercopithecidae, Bovidae + Moschidae, Phasianoidea, and Cryptobranchoidea, while decreasing in Ranoidea and remaining relatively static in Lacertoidea (supplementary fig. S7b, Supplementary Material online). Based on k-value likelihood ratio tests, selection was significantly relaxed at the smallest ranges in Bovidae + Moschidae, Phasianoidea, and Cryptobranchoidea, while in Lacertoidea selection was significantly relaxed in some bins across all range sizes. Selection was significantly intensified at the largest ranges in Cercopithecidae, Phasianoidea, and Cryptobranchoidea, while in Ranoidea responses were highly idiosyncratic: selection was significantly intensified at the smallest range sizes, and both intensified and relaxed in different bins of the largest range sizes.

Discussion

d_N/d_S Ratios Tend to Increase With Elevation

We detected a modest but widespread signature of elevated d_N/d_S ratios in the mtDNA of high-elevation species relative to their low-elevation sister taxa across several comparisons. This effect holds for almost all taxonomic groups, across mitochondrial genes, and is consistent in magnitude and direction under different inclusion criteria (Fig. 1c). The effect is also detectable in most large clades examined (Fig. 5a). This indicates a weak but pervasive bias in nonsynonymous substitution rates in the mitochondrial CDS of high-elevation species, echoing previous results describing this pattern in specific species, within certain families, and among vertebrates as a whole (Xu et al. 2005; Hassanin et al. 2009; Wang et al. 2021) and suggests that it may be a general pattern at least for terrestrial animals.

Importantly, this modest increase in d_N/d_S could reflect either positive or relaxed selection. Under positive selection, we would expect d_N/d_S to increase more at higher elevations. While relative elevation difference weakly predicted relative d_N/d_S across the triads (Fig. 3a), constraining the dataset to stricter cases (such as large differences in midpoint elevations) barely changed median d_N/d_S ratios or the difference between high and low medians (Fig. 2c, supplementary table S9, Supplementary Material online). In the clade-based analysis, only the Ranoidea displayed increasing k-values and significantly intensified selection with elevation (Fig. 5b). However, the particular bins that had significantly intensified selection were those with d_N/d_S ratios below neutrality (Fig. 5a). This indicates that intensified purifying, not positive selection had taken place in high-elevation species in Ranoidea, a potential signature of selection for energetic efficiency and has been documented in some groups (Graham et al. 2024). By contrast, all other groups showed flat or decreasing trends in k with elevation (Fig. 5b), and significantly relaxed k-values were often found in bins with relatively high d_N/d_S ratios.

Adaptation might be expected to cause similar trends for species living at high latitudes as for those living at high elevations, if temperature is an important environmental variable. However, there was no effect of absolute or relative latitude on d_N/d_S ratios across triads, and relationships were mixed among the clade analyses (supplementary fig. S5, Supplementary Material online). Among clades, no group showed unambiguous evidence for positive selection (elevated d_N/d_S and significantly intensified k-values) with midpoint latitude or CI. Responses of d_N/d_S among clades to the CI, which incorporates both elevation and latitude, tended to resemble responses to midpoint elevation more than latitude (supplementary figs. S5a and S6a, Supplementary Material online), perhaps suggesting a greater impact on d_N/d_S of environmental variation attributable solely to elevation over temperature. This could mean that hypoxia is a more important selective pressure on mtDNA than temperature. In support of this, freshwater fishes are the only group of triads that trended toward higher d_N/d_S ratios among low-elevation, rather than high-elevation, species. While high-elevation waters still tend to be less well-oxygenated than those at low elevation, the increase in solubility of gasses in cold water partially compensates for the decrease in atmospheric oxygen pressure (Jacobsen and Dangles 2017a, 2017b). Tumbling mountain rivers experienced by some high-elevation species may tend to be well-oxygenated, while lowland waters experienced by their relatives are sluggish and prone to hypoxia. In support of this, mt adaptation for greater oxygen affinity appears to play a role in adapting fishes to warmer waters (Chung et al. 2017). Therefore, even intermittent low-oxygen conditions could be driving a signal of mt adaptation in low- rather than high-elevation fishes. Unfortunately, the lack of distributional information for freshwater fish and empirical oxygen series from water at different elevations limited our ability to assess this hypothesis in a clade-based framework.

Convergent adaptation in response to elevation might be expected to result in elevated d_N/d_S in one or a few specific genes across many taxa. However, relaxed selection in high-elevation species due to lower N_e should be a genome-wide phenomenon affecting all genes (including nuclear-encoded genes). Our finding of similar patterns in high versus low elevation d_N/d_S values across all mt genes (Fig. 1b) might therefore better reflect relaxed selection. However, because the mt genome is one non-recombining linkage group, this result could also reflect genetic draft (hitchhiking) in mt genes rather than genome-wide relaxed selection (Hill 2020). Within genes, there was high variability between species and vice versa, potentially masking important differences in the strength of selection for individual cases. Moreover, our method examining d_N/d_S across whole genes and mitogenomes prevents the identification of positive selection on one or a few critical amino acid residues that may drive adaptation to high elevation. Relatively few studies have pinpointed individual amino acid substitutions in mtDNA and supported them with plausible adaptive scenarios at high elevations (see Scott et al. 2011; Ji et al. 2012; Kostin and Lavrenchenko 2018) and adaptation may often occur via different substitutions in the genes of different species (Natarajan et al. 2016; Storz 2017). For instance, a recent analysis showed almost no repeatability of substitutions in relation to high-elevation adaptation in the mtDNA of birds, nor in relation to several other selective pressures (Burskaia et al. 2021). The low-repeatability of molecular adaptation to the environment probably reflects the diversity of functional alterations that can be achieved within even one enzyme (Natarajan et al. 2016).

Relaxed Selection as an Alternative Hypothesis for Elevated d_N/d_S Ratios

We explicitly investigated relaxed selection within triads by calculating k-values and examining the relationship between d_N/d_S and k. While k-values were generally centered around 1, relatively high d_N/d_S values tended to reflect relatively lower k-values (Fig. 3d). This suggests that greater relative d_N/d_S ratios among high-elevation species are associated with more relaxed selection, supporting the results discussed above from our clade-based analyses. This evidence for relaxed selection is supplemented by indirect evidence from the relationship between d_N/d_S and several ecological variables which are reflective of N_e. A positive relationship between N_e and the effectiveness of selection is predicted from theory (Otto and Whitlock 1997), and relaxed selection has been detected in lineages undergoing reductions in effective population size, such as arthropods transitioning to eusociality (Chak et al. 2021; Weyna and Romiguier 2021). In our dataset, species at higher elevations had smaller ranges (Fig. 3c) while both higher relative elevation and smaller relative ranges predicted higher relative d_N/d_S (Fig. 3, a and b). The only subsets of our full dataset in which the high–low d_N/d_S difference was not significant were among the ectotherms and in taxa where the high-elevation species had a larger range (Fig. 2c; supplementary table S9, Supplementary Material online). Many of these species live on the vast Qinghai-Tibetan Plateau, where adjacent lower-elevation habitats occupy smaller areas (see supplementary table S1, Supplementary Material online).

While relative range size did not predict k-values across triads, results from the clade-based analysis lent support to a relationship between the two. In the Bovidae + Moschidae, Phasianoidea, and Lacertoidea, bins with the smallest range sizes had elevated d_N/d_S and significantly relaxed k-values (supplementary fig. S7b, Supplementary Material online). In two other groups, Cercopithecidae and Cryptobranchoidea, bins with the largest range sizes had evidence for positive selection, part of the same hypothetical continuum. Thus, there is support from 5/6 groups for a positive relationship between range size and the effectiveness of selection. Together, these findings suggest that faster mtDNA evolution in high-elevation species may be caused more by relaxed selection than adaptive evolution.

Results from the analysis of body size also suggest that demographic factors are more important than adaptation for overall substitution rates. Across our triads, endotherms exhibited a strong, positive relationship between body size and d_N/d_S, which may be a signal of relaxed selection due to low effective population sizes (Fig. 4a). Body size has a well-supported inverse relationship with population size due to metabolic scaling relationships and habitat-carrying capacities (Cotgreave 1993; Savage et al. 2004). Additionally, d_S values decreased significantly with body size in endotherms (supplementary fig. S4, Supplementary Material online), which likely reflects the effect of generation time. However, among ectotherms, the relationship between body size and both d_N/d_S and d_S was non-significant (Fig. 4c). In ectotherms, the relationship between body size and population size is apparent at warm temperatures, but less apparent at cool temperatures where ectotherm energetic needs are lower (Buckley et al. 2008). The metabolic flexibility of ectotherms and the intentional inclusion of many high-elevation and high-latitude species in our dataset may explain why only endotherms show the expected relationship between body size and d_N/d_S. If environmental adaptation rather than low population size had been driving d_N/d_S ratios among high-elevation species, we might have expected smaller endotherms, which are more susceptible to heat loss, to exhibit the strongest signals of selection.

Overall, we suggest that relaxation of purifying selection at high elevation is the most likely cause of modestly increased rates of overall mtDNA evolution. Small range sizes and low N_e are likely when transitioning upward in elevation range, and might thus provide a better explanation for mitogenome-wide elevated d_N/d_S than pervasive positive selection, which would be less efficient on all loci. However, this does not preclude the possibility that a number of important substitutions have also been fixed by positive selection, as has been functionally demonstrated in a number of cases (Scott et al. 2011; Ji et al. 2012; Kostin and Lavrenchenko 2018).

The Importance of Distinguishing Positive From Relaxed Selection

Tests of relaxed selection paired with d_N/d_S estimates can distinguish between positive, relaxed, and purifying selection and neutrality, allowing discrimination among competing hypotheses. Without this framework, alternative stories can be told from the same data and opposing datasets can be used to support the same hypothesis (Wertheim et al. 2015; Zwonitzer et al. 2023). For example, island species are thought to have severely reduced N_e relative to their mainland ancestors, and studies have detected a strong signal of elevated d_N/d_S ratios among island species relative to paired mainland taxa, arguing that this was a consequence of increased fixation of slightly deleterious substitutions due to low N_e (Johnson and Seger 2001; Woolfit and Bromham 2005). However, another study with a similar paired design found the opposite: birds from large landmasses had higher overall rates of substitution than their island counterparts, with a trend toward higher d_N/d_S (Wright et al. 2009). In this case, the authors suggested that higher d_N/d_S values in mainland species resulted from more effective natural selection in large populations (intensified positive selection). This illustrates how opposite results can theoretically fit the same hypothesis, making such hypotheses hard to falsify based solely on d_N/d_S ratios and similar tests. Our dataset contains examples of higher d_N/d_S being associated with both smaller and larger range sizes (supplementary fig. S7, Supplementary Material online), and RELAX analyses allowed us to distinguish between relaxed and positive selection.

A bias to interpret an adaptive signal in d_N/d_S ratios may be especially prevalent with animal mtDNA, which is increasingly available as whole mitochondrial genomes (Zwonitzer et al. 2023). Because mitochondrial DNA influences animal energetics (Hill et al. 2019; Norin and Metcalfe 2019), adaptation may be invoked to explain any signal that correlates with an interesting environmental or physiological niche. For example, after reanalyzing energetic hypotheses based on elevated d_N/d_S ratios in animal mtDNA, we found that many previous cases hypothesized to demonstrate adaptive mtDNA evolution may be driven by relaxed selection instead (Zwonitzer et al. 2023). Other recent studies have also called into question the generality and adaptive significance of mtDNA substitution bias in relation to energetics (Burskaia et al. 2021; Claramunt and Haddrath 2023). A growing body of work takes elevated mtDNA substitution rates as a de facto signal of adaptation to elevation (see Introduction). However, relaxed selection on mtDNA should be explicitly tested in a phylogenetic framework, as some studies have begun to do (i.e. Gutiérrez et al. 2023).

Here, we found that while d_N/d_S ratios in mtDNA do tend to be slightly elevated across high-elevation animals (Fig. 1), there is ample variation among lineages that remains to be explained and the causes of this pattern are not straight-forward. Our work does not exclude the possibility of mtDNA adaptation to elevation or cast doubt on specific studies that have highlighted particular adaptive substitutions; however, our triad- and clade-based explicit tests tended to show more evidence for relaxed selection than adaptation shaping overall mtDNA evolution at high elevation. When combined with ecological data on species ranges, our results lend support to an alternative explanation for elevated d_N/d_S ratios: high-elevation species tend to inhabit smaller ranges and likely live at lower effective population sizes, causing an uptick in slightly deleterious nonsynonymous substitutions as purifying selection is relaxed. This demographic scenario explains the mitogenome-wide signature of higher d_N/d_S (and predicts a similar signature in nuclear loci, a testable hypothesis), without precluding adaptation that might occur at specific loci in an idiosyncratic fashion. Overall, we caution against over-interpretation of simple molecular signatures like d_N/d_S and reaffirm the need to rigorously test adaptive hypotheses and support them with thorough mechanistic investigation.

Supplementary Material

Supplementary material is available at Molecular Biology and Evolution online.

Acknowledgments

We thank Allyson Gunderson, Nat Jennings, and Adam Zambie for their help with data collection. We thank Mike Ryan, Misha Matz, Geoff Hill, Molly Schumer, and members of the Havird Lab for their feedback. We thank Kendra Zwonitzer for the script for concatenating mitochondrial genomes. This work was funded by the National Institutes of Health (1R35GM142836) and the Stengl-Wyer Endowment at the University of Texas at Austin. This article was written without the use of artificial intelligence.

Data Availability

The Python script used for retrieving, aligning, and concatenating mt CDSs is available at https://github.com/thekzwon/mito_accessions_to_full_alignment. The Python script used for extracting geospatial data from IUCN shapefiles is available at https://github.com/abbycriswell/spatial-data-processing.

Sequences generated for this study can be found on FigShare: doi.org/10.6084/m9.figshare.28459061. Supplementary Files S1, S2, S3, and S4 can also be found on FigShare: doi.org/10.6084/m9.figshare.28590725.

Author notes

Conflict of Interest: The authors declare no conflict of interest.