Evolutionary impacts of purine metabolism genes on mammalian oxidative stress adaptation

INTRODUCTION

There are more than 6 000 species of extant mammals (Burgin et al., 2018), a taxonomic class that manifests remarkable biological diversity. Habitat conditions range from chronic hypoxia in tissues of highland and subterranean animals to repeated cycles of normoxia and hypoxia in diving mammals (i.e., cetaceans, pinnipeds, and sirenians), to flight and hibernation in bats (Hermes-Lima & Zenteno-Savı́n, 2002; Ramirez et al., 2007). These are all conditions of oxidative stress. Marine mammals are fascinating due to their daily routine of breath-hold diving, which is characterized by reduced blood flow (ischemia) and decreased tissue oxygen content (hypoxia) until resurfacing. Oxidative stress is an umbrella term referring to an imbalance between reactive oxygen species (ROS) production and their elimination. Major cellular components – lipids (peroxidation of unsaturated fatty acids in membranes), proteins (denaturation), carbohydrates, and nucleic acids – are susceptible to damage by ROS (Blokhina et al., 2003).

Numerous studies have revealed that mammals in extreme environments exhibit obvious oxidative stress adaptations that manifest as changes in physiological features, oxidants, antioxidant enzyme activity, genetic damage, and even disease status (Ramirez et al., 2007; Yu et al., 2020). For example, diving mammals (Wilhelm Filho et al., 2002), hibernating bats (Yin et al., 2016), and underground rodents (Yu et al., 2020) exhibit enhanced glutathione and higher activity and expression of antioxidant enzymes, including superoxide dismutase, catalase, and glutathione reductase. Furthermore, genome-wide studies found that genes related to oxygen storage and delivery, and energy metabolism are subjected to strong and consistent selective forces in species under oxidative stress (including highland species, diving mammals, and flying bats) (Nery et al., 2013; Qiu et al., 2012; Shen et al., 2010; Storz & Moriyama, 2008; Tian et al., 2016, 2017). Despite these advances, our understanding of genetic changes for adaptations to oxidative stress remains incomplete.

Purines are one of the most abundant metabolic substrates. In addition to serving as building blocks for DNA and RNA, purines are necessary for cell survival and proliferation (Pedley & Benkovic, 2017). Two pathways regulate cellular purine levels: purine salvage and de novo biosynthesis (Pedley & Benkovic, 2017; Yin et al., 2018). Most cellular purines stem from purine salvage, catalyzed by hypoxanthine-guanine phosphoribosyltransferase (HGPRT) and adenine phosphoribosyltransferase (APRT). HGPRT recycles hypoxanthine (HX) and guanine to inosine monophosphate (inosine 5'-monophosphate; IMP) and guanine monophosphate (GMP), respectively. APRT catalyzes adenine and 5'-phosphoribosyl 1-pyrophosphate (PRPP) to generate adenosine monophosphate (AMP). The de novo biosynthetic pathway is energy-intensive and utilizes numerous amino acid substrates and cofactors (e.g., ATP, glutamine, and aspartate) to generate IMP from PRPP in a nine-enzyme complex (Pareek et al., 2020). IMP is thus the main product of both purine pathways and essential for the interconversion to various intermediates (e.g., AMP, GMP, adenosine, guanosine, and inosine). Purine nucleoside phosphorylase (PNP) converts inosine and guanosine to the purine base HX and guanine. HX and guanine are catalyzed in separate reactions to form xanthine by xanthine oxidase (XO) and guanine deaminase. Notably, XO uses molecular oxygen as an electron acceptor and generates superoxide anion and other reactive oxygen products. XO can further oxidize xanthine to form a final product, uric acid, in hominoids, birds, and reptiles (Maiuolo et al., 2016). In other species (some teleost fish and amphibians), uric acid can be further degraded to allantoin by uricase, to allantoic acid by allantoinase, to urea by allantoicase, and to ammonia by urease (Hayashi et al., 2000). These enzymatic processes represent the terminal steps of purine degradation. It should be noted that uric acid is an effective antioxidant. Increased serum uric acid levels could protect against oxidative stress-provoked aging and cancer (Glantzounis et al., 2005).

All cells require a balanced purine metabolism. Dysregulated purine metabolism is related to oxidative stress that results in numerous human diseases and disorders. For example, HGPRT deficiency results in hyperuricemia (i.e., high uric acid levels in the blood), triggering a series of clinical consequences, such as gout and renal failure (Nyhan, 2005). Mutations or deficiency of adenosine deaminase (ADA), a key enzyme for purine degradation or salvage into the purine pool, increases susceptibility to infections and autoimmunity (Yin et al., 2018). Deleting erythrocyte equilibrative nucleoside transporter 1 (eENT1), an enzyme eliminating extracellular adenosine, allows rapid accumulation of plasma adenosine to counteract hypoxic tissue damage in mice (Song et al., 2017). Many studies have focused on the adenosine-signaling network in stress or hypoxic tolerance. Adenosine is a potent indicator of tissue and cell damage from hypoxia. Elevated adenosine signaling coupled with 5' AMP-activated protein kinase (AMPK) and adenosine receptors (e.g., ADORA2B) contributes to human adaptation to high-altitude hypoxia by inducing 2,3-BPG production and triggering O₂ release (Liu et al., 2016; Sun et al., 2017). Moreover, Nemkov et al. (2018) revealed that hypoxia decreases purine oxidation and enhances purine salvage reactions in human and mouse red blood cells (Nemkov et al., 2018). In addition, under conditions of reduced oxygen availability (i.e., ischemia and hypoxia), hydrolyzation of depleted ATP generates the XO substrate HX, providing free oxygen radicals for ROS generation during reperfusion (Robert & Robert, 2014). Apparent enhanced purine metabolism by marine mammals has been reported. The bottlenose dolphin (Tursiops truncates) shows enhanced purine salvage compared to humans (López-Cruz et al., 2016). Similarly, ringed seal (Pusa hispida) heart and kidney tissues can maintain HX levels under ischemic conditions (Elsner et al., 1998). However, whether and how molecular changes in purine metabolism genes enable mammals to resist oxidative damage and purine disorders throughout their relatively long lifespans remains an open question.

In this study, we investigated the molecular adaptations of genes involved in purine metabolism pathways in mammals. We performed codon-based selection analyses to identify if positive selection took place in oxidative stress-tolerant species. We provide evidence that convergent molecular adaptations in this pathway underlie the evolution of oxidative stress tolerance across different mammal lineages. Using a cell-based functional assay, we further investigated the consequences of unique amino acid changes in cetacean XDH (xanthine dehydrogenase), a key gene in the purine metabolism pathway, which can be converted to XO by oxidation of sulfhydryl groups or limited proteolysis. Our study contributes to a better understanding of molecular and genetic adaptations associated with oxidative stress defenses of mammals.

MATERIALS AND METHODS

Taxon coverage, dataset preparation, and sequence alignments

Our data set encompassed seven mammal lineages that independently evolved oxidative stress tolerance: cetaceans (n=15), pinnipeds (n=3), sirenians (n=1; manatee), highland species (n=2), subterranean species (n=1), flying bats (n=4). Additional mammals (n=14) were included as control (Supplementary Table S1). The 117 genes (Supplementary Table S2) involved in purine metabolism were obtained from the Kyoto Encyclopedia of Genes and Genomes ( KEGG ) database (map00230) (Kanehisa et al., 2007). KEGG gene IDs were converted to gene symbol using DAVID v.6.8 (Dennis et al., 2003). For species with a publicly available genome and gene set in Ensembl, the coding sequences (CDS) of one-to-one orthologous coding genes were downloaded from Ensembl BioMart (Kinsella et al., 2011). Otherwise, the candidate CDSs were obtained using BLASTn v2.4.0+ (Camacho et al., 2009) (E-value cutoff of 10^-5) and human genes as queries. Coding sequences were aligned with MEGA v7.0.62 (Kumar et al., 2016), followed by removal of gaps and non-homologous regions using Gblocks v.0.91b (Talavera & Castresana, 2007).

Molecular evolutionary analyses

Multiple sequence alignments were analyzed, with the species tree as the guide tree, using “codeml” program of PAML v.4.7 (Yang, 2007). For the all-mammal dataset (40 mammals), we ran the branch-site model along cetaceans and other lineages with oxidative stress resistance (branches A–G in Figure 1). We also repeated the tests on closely-related sister taxa with no apparent oxidative stress resistance (e.g., the cattle, Bos taurus as a sister species to cetaceans) (branches H–N in Figure 1). The modified model A (Ma) estimates separate site-wise variation in omega (ratio of nonsynonymous (d_N) to synonymous (d_S) rates) for the foreground branch (ω2>1) and the remaining background branches. This model was compared with a null model (Ma0), which is similar to Ma but allow the fixed ω2=1 in the foreground lineage.

Figure  1.  Phylogenetic tree of mammals

The topology of the tree is based on previous studies (McGowen et al., 2020; Scornavacca et al., 2019). Shaded boxes (A: pinnipeds, B: cetaceans, C: Tibetan yak, D: Tibetan antelope, E: bats, F: naked mole rat, and G: manatee) represent oxidative stress-tolerant taxa, while branch H–N represents sister taxa. External bars represent the diving abilities (maximum dive time, average dive time, maximum dive depth, average dive depth) of cetaceans obtained from published data (Supplementary Table S3).

DownLoad: Full-Size Img PowerPoint

Clade model C (CmC) (Weadick & Chang, 2012), which allows a class of sites to evolve differently between two or more partitions along a phylogeny, was used to test for a shift in selection pressures along particular branches. We partitioned the 15 cetacean species into deep/long diving (DLD) and shallow/short diving (SSD) groups according to their diving capacity (Figure 1; Supplementary Table S3). Analyses were conducted with the branch leading to all cetaceans, DLD and SSD, designated as the foreground. The null model M2a_rel does not allow for divergence in the foreground clade but allows ω to be unconstrained.

Nested models were compared using a likelihood ratio tests (LRT) with a χ² distribution. P-values <0.05 indicates a significantly better fit by the alternative model. P-values were adjusted using the Benjamini-Hochberg method (Benjamini & Hochberg, 1995), with a significance cut-off of 0.05. Amino acid sites under positive selection were detected by Bayesian empirical bayes (BEB) inference, with an 80% posterior probability cut-off. Overlap in positively selected genes was visualized using UpSet plots (Lex et al., 2014). The significance of observed intersections (with all genes tested as background) was calculated using the supertest function of SuperExactTest v.0.99.4 (Wang et al., 2015).

Gene–phenotype association analysis

The free-ratios models in PAML were used to evaluate the average d_N/d_S ratio from the ancestral cetaceans to each terminal species. The log₁₀-root-to-tip d_N/d_S were then set as species data and used in phylogenetic generalized least squares (PGLS) regression (Martins & Hansen, 1997) against log₁₀-transformed diving ability variables. A well-supported species tree from McGowen et al. (2020) and divergence times from TimeTree (Kumar et al., 2017) was employed. The diving ability variables, obtained from previously published data (Supplementary Table S3), were maximum dive time (min), average dive time (min), maximum dive depth (m), and average dive depth (m). We performed PGLS using R v3.2.5 and the R packages “ape” (Popescu et al., 2012), “nlme” (Pinheiro et al., 2013), and “phytools” (Revell, 2012).

Evaluation of gene loss and its evolution

We searched for genes that possess inactivating mutations in cetaceans by querying multiple sequence alignment of the 117 mammalian purine metabolism-related genes. We assessed stop codons, frameshifting insertions, and deletions using the human as the reference. This analysis revealed putative loss of allantoicase (ALLC). Loss of ALCC was validated by polymerase chain reactions (PCRs) on genomic DNA from skeletal muscle of dolphin, finless porpoise, pygmy sperm whale, and beluga whale. PCR primers listed in Supplementary Table S4. We employed a 2X EasyTaq PCR SuperMix (Promega, USA), purified amplicons using a QIAquick Gel Extraction Kit Protocol (QIAGEN, Germany), and sequenced amplicons in both direction by Sanger sequencing.

We used RELAX (Wertheim et al., 2015) implemented in Datamonkey (http://www.datamonkey.org) (Weaver et al., 2018), to determines if ALLC evolves under relaxed selection in the cetaceans relative to other mammals. RELAX compares a test branch (here: cetaceans) is under relaxed or intensified selective strength relative to reference branch/branches (here, all other lineages) by allowing variation of evolutionary rates both between sites and across branches (Wertheim et al., 2015). Relaxed selection strength on pseudogene branches can thus be assessed. We also applied “codeml” program of PAML to evaluate whether ALLC has evolved neutrally. Codeml model A assumes that all branches have a single ω value, while model B allow ω=1 in all branches. Purifying selection was detected by comparing model A and B. Model C was used to estimate ω2 for ALLC-pseudogenized cetaceans and ω1 for other mammals, separately. A comparison between A and C models indicates that selective pressure is relaxed. To investigate whether the selective pressure is completely removed from the branches with a pseudogenized ALLC, we tested model D which allows a fixed ω2=1 for pseudogenized species and a ω1 for all remaining branches. We also tested a model E, where ω is allowed to vary among branches.

To estimate how long ALCC has been evolving neutrally in cetaceans, we employed a previously described gene loss dating procedure using the formula (Sharma et al., 2018):

We estimated the K_a/K_s value (K) for the entire branch, where the K_s indicated the gene under selection. K_n was set to 1, suggesting the gene evolved neutrally. The “T” denotes time in million years (Ma) since cetaceans split from the last common ancestor (53~59 Ma). Using the lower (53 Ma) and upper (59 Ma) bound of the confidence interval for the species divergence time T obtained from TimeTree (Kumar et al., 2017) and the K_a/K_s value of mammals with a functional ALLC gene (K_s=0.187), we estimated a lower and upper bound for T_n (estimate of how long ALLC evolved neutrally) as T (K–K_s)/(1–K_s).

Detection of unique or convergent amino acid residues

To examine the convergent amino acids among species under oxidative stress, we reconstructed the ancestral amino acid sequences using “codeml” program of PAML (Yang, 2007). Parallel and convergent sites were identified using the JTT-f_gene amino acid substitution model (Zou & Zhang, 2015). We estimated the branch length and substitution rate of each site relative to the average rate of the entire protein to infer the expected number of convergent sites in each gene. A Poisson test was used to verify whether the observed number of convergent sites of each gene was significantly more than the expected number (P<0.05). We identified genes under convergence with the following criteria: (1) identical amino acid residues in any two of seven groups (branch A–G; Figure 1) leading to species under oxidative stress; (2) an amino acid change has occurred between the species and its most recent common ancestor; (3) the changes in amino acids differ from the homologous positions in the non-oxidative stress species and its ancestral pair. We removed gaps and ambiguous sites. Substitutions were considered “parallel” if the independent changes have occurred from the identical ancestral amino acid, while changes from different ancestral amino acids were defined as “convergent” (Zhang & Kumar, 1997). Genes subjected to positive selection and with possessed non-random convergent amino acid sites between oxidative stress-tolerant species were considered adaptively convergent genes. Specific amino acid changes in oxidative stress-tolerance species were identified from multiple sequence alignments using FasParser v2.10.0.0 (Sun, 2017).

Protein expression and functional assays

Wild-type XDH (XDH_cow) for cow was synthesized using with 5' (NheI) and 3' (BamHI) restriction sites for insertion into the pCDNA3.1(+) expression vector. The cow was chosen for this experiment as it was a close relative of cetaceans, and XDH_cow had high amino acid sequence identity with cetacean XDH except for six cetacean-unique substitutions (Figure 3A). The T241A, E332Q, S495P, M585L, T750A, K949R mutations in cow wild-type XDH were generated via QuickMutation Site-Directed Mutagenesis Kit (Beyotime Institute of Biotechnology, China). These mutants were verified by Sanger sequencing. Expression vectors containing wild-type and mutant XDH (XDH_cetacean) were transiently transfected into HEK293T cells and harvested after 48 h.

Figure  3.  Unique amino acid changes and functional assessment of cetacean XDH

A: Protein domain structure (obtained from UniProt) and overview of mammalian XDH. B: Thermal stability of wild-type (cow) and mutant-type XDH (cow sequence with six amino acid substitutions unique to cetaceans). The enzyme was incubated at 25–60 °C for 30 minutes. Top, immunoblotting of FLAG-tagged wild-type and mutant-type XDH. Bottom, densitometry of FLAG-tagged XDH immunoreactivity, normalized to β-tubulin. C: Analysis of XDH to XO conversion in response to trypsin. Top, proteins subjected to anti-XDH immunoblotting. β-actin served as protein loading control. Bottom, percentages of XO compared with total labelled XO plus XDH. D, E: Xanthine-O₂ activity (D) and xanthine-NAD activity (E) were assayed by adding various concentration of xanthine with/without NAD⁺, respectively. Mean±SD of three separate experiments. Differences between wild-type and mutant-type were performed by unpaired, two-tailed Student’s t-test. Data were analyzed using GraphPad Prism or Origin software. P≤0.05 were considered statistically significant (^***: P<0.001; ^**: P<0.01; ^*: P ≤ 0.05).

DownLoad: Full-Size Img PowerPoint

To determine thermal stability of XDH protein, cells were immediately lysed using RIPA lysis buffer (Beyotime Institute of Biotechnology, China) supplemented with protease inhibitor PMSF. Cellular protein lysates were incubated at different temperatures (25 °C, 37 °C, 49 °C, and 61 °C) for 30min. After centrifugation, soluble proteins were used in Western blot assay. The remaining protein contents after temperature exposures were calculated from the control (unheated protein, incubated in the ice-water mixture for 30min) taking the protein contents as 100%.

To investigate the conversion rate of XDH to XO in response to trypsin treatment, collected cells were resuspended in the extraction buffer (Solarbio, China) and disrupted by sonication on ice. The soluble supernatants were collected and treated with 0.05% w/w trypsin (KeyGen Biotech, China) from 5 to 60 min (i.e., 5, 10, 15, 40, 50, and 60 min). The XDH and XO protein expression after trypsin exposure were measured by Western blot assay. Protein band density from three biological replicates were acquired using Image J software and data were plotted using GraphPad Prism software.

Enzyme activity assays of xanthine oxidoreductase (XOR: XDH/XO) was performed as previously described (Saito & Nishino, 1989), with some modification. XDH shows a preference for NAD⁺ reduction at the FAD reaction, whereas XO exclusively uses dioxygen as its substrate, generating superoxide and H₂O₂. For the steady-state kinetics analysis of xanthine-O₂ activity, the reaction solution was 250 μL of potassium phosphate buffer (50 mmol/L, pH 7.8), containing 0.4 mmol/L EDTA, 6.25 μL XDH crude enzyme and various concentrations of xanthine (Sigma-Aldrich, USA) substrate ranged from 0.5 to 2.5 μmol/L. In steady-state kinetics of xanthine-NAD activity, the same reaction system was used with 2.5 μmol/L xanthine and various concentrations of β-NAD (Sigma-Aldrich, USA) substrate ranged from 0.1 to 20 μmol/L. Before reactions were started, each enzyme and substrates were incubated separately for 10 min at 25 °C. The absorbance change of at 295 nm (due to urate formation) was recorded for 1 min with UV-1800PC spectrophotometer (Mapada, China). Controls were run using bovine milk xanthine oxidase (Sigma-Aldrich, USA) and a reaction medium that contained no or heat-treated (99°C, 10 min) enzyme. The kinetics constant values were calculated by non-linear fitting (GraphPad Prism 8.0) according to the Michaelis-Menten equation.

RESULTS

Unique signature of selection in purine metabolism genes of oxidative stress-tolerant mammals

We obtained 117 genes related to purine metabolism in 40 mammalian species. To better understand the selection pressures acting on mammal lineages with improved tolerance to oxidative stress, we analyzed the annotated coding sequences of purine metabolism genes in cetaceans, pinnipeds, bats, highland species (Tibetan yak (Bos mutus) and Tibetan antelope (Pantholops hodgsonii)), sirenians (the manatee (Trichechus manatus latirostris)), and a subterranean species (the naked mole rat (Heterocephalus glaber)) using the PAML branch-site model. We identified 72 genes (43 after multiple testing correction) with signals of positive selection across branches (described in Figure 1) from environments with elevated oxidative stress: nine in pinnipeds (clade A), 32 in cetaceans (clade B), eight in Tibetan yak (clade C), six in Tibetan antelope (clade D), two in bats (clade E), 13 in naked mole rat (clade F), four in the manatee (clade G) (Table 1; Supplementary Tables S5, 6). This number of positively selected genes (PSGs) is comparable to that of control lineages (branches H–N) (Figure 1; Supplementary Tables S5, 6); however, there was little overlap (Table 1; Supplementary Table S7). For example, of the 32 genes with evidence of positive selection along the ancestral cetacean branch (Supplementary Table S5), 25 genes were significant after FDR correction (Supplementary Table S6). Six out of eight PSGs identified in the terrestrial counterpart (e.g., cattle) branch (Supplementary Table S5, branch I) were significant after FDR correction (Supplementary Table S6). There is minimal overlap among the identified PSGs between cetaceans and artiodactyls (PDE3B, ADCY5, and AK9). We next compared the PSGs of species under higher oxidative stress. The seven oxidative stress-tolerant lineages shared no PSGs, and there were no significant intersections. Eighteen PSGs were shared by at least two branches (Figure 2A; Supplementary Table S8). Among these, PDE8B, NME3, and ADCY5 are PSGs in three lineages (Figure 2A; Supplementary Table S8).

Table  1.  Conclusion of positively selected genes detected by PAML in lineages with or not with oxidative stress lifestyle

Focal branch name	PSGs in focal branches	Control branch name	PSGs in control branches	Overlapped PSG
Branch A	9	Branch H	2	1
Branch B	32	Branch I	8	3
Branch C	8	Branch J	4	1
Branch D	6	Branch K	8	1
Branch E	2	Branch L	12	0
Branch F	13	Branch M	6	1
Branch G	4	Branch N	7	0
Total	72		47	7
Oxidative stress-resistant branches are branch A (pinnipeds), B (cetaceans), C (Tibetan yak), D (Tibetan antelope), E (bats), F (naked mole rat), and G (sirenians; the manatee). H denotes sea otter and European polecat; I, ancestor of artiodactylas; J, cow; K, sheep; L, ancestor of Ferungulata; M, rat and mouse; N, African savanna elephant. Clades label in Figure 1.

 | Show Table

DownLoad: CSV

Figure  2.  Evolution of purine metabolism genes in oxidative stress tolerance mammals

Number and overlap of positively selected genes (A), parallel/convergent amino acid substitutions (B) and genes with unique amino acid substitutions (C) in each mammalian lineage with oxidative stress tolerance, visualized as an UpSet plot. (D) represents the overview of adaptive amino acid changes in ATIC exclusive to oxidative stress tolerant species.

DownLoad: Full-Size Img PowerPoint

Convergent/parallel substitutions are specific to oxidative stress-tolerant mammals

To investigate molecular convergence, convergent or parallel amino acid residue changes were identified in the lineages of species under higher oxidative stress. Ancestral protein-coding sequences of each node in the phylogenetic tree (Figure 1) were reconstructed. In total, we identified 32 genes with 54 non-random parallel or convergent amino acid changes in the high oxidative stress groups (Figure 2B; Supplementary Table S9). Of these genes, 26 genes (81.2%) were shared between any two of the seven with oxidative stress tolerance (cetaceans vs. bats: seven; cetaceans vs. pinnipeds: eight; cetaceans vs. highland: six; cetaceans vs. manatee: four; highland vs. pinnipeds: one). Of these, 18 parallel amino acid residue substitutions in six genes (18.7%) (ENPP3, GUCY2C, NPR2, NUDT5, PDE6A, and PDE7A) were shared among three of the seven oxidative stress-tolerant species (Figure 2B). Notably, 19 of 32 genes have undergone positive selection (Supplementary Table S9).

Unique amino acid residue substitutions in species under higher oxidative stress

We next identified unique amino acid residue substitutions in species from a high oxidative stress environment. We found 37 genes with 67 unique sites in pinnipeds, 53 genes with 108 unique sites in cetaceans, 61 genes with 152 unique sites in highland species, 42 genes with 65 unique sites in bats, 37 genes with 133 unique sites in naked mole rat, and 51 genes with 145 unique sites in the manatee (Figure 2C; Supplementary Table S10), respectively. Six lineage-specific amino acid replacements (sites 241, 332, 495, 585, 750, 949) were detected in XDH of cetaceans (Figure 3A). For instance, Site 241 of XDH is threonine (Thr) in cetaceans and alanine (Ala) in background species (Figure 3A). Notably, of these loci, seven sites were specific across two lineages although the substitutions were not shared (Supplementary Table S10) (bats (R/K) and naked mole rat (C): site 131 of AK2; bats (S) and manatee (I): site 112 of AK8; pinnipeds (S) and naked mole rat (N): site 253 of CANT1; naked mole rat (A) and manatee (R): site 611 of GUCY1B3, bats (A) and cetaceans (V): site 63 of PDE5A; Tibetan Antelope (G) and naked mole rat (K): site 488 of XDH; cetaceans (R) and naked mole rat (Q): site 949 of XDH). XDH site 949 is lysine/K in all mammals except for cetaceans (arginine/R) and the naked mole rat (glutamine/Q) (Figure 3A). In addition, we also identified non-synonymous amino acid changes at different sites in six genes (XDH, NTPCR, ADK, PDE7B, ADCY3, and ATIC) in at least five oxidative stress-tolerant groups (Figure 2C, D). These results indicate that purine metabolism genes might show evidence of convergent evolution in mammals under higher oxidative stress.

Functional effects of cetacean specific unique amino acid substitutions

To investigate whether unique amino acid substitutions are associated with functional changes, we experimentally assayed the XDH of cetaceans (six specific amino acid changes), representatives of oxidative stress-tolerant species (Figure 3A; Supplementary Table S11). We mutated cow XDH to obtain a mutant-type enzyme (XDH_cetacean). We then expressed wild-type (XDH_cow) and mutant-type (XDH_cetacean) XDH in vitro. Although wild-type and mutant-type XDH protein levels both decreased after heat treatment compared to the unheated control, the mutant-type protein was slightly more abundant at 25 °C (two-tailed Student’s t-test, P=0.001) and 37 °C (P=0.002) (Figure 3B), suggesting that it exhibits higher thermal stability at body temperature. XDH protein was exposed to trypsin to assess XDH/XO conversion. After 5 to 60 minutes, wild-type and mutant-type XDH levels decreased, whereas XO levels increased (Figure 3C). Mutant-type XDH generated significantly less XO at 40 min (two-tailed Student’s t-test, P=0.047), 50 min (P=0.012), and 60 min (P=0.001) (Figure 3C). A catalase activity assay showed that wild-type and mutant-type XDH display typical Michaelis-Menten kinetics with xanthine substrates or NAD⁺ as a coenzyme. Mutant-type enzyme incubated with various concentrations of xanthine showed a significant increase in xanthine-O₂ activity (Figure 3D). Similarly, significantly higher xanthine-NAD activity was also observed for the mutant-type at various concentrations of NAD⁺. Notably, the mutant-type Michaelis constant (K_m) value for xanthine-O₂ activity (K_m=2.572±0.154 µmol) was lower than that of the wild-type (K_m=3.521±0.326 µmol) (Figure 3D; Table 2). The K_m of NAD⁺ for the mutant enzyme (K_m=0.195±0.012 µmol) was also lower than that of wild-type (K_m=0.420±0.027 µmol) (Figure 3E; Table 2). The V_max values of the mutant-type xanthine-O₂ and xanthine-NAD activity were close to that of the wild-type (Table 2). These results reflect a higher affinity of the mutant-type enzyme for xanthine and NAD⁺.

Table  2.  Kinetic parameters of recombinant XDH

Substrate	XDHwild			XDHcetacean
	Km (mmol/L)	Vmax (μmol/min/mg/)		Km (mmol/L)	Vmax (μmol/min/mg/)
Xanthine	3.521±0.326	0.0130±0.0004		2.572±0.154	0.0269±0.0005
NAD+	0.420±0.027	0.0262±0.0004		0.195±0.012	0.0315±0.0004
All assays were performed at 25 °C in 50 mmol/L Potassium phosphate buffer, pH 7.8. Kinetic parameters were determined by monitoring the increase in uric acid at 295 nm.

 | Show Table

DownLoad: CSV

Loss of the ALLC in cetaceans

During our mammalian purine metabolism gene survey, we failed to obtain ALLC from the genomes of six cetaceans (Figure 4). An analysis of the ALLC sequences that could be recovered revealed that cetacean ALLC genes share point mutations that result in numerous premature stop codons (Figure 4). Four cetaceans (pygmy sperm whale (Kogia breviceps), humpback whale (Megaptera novaeangliae), gray whale (Eschrichtius robustus), and North Pacific right whale (Eubalaena japonica)) have additional one-bp and five-bp deletions, leading to frameshifts at the end of ALLC exon 6 (Figure 4). The shared nonsense mutation in the middle of ALLC exon 5 was validated by PCR and Sanger sequencing of genomic DNA from skeletal muscle of bottlenose dolphin, finless porpoise (Neophocaena phocaenoides), beluga whale (Delphinapterus leucas), and pygmy sperm whale (Figure 4).

Figure  4.  Inactivating mutations of cetacean ALLC

Schematic showing locations of inactivating mutations. Exons 1–11 are shown as yellow boxes (drawn to scale), introns are shown as black lines (not drawn to scale). The black dotted line surrounding exon 5 indicates that it was PCR-amplified and Sanger sequenced in four species: bottlenose dolphin (Tursiops truncates), finless porpoise (Neophocaena phocaenoides), beluga whale (Delphinapterus leucas), and pygmy sperm whale (Kogia breviceps). Right, overview of Sanger sequencing-validated premature stop codons.

DownLoad: Full-Size Img PowerPoint

Evolutionary analyses of ALLC showed that the average ω across the tree is 0.231 (model A), significantly lower than 1.0 (model B, Supplementary Table S12), indicative of strong purifying selection on ALLC in mammalian evolution. In addition, model D (Supplementary Table S12), which allows two different ω values for pseudogenized (ω2=0.926) and functional (ω1=0.187) lineages, fit the data significantly better than the null model A that assumes a single ω for all branches (Supplementary Table S12). The selective purifying pressure is markedly relaxed (increased ω) in cetacean branches with a pseudogenized ALLC. Furthermore, the ω of cetaceans (ω2=0.926) (model D vs. model C, P > 0.05; Supplementary Table S12), indicates that the selective pressure of cetacean ALLC has shifted toward neutrality (i.e., ω=1.0). Finally, model E, which allows ω variation in all branches, fit significantly better than model C, suggesting that ω varies significantly among non-cetacean lineages. Similarly, we also found that ALLC evolves under relaxed selection (K<1 at P<0.001) in all cetaceans using the RELAX model (Supplementary Table S13). For cetaceans, we estimate that loss of ALLC occurred around 48.2~53.6 Ma. That is, after the split of the cetacean and artiodactyl lineages 53~59 Ma (Kumar et al., 2017) (Supplementary Table S14).

Selective constraint shift and correlation analysis in purine metabolism-related gene of cetaceans with different diving ability

A previous study reported that oxidative stress levels are related to diving abilities, such as diving depth and durations (Cantú-Medellín et al., 2011). Considering that cetaceans are a diverse group of marine mammals with a range of diving capacities (Perrin et al., 2009), we collected the phenotype data of cetaceans (Supplementary Table S3) to test whether shifts in evolutionary constraints on purine metabolism-related gene exist between cetaceans with different diving capacities. Two independent rounds of Clade model C analyses showed significantly divergent selection in 26 genes between cetacean and other mammals partitions, as well as deep/long diving (DLD) and shallow/short diving (SSD) partitions (Supplementary Table S15). Among them, seven genes (GUCY2C, PDE8A, PDE10A, NT5C3A, PDE3B, ENPP3, and HPRT1) have been both subject to divergent selection and positive selection in cetaceans (Supplementary Table S15). Moreover, several genes (GMPR2, PDE8A, NT5C3A, PDE1B, PAICS, AMPD3, and HPRT1) showed significantly elevated ω in the DLD group (Supplementary Table S15), suggesting that deep/long diving cetaceans experienced an accelerated selection. We next performed phylogenetic generalized least squares (PGLS) analysis to assess the molecular evolution of purine metabolism-related genes and the diving ability of cetaceans (Figure 1; Supplementary Table S16). We found a significant association between the molecular rate (root-to-tip d_N/d_S) of five genes and average dive depth, two genes with average dive time, four genes with maximum dive depth, and four genes with maximum dive time (Supplementary Table S16). Interestedly, four genes (ADCY6, PAICS, PED6A, and PRPS2) presented significant associations with at least two diving associated variables (Supplementary Figure S1). Among these, the de novo purine biosynthesis enzyme phosphoribosylaminoimidazole carboxylase (PAICS) showed association with three variables (Supplementary Table S16 and Figure S1).

DISCUSSION

Oxidative stress drives shared selection signatures of purine biosynthesis genes across mammals

Although species under high oxidative stress live in different natural environment, they all experience ischemia/reperfusion that is damaging to most mammals (Hermes-Lima & Zenteno-Savı́n, 2002; Ramirez et al., 2007). Therefore, these divergent species may have arrived at broadly similar ecologies and phenotypes through convergent or independent genetic routes (Hao et al., 2019). In recent years, studies on convergent molecular mechanisms have illuminated that adaptive convergence may occur at different hierarchical levels, including protein-coding genes, amino acid substitutions, and non-coding regions. For example, convergent amino acid substitution in marine mammals possibly contribute to their musculoskeletal system and thermoregulation (Yuan et al., 2021). In addition, convergent accelerated evolution of conserved non-coding elements may underlie the loss of flight in the bird clade Palaeognathae (Sackton et al., 2019). In line with these studies, we found evidence of positive selection in 72 out of 117 purine metabolism genes in oxidative stress-tolerant mammals. A quarter (25%) (18/72) of these PSGs (i.e., NT5E, P2RX6, P2RX5) are shared by at least one oxidative stress tolerant species, and most of them (12/18) are involved in de novo purine biosynthesis pathway. Moreover, three genes (NME3, ADCY5, and PDE8B) are under selection in three lineages. A similar trend was observed for convergent analyses. Out of 32 genes with putative convergent amino acid changes, a fraction of genes (e.g., ADA, ENPP3, NPR2, GUCY2C, PDE6A, and PDE7A) containing putative convergent amino acids along lineages under high oxidative stress, and 19 genes (59.3%) were also under positive selection (Supplementary Table S8). These convergent genes may be related to resistance to oxidative stress of mammals.

NT5E encodes CD73, which converts AMP to adenosine. Elevated CD73 is associated with high-altitude hypoxia. CD73-deficient mice showed severe tissue hypoxia, inflammation, and lung injury (Liu et al., 2016). Convergent signals of NT5E detected in oxidative stress-tolerant mammals (e.g., cetaceans and naked mole rat) may suggest that selection acted on the same gene to produce similar adaptations to oxidative stress. NME3 encodes nucleoside diphosphate kinase, an enzyme necessary for generating various NTPs, including the signaling molecule GTP, and thus essential for purine biosynthesis (Yu et al., 2017). One convergent positively selected gene, ADA (adenosine deaminase), encodes a key purine salvage enzyme. Its deficiency leads to an accumulation of toxic purine degradation by-products (Flinn & Gennery, 2018). Intracellular cAMP and cGMP are ubiquitous intracellular second messengers that regulate important cellular functions (Wu et al., 2013). Intracellular levels of cGMP and cAMP are dynamically regulated by the rate of their synthesis by cyclases and their hydrolysis by phosphodiesterases (PDEs) (Lucas et al., 2000; Sutherland et al., 1968). Interestedly, we observed that several genes with cGMP and cAMP regulation roles exhibited convergent signals in multiple species adapted to high oxidative stress environments (Figure 5). Of these genes, ADCY3 (adenylate cyclase 3) and ADCY5 (adenylate cyclase 5) convert ATP/GTP into cyclic AMP/GMP (Furuhashi, 2020). NPR2 and GUCY2C are guanylyl cyclases that catalyze the conversion of GTP to cGMP. Stimulation of guanylyl cyclases and the resultant accumulation of cGMP regulates complex signaling cascades (Lucas et al., 2000). Moreover, GUCY2C mutations in mice or humans lead to dysregulation of the GUCY2C-cGMP signaling axis and disease (Rappaport & Waldman, 2018). PDE6A and PDE7A, which encodes enzymes that hydrolyze cyclic nucleotide second messengers (cAMP and cGMP) (Kokkonen & Kass, 2017), have convergent sites shared by cetaceans, seals, and the black flying fox (Supplementary Table S9). It should also be noted that cAMP and cGMP levels are related to oxidative stress. For example, increasing the levels of cAMP in spontaneously hypertensive rats may have a protective effect by reducing oxidative stress (Saha et al., 2008). Additionally, extracellular cGMP and its metabolite GMP have been reported to be protective compounds in the paradigm of oxidative glutamate toxicity (Albrecht et al., 2013). We, therefore, hypothesize that these genetic innovations contribute to oxidative stress adaptation of oxidative increasing the levels of cAMP and cGMP.

Figure  5.  Simplified diagram of purine metabolism pathways

Red ellipses indicate convergent genes involved in oxidative stress adaptation, orange ellipses represent genes under positive selection or with unique amino acid changes in cetaceans, blue ellipses genes inactivated or lost in cetaceans. The composition of the purinosome was obtained from (Pedley & Benkovic, 2017). Abbreviations are as follows: 5′-phosphoribosyl 1-pyrophosphate (PRPP), Phosphodiesterases (PDEs) phosphoribosylpyrophosphate amidotransferase (PPAT), phosphoribosylglycinamide synthetase/phosphoribosylglycinamide formyltransferase/phosphoribosylaminoimidazole synthetase (GART), phosphoribosyl formylglycinamidine synthase (FGAMS or PFAS), phosphoribosyl aminoimidazole carboxylase (CAIRS)/phosphoribosyl aminoimidazole succinocarboxamide synthetase (SAICARS) (PAICS), adenylosuccinate lyase (ADSL), 5-aminoimidazole-4-carboxamide ribonucleotide formyltransferase/inosine 5'-monophosphate cyclohydrolase (ATIC), inosine 5'-monophosphate (IMP), IMP dehydrogenase (IMPDH), xanthosine monophosphate (XMP), guanosine monophosphate (GMP), adenosine monophosphate (AMP), adenylosuccinate (AMPS), adenylosuccinate synthase (ADSS1, ADSS2), NME/NM23 nucleoside diphosphate kinase 3 (NME3), adenosine kinase (ADK), adenosine deaminase (ADA), 5'-Nucleotidase Ecto (NT5E), ectonucleotide pyrophosphatase/phosphodiesterase 3 (ENPP3), nucleoside-triphosphatase cancer-related (NTPCR), guanylate cyclase (GUCY2C, NPR2), adenylate cyclase (ADCY3, ADCY5), GMP synthase (GMPS), purine nucleoside phosphorylase (PNP), hypoxanthine-guanine phosphoribosyl transferase (HGPRT), xanthine dehydrogenase (XDH), xanthine oxidase (XO), urate oxidase (UOX), and allantoicase (ALLC).

DownLoad: Full-Size Img PowerPoint

Despite convergent phenotypes across oxidative stress-tolerant species, examples of identical molecular convergence are rare (Hindle, 2020). It should be noted that molecular convergence in protein function driven by nonparallel amino acid substitutions has been reported across a wide range of proteins. For instance, increases in hemoglobin-oxygen binding affinity occurred in high-altitude species, but the molecular changes underlying the hemoglobin amino acid substitutions are variable (Natarajan et al., 2016). Furthermore, substitutions at different sites of rhodopsin in high-altitude catfishes convergently accelerated the kinetic rates of rhodopsin, destabilizing its ligand-bound forms (Castiglione et al., 2018). These findings suggest that while an adaptive phenotypic change is, to some extent, predictable, the molecular pathways to these changes may not be because of differences in the genetic background of species. In line with this, we also found nonparallel amino acid changes unique to purine metabolism-related genes of oxidative stress-tolerant species. Enzymes involved in de novo purine biosynthesis pathway activation form dynamic multienzyme complexes, termed purinosomes (Pareek et al., 2020; Pedley & Benkovic, 2017). ATIC, located in purinosomes, catalyzes the final two steps of purine de novo biosynthesis (Pedley & Benkovic, 2017). Interestingly, nonoverlapping amino acid changes exclusive to ATIC were identified in oxidative stress-tolerant species (Figure 2D), suggesting that similar high oxidative stress environments target different sites of ATIC in these species by nonparallel molecular mechanisms. Strikingly, specific substitutions of cetaceans (G106A and G113A), Tibetan antelope (V17E), and manatee (A110S) are found in the ATIC MGS-like domain. ATIC mutations affect the assembly of functional purinosomes (Baresova et al., 2012; Ramond et al., 2020). A similar trend was also found in XDH. As we describe in detail below, thermal stability and conversion rate, combined with XDH enzyme activity, imply that natural selection has targeted these unique sites due to their potential capacity to modulate functions of these convergent genes at high oxidative stress environments.

Taken together, although further studies on the genes identified in our study are warranted, we propose that our results reflect that positive and convergent selection on purine metabolism-related genes, especially in purine de novo biosynthesis, contribute to the convergent phenotypes between species that tolerate high oxidative stress environments.

Purine salvage genes contributing to oxidative stress adaptation in cetaceans

Cetaceans are excellent models for studying the molecular basis of oxidative stress tolerance. They are exposed to reperfusion/reoxygenation during diving. During this process of oxygen deprivation, purine metabolites accumulate due to ATP breakdown, which is connected with many pathological conditions in humans. The physiological adaptations to breath-hold diving in cetaceans have been well described (López-Cruz et al., 2016; Ramirez et al., 2007). However, the precise genetic mechanism that facilitated oxidative stress adaptation in cetaceans remains unknown. In the present study, we discovered evidence of positive selection on genes associated with the purine salvage pathway in cetaceans (Figure 5). The protein product of HPRT1 is hypoxanthine-guanine phosphoribosyltransferase (HGPRT), an important enzyme for the re-utilization of HX and guanine (the conversion of hypoxanthine/guanine and IMP/GMP, respectively). HGPRT-deficiency manifests as an accumulation of DNA and nucleotide lesions and, hence, replication errors in genetic disorders such as Lesch-Nyhan disease and chronic gout (Agrahari et al., 2019; Nguyen & Nyhan, 2016; Schretlen et al., 2016). The bottlenose dolphin has higher plasma HX and IMP concentrations (i.e., HGPRT substrate and product, respectively) but a lower intraerythrocytic HGPRT activity than humans, suggesting an elevated purine recycling rate (López-Cruz et al., 2016). Additionally, we detected positive selection and specific substitutions in cetacean DGUOK, an enzyme that mediates the phosphorylation of purine deoxyribonucleosides in the first step of the mitochondrial purine salvage pathway (El-Hattab et al., 2017). We propose that selection acting on HPRT1 and DGUOK may be attributable to adaptive evolution to limit the accumulation of non-recyclable purines (xanthine and uric acid) and to deliver nucleotide precursors to ischemic tissues during diving.

An alternative adaptation to oxidative stress is evident from the functional changes of XDH. We found six unique amino acid changes in cetacean XDH, despite no evidence of positive selection acting on this gene, and assessed its function in vitro. XDH is a rate-limiting purine degradation enzyme, catalyzing the final two-step oxidation of HX and xanthine to uric acid (Figure 5) (Battelli et al., 2016). XDH can be converted into XO to promote ROS generation during ischemic/hypoxic conditions. Since XO increases ROS by generating superoxide and hydrogen peroxide (Furuhashi, 2020), the decreased conversion rate from XDH to XO and higher thermal stability of cetacean XDH may reduce oxidative stress. We show that the unique substitutions in cetacean XDH enhance enzyme activity and mediate a higher affinity for xanthine and NAD⁺, possibly resulting in elevated levels of uric acid in cetaceans. Although sustained elevated uric acid levels come at the cost of high risk for gout (Maiuolo et al., 2016), uric acid is also considered an antioxidant and neuroprotectant (El Ridi & Tallima, 2017). Therefore, there might be a balance between the deleterious and protective properties of uric acid in cetaceans.

Inactivation of the cetacean ALLC gene and its implication for mammalian uric acid metabolism

It is important to note that it is not well-established how most mammals metabolize uric acid. For example, the elephant and manatee have lost urate oxidase (UOX, which converts uric acid to allantoin) yet can maintain low uric acid levels via unknown molecular mechanisms (Sharma & Hiller, 2020). In line with this idea, we found that ALLC (which converts allantoin to allantoic acid) is “missing” or inactivated in all cetaceans, suggesting that the enzyme became dispensable upon migration to an aquatic environment. We show an increased ALLC d_N/d_S value in cetaceans, implying relaxation of selective constraints at this locus. A d_N/d_S value of 0.926 on the cetacean branch indicates that cetacean ALLC shifted from purifying selection to neutral evolution in early evolutionary time (inactivation soon after cetaceans diverged from their sister group). In agreement, pseudogene dating suggests that ALLC was inactivated 48.2~53.6 Ma, after the split of the cetacean and artiodactyl lineages.

This finding begs the question of how cetaceans can “afford” to lose the ALLC. The formation of uric acid from purine metabolism is common to all vertebrates; however, uric acid degradation varies from species to species. Reptiles, birds, and hominoids among mammals excrete uric acid, whereas other mammals (like cetaceans) possess uricase and excrete allantoin (Álvarez-Lario & Macarrón-Vicente, 2010). Although an intact ALLC gene or transcripts are found in some species (e.g., humans and mice), it is widely accepted that allantoicase activity is lost in mammals (Vigetti et al., 2002a, 2002b). Therefore, it is plausible that ALLC has been rendered obsolete in cetaceans. However, for species with an intact allantoicase gene, we hypothesize that ALLC of the mammal ancestor acquired one or more mutations that led to the loss of enzymatic activity associated with purine catabolism but escaped the fate of becoming a pseudogene – similar to the loss of UOX in vertebrates (Keebaugh & Thomas, 2010; Sharma & Hiller, 2020). Nevertheless, we cannot rule out that the retained ALLC gene in other mammals has critical functions beyond purine metabolism.

Divergent selective pressures are associated with various diving phenotypes in cetaceans

A previous study found that shallow/short and deep/long divers among marine mammals have distinct defence strategies against ROS formation and oxidative stress (Cantú-Medellín et al., 2011). For example, the deep/long-diving Kogia spp. have higher superoxide radical production than the shallow/short-diving bottlenose dolphin (Cantú-Medellín et al., 2011). In addition, oxidative stress biomarkers (e.g., superoxide dismutase, catalase, and glutathione peroxidase activity) in the deep/long-diving Southern elephant seal (Mirounga leonine) were significantly higher than the shallow/short South American fur seal (Arctocephalus australis) (Righetti et al., 2014). These findings suggest that a higher antioxidant content is needed to counteract the higher ROS production. Here, we found that the divergent selective pressures of purine metabolism-related genes, such as AMPD3, among cetacean groups with different diving abilities, suggest that distinct responses to ischemia (different oxidative stress) underly this trait. AMPD3 is the red blood cell-specific isoform of AMP deaminase stimulated by oxidative stress. Nemkov et al. (2018) reported that storage-dependent hypoxanthine accumulation could be ameliorated by hypoxia-induced decreases in AMPD3 activity. The evolutionary rate of the AMPD3 gene is significantly higher in cetaceans with excellent diving ability, suggesting that these cetaceans have an enhanced ability to tolerant higher oxidative stress (Cantú-Medellín et al., 2011). Coincident with divergent selective pressure between long and short diving cetaceans, there was also a significant positive association between average d_N/d_S of ten genes and diving phenotypes, indicating a high rate of evolution of purine metabolism-related genes in long/deep-diving cetaceans. The increase in evolutionary rates of purine-metabolism genes could be related to an accumulation of non-synonymous mutations, which could drive molecular innovations of oxidative stress defences of diving cetaceans.

SUPPLEMENTARY DATA

Supplementary data to this article can be found online.

COMPETING INTERESTS

The authors declare that they have no competing interests.

AUTHORS’ CONTRIBUTIONS

G.Y. and I.S. conceived and supervised the study. R.T., S.M.C., and H.G. collected the data and conducted the genomic analyses. C.Y. performed cell-based experiments. R.T. wrote the manuscript, and G.Y. and I.S. revised and directed it. All authors read and approved the final version of the manuscript.