Skip to main content

Comparative transcriptome reveal the potential adaptive evolutionary genes in Andrias davidianus


To search the evidence of molecular evolution mechanism for aquatic and cave habitat in Andrias davidianus, the evolution analysis was carried out among several species transcriptome data. The transcriptome data of Notophthalmus viridescens, Xenopus tropicalis, Cynops pyrrhogaster, Hynobius chinensis and A. davidianus were obtained from the Genbank and reassembled except Xenopus tropicalis. The BLAST search of transcriptome data obtained 1244 single-copy orthologous genes among five species. A phylogenetic tree showed A. davidianus to have the closest relationship to H. chinensis. Fourteen positively selected genes were detected in A. davidianus and N. vridescens group and fifteen in A. davidianus and H. chinensis group. Five genes were shared in the both groups which involved in the immune system, suggesting that A. davidianus adaptation to an aquatic and cave environment required rapid evolution of the immune system compared to N. viridescens and H. chinensis.


Amphibians played an important role as a transitional group linking aquatic to terrestrial in the evolution of vertebrates [1]. To elucidate evolutionary history, the genome and mitochondrial DNA are traditionally used to estimate divergence time [2]. Transcriptome sequencing has become a viable alternative to provide rapid developing genomic resources in non-model organisms [3, 4]. Comparative transcriptome analysis is used to estimate the non-synonymous substitution (Ka) and synonymous substitution (Ks) rates to calculate the evolutionary rate [5, 6] and hence, to identify genes involved in environmental adaptation. Distribution of synonymous substitutions can be used to calculate the divergent time based on the coding sequence [2, 7].

The Chinese giant salamander Andrias davidianus is a typical urodele, and an important species both as a biological resource and with respect to its value as a living fossil [8]. The species was historically widespread in China, but environmental degradation and human killing have led to its severe decline in the wild. From 1980s, it is classified as endangered by the International Union for Conservation of Nature and Nature Resources. Because of its irreplaceable protection status and good taste, artificial propagation technology was studied and succeeded at the end of 1990s. Success of artificial propagation technology provided a value way to protect the wild resources. In wild, it is aquatic in all life stages and typically inhabits rocky crevices in banks of streams and lakes, as well as subterranean rivers. To identify genes possibly related to A. davidianus adaptation to its aquatic life history and to a cave habitat, transcriptome data of other amphibian species were obtained from GenBank, and comparative transcriptome analysis was carried out to detect genes positively selected for in evolution.


RNA extraction and sequencing

Total RNA was extracted from five ovaries and testes using Trizol reagent (Invitrogen, USA) according manufacturer’s instructions and treated with RNase-free DNase I (Takara, China) to remove the genomic DNA, respectively. After RNA quality and quantity test, RNA was broken into short fragment, and first-strand cDNA was synthesized, and then the sequencing adapter was added. The cDNA libraries were constructed and sequenced on the Illumina sequencing platform (Illumina HiSeq™ 2500). All raw reads, low quality sequences, and reads containing adaptor sequences were removed, and the clean reads were obtained.

Identification of orthologues genes and phylogenetic analysis

Two gonad transcriptome data (SRR3308418 and SRR3308420) of A. davidianus were provided by my lab. To expand data of A. davidianus, transcriptome data of skin (SRX729810) and spleen (SRX729743) were obtained from NCBI database and reassembled with the gonad transcriptome data. Transcriptome data of N. viridescens from heart, lens, brain, eye, liver, lung, spleen, kidney, testis, and ovary (ERR108189), C. pyrrhogaster lens and neural retina (SRR1051839), H. chinensis whole body (SRR1042328) and X. tropicalis from genome sequencing (GCA_000004195) were also obtained. The unigenes were reassembled from the downloaded raw reads, except for X. tropicalis. The numbers of unigenes for each species is given in Table 1. BLASTN software was used to align sequences, with the cutoff E-value set at 1e-7 [9]. OrthoMCL software was applied to classify the gene family [10]. Orthologous genes were obtained, and Venn diagrams were used to obtain the gene number [11]. The orthologous genes were used to construct the phylogenetic tree by the NJ method with 1000 bootstrap replications.

Table 1 Results of the assembly for each study species

Estimate of substitution rates among species

