A list of human energy metabolism-related genes was retrieved from the KEGG (Kyoto Encyclopedia of Genes and Genomes) database (Kanehisa and Goto, 2000). The following pathways were retrieved: Glycolysis/Gluconeogenesis (EMP/GNG): ko00010; Citrate cycle: ko00020 (TCA); Pyruvate metabolism (PM): ko00620; and Oxidative phosphorylation (OXPHOS): ko00190 (Figure 1). The resulting gene list was used to query Ensembl BioMart and retrieve non-human Ensembl gene identificators (Kinsella et al., 2011). The gene IDs were used to search annotated coding sequences (CDS) of mammals in the NCBI¹ and Ensembl databases² (Supplementary Table S1). Additionally, BLASTn and tBLASTn searches (Johnson et al., 2008) were performed to identify one-to-one orthologs in the genomes of 12 hypoxia-tolerant species. This included 9 marine mammals, cetaceans: bottlenose dolphin (Tursiops truncatus), baiji (Lipotes vexillifer), killer whale (Orcinus orca), sperm whale (Physeter macrocephalus), bowhead whale (Balaena mysticetus), minke whale (Balaenoptera acutorostrata); pinnipeds: Weddell seal (Leptonychotes weddellii), Pacific walrus (Odobenus rosmarus divergens), sirenians: West Indian manatee (Trichechus manatus latirostris); 2 highland mammals: Tibetan yak (Bos mutus), Tibetan antelope (Pantholops hodgsonii); and 1 subterranean mammals: naked mole rat (Heterocephalus glaber) (Supplementary Figure S1). For genes with multiple splice variants, the longest transcript was retained. The resulting dataset, derived from 54 mammalian genomes, included 181 nuclear-encoded energy metabolism-related genes (Supplementary Table S1). In addition, nucleotide sequences of 13 mitochondrial protein-coding genes from 101 species were downloaded from NCBI and the MitoZoa database³ (Supplementary Figure S2). A total of 194 energy metabolism-related (nuclear and mitochondrial) genes were obtained. The EMP/GNG, PM, TCA, and OXPHOS metabolic pathways were composed of 49, 34, 29, and 111 genes, respectively (Supplementary Figure S3). Genes shared by different pathways are depicted in Supplementary Figure S3.

Coding (amino acid) sequence for each gene was aligned using MUSCLE v3.2 (Edgar, 2004) with default settings, manually corrected, and used to align nucleotide sequences in MEGA 6 (Tamura et al., 2013). Because the quality of sequence alignments dramatically affects the estimation of evolutionary parameters, we deleted all gaps and undetermined bases ‘N’ in alignments to reduce the rate of false-positive prediction, and only high-quality sequences were used in downstream analyses. The species used for each gene, and associated accession numbers, are presented in Supplementary Table S1.

The impact of natural selection on energy metabolism related genes was determined by estimating the ratio of non-synonymous (d_N) and synonymous (d_S) mutations (ω = d_N/d_S) implemented in CodeML (implemented in PAML v4.8) (Yang, 2007). Briefly, ω > 1, ω = 1, and ω < 1 indicate positive selection, neutral, and negative selection, respectively. To determine whether codon positions in hypoxia-tolerant branches were under positive selection, we employed the PAML branch-site model, which assesses selective pressure (ω) at sites in foreground lineages of interest, specified a priori, and remaining background lineages (Zhang et al., 2005). In separate tests, branches leading to cetaceans, pinnipeds, sirenians, highland (high-altitude) species, or subterranean (burrowing) mammals were specified as foreground branches, with the remaining of the tree as the background branch (Supplementary Figures S1, S2). Additionally, we also analyzed the combined branches, all hypoxia-tolerant mammals, independently. Compared with the corresponding null hypothesis of neutral evolution, foreground branches in the modified model A (test 2) that have a class of sites with the ratio ω₂ > 1 are candidates for positive selection. Likelihood ratio tests (LRTs) with a chi-squared test (df = 2) were conducted on nested models. Such sites under positive selection were detected by the Bayes empirical Bayes (BEB) method, with posterior probabilities of ≥0.80 (Yang et al., 2005). Resulting P-values were corrected by employing a false discovery rate (FDR) cutoff of 0.05 (Storey and Tibshirani, 2003). A well-accepted species tree among mammals (Supplementary Figures S1, S2) was used as the input tree in all analyses (Murphy et al., 2001).

