Comparative transcriptomic analysis revealed adaptation mechanism of Phrynocephalus erythrurus, the highest altitude Lizard living in the Qinghai-Tibet Plateau
© Yang et al.; licensee BioMed Central. 2015
Received: 8 January 2015
Accepted: 29 April 2015
Published: 2 June 2015
Organisms living at high altitudes must overcome three major environmental challenges: hypoxia, cold, and intense UV radiation. The molecular mechanisms that enable these challenges to be overcome have mainly been studied in endothermic organisms; relatively little attention has been paid to poikilothermic species. Here, we present deep transcriptome sequencing in two closely related lizards, the high altitude-dwelling Phrynocephalus erythrurus and the lowland-dwelling P. putjatia, to identify candidate genes under positive selection and to explore the convergent evolutionary adaptation of poikilothermic animals to high altitude life.
More than 70 million sequence reads were generated for each species via Illumina sequencing. De novo assembly produced 56,845 and 63,140 transcripts for P. erythrurus and P. putjatia, respectively. P. erythrurus had higher Ka/Ks ratios than P. putjatia, implying an accelerated evolutionary rate in the high altitude lizard lineage. 206 gene ontology (GO) categories with accelerated evolutionary rates and 43 candidate positively selected genes were detected along the P. erythrurus lineage. Some of these GO categories have functions associated with responses to hypoxia, energy metabolism and responses to UV damage. We also found that the high-altitude ranid frog R. kukunoris had higher Ka/Ks ratios than the closely related low-altitude frog R. chensinensis, and that the functional categories with accelerated evolutionary rates in R. kukunoris overlapped extensively with those detected along the P. erythrurus lineage.
The mechanisms of high altitude adaptation in P. erythrurus were tentatively inferred. By comparing two pairs of low- and high-altitude poikilothermic species, we found that similar functional categories had undergone positive selection in high altitude-dwelling Phrynocephalus and Rana lineages, indicating that similar mechanisms of adaptation to high altitude might have evolved in both genera. Our findings provide important guidance for future functional studies on high altitude adaptation in poikilothermic animals.
Survival at high altitudes is very challenging. Nevertheless, many native peoples and animals can thrive under the cold and hypoxic conditions associated with high altitude environments . The molecular responses to high-altitude stress have been studied for over a century, but most of the previous studies focused on a single or a few candidate genes in model systems, which has limited our understanding of the molecular basis of adaptation in non-model systems [2–4]. The recent development of next-generation sequencing (NGS) technologies has significantly accelerated the speed of gene discovery and genomics studies, greatly facilitating evolutionary and ecological research on non-model organisms [5, 6]. Numerous studies on the adaptation of endothermic animals to high altitude have been published; species investigated at the genome scale include the yak (descendants of wild yak) , Tibetan wild boar , ground tit , snow leopard  and Tibetan mastiff , while transcriptome-scale studies have examined the blind subterranean mole rat , deer mice [13, 14] and rufous-collared sparrows . These studies resulted in the identification of many genes associated with hypoxia tolerance that have undergone natural selection in high altitude-dwelling species or have different expression patterns in high altitude species between closely related low altitude dwellers. However, there have been few investigations into the mechanisms of high altitude adaptation in poikilothermic animals.
Toad-headed sand lizards of the genus Phrynocephalus are widespread in the Qinghai-Tibet Plateau (QTP) . They comprise over 40 species and represent an important component of low and high-altitude desert ecosystems . Despite increasing awareness of their ecological importance, the ecological niches and roles of Phrynocephalus species in lizard communities of extreme environments are poorly studied. Phrynocephalus erythrurus, which typically lives at altitudes of 4,500 m above sea level (a.s.l) or more, is considered to be the most high altitude-adapted lizard in the world . It has a congener, P. putjatia, which is widely distributed across the Gobi desert and the semi-desert areas of the Qinghai Lake Basin, and commonly lives at lower altitudes of around 2,500 m.a.s.l . The adaptations of these two Phrynocephalus species to their respective habitats are likely to have contributed significantly to species diversification within the genus, making Phrynocephalus an ideal system for studying high altitude adaptation in poikilothermic animals. Previous investigations have examined differences in the efficiency of ATP production in the mitochondria of sand lizard species dwelling at different altitudes, as well as differences in their metabolic activity and physiology [20, 21]. However, little is known about the molecular mechanisms that underpin adaptation to life at very high altitudes in this genus. Such studies have been hindered by the absence of genomic and transcriptomic resources for these non-model organisms.
Here we present a large scale, multi-tissue, multi-individual transcriptome for P. erythrurus and P. putjatia that was derived using RNA-seq. This transcriptome was created for two main reasons. First, to develop a molecular resource to support studies on sand lizards and make it available to the scientific community. Second, to identify adaptive evolutionary patterns in sand lizards, which will hopefully provide important insights into the genetic basis of high altitude adaptation and speciation. An additional goal of this study was to analyse transcriptomic data for ranid species dwelling at high and low altitudes , and to determine whether the evolutionary adaptations to life at high altitude in these cold-blooded organisms were functionally convergent with those seen in the sand lizards.
Results and Discussion
Sequencing and De novo assembly
Summary of transcriptome data for P. erythrurus and P. putjatia
Total number of reads
Read length (bp)
Total length of the final assembly (Mb)
Total number of transcripts
N50 length of assembly (bp)
Mean length of assembly(bp)
Median length of assembly(bp)
Total number of predicted proteins
N50 length of predicted proteins (AA)
Mean length of predicted proteins (AA)
bp = base pair; Mb = mega base pairs.
Functional annotation of unigenes
The GO category distributions of the unigenes for both species were highly similar (Fig. 1c). The three main functional categories in each case were ‘Biological process’ (represented by 3,799 unigenes from P. erythrurus and 4,235 from P. putjatia), ‘Cellular component’ (4,037 from P. erythrurus and 4,552 from P. putjatia), and ‘Molecular function’ (3,854 from P. erythrurus and 4,218 from P. putjatia). Large groups of unigenes from both species (13,697 and 15,219 from P. erythrurus and P. putjatia, respectively) were assigned to all three categories. The COG database was used to define the orthologous functions of the identified unigenes. The distributions of unigene function categories for the two sand lizards were also very similar (Fig. 1d).
Lizard phylogenetic trees presented in previous publications were generally constructed on the basis of mitochondrial DNA sequence data and/or a few nuclear genes [24–26]. The development of transcriptome sequencing technologies in recent years has made it possible to determine phylogenies more efficiently than using traditional PCR- and EST-based methods . Several publications have demonstrated the utility of transcriptome data for resolving the relationships of Spalacidae , Passerida , and the large tetrapod group consisting of turtles, birds and crocodiles .
Hypotheses concerning mechanisms of high altitude adaptation in P. erythrurus
By combining our results with those of earlier studies, we were able to tentatively infer the mechanisms that underpin adaptation to high altitude in P. erythrurus. The partial pressure of oxygen decreases with increasing altitude, so organisms living at high altitude must be tolerant of chronic hypoxia . This tolerance is established in two ways, the first of which is based on increasing the capacity to take up and transport oxygen . Genes associated with the “response to hypoxia” promote oxygen binding and transport to facilitate the capture of oxygen from the air. In particular, the HIF-1α pathway contributes to hypoxia tolerance by activating the expression of vascular endothelial growth factor to promote the formation of new blood and thereby maintain an adequate supply of oxygen . The capillary density and blood hemoglobin concentration of the high-altitude lizards are both significantly greater than those of their low-altitude counterparts, further increasing their capacity for oxygen transport . The second way of adapting to the low-oxygen conditions associated with high altitudes is to reduce oxygen consumption. This typically involves reducing the rates of metabolic processes and ATP demand (hypometabolism) . Metabolic suppression entails reducing the mitochondrial rate while maintaining a core temperature high enough to permit survival and supplying enough ATP to sustain core physiological functions. This necessitates an increase in mitochondrial efficiency such that more ATP and heat are generated while consuming less oxygen. ATP production in such cases is largely driven by aerobic respiration, and lipid oxidation accounts for a relatively high proportion of the organism’s energy supply . The finding that genes associated with energy metabolism and lipid metabolic process show accelerated divergence in P. erythrurus may thus reflect adaptation to high altitude. In addition to needing an increased capacity for oxygen uptake and transport, and a reduced metabolic rate with increased mitochondrial efficiency, high altitude dwelling organisms must cope with higher levels of UV radiation than their lowland counterparts. UV exposure can damage DNA molecules by generating highly reactive chemical intermediates such as oxygen radicals . This presents a particularly acute challenge to reptiles, which often rely on basking for thermoregulation . This increased level of DNA damage appears to have caused accelerated divergence of genes associated with DNA damage checking and DNA repair in P. erythrurus. While these findings are intriguing, it should be noted that analyses based on GO enrichment and computational approaches for detecting genes with accelerated evolutionary rates can yield high false-positive rates. Consequently, additional functional and physiological experiments are needed to verify the contributions of the identified genes to high altitude adaptation .
Detecting candidate genes under positive selection in P. erythrurus
To identify genes that are likely to be important in the high altitude adaptation of P. erythrurus, we used rigorous criteria to identify candidate positively selected genes (PSGs) in the lineage. Previous studies that attempted to detect genes subject to positive selection in specific lineages often had worryingly high false positive rates due to the use of overly loose criteria [36–38], so rigorous criteria were employed in this work to obtain more robust results. We initially used the PRANK (Probabilistic Alignment Kit) to align the orthologs, which is considered to be more accurate than alternative methods and to provide a lower level of false positives . Rather than performing simple pairwise comparisons of Ka/Ks ratios for P. erythrurus and P. putjatia, we used the branch-site likelihood ratio test to generate a list of candidate PSGs . This approach was shown to have reasonable power and an acceptable false positive rate under a variety of selection schemes. We then manually removed all candidate PSGs with potential errors in their alignments to minimize the false discovery rate. After applying these strict criteria, 43 candidate high altitude responsive genes remained on the list for P. erythrurus (Additional file 5: Table S3).
Convergent evolution in ectothermic animals adapted to life at high altitudes
Convergence is the independent evolution of similar features in different species. Its occurrence supports the hypothesis that specific environmental challenges can induce species to evolve in predictable and repeatable ways . Studies on convergence between native animals living at high altitudes was very important to understanding the mechanisms of adaptation to high altitude. Genes associated with three distinct functions (response to hypoxia, energy metabolism and response to UV damage) were found to have been subject to adaptive selection in P. erythrurus during its adaptation to life at high altitude. We therefore sought to determine whether adaptation based on these functions is common in poikilothermic organisms living at high altitudes.
One notable early study on high altitude adaptation involved a pair of ranid frog species, the high altitude dwelling Rana kukunoris and the low altitude dwelling R. chensinensis . We reanalyzed the sequence data for these two species to identify GO categories exhibiting accelerated divergence using the protocol developed for use with the toad-headed lizards. R. kukunoris was found to have a higher branch Ka/Ks ratios than R. chensinensis (P < 2.2 × 10−22, binomial test, Fig. 2b). 255 GO categories showing putatively accelerated evolutionary rates (P < 0.05, binomial test, Additional file 4: Table S2) were detected in R. kukunoris, with functional annotations similar to those identified in P. erythrurus. Thus, GO categories associated with the functions response to hypoxia, energy metabolism, and response to UV damage were all found to exhibit accelerated divergence in the high altitude dwelling R. chensinensis. Specifically, the GO categories “response to hypoxia” and “response to oxidative stress” were associated with the function of response to hypoxia; the categories “respiratory electron transport chain”, “regulation of lipid metabolic process”, “glycerophospholipid metabolic process” and “lipid biosynthetic process” were associated with energy metabolism; and the categories “response to UV”, “DNA integrity checkpoint”, “DNA damage checkpoint”, and “single-stranded DNA binding” were associated with the function of response to UV damage. We also detected 49 candidate PSGs along the R. kukunoris lineage (Additional file 5: Table S3), two of which (UBE2D1 and NBN ) overlap with candidate PSGs for the P. erythrurus lineage. The overlap presumably reflects convergent evolution relating to the response to hypoxia and DNA damage repair. This finding is consistent with previous studies that have demonstrated convergent evolution associated with hypoxia adaptation in different high-altitude Andean hummingbirds  and in Tibetan Mastiffs and Tibetan people .
Studies on two separate pairs of high- and low-altitude poikilothermic species from the QTP revealed significantly higer Ka/Ks ratios for high altitude species than for their low altitude counterparts, which suggests that the high altitude species may have experienced adaptive evolution that allows them to cope with their extremely inhospitable environment. However, the elevated Ka/Ks ratios could also be due to a relaxation of selective pressure; further population genomic analyses will be required to test this alternative hypothesis. We identified three main function groups showing evidence of accelerated evolution, suggesting that genes from these categories may play key roles in adaptation to high altitude life among poikilotherms. While we were writing this article, another paper on high altitude adaptation in lizards was published in which the highland dwelling P. vlangalii was compared to its lowland relative P. przewalskii . That work similarly indicated that genes associated with responses to hypoxia and UV damage had undergone adaptive evolution in the high altitude species, which is consistent with our findings and suggests that the mechanisms of adaptation identified for P. erythrurus may be common to many highland-dwelling cold-blooded organisms. To properly test this hypothesis, it will be necessary to perform additional functional and physiological experiments that should be integrated with genomic/transcriptomic analyses of a wider range of high altitude-dwelling poikilothermic species [13, 14].
We have successfully sequenced and annotated large-scale, multi-tissue, multi-individual transcriptomes for two ectothermic vertebrate species and enriched the expressed sequence data available for the genus Phrynocephalus. The transcriptome data sets reported here represent a valuable resource that will support further expressional and functional studies that will help interested researchers to address ecological and evolutionary questions concerning sand lizards and other lizard species. High altitude adaptation is a complex process that involves numerous genes. We identified several functional groups and genes that may have undergone accelerated evolution and positive selection in the highland-dwelling P. erythrurus and R. kukunoris lineages. Similar GO categories with accelerated evolutionary rates and candidate positively selected genes were identified in both species. In addition, three main gene functions that might contribute to high altitude adaptation were identified: responses to hypoxia, energy metabolism, and responses to UV damage. We hypothesise that these three functions are generally important in adaptation to high altitude life among poikilotherms. The findings presented herein increase our understanding of the mechanisms by which cold blooded animals adapt to highland life and will support further studies on high altitude adaptation among lizards, other reptiles, and poikilothermic species in general.
All animal specimens were collected legally. Animal collection and utility protocols were approved by the Ethics Committee of Animal Experiments at Lanzhou University and in accordance with guidelines from the China Council on Animal Care.
Adult P. putjatia individuals were captured in June, 2013, from Gui’de County at an altitude of 2,312 m above sea level (m.a.s.l). Adult P. erythrurus individuals were captured in July, 2013 from the Tuotuo River at an altitude of 3,557 m.a.s.l. To include as many expressed genes as possible, we sampled seven different tissues from one male and one female of each species, including brain, heart, liver, kidney, testes/ootheca, lung and skeletal muscle. Tissue samples for RNA extraction were stored at −70 °C in a cryogenic refrigerator following euthanasia and dissection.
Total RNA was isolated using the TRIzol reagent (Invitrogen, Carlsbad, CA, USA) and the RNeasy Kit (Qiagen, Hilden, Germany). A single pooled RNA sample was prepared for each species by mixing equal volumes of the RNA extracted from each tissue sample. These pooled samples were then used for cDNA preparation and RNA-Seq. cDNA library construction and sequencing were performed by Encode Genomics (Suzhou, China). RNAs with poly(A) tails were purified from total RNA using oligonucleotide (dT) magnetic beads and fragmented into short sequences that were used to template cDNA synthesis by PCR. The resulting cDNA libraries were purified using a PCR extraction kit (QiaQuick) and had an insert size of 200–700 base pairs (bp). The paired ends of the library were sequenced with a read length of 90 bp using the Illumina HiSeq 2000 sequencing platform. The image data generated by the sequencer were converted into raw sequence data by base calling to yield the raw reads, which were stored in the FastQ format. The entire process followed a standardized procedure and was monitored using a standard Quality Control System. All reads were deposited in the NCBI Short Read Archive (SRA) under accession number SRP050887.
Data filtration and De novo assembly
We first cleaned the raw sequence reads by removing the exact duplicates from both sequencing directions with the FASTX Toolkit . We further trimmed reads by removing adapter sequences, reads for which unknown base calls (N) represented more than 5 % of all bases, low complexity sequences, and low quality sequences (i.e. sequences for which >45 % of the bases were called with a quality score of ≤7). De novo assembly of the clean reads was performed using the Trinity software package, which efficiently constructs and analyses sets of de Bruijn graphs and then fully reconstructs a large fraction of transcripts, including alternatively spliced isoforms and transcripts from recently duplicated genes . The calculations were performed using the –CuffFly option, which enforces the use of the cufflinks-like algorithm to report minimum transcripts. The transcripts were then clustered by CD-HIT-EST , with a sequence similarity cutoff of 95 %. The RSEM package  within Trinity was used to obtain abundance estimates for the transcriptome assemblies. The final assemblies were produced after deleting transcripts with FPKM values below 1 to ensure that all of the included transcripts could be detected. Partial and complete open reading frames (ORFs) were predicted using the TransDecoder script present in the Trinity package, with a minimum length of 100 amino acids.
Unigenes and predicted protein sequences were used as queries to search protein databases using the BLAST program . Queries were performed against the NCBI (National Center for Biotechnology Information) NR (non-redundant database), Swiss-Prot, COG (cluster of orthologous groups databases) and KEGG (Kyoto Encyclopedia of Genes and Genomes) databases. Homology searches were conducted using BLASTx with an e-value cut-off of 1E-05 (in the case of the COG database, 1E-20 was used for increased stringency). Gene ontology (GO) terms were obtained from NR hits using the Blast2GO program  with default parameters for the mapping and annotation steps, except that an e-value cutoff of 1E-6 was applied and only the first 20 hits from the BLAST process were considered. GO functional classifications were obtained for all unigenes using the WEGO (Web Gene Ontology Annotation Plot) software . COG function classification analysis of all unigenes sequences was performed using in-house Perl scripts. The metabolic pathways were mapped using the KAAS (KEGG Automatic Annotation Server)  with a bi-directional best-hit strategy to assign KEGG orthology terms (KO) to unigenes. The identified pathways were settled using their respective KO assignments.
Ortholog identification and sequence alignment
To explore the phylogenetic relationships of the major lizards, characterize the mechanisms of adaptation to high altitude in P. erythrurus, and determine whether there was any evidence of functional convergent evolution among highland-dwelling poikilotherms, we downloaded sequence data for five lizard species, three amphibians, and the zebra fish as an outgroup. For the three genome-sequenced species (Anolis carolinensis, Xenopus tropicalis and Danio rerio), the CDS and protein sequences from the Ensemble database were downloaded (release 74). For the four transcriptome-sequenced lizards (Pogona vitticeps, Chamaeleo chamaeleon, Chalcides ocellatus, and Sphenodon punctatus), the unigenes were downloaded from the supplementary data for the corresponding publications. For the two transcriptome-sequenced ranid frogs (Rana chensinensis and Rana kukunoris), sequencing reads were downloaded and processed using Trinity to assemble the transcripts, after which Transdecode was used to extract CDS and protein sequences. Additional file 6: Table S4 provides details of websites from which all of these datasets can be downloaded.
Orthologs were identified using InParanoid  and MultiParanoid  by considering protein sequences from 11 species (A. carolinensis, X. tropicalis, D. rerio, P. vitticeps, C. chamaeleon, C. ocellatus, S. punctatus, R. chensinensis, R. kukunoris, P. erythrurus and P. putjatia). In-house Perl scripts were used to identify high-confidence 1:1 orthologs based on rigorous criteria: genes in each cluster with confidence scores of less than 1 were removed and only clusters containing a single gene from each species were selected. All orthologs were aligned using PRANK  with the codon option, and alignments shorter than 150 bp after removing sites with ambiguous data were discarded. Finally, we obtained two multiple sequence alignments: (1) one alignment (9,175 bp, 10.1 % of total length of the 206 genes ) for the concatenated 4D-sites of the 206 high-confidence single-copy genes among 11 species, which was used to support the phylogenetic analysis and estimate divergence times; and (2) 4,094 individual alignments for each high-confidence single-copy gene within six of the studied species (Xenopus tropicalis, Danio rerio, P. erythrurus, P. putjatia, Rana chensinensis and Rana kukunoris) and at least one of the other five species (Anolis carolinensis, Pogona vitticeps, Chamaeleo chamaeleon, Chalcides ocellatus and Sphenodon punctatus), which was used to investigate the mechanisms of adaptation to high altitude.
The jModelTest program  was used to select the best fitting substitution model according to the Akaike information criterion based on the 9,175 bp concatenated 4D-sites. The GTR + I + G model was determined to be the best fitting, and the Phyml program  was employed to build the ML tree. The bootstrap support for all nodes was 100 %. Based on this phylogenetic tree, we estimated the divergence time of each node by the MCMCtree program in the PAML package  and the BEAST program . For the MCMCtree program, we set the clock parameter to 1 to ensure the use of a strict molecular clock. Two-round analyses were performed in which the first 5,000 iterations were discarded as a burn-in followed by sampling every 50 iterations until 20,000 samples had been gathered. These two-round analyses consistently yielded convergence; one representative set of obtained results is presented in this work. For the BEAST program, a strict clock was assumed and the analyses were run for 10,000,000 generations, sampling 1 generation in every 1,000; the first 1,000,000 generations were discarded as a burn-in. For both programs, the divergence times were calibrated (on the assumption of a constant molecular clock) against the timing of the divergences of the Zebra fish from the lineage leading to toads/lizards (421.75- 416 Mya), the toads from the lineage leading to lizards (350.1-330.4 Mya), and Sphenodon punctatus from the lineage leading to other lizards (275–285 Mya). The timings of the first two divergences were obtained from http://www.fossilrecord.net/ and that of the third from http://www.timetree.org/. The results obtained using the two program were similar; see Additional file 3: Table S1 for more details.
Analysis of accelerated evolutionary rates along each lineage
We used the Ka/Ks ratios to measure the evolutionary rate along a lineage. The values of Ka and Ks and the Ka/Ks ratios were estimated for each of the 4,094 single-copy orthologs using the Codeml program with the free-ratio model for each branch . 10000 concatenated alignments constructed from 150 randomly chosen orthologs were used to estimate lineage-specific mean values. Based on the GO classifications, the unique orthologs were clustered into different functional terms and the Ka, Ks and Ka/Ks ratios for each term was calculated. A comparison of evolutionary rates based on non-synonymous substitution between P. erythrurus and P. putjatia or R. chensinensis and R. kukunoris was conducted using a binomial test (see  for details of the method used). Only GO categories containing more than 20 genes were included in this analyses.
Detecting candidate genes under positive selection
To clarify the mechanism of adaptation to high elevations, we used the Codeml program from the PAML package with a branch-site model  to detect positively selected genes in any of the P. erythrurus, P. putjatia, R. chensinensis or R. kukunoris lineages, with each lineage being specified as the foreground branch sequentially. Very rigorous criteria were used to identify the PSGs. In the branch-site model, a PSG was required to have a Ka/Ks ratio greater than 1 in the model with positive selection on the foreground branch and to have positive sites for the foreground branch with a posterior probability in excess of 0.95 based on the Bayes empirical Bayes (BEB) results . A likelihood ratio test was also conducted to compare the model with positive selection on the foreground branch to a null model with no positive selection on the foreground branch for each orthologous gene, and only those genes with P values of less than 0.05 were selected [38, 40]. We then manually removed all PSGs with potential errors in their alignments to minimize the false discovery rate.
Protein structure homology modeling
We investigated the relationships between the amino acid sequence mutations and the structure modeling within 14 candidate PSGs, which were associated with response to hypoxia, response to UV and energy metabolism. These candidate PSGs’ functional domains were obtained from the Pfam database  and the protein structures were constructed with the Phyre2 server, which uses the alignment of hidden Markov models via HHsearch to significantly improve the accuracy of alignment and detection rates . We downloaded all of the high confidence results from the Phyre2 server and constructed ribbon representation of the structures with UCSF Chimera version 1.10.1 . The electrostatic potentials at the surfaces of the predicted protein structures were rendered using PyMol  with ABPS and PDB2PQR .
Availability of supporting data
The transcriptome reads of P. erythrurus and P. putjatia have been deposited in the NCBI Short Read Archive (SRA) database under the accession number SRP050887, while the links that can be used to download information concerning the other 9 species considered in this work were list in the Additional file 6: Table S4. The data of the phylogenetic analysis are available in the Dryad Digital Repository: http://0-dx.doi.org.brum.beds.ac.uk/10.5061/dryad.f587b .
This study was supported by grants from the National Natural Science Foundation of China (31200946 and 31322052), the National High Technology Research and Development Program of China (2013AA102505 3–2), the Fundamental Research Funds for the Central Universities (lzujbky-2013-205), and the 985 and 211 Projects of Lanzhou University.
The authors declare that they have no competing interests.
QQ conceived and designed the experiments. QC, XT and MM collected the samples. YY, XZ, QR and JH analyzed the data. YY, LW, KW and QQ wrote the paper. All authors read and approved the final manuscript.
- Storz JF, Scott GR, Cheviron ZA. Phenotypic plasticity and genetic adaptation to high-altitude hypoxia in vertebrates. J Exp Biol. 2010;213(Pt 24):4125–36.View ArticlePubMed CentralPubMedGoogle Scholar
- Pugh CW, Ratcliffe PJ. Regulation of angiogenesis by hypoxia: role of the HIF system. Nat Med. 2003;9(6):677–84.View ArticlePubMedGoogle Scholar
- To KK, Huang LE. Suppression of hypoxia-inducible factor 1alpha (HIF-1alpha) transcriptional activity by the HIF prolyl hydroxylase EGLN1. J Biol Chem. 2005;280(45):38102–7.View ArticlePubMed CentralPubMedGoogle Scholar
- Salnikow K, Kluz T, Costa M, Piquemal D, Demidenko ZN, Xie K, et al. The regulation of hypoxic genes by calcium involves c-Jun/AP-1, which cooperates with hypoxia-inducible factor 1 in response to hypoxia. Mol Cell Biol. 2002;22(6):1734–41.View ArticlePubMed CentralPubMedGoogle Scholar
- Stapley J, Reger J, Feulner PG, Smadja C, Galindo J, Ekblom R, et al. Adaptation genomics: the next generation. Trends Ecol Evol. 2010;25(12):705–12.View ArticlePubMedGoogle Scholar
- Ekblom R, Galindo J. Applications of next generation sequencing in molecular ecology of non-model organisms. Heredity (Edinb). 2011;107(1):1–15.View ArticleGoogle Scholar
- Qiu Q, Zhang G, Ma T, Qian W, Wang J, Ye Z, et al. The yak genome and adaptation to life at high altitude. Nat Genet. 2012;44(8):946–9.View ArticlePubMedGoogle Scholar
- Li M, Tian S, Jin L, Zhou G, Li Y, Zhang Y, et al. Genomic analyses identify distinct patterns of selection in domesticated pigs and Tibetan wild boars. Nat Genet. 2013;45(12):1431–8.View ArticlePubMedGoogle Scholar
- Qu Y, Zhao H, Han N, Zhou G, Song G, Gao B, et al. Ground tit genome reveals avian adaptation to living at high altitudes in the Tibetan plateau. Nat Commun. 2013;4:2071.PubMedGoogle Scholar
- Cho YS, Hu L, Hou H, Lee H, Xu J, Kwon S, et al. The tiger genome and comparative analysis with lion and snow leopard genomes. Nat Commun. 2013;4:2433.PubMed CentralPubMedGoogle Scholar
- Wang GD, Fan RX, Zhai W, Liu F, Wang L, Zhong L, et al. Genetic convergence in the adaptation of dogs and humans to the high-altitude environment of the tibetan plateau. Genome Biol Evol. 2014;6(8):2122–8.View ArticlePubMedGoogle Scholar
- Malik A, Korol A, Weber M, Hankeln T, Avivi A, Band M. Transcriptome analysis of the spalax hypoxia survival response includes suppression of apoptosis and tight control of angiogenesis. BMC Genomics. 2012;13:615.View ArticlePubMed CentralPubMedGoogle Scholar
- Cheviron ZA, Connaty AD, McClelland GB, Storz JF. Functional genomics of adaptation to hypoxic cold-stress in high-altitude deer mice: transcriptomic plasticity and thermogenic performance. Evolution. 2014;68(1):48–62.View ArticlePubMed CentralPubMedGoogle Scholar
- Cheviron ZA, Bachman GC, Connaty AD, McClelland GB, Storz JF. Regulatory changes contribute to the adaptive enhancement of thermogenic capacity in high-altitude deer mice. Proc Natl Acad Sci U S A. 2012;109(22):8635–40.View ArticlePubMed CentralPubMedGoogle Scholar
- Cheviron ZA, Whitehead A, Brumfield RT. Transcriptomic variation and plasticity in rufous-collared sparrows (Zonotrichia capensis) along an altitudinal gradient. Mol Ecol. 2008;17(20):4556–69.View ArticlePubMedGoogle Scholar
- Jin YT, Brown RP. Species history and divergence times of viviparous and oviparous Chinese toad-headed sand lizards (Phrynocephalus) on the Qinghai-Tibetan Plateau. Mol Phylogenet Evol. 2013;68(2):259–68.View ArticlePubMedGoogle Scholar
- Guo X, Wang Y. Partitioned Bayesian analyses, dispersal-vicariance analysis, and the biogeography of Chinese toad-headed lizards (Agamidae: Phrynocephalus): a re-evaluation. Mol Phylogenet Evol. 2007;45(2):643–62.View ArticlePubMedGoogle Scholar
- Zhu L, Liao P, Tong H, Jin Y. The complete mitochondrial genome of the subspecies, Phrynocephalus erythrurus parva (Reptilia, Squamata, Agamidae), a toad-headed lizard dwell at highest elevations of any reptile in the world. Mitochondrial DNA 2014.Google Scholar
- Jin Y. Evolutionary studies of Phrynocephalus (Agamidae) on the Qinghai–Xizang (Tibetan) Plateau. Ph. D. Thesis. Lanzhou University, Lanzhou (in Chinese with English abstract); 2008.Google Scholar
- Tang X, Xin Y, Wang H, Li W, Zhang Y, Liang S, et al. Metabolic characteristics and response to high altitude in Phrynocephalus erythrurus (Lacertilia: Agamidae), a lizard dwell at altitudes higher than any other living lizards in the world. PLoS One. 2013;8(8), e71976.View ArticlePubMed CentralPubMedGoogle Scholar
- He J, Xiu M, Tang X, Wang N, Xin Y, Li W, et al. Thermoregulatory and metabolic responses to hypoxia in the oviparous lizard, Phrynocephalus przewalskii. Comp Biochem Physiol A Mol Integr Physiol. 2013;165(2):207–13.View ArticlePubMedGoogle Scholar
- Yang W, Qi Y, Bi K, Fu J. Toward understanding the genetic basis of adaptation to high-elevation life in poikilothermic species: a comparative transcriptomic analysis of two ranid frogs. Rana chensinensis and R kukunoris BMC Genomics. 2012;13:588.View ArticleGoogle Scholar
- Chiari Y, Cahais V, Galtier N, Delsuc F. Phylogenomic analyses support the position of turtles as the sister group of birds and crocodiles (Archosauria). BMC Biol. 2012;10:65.View ArticlePubMed CentralPubMedGoogle Scholar
- Kumazawa Y. Mitochondrial genomes from major lizard families suggest their phylogenetic relationships and ancient radiations. Gene. 2007;388(1–2):19–26.View ArticlePubMedGoogle Scholar
- Vidal N, Hedges SB. The phylogeny of squamate reptiles (lizards, snakes, and amphisbaenians) inferred from nine nuclear protein-coding genes. C R Biol. 2005;328(10–11):1000–8.View ArticlePubMedGoogle Scholar
- Pyron RA, Burbrink FT, Wiens JJ. A phylogeny and revised classification of Squamata, including 4161 species of lizards and snakes. BMC Evol Biol. 2013;13:93.View ArticlePubMed CentralPubMedGoogle Scholar
- Hittinger CT, Johnston M, Tossberg JT, Rokas A. Leveraging skewed transcript abundance by RNA-Seq to increase the genomic depth of the tree of life. Proc Natl Acad Sci U S A. 2010;107(4):1476–81.View ArticlePubMed CentralPubMedGoogle Scholar
- Lin GH, Wang K, Deng XG, Nevo E, Zhao F, Su JP, et al. Transcriptome sequencing and phylogenomic resolution within Spalacidae (Rodentia). BMC Genomics. 2014;15:32.View ArticlePubMed CentralPubMedGoogle Scholar
- Nabholz B, Jarvis ED, Ellegren H. Obtaining mtDNA genomes from next-generation transcriptome sequencing: a case study on the basal Passerida (Aves: Passeriformes) phylogeny. Mol Phylogenet Evol. 2010;57(1):466–70.View ArticlePubMedGoogle Scholar
- Bickler PE, Buck LT. Hypoxia tolerance in reptiles, amphibians, and fishes: life with variable oxygen availability. Annu Rev Physiol. 2007;69:145–70.View ArticlePubMedGoogle Scholar
- Yee Koh M, Spivak-Kroizman TR, Powis G. HIF-1 regulation: not so easy come, easy go. Trends Biochem Sci. 2008;33(11):526–34.View ArticlePubMedGoogle Scholar
- He J, Xiu M, Tang X, Yue F, Wang N, Yang S, et al. The different mechanisms of hypoxic acclimatization and adaptation in Lizard Phrynocephalus vlangalii living on Qinghai-Tibet Plateau. J Exp Zool A Ecol Genet Physiol. 2013;319(3):117–23.View ArticlePubMedGoogle Scholar
- Svobodova AR, Galandakova A, Sianska J, Dolezal D, Lichnovska R, Ulrichova J, et al. DNA damage after acute exposure of mice skin to physiological doses of UVB and UVA light. Arch Dermatol Res. 2012;304(5):407–12.View ArticlePubMedGoogle Scholar
- Huey RB. Temperature, physiology, and the ecology of reptiles. Biology of the Reptilia. 1982;12:25–91.Google Scholar
- Pavlidis P, Jensen JD, Stephan W, Stamatakis A. A critical assessment of storytelling: gene ontology categories and the importance of validating genomic scans. Mol Biol Evol. 2012;29(10):3237–48.View ArticlePubMedGoogle Scholar
- Messier W, Stewart CB. Episodic adaptive evolution of primate lysozymes. Nature. 1997;385(6612):151–4.View ArticlePubMedGoogle Scholar
- Zhang J, Kumar S. Detection of convergent and parallel evolution at the amino acid sequence level. Mol Biol Evol. 1997;14(5):527–36.View ArticlePubMedGoogle Scholar
- Yang Z. Likelihood ratio tests for detecting positive selection and application to primate lysozyme evolution. Mol Biol Evol. 1998;15(5):568–73.View ArticlePubMedGoogle Scholar
- Fletcher W, Yang Z. The effect of insertions, deletions, and alignment errors on the branch-site test of positive selection. Mol Biol Evol. 2010;27(10):2257–67.View ArticlePubMedGoogle Scholar
- Zhang J, Nielsen R, Yang Z. Evaluation of an improved branch-site likelihood method for detecting positive selection at the molecular level. Mol Biol Evol. 2005;22(12):2472–9.View ArticlePubMedGoogle Scholar
- Wang W, Huang Y, Zhou Z, Tang R, Zhao W, Zeng L, et al. Identification and characterization of AGTRAP, a human homolog of murine Angiotensin II Receptor-Associated Protein (Agtrap). Int J Biochem Cell Biol. 2002;34(1):93–102.View ArticlePubMedGoogle Scholar
- Camelo Jr JS, Martins AR, Rosa E, Ramos SG, Hehre D, Bancalari E, et al. Angiotensin II type 1 receptor blockade partially attenuates hypoxia-induced pulmonary hypertension in newborn piglets: relationship with the nitrergic system. Braz J Med Biol Res. 2012;45(2):163–71.View ArticlePubMed CentralGoogle Scholar
- Gehrke SG, Riedel HD, Herrmann T, Hadaschik B, Bents K, Veltkamp C, et al. UbcH5A, a member of human E2 ubiquitin-conjugating enzymes, is closely related to SFT, a stimulator of iron transport, and is up-regulated in hereditary hemochromatosis. Blood. 2003;101(8):3288–93.View ArticlePubMedGoogle Scholar
- Vanni E, Gatherer D, Tong L, Everett RD, Boutell C. Functional characterization of residues required for the herpes simplex virus 1 E3 ubiquitin ligase ICP0 to interact with the cellular E2 ubiquitin-conjugating enzyme UBE2D1 (UbcH5a). J Virol. 2012;86(11):6323–33.View ArticlePubMed CentralPubMedGoogle Scholar
- Loublier S, Bayot A, Rak M, El-Khoury R, Benit P, Rustin P. The NDUFB6 subunit of the mitochondrial respiratory chain complex I is required for electron transfer activity: a proof of principle study on stable and controlled RNA interference in human cell lines. Biochem Biophys Res Commun. 2011;414(2):367–72.View ArticlePubMedGoogle Scholar
- Whitfield AJ, Barrett PH, van Bockxmeer FM, Burnett JR. Lipid disorders and mutations in the APOB gene. Clin Chem. 2004;50(10):1725–32.View ArticlePubMedGoogle Scholar
- Kadyrov FA, Dzantiev L, Constantin N, Modrich P. Endonucleolytic function of MutLalpha in human mismatch repair. Cell. 2006;126(2):297–308.View ArticlePubMedGoogle Scholar
- Merkle CJ, Karnitz LM, Henry-Sanchez JT, Chen J. Cloning and characterization of hCTF18, hCTF8, and hDCC1. Human homologs of a Saccharomyces cerevisiae complex involved in sister chromatid cohesion establishment J Biol Chem. 2003;278(32):30051–6.Google Scholar
- Liu Y, Zhou K, Zhang H, Shugart YY, Chen L, Xu Z, et al. Polymorphisms of LIG4 and XRCC4 involved in the NHEJ pathway interact to modify risk of glioma. Hum Mutat. 2008;29(3):381–9.View ArticlePubMedGoogle Scholar
- Chiruvella KK, Renard BM, Birkeland SR, Sunder S, Liang Z, Wilson TE. Yeast DNA ligase IV mutations reveal a nonhomologous end joining function of BRCT1 distinct from XRCC4/Lif1 binding. DNA Repair (Amst). 2014;24:37–45.View ArticleGoogle Scholar
- Gatei M, Young D, Cerosaletti KM, Desai-Mehta A, Spring K, Kozlov S, et al. ATM-dependent phosphorylation of nibrin in response to radiation exposure. Nat Genet. 2000;25(1):115–19.View ArticlePubMedGoogle Scholar
- Varon R, Vissinga C, Platzer M, Cerosaletti KM, Chrzanowska KH, Saar K, et al. Nibrin, a novel DNA double-strand break repair protein, is mutated in Nijmegen breakage syndrome. Cell. 1998;93(3):467–76.View ArticlePubMedGoogle Scholar
- Sevcik J, Falk M, Macurek L, Kleiblova P, Lhota F, Hojny J, et al. Expression of human BRCA1Delta17-19 alternative splicing variant with a truncated BRCT domain in MCF-7 cells results in impaired assembly of DNA repair complexes and aberrant DNA damage response. Cell Signal. 2013;25(5):1186–93.View ArticlePubMedGoogle Scholar
- Somyajit K, Basavaraju S, Scully R, Nagaraju G. ATM- and ATR-mediated phosphorylation of XRCC3 regulates DNA double-strand break-induced checkpoint activation and repair. Mol Cell Biol. 2013;33(9):1830–44.View ArticlePubMed CentralPubMedGoogle Scholar
- Murray JM, Carr AM. Smc5/6: a link between DNA repair and unidirectional replication? Nat Rev Mol Cell Biol. 2008;9(2):177–82.View ArticlePubMedGoogle Scholar
- Piwko W, Buser R, Peter M. Rescuing stalled replication forks: MMS22L-TONSL, a novel complex for DNA replication fork repair in human cells. Cell Cycle. 2011;10(11):1703–5.View ArticlePubMedGoogle Scholar
- Wu S, Shi Y, Mulligan P, Gay F, Landry J, Liu H, et al. A YY1-INO80 complex regulates genomic stability through homologous recombination-based repair. Nat Struct Mol Biol. 2007;14(12):1165–72.View ArticlePubMed CentralPubMedGoogle Scholar
- Stern DL. The genetic causes of convergent evolution. Nat Rev Genet. 2013;14(11):751–64.View ArticlePubMedGoogle Scholar
- Projecto-Garcia J, Natarajan C, Moriyama H, Weber RE, Fago A, Cheviron ZA, et al. Repeated elevational transitions in hemoglobin function during the evolution of Andean hummingbirds. Proc Natl Acad Sci U S A. 2013;110(51):20669–74.View ArticlePubMed CentralPubMedGoogle Scholar
- Yang W, Qi Y, Fu J. Exploring the genetic basis of adaptation to high elevations in reptiles: a comparative transcriptome analysis of two toad-headed agamas (genus Phrynocephalus). PLoS One. 2014;9(11), e112218.View ArticlePubMed CentralPubMedGoogle Scholar
- Gordon A, Hannon G: FASTQ/A short-reads pre-processing tools. [http://hannonlab.cshl.edu/fastx_toolkit]
- Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, et al. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat Biotechnol. 2011;29(7):644–52.View ArticlePubMed CentralPubMedGoogle Scholar
- Fu L, Niu B, Zhu Z, Wu S, Li W. CD-HIT: accelerated for clustering the next-generation sequencing data. Bioinformatics. 2012;28(23):3150–2.View ArticlePubMed CentralPubMedGoogle Scholar
- Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics. 2011;12:323.View ArticlePubMed CentralPubMedGoogle Scholar
- Altschul SF, Madden TL, Schaffer AA, Zhang J, Zhang Z, Miller W, et al. Gapped BLAST and PSI-BLAST: a new generation of protein database search programs. Nucleic Acids Res. 1997;25(17):3389–402.View ArticlePubMed CentralPubMedGoogle Scholar
- Conesa A, Gotz S, Garcia-Gomez JM, Terol J, Talon M, Robles M. Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research. Bioinformatics. 2005;21(18):3674–6.View ArticlePubMedGoogle Scholar
- Ye J, Fang L, Zheng H, Zhang Y, Chen J, Zhang Z. WEGO: a web tool for plotting GO annotations. Nucleic Acids Res. 2006;34((Web Server issue):W293–7.View ArticlePubMed CentralPubMedGoogle Scholar
- Moriya Y, Itoh M, Okuda S, Yoshizawa AC, Kanehisa M. KAAS: an automatic genome annotation and pathway reconstruction server. Nucleic Acids Res. 2007;35(Web Server issue):W182–5.View ArticlePubMed CentralPubMedGoogle Scholar
- Remm M, Storm CE, Sonnhammer EL. Automatic clustering of orthologs and in-paralogs from pairwise species comparisons. J Mol Biol. 2001;314(5):1041–52.View ArticlePubMedGoogle Scholar
- Alexeyenko A, Tamas I, Liu G, Sonnhammer EL. Automatic clustering of orthologs and inparalogs shared by multiple proteomes. Bioinformatics. 2006;22(14):e9–15.View ArticlePubMedGoogle Scholar
- Loytynoja A, Goldman N. An algorithm for progressive multiple alignment of sequences with insertions. Proc Natl Acad Sci U S A. 2005;102(30):10557–62.View ArticlePubMed CentralPubMedGoogle Scholar
- Darriba D, Taboada GL, Doallo R, Posada D. jModelTest 2: more models, new heuristics and parallel computing. Nat Methods. 2012;9(8):772.View ArticlePubMedGoogle Scholar
- Guindon S, Dufayard JF, Lefort V, Anisimova M, Hordijk W, Gascuel O. New algorithms and methods to estimate maximum-likelihood phylogenies: assessing the performance of PhyML 3.0. Syst Biol. 2010;59(3):307–21.View ArticlePubMedGoogle Scholar
- Yang Z. PAML 4: phylogenetic analysis by maximum likelihood. Mol Biol Evol. 2007;24(8):1586–91.View ArticlePubMedGoogle Scholar
- Drummond AJ, Rambaut A. BEAST: Bayesian evolutionary analysis by sampling trees. BMC Evol Biol. 2007;7:214.View ArticlePubMed CentralPubMedGoogle Scholar
- Yang Z, Nielsen R. Codon-substitution models for detecting molecular adaptation at individual sites along specific lineages. Mol Biol Evol. 2002;19(6):908–17.View ArticlePubMedGoogle Scholar
- Chimpanzee S, Analysis C. Initial sequence of the chimpanzee genome and comparison with the human genome. Nature. 2005;437(7055):69–87.View ArticleGoogle Scholar
- Yang Z, Wong WS, Nielsen R. Bayes empirical bayes inference of amino acid sites under positive selection. Mol Biol Evol. 2005;22(4):1107–18.View ArticlePubMedGoogle Scholar
- Punta M, Coggill PC, Eberhardt RY, Mistry J, Tate J, Boursnell C, et al. The Pfam protein families database. Nucleic Acids Res. 2012;40(Database issue):D290–301.View ArticlePubMed CentralPubMedGoogle Scholar
- Kelley LA, Sternberg MJ. Protein structure prediction on the Web: a case study using the Phyre server. Nat Protoc. 2009;4(3):363–71.View ArticlePubMedGoogle Scholar
- Pettersen EF, Goddard TD, Huang CC, Couch GS, Greenblatt DM, Meng EC, et al. UCSF Chimera–a visualization system for exploratory research and analysis. J Comput Chem. 2004;25(13):1605–12.View ArticlePubMedGoogle Scholar
- DeLano, WL. The PyMOL Molecular Graphics System. [www.pymol.org].
- Unni S, Huang Y, Hanson RM, Tobias M, Krishnan S, Li WW, et al. Web servers and services for electrostatics calculations with APBS and PDB2PQR. J Comput Chem. 2011;32(7):1488–91.View ArticlePubMed CentralPubMedGoogle Scholar
- Yang Y, Wang L, Han J, Tang X, Ma M, Wang K, Zhang X, Ren Q, Chen Q, Qiu Q. Data from: Comparative transcriptomic analysis revealed adaptation mechanism of Phrynocephalus erythrurus, the highest altitude Lizard living in the Qinghai-Tibet Plateau. BMC Evol Biol 2015, http://0-dx.doi.org.brum.beds.ac.uk/10.5061/dryad.f587b.
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly credited. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.