Form the orthologous gene, only one orthologous gene in other species was classed as single-copy orthologous by PERL package [12, 13]. The single-copy orthologous genes were identified to calculate the synonymous substitution rates (Ks) and non-synonymous rate (Ka). The amino acid sequences were aligned by Muscle software [14]. The aligned sequences were converted to corresponding nucleotide sequences. Synonymous substitution rates (Ks) and non-synonymous rates (Ka) were estimated between species pairs by sit model under Codeml program in PAML package [15]. The best threshold was set at 0.5 based on the Ka/Ks value according to previous reports [5, 6]. Value of two fold log-likelihood difference was used to perform a Chi-squared test and the difference of the parameter number was set as the degree in the Chi-squared. Positively selected sites were allowed when P was < 0.05 and posterior probability was > 0.95 [16]. A Ka/Ks value > 1 indicated strong positive selection, from 0.5 to 1 indicated weak positive selection, and a value < 0.1 indicated negative selection.


Orthologue identification and phylogenetic analysis

To identify the phylogenetic relationship among the species, large-scale transcriptome characterizations were carried out for N. viridescens, X. tropicalis, C. pyrrhogaster, H. chinensis, and A. davidianus, and transcriptome data were downloaded and reassembled (Table 1). Comparative analysis yielded 4279 gene families and 34,246 putative orthologous genes (Fig.1). To construct the phylogenetic tree with X. tropicalis as out-group, 1244 single-copy orthologous genes were identified. The phylogenetic tree showed A. davidianus to have the closest relationship to H. chinensis, with N. viridescens and C. pyrrhogaster clustered on one separate branch (Additional file 1: Figure S1).

Fig. 1
figure 1

Venn diagrams showing the unigenes for comparative transcriptomes. The superscript indicates the protein family and the subscript indicates the unigenes

Evolutionary profile of Andrias davidianus genes

We analyzed the evolutionary pattern of 1244 single-copy orthologous genes in A. davidianus, H. chinensis, and N. viridescens. Synonymous (Ks) and non-synonymous (Ka) substitutions per site were observed (Fig. 2). A majority of sequence pairs showed a Ka /Ks < 0.5, implying that these genes involved negative selection. Fifteen rapidly evolving sequences were identified with Ka/Ks > 0.5 between A. davidianus and H. chinensis, and 14 such sites were observed between A. davidianus and N. viridescens (Additional file 2: Table S1).

Fig. 2
figure 2

Ka/Ks ratio of 1244 single-copy orthologous genes. a. Ka/Ks distribution in Andrias davidianus and Hynobius chinensis, b. Ka/Ks distribution in Andrias davidianus and Notophthalmus viridescens. The solid line shows the threshold of Ka/Ks = 1, the dashed line marked the weak positive selection threshold of Ka/Ks = 0.5, and the short dashed line represented threshold of Ka/Ks = 0.1


Next-generation sequencing technology yielded a large number of sequences at the low cost and provides more sequences compared to traditional sequencing methods [17, 18]. Due to the cost and the throughput, genome-wide detection of the adaptive evolution gene was performed in many species by next-generation sequencing [17, 19]. Comparative phylogenetic analysis at the genome level improved the precision of evolutionary inference compared to single gene [20]. However, because of the large genome of the A. davidianus, evolutionary analysis by comparative genome was hard to carry out. Transcriptome sequencing was a valuable way to obtain large-scale sequences without reference genome [21, 22]. Phylogenetic analysis of transcriptome sequence data exhibited high supported tree topologies in many species [23, 24].

To elucidate the phylogenetic evolution of A. davidianus, comparative transcriptome analysis was conducted to construct the phylogenetic tree with X. tropicalis as out-group. To search adaptive gene for aquatic and cave life, molecular evolution was analyzed among the related species. Synonymous substitution rates (Ks) and non-synonymous substitution rates (Ka) were calculated according to the phylogenetic tree by PAML software [15, 25], with the optimal threshold for selecting the positively expressed sequence tag (EST) of 0.5 based on previous study [25]. Several positively selected genes were detected. Similar results were found in topmouth culter Erythroculter ilishaeformis and zebrafish Danio rerio, in which 38 candidate genes exhibited signs of positive selection with dN/dS ratios > 0.5 [6]. Five genes related to the immune system [26,27,28,29] [cystatin-like, oncostatin-M-specific receptor subunit beta isoform X1(OSMF), exonuclease, cell death regulator Aven, and centromere protein H] were detected in the A. davidianus/H. chinensis and A. davidianus/N. viridescens groups. Andrias davidianus is aquatic and inhabiting subterranean rivers and caves while N. viridescens and H. chinensis are mainly terrestrial and only special stage in water.

