High‐altitude adaptation in vertebrates as revealed by mitochondrial genome analyses

Abstract The high‐altitude environment may drive vertebrate evolution in a certain way, and vertebrates living in different altitude environments might have different energy requirements. We hypothesized that the high‐altitude environment might impose different influences on vertebrate mitochondrial genomes (mtDNA). We used selection pressure analyses and PIC (phylogenetic independent contrasts) analysis to detect the evolutionary rate of vertebrate mtDNA protein‐coding genes (PCGs) from different altitudes. The results showed that the ratio of nonsynonymous/synonymous substitutions (dN/dS) in the mtDNA PCGs was significantly higher in high‐altitude vertebrates than in low‐altitude vertebrates. The seven rapidly evolving genes were shared by the high‐altitude vertebrates, and only one positive selection gene (ND5 gene) was detected in the high‐altitude vertebrates. Our results suggest the mtDNA evolutionary rate in high‐altitude vertebrates was higher than in low‐altitude vertebrates as their evolution requires more energy in a high‐altitude environment. Our study demonstrates the high‐altitude environment (low atmospheric O2 levels) drives vertebrate evolution in mtDNA PCGs.

(mtDNA PCG) encoding thirteen proteins involved in aerobic respiration (Chong & Mueller, 2013). To adapt to hypoxia environments, aerobic respiration in high-altitude vertebrates may undergo natural selection to increase efficiency. Therefore, the mitochondrial evolutionary rate may be affected by the low atmospheric O2 levels at high altitudes (Scott et al., 2009).
At present, several studies using mitochondrial genome analyses have demonstrated that the thirteen animal mtDNA PCGs have different environmental adaptations based on different evolutionary rates and selection constraints, including mammals (Gu et al., 2012;Luo et al., 2008), aves (Scott et al., 2009;Zhou et al., 2014), and actinopterygii (Li et al., 2013;Wang et al., 2016). For example, the NS/S values of the ATP6, ATP8, and Cytb genes were larger (>1) in the Tibetan than the Han population (Gu et al., 2012). Similarly, the bar-headed geese had a striking alteration in the kinetics of the cytochrome c oxidase (COX) in high-altitude environments (Scott et al., 2009). Wang et al. (2016) discovered that Tibetan loaches accumulated more nonsynonymous mutations and exhibited rapid evolution when compared to non-Tibetan loaches. It is possible that the mtDNA of low-altitude vertebrates might have a lower evolutionary rate. In contrast, the mtDNA of high-altitude vertebrates may have a faster evolutionary rate as they adapt to the low oxygen environment to maintain an efficient energy metabolism. Therefore, we tested the rate of mtDNA evolution between vertebrates living at different altitudes using the rate ratio of nonsynonymous/synonymous nucleotide substitutions (dN/dS (ω)). We aimed to determine whether there is a molecular basis for high-altitude adaptation that is related to vertebrate mitochondrial evolutionary rates. This adaptation to low atmospheric O 2 levels shows the mtDNA evolution is a critical molecular process enabling the survival of vertebrates at high altitudes and hypoxic environments driving the evolution of vertebrate mitochondrial genes.

| Species sample and mtDNA sequence data
Based on previous studies and the existing mitochondrial genome data of NCBI, we selected 104 vertebrate species (with habitats at varying altitudes; high-altitude vertebrate, 52; low-altitude vertebrate, 52). We downloaded mitochondrial genomes from the NCBI GenBank database (http://www.ncbi.nlm.nih.gov/). Our initial criterion was to only select species with clear information available on altitude. If there was no clear altitude information, our second criteria were to include species specific to the plateau. The third criteria were to select the low-altitude species related to the high-altitude species. The data set covered five taxa (mammal, aves, reptile, amphibian, and actinopterygii), and the data set of each taxon was greater than or equal to 20 mitochondrial genomes. For most families investigated within each class (mammal, aves, reptile, amphibian, and actinopterygii) in this study, more than one representative of each species was selected. The species and mitochondrial genome accession numbers are listed in Appendix S1.

| Phylogenetic construction
We retrieved the 13 mtDNA PCGs of the mitochondrial genomes from the mtDNA sequence data and aligned the mtDNA PCGs using MUSCLE v3.8.31 (Edgar, 2004). We obtained fourteen sequence datasets (Appendix S1). The phylogenetic relationships of the 104 vertebrate species were determined using 13 mtDNA PCGs datasets and their Bayesian inference (BI). We used Mega 7.0 to assess the base substitution saturation in the 13 mtDNA PCGs datasets and used DAMBE 6.4 (Xia, 2001) to calculate the I ss and I ss.c for testing the substitution saturation. The optimal model (GTI +T + G) was selected in the modelfinder function of PhyloSuite software  based on the greedy search algorithm and the Bayesian information criterion. The concatenated matrix was performed with four independent Markov Chain Monte Carlo (MCMC) chains for 2,000,000 generations and sampling one tree every 1000 generations.

| Selection pressure analyses
We used the one ratio (M0) model and the free ratio model (model = 1) to estimate the ratio of nonsynonymous (dN) to synonymous (dS) substitutions rates (ω = dN/dS) using the CodeML program in PAML version 4 (Yang, 2007). The likelihood ratio tests (LRTs) were used to determine which mtDNA PCG free ratio model was not effective models. The free ratio model allowed different ω values for each branch on a tree. The ω values of the 13 mtDNA PCGs from the 104 species and each mtDNA PCG with their mitochondrial genomes were computed separately for the terminal branches to evaluate the selective pressure. If the dN or dS values were equal to zero, the ω values could be smaller or larger and we did not use these ω values in our analysis, and we assigned these ω values as NA data in the ω dataset. The ω values from the 5 high-altitude groups (MH, BH, RH, AH, and FH) and 5 low-altitude groups (ML, BL, RL, AL, and FL) were compared using one-way ANOVA and multiple comparisons analysis. We used false discovery rates (FDR) to correct the p values. We used a Wilcoxon test to evaluate the statistical significance in the differences in ω values of 13 mtPCGs and each mtPSG between the species at different altitudes (high-altitude vertebrate, HV; low-altitude vertebrate, LV).
To identify the positive selection gene or rapidly evolving gene for high-altitude adaptation, we used a branch model (one ratio (M0) model, two ratio (M2) model, and NSsites = 0) to detect each mtDNA PCGs in the 104 species. In our research, the branches with highaltitude vertebrates were used as foreground branches, and the lowaltitude vertebrates were the background branches. Furthermore, we used a branch-site model (one ratio (M0) model, two ratio (M2) model, and NSsites = 2) in PAML to detect each mtDNA PCG of the 104 species (foreground branches, HV; background branches, LV).
We used LRTs to determine which mtDNA PCGs were positive selection genes or rapidly evolving genes. These analyses were based on the BI tree ( Figure 1).

| Phylogenetic independent contrasts analysis
The closely related species (shared inheritance) may affect the species comparative analysis. Therefore, we used the phylogenetic independent contrasts (PIC) analysis (Felsenstein Joseph, 1985) to remove the closely related species influences and explore the relationship between altitude and the ω values of the 13 mtDNA PCGs using the ape package in R software. Firstly, we used FigTree (v1.4.3) to transform the BI tree into a binary tree. The binary tree was used as the input file in the analysis. Secondly, we classified the 104 species into two groups (high-altitude and low-altitude). The high-altitude group and low-altitude groups were coded 0 and 1, respectively. We used the ω values, 0 and 1, as the character data in the PIC analysis.

| Phylogenetic analyses
We used 13 mtDNA PCGs datasets of the 104 vertebrates to conduct the phylogenetic analysis. The BI phylogenetic analyses of the concatenated datasets yielded consistent topological relationships between vertebrates with high bootstrap support values and Bayesian posterior probabilities. The BI tree is divided into 5 large branches: the mammal, aves, reptile, amphibian, and actinopterygii ( Figure 1).
Our results are consistent with previous studies (Townsend et al., 2008), and this was credible for the following analyses.

| Free ratio model analysis
The ω values (13 mtDNA PCGs and each mtDNA PCG of 104 vertebrates) were estimated using CodeML within the PAML package. The

F I G U R E 1
The BI phylogenetic tree of 104 vertebrates based on 13 PCGs of mitochondrial genomes (red background, high-altitude aves; green background, high-altitude mammal; yellow background, reptile; blue background, high-altitude actinopterygii; brown background, high-altitude amphibian) ω values <1 indicate a purifying selection, the ω values = 1 indicate a neutral selection, and the ω values >1 indicate a positive selection.
We found that the ω values were lower than 1 for the 13 mtDNA

| Positive selection analysis
We used the branch model and branch-site model analyses to test our hypothesis. The branch model analysis identified seven rapid evolutionary genes (LRTs, p < .05) in the HV, that is, the ND1 gene  Table 1). The results confirmed the effectiveness of the Wilcoxon test results. We used a branch-site model to detect the positive site of each mtDNA PCG in the HV. We found that only one positive selection gene (ND5 gene, corresponding to site 247 and 524) was detected in the HV using the branch-site model analysis (Table 2).

| Phylogenetic independent contrasts analysis
We used a PIC analysis to overcome the potential phylogenetic biases, for example, variation in the branch lengths and any phylogenetic inertia. The PIC showed a significant decreasing trend in the F I G U R E 2 Comparisons of ω values among vertebrates of different altitudes based on 13 mtDNA PCGs (*p < .05; **p < .01) dN/dS ratio with elevation reduction (R 2 = 0.057, p = .014; Figure 3).
We identified the evolutionary rate of the 13 mtDNA PCGs were faster in the high-altitude vertebrates when compared with lowaltitude vertebrates.

| D ISCUSS I ON
This research focused on the adaptation of vertebrate mtDNA PCGs to high-altitude environments. In high-altitude environments, vertebrates must maintain normal energy production under hypoxic pressure (Qiu et al., 2012;Wang et al., 2011). Therefore, high-altitude vertebrates require more energy than low-altitude vertebrates.
Mitochondria (as the energy metabolism center of cells) provide energy for the animals through oxidative phosphorylation (Maes et al., 2004). Therefore, vertebrate mitochondria play important roles in high-altitude adaptation. Previous research has mainly focused on one kind or family of vertebrates, such as Tibetans (Chen et al., 2020), Tibetan pigs (Ma et al., 2019), Antilopinae (Hassanin et al., 2009), Galliform birds (Zhou et al., 2014), and Schizothoracine fishes (Li et al., 2013). However, previous research has not focused on different altitudes when assessing the five vertebrate taxa. Therefore, we studied the evolutionary rate of 13 mtDNA PCGs within vertebrates at different altitudes. Our research revealed a consistent role of vertebrate mitochondria in high-altitude adaptation.

| Ratio of nonsynonymous/synonymous nucleotide substitutions
The ratio of nonsynonymous/synonymous nucleotide substitutions (i.e., dN/dS ratio or ω value) has been widely used to measure the selection pressure intensity of the protein-coding genes (Andrieux & Arenales, 2014;Xia et al., 2019). In our research, the ω values The ω values of high-altitude vertebrates (HV) The Background ω 0.05720 1.00000 0.05720 1.00000 Foreground ω 0.05720 1.00000 20.80002 20.80002 Note: *represents the significant level; PPs of Bayes Empirical Bayes (BEB) analysis with p > .95 was regarded as candidates for selection (*> 0.95, **> 0.99) 2020) were under purifying selection during their evolution. This relates to the highly conserved mitochondrial DNA protein-coding genes (Wolstenholme, 1992). Purifying selection can delete deleterious substitutions to maintain the normal function of the proteincoding genes.
The ω values of 13 mtDNA PCGs and each mtDNA PCG were generally high among vertebrates in high altitudes. This result was consistent with previous research (Wang et al., 2016;Xu et al., 2007) and indicates that high-altitude environments affect the evolutionary rate of vertebrate mitochondrial DNA proteincoding genes.

| Selection pressure comparison
We used the Wilcoxon test, a branch model, and a branch-site model to detect any rapidly evolving genes and positive selection genes shared among the five high-altitude taxa. We used the Wilcoxon test and one-way ANOVA analysis (ω values of 13 mtDNA PCGs and each mtDNA PCG in the terminal branches) to detect the evolutionary rate variations between the vertebrates at different altitudes and the evolutionary rate differences between different taxonomies at the same altitude (i.e., the evolutionary rate differences among mammal, aves, reptile, amphibian, and actinopterygii at high or low altitudes), respectively. The ω values (13 mtPCGs and each PCG) showed no significant difference between the five taxonomies at high-or low-altitude environments. This may be due to the similar selection pressures (oxygen levels).
The ω values of the COX2 gene, COX3 gene, ND2 gene, and Cytb gene suggested they were all rapidly evolving in the highaltitude vertebrates. These results show the high-altitude environment (hypoxia) affects the evolutionary rate of the mitochondria in vertebrates, and the partial gene evolutionary rate in HV was faster than LV (Cheviron & Brumfield, 2012;Storz & Scott, 2019;Storz et al., 2010).
The rapidly evolving genes of the ND1 gene, ND2 gene, ND4 gene, We detected only one positive selection gene (ND5 gene) using the branch-site model. The ND5 gene is involved in the regulation of electron transfer in the respiratory chain, which provides necessary energy for vertebrate activity (Brandt, 2006;Javadov et al., 2021;Sousa et al., 2018). We suggest that the ND5 gene is responsible for high-altitude adaptation in vertebrate. These results suggested that high-altitude vertebrates probably employ the same genic toolkit to adapt to the extreme environment.

| CON CLUS ION
In summary, we investigated the evolution of vertebrate mtDNA PCGs at different altitudes. The ω value for the 13 mtDNA PCGs was higher in high-altitude vertebrates than in the low-altitude vertebrates as analyzed by the Wilcoxon test, branch model, branchsite model, and PIC analysis. Although mitochondrial evolution

CO N FLI C T O F I NTE R E S T
No potential conflict of interest was reported by the authors.

DATA AVA I L A B I L I T Y S TAT E M E N T
All the mitochondria genome sequences used in this study were accessed through the GenBank database using the accession numbers and DOI accession numbers in Appendix S1.