In order to identify parallel/convergent amino acid changes in hypoxia-tolerant mammals, the sequences of the ancestral nodes were reconstructed with the help of the CodeML program implemented in PAML (Yang, 2007). Convergent/parallel amino acid replacements, from the last common ancestor-state to present-species along each hypoxia-tolerant branch, were detected by an in-house Python script (Supplementary Data Sheet 1). Amino acid residue changes were classified into convergent or parallel changes. Convergent changes: any amino acid at the same position changes to the same amino acid in two or more distinct lineages (for example, Ala10Glu in cetaceans; His10Glu in pinnipeds). Parallel changes: an identical amino acid at the same position changes to the same amino acid in two or more distinct lineages (for example, Ala10Glu in cetaceans; Ala10Glu in pinnipeds). To test whether observed parallel/convergent substitutions in focal hypoxia-tolerant branches were fixed by random chance or by natural selection, the statistical significance of amino acid changes was tested by the method developed by Zhang and Kumar (). To carry out in-depth investigation of convergence between hypoxia-tolerant species, we measured the fit (site-wise log-likelihood support; SSLS) (Parker et al., 2013) of each amino acid along the orthologous CDS alignment to the generally accepted species tree (H₀) and a convergent tree (H_A) that constrained hypoxia-tolerant species into monophyletic clades. The difference in SSLS for a single site under different tree topologies was calculated as follows: ΔSSLS = SSLS (H_A) - SSLS (H₀) (Supplementary Figure S4), where a negative ΔSSLS implies support for convergence. SSLS for each gene alignment was calculated by RAxML v8.0.26 (Stamatakis, 2014).

The expression of the 194 energy metabolism-related genes were examined in de novo assembled transcriptomes (Seim et al., 2014). We contrasted cetacean (bowhead whale and minke whale) and terrestrial (naked mole rat, domestic yak, Brandt’s bat, Chinese tree shrew, rhesus monkey, rat, and mouse) liver transcriptomes. Briefly, library size-normalized read counts were subjected to the voom function (variance modeling at the observation-level) function in the R package ‘limma’ v3.18.13 (Law et al., 2014) with trend = TRUE for the eBayes function and correction for multiple testing (Benjamini–Hochberg false discovery P ≤ 0.05). Following limma analysis, genes with a log₂ 1.2-fold-change differences between cetaceans (bowhead whale, and minke whale) and other mammals were considered differentially expressed. RNA-seq gives very accurate measurements of relative expression across a broad dynamic range (Wang et al., 2008). However, in cross-species analyses (here, two cetaceans versus seven terrestrial mammals), one would expect more heterogeneity than when comparing within a species (where typically a log₂ 2.0-fold-change cutoff is employed). Thus, a gene could be modestly upregulated, but still be biologically relevant.

To estimate the selective pressure acting on genes in energy metabolism pathways in mammals, we ran separate branch-site tests on the ancestral and terminal branches leading to cetaceans (bottlenose dolphin, baiji, killer whale, sperm whale, minke whale, and bowhead whale), pinnipeds (Weddell seal and Pacific walrus), sirenian (West Indian manatee), highland species (Tibetan yak, and Tibetan antelope), and a subterranean species (naked mole rat) (Figure 2 and Supplementary Table S2). In cetaceans, eight genes showed evidence of positive selection (positively selected genes; PSGs) (test 2, Supplementary Table S2) along the lineages leading to bottlenose dolphin (Figure 2: branch i, ALDOA, CS, and PC), killer whale (Figure 2: branch j, NDUFV3), sperm whale (Figure 2: branch l, LHPP, and NDUFA9), minke whale (Figure 2: branch m, LHPP), bowhead whale (Figure 2: branch n, ATP6V0A4, ENO2, and PC), and the last common ancestor of baleen whales (mysticetes; Figure 2: branch e, ENO2) (Table 1) after correcting for multiple testing. For the remaining aquatic mammals, the LRT test showed evidence for positive selection on the genes PC and OGDH, respectively, in the lineages leading to Weddell seal (Figure 2: branch t) and West Indian manatee (Figure 2: branch x). When considering highland species, one gene (NDUFA10) was under positive selection in the Tibetan antelope (Figure 2: branch q), whereas in Tibetan yak four genes (PGM1, PKLR, ATP5F1, and NDUFA10) were under positive selection (Figure 2: branch o). When we set the subterranean (naked mole rat) as the foreground branch (Figure 2: branch v), positive selection was detected in two mitochondrial genes: ATP8 and ND5. Notably, the PSGs identified in hypoxia-tolerant species were not found in their sister taxa, although some genes in the sister taxa were also positively selected (Supplementary Table S2). These finding suggests that the PSGs identified are robust candidate hypoxia adaptation genes.

To assess the potential function of PSGs in hypoxia-tolerant mammals, we categorized them at the pathway level. In cetaceans, PSGs were found in all four energy metabolism pathways of interest: 6.9% (2/29 genes) in TCA, 4.1% (2/49 genes) in EMP/GNG, 3.6% (4/111 genes) in OXPHOS, and 2.9% (1/34 genes) in PM (Figure 3 and Table 1). Notably, PSGs were enriched in TCA signaling pathway in the cetacean lineages, with approximately twice as many PSGs as the PM pathway (Figure 3 and Table 1). Similarly, PSGs were also enriched in TCA signaling pathway in the pinniped and sirenian branches (Figure 3). In contrast, PSGs in the Tibetan antelope (0.9%, 1/111 genes) and naked mole rat (1.8%, 2/111) were exclusively found in the OXPHOS pathway. Tibetan yak showed positive selection in EMP/GNG (4.1%, 2/49 genes), PM (2.9%, 1/34 genes) pathways, and OXPHOS (1.8%, 2/111 genes) pathways.

To gather additional evidence on the adaptive evolution of energy metabolism related genes in cetaceans, we reconstructed ancestral sequences and mapped species-specific amino acid changes along cetacean lineages. Compared with other mammals, a total of 55 unique amino acid substitutions in 44 genes were observed in all cetaceans; three of which (ATP6V0A4, ENO2, and PC) were inferred to have evolved under positive selection (Supplementary Table S3). Given that differential gene expression could conceivably also contribute to hypoxia-tolerance (Böhne et al., 2013; Fang et al., 2014), we contrasted liver (a key metabolic organ) transcriptomes from two cetaceans (bowhead whale, and minke whale) and seven terrestrial mammals, a liberal cutoff (log₂ fold-change ≥ 1.2) revealed no significant differential expression of the 194 energy metabolism genes interrogated in this study (Supplementary Figure S5 and Supplementary Table S4).

To identify convergent evolution of energy metabolism related genes in hypoxia-tolerant mammals from different habitats, we reconstructed ancestral gene sequences for the internal nodes of the species tree and identified shared amino acid substitutions along lineages leading to hypoxia-tolerant species (Figure 4). In total, 7 parallel amino acid residue substitutions in 7 genes were identified in the marine group (branch a, g, x), as well as 2 parallel substitutions in 2 genes in high-altitude species (branch p, q). We identified 26, 3 and 2 parallel substitutions in 23, 3, and 1 genes in marine + subterranean mammals (branch v), marine + highland mammals, and highland + subterranean mammals, respectively (Table 2). Taken together, 38 parallel non-synonymous amino acid substitutions in 29 genes occurred along hypoxia-tolerant species. The number of substitutions was not random (P ≤ 0.05; data not shown). Three of these genes (NDUFA9, NDUFA10, and NDUFV3) were inferred to have evolved under positive selection in the hypoxia-tolerant lineages. Remarkably, parallel non-synonymous amino acid substitutions in 5 genes (CYC1, LDHD, NDUFA10, NDUFC2, and PGAM2) occurred along more than two hypoxia-tolerant species lineages (Table 2). Moreover, all 38 substitutions also ranked highly in the convergence signal distribution, and 92.1% (35/38) loci in the top 5% supporting convergent tree by ΔSSLS test (Table 2). Additionally, 76.3% (29/38) of substitutions were radical amino acid changes.

We found evidence of positive selection acting on energy metabolism related genes in diverse hypoxia-tolerant species. PSGs related to the OXPHOS pathway mainly occurred along hypoxia-tolerant lineages with a terrestrial lifestyle (such as Tibetan yak, Tibetan antelope, and naked mole rat), whereas PSGs associated with the TCA pathway appeared to have a stronger signal of positive selection in marine mammals. Species-specific PSGs were examined in the hypoxia-tolerant mammals, providing useful pointers to the evolution of unique traits. For example, ALDOA, ENO2, CS, ATP6V0A4, LHPP, NDUFA9, and NDUFV3 PSGs were unique to cetaceans.

Another PSG (PC) found in two cetaceans and the Weddell seal encodes pyruvate carboxylase, an enzyme which converts pyruvate to oxalacetic acid - the first control step of gluconeogenesis (Söling and Willms, 2013). Adaptive evolution of this gene suggests that a hypoxia adaptation of marine mammals involves the Cori cycle and improved lactic acid metabolism during recovery from dives (Cori, 1981; Castellini, 2012). Similarly, PKLR, which encodes pyruvate kinase, a rate-limiting enzyme in glycolysis, was unique to Tibetan yak. This finding is consistent with a previous study on this species (Ge et al., 2013), suggesting that glycolysis is essential for survival in high-altitude hypoxic environments. Two mitochondrial genes (ND5 and ATP8) in the OXPHOS pathway were exclusively identified to be under positive selection in the naked mole rat, suggesting that mitochondria play a role in subterranean hypoxia-tolerance. Various studies have investigated hypoxia-tolerance in the naked mole rat (Fang et al., 2014; Park et al., 2017), but its mitochondrial energetics remains largely unexplored, which are likely critical to its hypoxia-tolerance and extraordinary lifespan (Stoll et al., 2016).

The above observations are indicative of species- or lineage-specific adaptations to hypoxia. There are a several, non-mutually exclusive explanations for this deduction. It is well established that physiological features supporting a hypoxic lifestyle have evolved independently. Marine mammals routinely experience acute hypoxia, which causes them to depend on a finite supply of intrinsic oxygen stores (Ramirez et al., 2007). A dive is manifested by progressive asphyxia – concomitant with increased hypoxia, hypercapnia, and acidosis (Elsner and Gooden, 1983) – yet marine mammals exhibit an extraordinary ability to circumvent the deleterious effects of hypoxia. These include apnea (cease of breath), bradycardia, and subsequent decreased cardiac output, as well as peripheral vasoconstriction (selective shutdown of blood flow in non-essential, peripheral tissues) (Panneton, 2013). In contrast to marine mammals, terrestrial animals that live at high-altitudes or in burrows cope with a different form of hypoxic stress, chronic hypoxia due to a lower partial oxygen pressure (pO₂), by improving O₂ affinity and increasing ventilation and cardiac output (Ramirez et al., 2007). For example, the genomes of the Tibetan antelope and Tibetan yak show independent signals of adaptive evolution and gene-family expansion. Genes associated with oxygen transmission and mitochondrial membranes expanded in Tibetan antelope (Ge et al., 2013), whereas genes related to sensory perception and olfactory receptor activity expanded in Tibetan yak (Qiu et al., 2012), respectively. A series of genome scans of high-altitude mammals have identified genes under positive selection: EGLN1 and EPAS1 in Tibetan humans (Simonson et al., 2010; Yi et al., 2010); and ADAM17, ARG2, and MMP3 in Tibetan yak (Qiu et al., 2012). Moreover, candidate high-altitude gene loci in snub-nosed monkey populations are manifested by habitat-specific molecular evolution (Yu et al., 2016). Taken together, ecological and/or biogeographic differences, as well as distinct physiological traits, may have imposed divergent selective pressures on the evolution of hypoxia-tolerance.

Cetaceans, one of three extant marine mammal orders descended from terrestrial ancestors that recolonized the sea, offer an ideal case study. During the habitat transition of cetacean ancestors, dramatic changes in anatomy and physiology occurred, e.g., streamlined body, paddle-shaped fore-flippers, degenerated hindlimbs, and lack of hair (Berta et al., 2006; Muizon, 2009). Cetaceans are also further challenged by acute episodes of hypoxia during diving. Cetacean adaptations to hypoxia are manifested as increased oxygen storage in blood and muscle, a reduced metabolic rate, and selective vasoconstriction (Costa, 2007; Ramirez et al., 2007). Molecular mechanisms underlying these innovations have only begun to be explored (Mirceta et al., 2013; McGowen et al., 2014; Tian et al., 2016).

In this study, we also identified genes under positive selection in cetaceans and investigated their functional enrichment, at the energy metabolism pathway level. We found a signature of positive selection across cetaceans in a total of eight genes, and in all four energy metabolism pathways (EMP/GNG: 2; TCA: 2; PM: 1; and OXPHOS: 4). Interestingly, the highest proportion of PSGs was in the TCA cycle. The TCA cycle is a central part of aerobic metabolism and a chief source of electrons for ATP generation from the chemical breakdown of carbohydrates, fats, and proteins (Weitzman, 1987). Of note, this includes the citrate synthase (CS) gene, which encodes a key rate-limiting enzyme in energy production via the TCA cycle. We speculate that these changes mediate an enhanced capability for aerobic metabolism by cetaceans during dives (Castellini, 2012). Of particular interest, six genes (LDHA, LDHD, PC, PCK1, FBP1, and GPI) with cetacean-specific amino acid changes are rate-limiting enzymes in the gluconeogenesis pathway (GNG) (Cori, 1981) (Figure 1 and Supplementary Figure S6). Changes in these genes may result in improved conversion of lactate to glucose or glycogen; clearing lactate from the blood and restoring circulating levels of glucose. Two genes (ALDOA and ENO2) associated with anaerobic respiration adaptively evolved in cetaceans, further supporting the role of glycolysis (EMP/GNG pathway) as a critical energy supply during prolonged, deep dives (Butler and Jones, 1997).

It is appreciated, however, that stratifying hypoxia-resistance and other unique adaptations at the genetic level is not straightforward. Cetaceans have large, energy-demanding brains and previous genome-wide studies of mammals revealed that genes associated with energy metabolism are targets of natural selection in cetaceans (McGowen et al., 2012; Shen et al., 2012; Sun et al., 2013). It then follows that selection acting on energy metabolism-related genes may be in response to the evolution of hypoxia adaptation, large brains, as well as aquatic locomotion (reflecting the historic shift from running to swimming). Similarly, flying insects and mammals (bats) also require significant energy for locomotion. Indeed, studies on bats (Shen et al., 2010, 2012) revealed a critical role for genes in the OXPHOS pathway in flight. Previous molecular studies on primates also showed evidence of adaptive evolution in energy metabolism-related genes, suggesting that these genes played a major role in overcoming elevated energy demands during the expansion of the neocortex (Grossman et al., 2001). Taken together, positive selection of energy metabolism-related genes in non-hypoxia-tolerant species is expected, and not mutually exclusive to a role in hypoxia-tolerance, given that such genes had multiple advantageous phenotypic effects (pleiotropy) during mammalian evolution.

Conceptually, hypoxia-tolerance provides a classical example of convergent evolution. Even though there are clearly environmental differences, inherent similarities may drive recurrent evolution at the sequence level, manifested as convergent or parallel amino acid changes. Recent work, relying on whole-genome sequences of marine mammals (Foote et al., 2015), revealed that convergent amino acid substitutions in genes associated with phenotypic adaptations can be identified. It is appreciated that hypoxia-tolerance is a very complex phenomenon that involves a suite of genes. Although the morphologies and physiological characteristics of hypoxia-tolerance in marine and terrestrial mammals differ greatly, we reveal sequence convergence in energy metabolism-related genes. These genes could reflect core hypoxia-tolerance associated functions. Interestingly, one PSG (PFKP) was detected in the ‘combined’ hypoxia-tolerant mammal branches (Supplementary Table S1). Five genes (NDUFA9, NDUFA10, NDUFAB1, NDUFC2, and NDUFV3) encodes the nuclear-encoded subunits of complex I in the OXPHOS pathway (Figure 1). Parallel/convergent evolution of these gene may reflect the importance of oxidative phosphorylation for survival in hypoxic environments. Hypoxia also induces oxidative stress, manifested by lipid peroxidation, protein oxidation, and DNA damage (Ramirez et al., 2007). ALDH3A1 (encodes cytosolic aldehyde dehydrogenase 3A1; Figure 1) protects against oxidative stress in various ways, including enhanced DNA and cell repair (Estey et al., 2007). Another gene, PFKP (phosphofructokinase platelet type; Figure 1), catalyzes the glycolytic pathway by converting fructose-6-phosphate to fructose-1, 6-bisphosphate (Webb et al., 2015). Mutations inhibiting phosphofructokinase activity cause glycogen storage disease (Webb et al., 2015). Taken together, we reveal several convergent genes in hypoxia-tolerant species.