Aquatic and cave dwelling organisms generally encounter more bacteria than do terrestrial animals. Thus, the A. davidianus immune system should show more rapid mutations, as was confirmed in our investigation. Due to lack of full-length according to the transcriptome sequencing, many gene relevant to positive selection was omitted and Ka/Ks ratio was decreased from normal level [6]. Further study will be carried out to identify the genes under positive selection.



Expressed sequence tag


Non-synonymous substitution rates


Synonymous substitution rates


Ooncostatin-M-specific receptor subunit beta isoform X1


  1. Lowcock LA, Green DM, Sessions SK. Amphibian cytogenetics and evolution. Copeia. 1991;1992(2):588.

    Article  Google Scholar 

  2. Zhang P, Chen YQ, Zhou H, Liu YF, Wang XL, Papenfuss TJ, et al. Phylogeny, evolution, and biogeography of Asiatic salamanders (Hynobiidae). Proc Natl Acad Sci U S A. 2006;103(19):7360–5.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  3. Hiremath PJ, Farmer A, Cannon SB, Woodward J, Kudapa H, Tuteja R, et al. Large-scale transcriptome analysis in chickpea (Cicer arietinum L.), an orphan legume crop of the semi-arid tropics of Asia and Africa. Plant Biotechnol J. 2011;9(8):922–31.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  4. Fraser BA, Weadick CJ, Janowitz I, Rodd FH, Hughes KA. Sequencing and characterization of the guppy (Poecilia reticulata) transcriptome. BMC Genomics. 2011;12(1):202.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  5. Zhang L, Yan HF, Wu W, Yu H, Ge XJ. Comparative transcriptome analysis and marker development of two closely related primrose species (Primula poissonii and Primula wilsonii). BMC Genomics. 2013;14(1):329.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  6. Ren L, Tan XJ, Xiong YF, Xu K, Zhou Y, Zhong H, et al. Transcriptome analysis reveals positive selection on the divergent between topmouth culter and zebrafish. Gene. 2014;552(2):265–71.

    Article  CAS  PubMed  Google Scholar 

  7. Niu SH, Li ZX, Yuan HW, Chen XY, Li Y, Li W. Transcriptome characterisation of Pinus Tabuliformis and evolution of genes in the Pinus phylogeny. BMC Genomics. 2013;14(1):263.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  8. Gao KQ, Shubin NH. Earliest known crown-group salamanders. Nature. 2003;422(6930):424–8.

    Article  CAS  PubMed  Google Scholar 

  9. Zhang Z, Schwartz S, Wagner L, Miller W. A greedy algorithm for aligning DNA sequences. J Comput Biol. 2000;7(1-2):203–14.

    Article  CAS  PubMed  Google Scholar 

  10. Li L, Jr SC, Roos DS. OrthoMCL: identification of ortholog groups for eukaryotic genomes. Genome Res. 2003;13:2178–89.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  11. Chen H, Boutros PC. VennDiagram: a package for the generation of highly-customizable Venn and Euler diagrams in R. BMC Bioinformatics. 2011;12(1):35.

    Article  PubMed  PubMed Central  Google Scholar 

  12. Gabriel V, Francesc R, Gabriel C. A perl package and an alignment tool for phylogenetic networks. Bmc Bioinformatics. 2008;9(1):1–5.

    Article  Google Scholar 

  13. Wu F, Mueller LA, Crouzillat D, Pétiard V, Tanksley SD. Combining bioinformatics and phylogenetics to identify large sets of single-copy orthologous genes (COSII) for comparative, evolutionary and systematic studies: a test case in the euasterid plant clade. Genetics. 2006;174(3):1407.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  14. Edgar RC. MUSCLE: multiple sequence alignment with high accuracy and high throughput. Nucleic Acids Res. 2004;32(5):1792–7.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  15. Yang Z. PAML 4: Phylogenetic analysis by maximum likelihood. Mol Biol Evol. 2007;24(8):1586–91.

    Article  CAS  PubMed  Google Scholar 

  16. Yang Z, Bielawski JP. Statistical methods for detecting molecular adaptation. Trends Ecol Evol. 2000;15(12):496.

    Article  CAS  PubMed  Google Scholar 

  17. Tine M, Kuhl H, Gagnaire PA, Louro B, Desmarais E, Martins RS, et al. European sea bass genome and its variation provide insights into adaptation to euryhalinity and speciation. Nat Commun. 2014;5:5770.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  18. Liu H, Chen CL, Gao ZX, Min JM, Gu YM, Jian JB, et al. The draft genome of blunt snout bream (Megalobrama amblycephala) reveals the development of intermuscular bone and adaptation to herbivorous diet. Gigascience. 2017;6:1–13.

    Google Scholar 

  19. Shao CW, Bao BL, Xie ZY, Chen XY, Li B, Jia XD, et al. The genome and transcriptome of Japanese flounder provide insights into flatfish asymmetry. Nat Genet. 2017;49(1):119–24.

    Article  CAS  PubMed  Google Scholar 

  20. Clark AG, Eisen MB, Smith DR, Bergman CM, Oliver B, Markow TA, et al. Evolution of genes and genomes on the drosophila phylogeny. Nature. 2007;450(7167):203–18.

    Article  PubMed  Google Scholar 

  21. Che RB, Sun YN, Wang RX, Xu TJ. Transcriptomic analysis of endangered Chinese salamander: identification of immune, sex and reproduction-related genes and genetic markers. PLoS One. 2014;9:e87940.

    Article  PubMed  PubMed Central  Google Scholar 

  22. Cai JF, Yang W, Chen D, Zhang YZ, He Z, Zhang WM, et al. Transcriptomic analysis of the differentiating ovary of the protogynous ricefield eel Monopterus albus. BMC Genomics. 2017;18(1):573.

    Article  PubMed  PubMed Central  Google Scholar 

  23. Roeding F, Borner J, Kube M, Klages S, Reinhardt R, Burmester T. A 454 sequencing approach for large scale phylogenomic analysis of the common emperor scorpion (Pandinus imperator). Mol Phylogenet Evol. 2009;53(3):826–34.

    Article  CAS  PubMed  Google Scholar 

  24. Meusemann K, von Reumont BM, Simon S, Roeding F, Strauss S, Kuck P, et al. A phylogenomic approach to resolve the arthropod tree of life. Mol Biol Evol. 2010;27(11):2451–64.

    Article  CAS  PubMed  Google Scholar 

  25. Swanson WJ, Wong A, Wolfner MF, Aquadro CF. Evolutionary expressed sequence tag analysis of drosophila female reproductive tracts identifies genes subjected to positive selection. Genetics. 2004;168(3):1457–65.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  26. Li FM, Gai XM, Wang LL, Song LS, Zhang H, Qiu LM, et al. Identification and characterization of a cystatin gene from Chinese mitten crab Eriocheir sinensis. Fish Shellfish Immunol. 2010;29(3):521–9.

    Article  PubMed  Google Scholar 

  27. Kutuk O, Temel SG, Tolunay S, Basaga H. Aven blocks DNA damage-induced apoptosis by stabilising Bcl-xL. Eur J Cancer. 2010;46(13):2494–505.

    Article  CAS  PubMed  Google Scholar 

  28. He WL, Li YH, Yang DJ, Song W, Chen XL, Liu FK, et al. Combined evaluation of centromere protein H and Ki-67 as prognostic biomarker for patients with gastric carcinoma. Eur J Surg Oncol. 2013;39(2):141–9.

    Article  CAS  PubMed  Google Scholar 

  29. Larrea E, Echeverria I, Riezu-Boj JI, Aldabe R, Guembe L, Sola I, et al. Characterization of the CD40L/Oncostatin M/Oncostatin M receptor axis as an antiviral and immunostimulatory system disrupted in chronic HCV infection. J Hepatol. 2014;60(3):482–9.

    Article  CAS  PubMed  Google Scholar 

Download references


The authors would like to thank the Linen lab group for comments on data analysis.


This work was supported by National Nature Science Foundation of China (31502155). This work was supported by National Nonprofit Institute Research Grant (2017JBF0205). Key Laboratory of Freshwater Aquatic Biotechnology and Breeding, Ministry of Agriculture (KF-2017-06).

Availability of data and materials

The variation data reported in the paper was deposited in the GenBank.

Author information

Authors and Affiliations



QH analyze the data and write the manuscript. QW collected the data. YM and HT participated in the data analyzed and revised the manuscript. The study was conceived by HX. All the authors participated during the discussion and approved of its final manuscript.

Corresponding authors

Correspondence to Qiaomu Hu or Hanbing Xiao.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Additional files

Additional file 1:

Figure S1. Phylogenetic tree of selected species based on 1244 single-copy orthologous genes. (TIFF 212 kb)

Additional file 2:

Table S1. Orthologs gene under positive selection among species. (DOCX 18 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (, which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Hu, Q., Wang, Q., Meng, Y. et al. Comparative transcriptome reveal the potential adaptive evolutionary genes in Andrias davidianus. Hereditas 155, 18 (2018).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: