- Email: [email protected]
Divergent adaptation to Qinghai-Tibetan Plateau implicated from transciptome study of Gymnocypris dobula and Schizothorax nukiangensis
Biochemical Systematics and Ecology 71 (2017) 97e105
Contents lists available at ScienceDirect
Biochemical Systematics and Ecology journal homepage: www.elsevier.com/locate/biochemsyseco
Divergent adaptation to Qinghai-Tibetan Plateau implicated from transciptome study of Gymnocypris dobula and Schizothorax nukiangensis Mengchao Yu a, 1, Dongsheng Zhang a, 1, Peng Hu a, Sihua Peng a, Weiwen Li b, Shunping He c, Wanying Zhai a, Qianghua Xu b, **, Liangbiao Chen a, d, * a
Key Laboratory of Aquaculture Resources and Utilization, Ministry of Education, College of Fisheries and Life Science, Shanghai Ocean University, Shanghai, China Key Laboratory of Sustainable Exploitation of Oceanic Fisheries Resources, Ministry of Education, College of Marine Sciences, Shanghai Ocean University, Shanghai, China c The Key Laboratory of Aquatic Biodiversity and Conservation of Chinese Academy of Sciences, Institute of Hydrobiology, Chinese Academy of Sciences, Wuhan, China d Marine Biology and Biotechnology Laboratory, Qingdao National Laboratory for Marine Science and Technology, Qingdao, China b
a r t i c l e i n f o
a b s t r a c t
Article history: Received 29 December 2016 Received in revised form 2 February 2017 Accepted 4 February 2017
The Schizothoracine fishes are widely distributed in the Qinghai-Tibetan Plateau (QTP) area and its peripheral regions, which provide a prime example of adaptation in highland aquatic environments. Recent progresses have revealed various genetic adaptations of these fishes by comparing to distantly related lowerland species, however, comparative studies on closely-related species of different altitudes are still lacking. In this study, we sequenced and annotated a primitive Schizothoracine fish Schizothorax nukiangensis Tsao and a highly specialized one Gymnocypris dobula. We performed evolutionary analyses to investigate the candidate genes and signaling pathways involved QTP highland adaptation in both Schizothoracine fishes. Analysis of the 11,007 one-copy orthologs to the primitive cyprinid species, Danio rerio, revealed that both G. dobula and S. nukiangensis showed elevated evolutionary rates. A large number of genes related to hypoxia, including genes involved metabolic processes and cardiovascular system development, exhibited signatures of positive selection in both Schizothoracine fishes, but very few positively selected genes were found overlapping among these Schizothoracines. Our results indicated divergent genetic adaptation to highland environment for aquatic species living in QTP. © 2017 Elsevier Ltd. All rights reserved.
Keywords: Tibetan plateau Adaptive evolution Schizothoracine fish Hypoxia
1. Introduction The Qinghai-Tibet Plateau (QTP) is the highest and one of the biggest plateaus on earth, covering 2.5 106 square kilometers with an elevation of 3000e5000 m for most parts of the area. The QTP has been uplifting since approximately 45
* Corresponding author. ** Corresponding author. E-mail addresses: [email protected] (Q. Xu), [email protected] (L. Chen). 1 These authors contributed equally. http://dx.doi.org/10.1016/j.bse.2017.02.003 0305-1978/© 2017 Elsevier Ltd. All rights reserved.
98
M. Yu et al. / Biochemical Systematics and Ecology 71 (2017) 97e105
million years ago, resulted from collision of the India plate and the Eurasia plate (Li and Fang, 1999; Favre et al., 2014). The uplifting dramatically changed the environment conditions from an originally humid and warm climate to currently a dry and cold one (Wu et al., 2008). Characterized by factors including hypoxia, cold and strong ultraviolet radiation, the QTP environment posed harsh challenges to the endemic animals (Scheinfeldt and Tishkoff, 2010). Using comparative genomics and transcriptomics analysis, recent studies have identified various genes and signaling pathways that may be responsible for highland adaptation in both terrestrial and aquatic vertebrates. For example, genetic modifications to metabolism and cardiovascular system development are frequently found in the highland adapted animals (Qiu et al., 2012; Qu et al., 2013; Gou et al., 2014). The family Cyprinidae is the most diverse clade in freshwater fishes (Chen and Mayden, 2009). The Schizothoracine fishes (Teleostei: Cyprinidae) dominate the lakes and rivers throughout the Tibetan Plateau and its peripheral regions. More than 70 species in twelve genera in the subfamily Schizothoracinae are endemic to the Tibetan Plateau (He and Chen, 2006). According to the degree of specialization of their morphological traits, the Schizothoracine fishes are divided into three grades: Primitive, specialized, and highly specialized, with Primitive ones live in the peripheral regions at around 1000 m above sea level (a.s.l) and the other two groups are intermingled with each other in the central part of QTP at above 3000 m a.s.l (He and Chen, 2006; Qi et al., 2012). Recent transcriptome studies on Schizothoracines and other endemic Tibetan fishes have identified multiple biological processes and genes involved in highland adaptation, including genes that participate in metabolic processes and responses to hypoxia (Yang et al., 2014; Ma et al., 2015). Despite great progresses made in these investigations, most of current studies are limited to single highland species with comparisons made to distantly related teleosts. Such comparisons are limited in power to decipher the detailed molecular adaptations associated with Schizothoracine expansion in QTP. Therefore, genomewide investigations of adaptive signals among Schizothoracine species of differently graded specialization are highly demanded to further reveal the genetic mechanisms of adaptation in these fishes. The objective of this study was to carry out a comparative genome-wide screen for genes that might be involved in highland adaptation from species living at different altitudes. To achieve these aims, we generated transcriptomes of a primitive Schizothoracine fish Schizothorax nukiangensis Tsao, which lives at about 1,000 m a.s.l, and a highly specialized one Gymnocypris dobula which lives at 3000e4000 m a.s.l, and we performed evolutionary analyses on these data. 2. Materials and methods 2.1. Ethics statement This study was approved by the Ethics Committee for the Use of Animal Subjects of Shanghai Ocean University. 2.2. Sample collection, cDNA library construction and illumina sequencing G. dobula was collected from Lake Duoqing, situated at eastern QTP (28 03.370, 8917.830 ) with an altitude of 4509 m a.s.l, and a temperature at 11.0 C, the concentration of dissolved oxygen in the water was 1.9 mg/L; while S. nukiangensis was collected from Nujiang, Yunnan, situated on the eastern margin of QTP (25 41.24’, 98 53.22), with an altitude of 1201 m a.s.l and at a temperature of 16.0 C, the concentration of dissolved oxygen in the water was 7.8 mg/L. The collection loci were labelled in Fig. 1. For both species, five individuals were captured from the same sampling locations, and fishes with similar size were selected for tissue dissection. The fish were live trapped, anesthetized. Tissues including brain, gill, head kidney, kidney, liver, muscle, skin, spleen were quickly biopsied and placed in RNAlater (QIAGEN), and stored at 80 C on arrival at the laboratory till RNA extraction. Total RNAs from each tissue sample were extracted using TRIzol reagent (Invitrogen Corp., Carlsbad, CA). Library constructions from G. dobula and S. nukiangensis were made using Illumina Hiseq1500 RNA sample preparation kits (Illumina, San Diego, CA) following manufacturer's instructions. RNAseq libraries were quantified using an Agilent 2100 Bioanalyzer (Agilent Technologies). The libraries were then sequenced on an Illumina HiSeq 1500 platform (Illumina, Inc.) with 100 bps paired-end reads. 2.3. De novo assembly and annotation of transcripts Raw reads were cleaned by removing adapter sequences using sequence pre-processing tool Trimmomatic. Reads with Phred score 25 were kept for sequence assembly. The transcriptomes of the tissues were pooled to generate an improved de novo assembly using the Trinity package, and open reading frames (ORFs) were predicted using the transdecoder with a minimum length of 50 amino acids (Haas et al., 2013). The assembled unigenes were first annotated by searching for homologous sequences against National Center for Biotechnology Information (NCBI) nonredundant protein (Nr) database and then were searched against the KOG database using BLASTP (evalue<1e-5). We also used the KOBAS software to annotate the unigenes with Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways (Xie et al., 2011). As zebrafish is currently the most closely related species to the Schizothoracines which has a completed and well annotated genome, we further annotated the unigenes of the two Schizothoracine fishes by mapping them to zebrafish proteins using BLASTP bidirectional best hit (BBH) method. The protein dataset of zebrafish (Danio rerio) was downloaded from
M. Yu et al. / Biochemical Systematics and Ecology 71 (2017) 97e105
99
Fig. 1. Geographical locations of sample collection.
Ensembl online database (ftp://ftp.ensembl.org/pub/release-84/fasta/danio_rerio/pep/) and GO annotation data for zebrafish genome was downloaded from Ensembl BioMart.
2.4. Ortholog determination and sequence alignment Orthologs among the two Schizothoracine fishes and zebrafish were determined using OrthoMCL (version_2.0.9) with default settings (Li et al., 2003). Ortholog groups of single copy genes (one-copy orthologs) that were shared by all three species were retained for further analysis. These one-copy orthologs were aligned at the protein level using the Prank pro€ytynoja and Goldman, 2008) and the corresponding codon alignments were produced using PAL2NAL (Suyama et al., gram (Lo 2006) with the unreliable alignments removed using the Gblocks program (Castresana, 2000). The trimmed alignments with lengths longer than 150 bps were applied for further evolutionary analysis.
2.5. Calculation of evolutionary rates To estimate the evolutionary rates for the two Schizothoracine fishes and zebrafish, we calculated ratios of nonsynonymous (dN) to synonymous (dS) substitution rate (u ¼ dN/dS) using the free-ratio model in codeml program of the PAML package (Yang, 1997). The dN/dS ratios were estimated with both individual one-copy orthologs and concatenated alignments from all one-copy orthologs. The values of dN, dS, and dN/dS ratios were retrieved from output files of the program. Genes were discarded if N*dN or S*dS < 1 or dS >1 according to the method described in a previous study (Goodman et al., 2010). We drew the boxplot for the distributions of dN/dS ratios. We also estimated lineage-specific mean values of dN/ dS using concatenated alignments from all one-copy orthologs. To investigate what molecular processes and functions are under rapid evolution, we estimated evolutionary rates of genes based on their Gene Ontology (GO) terms. Genes were categorized by the GO terms they were annotated, and only GO categories with more than 30 genes were included in the analyses. Lineage-specific mean dN/dS values were estimated by the concatenated alignments of all one-copy orthologs in each GO category, and the average dN/dS values for each GO category were calculated with free-ratio model in codeml program, GO terms experiencing relatively accelerated evolution were identified using a binomial test following the method applied in a previous study (Sequencing and Consortium, 2005).
2.6. Identification of positively selected genes (PSGs) Since the Branch-site model in codeml program has the advantage of detecting positive selection that affects only a few sites on a pre-specified (foreground) branch of the species tree (Zhang et al., 2005), this model was used to detect positive selection in the single-copy orthologous genes. We used likelihood ratio test (LTR) to compare the alternative hypothesis of positive selection on the foreground branch to a null hypothesis with no positive selection on the foreground branches for each orthologous gene. Positively selected genes were inferred only if their P-values were less than 0.05 after FDR adjusting (Benjamini and Hochberg, 1995). GO and KEGG enrichment analyses of PSGs were performed using the Fisher's exact test (Mehta and Patel, 1983).
100
M. Yu et al. / Biochemical Systematics and Ecology 71 (2017) 97e105
3. Results 3.1. De novo assembly and annotation of transcriptomes RNA-sequencing generated 188.70 million raw reads in about 94 Gb of data from S. nukiangensis, and 198.82 million raw reads in about 98 Gb of data of G. dobula. The original sequencing data were deposited in the NCBI and can be accessed in the Short Read Archive (SRA) under the accession numbers (SRP094994 and SRP095655). De novo assembly of these reads generated 91,854 unigenes for G. dobula and 171,555 unigenes for S. nukiangensis. The N50s of the two datasets are 1665 bp and 1661 bp, respectively. A total of 97,648 ORFs in S. nukiangesis and 51,782 ORFs in G. dobula were predicted (Supplementary Table 1). About 94.5% (92,307) and 94.7% (49,302) ORFs in S. nukiangensis and G. dobula had significant matches to protein sequences in the NR database, respectively. There are 87% and 58.2% of the BLASTP top-hits matched zebrafish protein sequences in S. nukiangensis and G. dobula, respectively (Fig. 2). As the Schizothoracine fishes and zebrafish belong to the same family Cyprinidae, the transcripts were further annotated against the protein database of zebrafish (Danio.rerio) using BBH method. In total, 19,520 and 21,043 protein-coding genes of G. dobula and S. nukiangensis mapped to D. rerio, respectively. We also did functional annotation using KOG. A total of 26,706 ORFs (27.35% of 97,648 ORFs) in S. nukiangesis and 35,632 ORFs (68.81% of 51,782 ORFs) in G. dobula were assigned to 26 eukaryotic orthologous groups. The top two categories were “Signal transduction” and “General function prediction only”. 3.2. Identification of orthologous genes and dN/dS analyses A total of 15,623 orthologous groups (also known as orthogroups) were detected for G. dobula, S. nukiangensis and D. rerio. Among these, 13,330 orthogroups were shared by all species, including 11,141 one-copy ortholog groups. After sequence alignment and trimming, 11,007 one-copy orthologous genes were retained for estimating the evolutionary constraints acting on G. dobula, S. nukiangensis and D. rerio. The lineage-specific evolutionary rates were estimated by measuring dN/dS ratios using the free-ratio model in PAML, which allows a separate dN/dS ratio for each branch. The dN/dS ratios for each branch were calculated using a concatenated alignment of 11,007 one-copy orthologs, and the evolutionary rates for the species were G. dobula (0.3101) > S. nukiangensis (0.2167) > D. rerio (0.1456) (Fig. 3A). Furthermore, dN/dS ratios for each ortholog were evaluated at gene level using the free ratio model. There were 7980 genes included in the analysis as they met the criteria N*dN or S*dS > 1 and dS < 1. Both G. dobula and S. nukiangensis had significantly higher dN/dS ratios than D. rerio (Wilcoxon rank sum test, P-values <2.2e-16). There was no significant difference between dN/dS ratios of G. dobula and S. nukiangensis (Wilcoxon rank sum test, P > 0.05) (Fig. 3B). 3.3. Accelerated evolution of the G. dobula lineages and S. nukiangensis lineages To identify lineage-specific accelerated GO categories in G. dobula and S. nukiangensis, the average dN/dS ratios for each GO category were calculated for each branch using concatenated alignment of genes belong to the GO category. We analyzed 1276 GO categories with more than 30 one-copy orthologs (Supplementary Table 2). When comparing S. nukiangensis to D. rerio, there were much more GO categories with significantly higher average dN/dS ratios in S. nukiangensis than in zebrafish
Fig. 2. Top five hit species distributions for S. nukiangensis and G. dobula based on BLASTP against Nr databases.
M. Yu et al. / Biochemical Systematics and Ecology 71 (2017) 97e105
101
Fig. 3. The average dN/dS ratios for each species estimated using concatenated alignments for all one-copy orthologs (A) and each ortholog (B).
lineage (1077 vs. 14); and there were 1008 GO categories with significantly higher average dN/dS ratios in G. dobula while no GO categories showed significantly higher dN/dS ratios in zebrafish. Furthermore, when comparing G. dobula to S. nukiangensis, we found that 936 GO categories had significantly higher dN/dS ratios in G. dobula, including “immune response”, “metabolic process” and “DNA repair”. While only 16 GO categories in S. nukiangensis: including “hemostasis” and “blood coagulation” (Fig. 4). 3.4. Positively selected genes of G. dobula and S. nukiangensis The improved branch-site model in PAML was used to detect positively selected genes (PSGs) among the 11,007 one-copy orthologous genes (Zhang et al., 2005). We identified 315 and 239 PSGs in G. dobula and S. nukiangensis by setting respective branches as the foreground branch. There were 18 PSGs overlapping between the two species, including nup214, which is involved in regulation of glucokinase, and bcam, a member of the immunoglobulin superfamily (Supplementary Table 3). When comparing to a previous study on a specialized Schizothoracine fish Gymnodiptychus pachycheilus (Yang et al., 2014), we found that nup214 was the only one in common for the three species and that there were two PSGs overlapping between G. dobula and G. pachycheilus. We also conducted GO and KEGG functional classification for these PSGs in these two species (Fig. 5). PSGs in G. dobula and S. nukiangensis covered similar GO categories, including “Metabolic process” (84 and 67 PSGs, respectively), “response to stress” (18 and 12), “immune system process” (12 and 10), “cardiovascular system development” (9 and 7) and “cell adhesion” (9 and 6). Some GO categories were specifically enriched in G. dobula or S. nukiangensis. For example, “receptor activity” and “cell adhesion” were enriched in G. dobula, while “ion binding” and “erythrocyte differentiation” were specifically enriched in S. nukiangensis. Among KEGG classification, the class ‘Environmental Information Processing’ contained the largest number of PSGs, for example: there were 26 and 10 PSGs assigned to its subclass “signal transduction” for G. dobula and S. nukiangensis respectively. Another predominant KEGG class was “immune system”, which included 11 and 9 PSGs, respectively. Some genes related to DNA damage, for example, ddb2 in S. nukiangensis, and msh5 in G. dobula, were identified in our study, which may relate to high UV exposure in Schizothoracines. Since hypoxia is one of the most challenging environmental factors for QTP endemic species, we further examined hypoxia-related PSGs by comparing to hypoxiaDB, which is a database of hypoxia-regulated proteins extracted from previous publications (Khurana et al., 2013; Zhang et al., 2014). We identified 23 and 25 PSGs from G. dobula and S. nukiangensis in the hypoxiaDB, respectively. There were no common hypoxia-related PSGs in these two species (Supplementary Table 4). As key regulator in hypoxia response, HIF-1 signaling pathway is believed to play pivotal roles in highland adaptation (Qiu et al., 2012; Qu et al., 2013; Gou et al., 2014). In this study, we identified four PSGs among the 20 genes annotated in the KEGG HIF-1 signaling pathway, including three PSGs in G. dobula and one in S. nukiangesis (Fig. 6). 4. Discussion In this study, we performed multi-tissue transcriptome sequencing of two Schizothoracine fishes and found many genes and functions might be involved in highland adaptation for these fishes. As increased evolutionary rates of genes are regarded as signals of adaptive evolution at molecular level, we investigated the evolutionary rates of Schizothoracines by estimating dN/dS ratios. We found that dN/dS ratios of both Schizothoracines were significantly higher than that of zebrafish. Our results are consistent with some previous studies which found elevated dN/dS ratios in species endemic to Tibet comparing to those living at lower altitudes (Yang et al., 2014; Ma et al., 2015). These results indicated that Schizothoracines have undergone rapid evolution to adapt to the extreme environment of the Tibetan Plateau. On the other hand, as the chromosome numbers in the Schizothoracine fishes are frequently duplicated (Wang et al.,
102
M. Yu et al. / Biochemical Systematics and Ecology 71 (2017) 97e105
Fig. 4. Comparisons of average dN/dS ratios for each GO category with more than 30 orthologs. (A: G. dobula vs. D. rerio; B: S. nukiangensis vs. D. rerio, C: G. dobula vs. S. nukiangensis). GO categories with significantly higher dN/dS ratios for species on Y-axis (red) and X-axis (blue) are highlighted. Grey dots stand for GO categories without significant differences in dN/dS ratios between the two species. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
2016), the increased dN/dS ratios in the Schizothoracines may thus associate with the specific genome duplication events in the subfamily. We also compared dN/dS ratios between the two Schizothoracines. There is no significant difference between them when the comparison was made at individual gene level, but dN/dS ratios estimated by the concatenated alignments of orthologs indicated that G. dobula has increased average dN/dS ratio than S. nukiangensis. The incongruence of the results above may come from the datasets used for the estimation. As many as 3027 genes were not included in the gene level analysis, and thus the mean evolutionary rates estimated from individual genes may be biased comparing with the ones from concatenated alignments. Investigation on evolutionary rates of GO categories also suggested accelerated evolution in Schizothoracines. Most GO categories have higher dN/dS ratios in both Schizothoracines than in zebrafish. In addition, G. dobula have much more GO categories with elevated dN/dS ratios than S. nukiangensis, including many GO categories that were related to the adaptation of high altitude. These results may indicate that extensive genetic adaptation to highland occurred in these Schizothoracines, and that species living at higher altitude in QTP are under stronger selection pressure. The most striking feature of the QTP highland is hypoxia. Adaptation to hypoxia requires decreased oxygen consumption and/or increased oxygen intake. The former is mostly related to metabolism, while the latter involves cardiovascular system development and hemopoiesis. In this study, we identified as many as 84 and 67 PSGs that belong to the GO category
M. Yu et al. / Biochemical Systematics and Ecology 71 (2017) 97e105
103
Fig. 5. A: Distribution of GO classification of PSGs. B: Distribution of KEGG classification of PSGs. (A) Cellular Processes; (B) Environmental Information Processing; (C) Genetic Information Processing; (D) Metabolism; (E) Organismal Systems.
Fig. 6. HIF-1 signaling pathway in KEGG populated with one-copy ortholog genes annotated in this study. The annotated genes are shown in colored boxes: red ones are PSGs identified in G. dobula, blue ones in S. nukiangensis. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
“metabolic process” in G. dobula and S. nukiangensis, respectively. KEGG pathway inference of the PSGs again suggested that various metabolic processes in both Schizothoracines might have been evolved under selection pressure, such as lipid
104
M. Yu et al. / Biochemical Systematics and Ecology 71 (2017) 97e105
metabolism, energy metabolism, carbohydrate metabolism, amino acid metabolism and nucleotide metabolism. Nine and seven PSGs that associated with the GO category “cardiovascular system development” were identified in G. dobula and S. nukiangensis respectively. Three genes involved in “erythrocyte differentiation” was found under positive selection in S. nukiangensis, including erythropoietin, a principal regulator in erythropoiesis, suggesting possible developmental adaptation of this species in hemopoiesis. Furthermore, there were 23 and 25 PSGs listed in hypoxiaDB for G. dobula and S. nukiangensis, respectively, suggesting the two scizothoracines in this study are possibly under selection pressure caused by hypoxia. Hypoxia-inducible factor (HIF) is a master switch that mediates major changes in gene expression under low oxygen. Previous studies have shown that some genes involved in the HIF-1 signaling pathway are under positive selection in scizothoracine fishes living in QTP (Guan et al., 2014; Yang et al., 2014; Shao et al., 2015). In this study, three PSGs were identified to be in the HIF-1 signaling pathway in G. dobula, suggesting hypoxia adaptation in this species. Very few PSGs were found overlapping among different Schizothoracinae species. nup214, which is involved in regulation of glucokinase, was the only common PSG identified in these two Schizothoracines investigated in this study and another specialized schizothoracine fish, G. przewalskii (Yang et al., 2014), indicating divergent evolution during highland adaptation for Schizothoracines. Acknowledgement This work was supported in part by the National Natural Science Foundation of China (31130049), the Research Project of the Chinese Ministry of Education (no. 213013A), the National Natural Science Foundation of China (grant no. 31572598), the “Shuguang Program” supported by the Shanghai Education Development Foundation and the Shanghai Municipal Education Commission (grant no. 13SG51), and the Shanghai Municipal Project for First-Class Discipline of Fishery to Shanghai Ocean University. Appendix A. Supplementary data Supplementary data related to this article can be found at http://dx.doi.org/10.1016/j.bse.2017.02.003. References Benjamini, Y., Hochberg, Y., 1995. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J. R. Stat. Soc. B 57, 289e300. http://dx.doi.org/10.2307/2346101. Castresana, J., 2000. Selection of conserved blocks from multiple alignments for their use in phylogenetic analysis. Mol. Biol. Evol. 17, 540e552. http://dx. doi.org/10.1093/oxfordjournals.molbev.a026334. Chen, W.J., Mayden, R.L., 2009. Molecular systematics of the Cyprinoidea (Teleostei: cypriniformes), the world's largest clade of freshwater fishes: further evidence from six nuclear genes. Mol. Phylogenet. Evol. 52, 544e549. http://dx.doi.org/10.1016/j.ympev.2009.01.006. €ckert, M., Pauls, S.U., Ja €hnig, S.C., Uhl, D., Michalak, I., Muellner-Riehl, A.N., 2014. The role of the uplift of the Qinghai-Tibetan Plateau for the Favre, A., Pa evolution of Tibetan biotas. Biol. Rev. 90 http://dx.doi.org/10.1111/brv.12107 na/na. Goodman, M., Sterner, K.N., Islam, M., Sherwood, C.C., Hof, P.R., Hou, Z., Lipovich, L., Jia, H., Grossman, L.I., Derek, E., Pease, L.R., Goodman, M., Sterner, K.N., Islam, M., Uddin, M., Sherwood, C.C., Hof, P.R., 2010. Phylogenomic analyses reveal convergent patterns of adaptive evolution in elephant and human ancestries. Correction for Goodman et al. Proc. Natl. Acad. Sci. 107, 8498. http://dx.doi.org/10.1073/pnas.1003435107. Gou, X., Wang, Z., Li, N., Qiu, F., Xu, Z., Yan, D., Yang, S., Jia Konga, X., Wei, Z., Lu, S., Lian, L., Wu, C., Wang, X., Li, G., Teng, M., Jiang, Q., Zhao, X., Yang, J., Liu, B., Wei, D., Li, H., Yang, J., Yan, Y., Zhao, G., Dong, X., Li, M., Deng, W., Leng, J., Wei, C., Wang, C., Mao, H., Zhang, H., Ding, G., Li, Y., 2014. Whole-genome sequencing of six dog breeds from continuous altitudes reveals adaptation to high-altitude hypoxia. Genome Res. 24, 1308e1315. http://dx.doi.org/10. 1101/gr.171876.113. Guan, L., Chi, W., Xiao, W., Chen, L., He, S., 2014. Analysis of hypoxia-inducible factor alpha polyploidization reveals adaptation to Tibetan plateau in the evolution of schizothoracine fish. BMC Evol. Biol. 14, 192. http://dx.doi.org/10.1186/s12862-014-0192-1. Haas, B.J., Papanicolaou, A., Yassour, M., Grabherr, M., Blood, P.D., Bowden, J., Couger, M.B., Eccles, D., Li, B., Lieber, M., Macmanes, M.D., Ott, M., Orvis, J., Pochet, N., Strozzi, F., Weeks, N., Westerman, R., William, T., Dewey, C.N., Henschel, R., Leduc, R.D., Friedman, N., Regev, A., 2013. De novo transcript sequence reconstruction from RNA-seq using the Trinity platform for reference generation and analysis. http://dx.doi.org/10.1038/nprot.2013.084. He, D., Chen, Y., 2006. Biogeography and molecular phylogeny of the genus Schizothorax (Teleostei: Cyprinidae) in China inferred from cytochrome b sequences. J. Biogeogr. 33, 1448e1460. http://dx.doi.org/10.1111/j.1365-2699.2006.01510.x. Khurana, P., Sugadev, R., Jain, J., Singh, S.B., 2013. HypoxiaDB: a database of hypoxia-regulated proteins. Database 2013, 1e12. http://dx.doi.org/10.1093/ database/bat074. Li, J., Fang, X., 1999. Uplift of the Tibetan Plateau and environmental changes. Chin. Sci. Bull. 44, 2117e2124. http://dx.doi.org/10.1007/BF03182692. Li, L., Stoeckert Jr., C.J., Roos, D.S., 2003. OrthoMCL: Identification of Ortholog Groups for Eukaryotic Genomes, pp. 2178e2189. http://dx.doi.org/10.1101/gr. 1224503.candidates. €ytynoja, A., Goldman, N., 2008. Phylogeny-aware gap placement prevents errors in sequence alignment and evolutionary analysis. Science 320, Lo 1632e1635. http://dx.doi.org/10.1126/science.1158395. Ma, X., Dai, W., Kang, J., Yang, L., He, S., 2015. Comprehensive Transcriptome Analysis of Six Catfish Species from an Altitude Gradient Reveals Adaptive Evolution in Tibetan Fishes, pp. 1e33. http://dx.doi.org/10.1534/g3.115.024448. Mehta, C.R., Patel, N.R., 1983. A network algorithm for performing Fisher's exact test in r c contingency tables. J. Am. Stat. Assoc. 78, 427e434. http://dx. doi.org/10.2307/2288652. Qi, D., Chao, Y., Guo, S., Zhao, L., Li, T., Wei, F., Zhao, X., 2012. Convergent, parallel and correlated evolution of trophic morphologies in the subfamily schizothoracinae from the qinghai-tibetan plateau. PLoS One 7, 1e10. http://dx.doi.org/10.1371/journal.pone.0034070. Qiu, Q., Zhang, G., Ma, T., Qian, W., Wang, J., Ye, Z., Cao, C., Hu, Q., Kim, J., Larkin, D.M., Auvil, L., Capitanu, B., Ma, J., Lewin, H.A., Qian, X., Lang, Y., Zhou, R., Wang, L., Wang, K., Xia, J., Liao, S., Pan, S., Lu, X., Hou, H., Wang, Y., Zang, X., Yin, Y., Ma, H., Zhang, J., Wang, Z., Zhang, Y., Zhang, D., Yonezawa, T., Hasegawa, M., Zhong, Y., Liu, W., Zhang, Y., Huang, Z., Zhang, S., Long, R., Yang, H., Wang, J., Lenstra, J.A., Cooper, D.N., Wu, Y., Wang, J., Shi, P., Liu, J., 2012. The yak genome and adaptation to life at high altitude. Nat. Genet. 44, 946e949. http://dx.doi.org/10.1038/ng.2343. Qu, Y., Zhao, H., Han, N., Zhou, G., Song, G., Gao, B., Tian, S., Zhang, J., Zhang, R., Meng, X., Zhang, Y., Zhang, Y., Zhu, X., Wang, W., Lambert, D., Ericson, P.G.P., Subramanian, S., Yeung, C., Zhu, H., Jiang, Z., Li, R., Lei, F., 2013. Ground tit genome reveals avian adaptation to living at high altitudes in the Tibetan plateau. Nat. Commun. 4, 2071. http://dx.doi.org/10.1038/ncomms3071.
M. Yu et al. / Biochemical Systematics and Ecology 71 (2017) 97e105
105
Scheinfeldt, L.B., Tishkoff, S.A., 2010. Living the high life: high-altitude adaptation. Genome Biol. 11, 133. http://dx.doi.org/10.1186/gb-2010-11-9-133yrgb2010-11-9-133 [pii]. Sequencing, T.C., Consortium, A., 2005. Initial sequence of the chimpanzee genome and comparison with the human genome. Nature 437, 69e87. http://dx. doi.org/10.1038/nature04072. Shao, Y., Li, J.-X., Ge, R.-L., Zhong, L., Irwin, D.M., Murphy, R.W., Zhang, Y.-P., 2015. Genetic adaptations of the plateau zokor in high-elevation burrows. Sci. Rep. 5, 17262. http://dx.doi.org/10.1038/srep17262. Suyama, M., Torrents, D., Bork, P., 2006. PAL2NAL: robust conversion of protein sequence alignments into the corresponding codon alignments. Nucleic Acids Res. 34 http://dx.doi.org/10.1093/nar/gkl315. Wang, X., Gan, X., Li, J., Chen, Y., He, S., 2016. Cyprininae phylogeny revealed independent origins of the Tibetan Plateau endemic polyploid cyprinids and their diversifications related to the neogene uplift of the plateau. Sci. China Life Sci. 59, 1149e1165. http://dx.doi.org/10.1007/s11427-016-0007-7. Wu, Zhenhan, Barosh, P.J., Zhonghai, W., Daogong, H., Xun, Z., Peisheng, Y., 2008. Vast early Miocene lakes of the central Tibetan plateau. Bull. Geol. Soc. Am. 120, 1326e1337. http://dx.doi.org/10.1130/B26043.1. Xie, C., Mao, X., Huang, J., Ding, Y., Wu, J., Dong, S., Kong, L., Gao, G., Li, C.Y., Wei, L., 2011. KOBAS 2.0: a web server for annotation and identification of enriched pathways and diseases. Nucleic Acids Res. 39 http://dx.doi.org/10.1093/nar/gkr483. Yang, Z., 1997. User Guide PAML: Phylogenetic Analysis by Maximum Likelihood, vol. 6, pp. 1e70. http://dx.doi.org/10.1093/molbev/msm088. Yang, L., Wang, Y., Zhang, Z., He, S., 2014. Comprehensive transcriptome analysis reveals accelerated genic evolution in a Tibet fish, Gymnodiptychus pachycheilus. Genome Biol. Evol. 7, 251e261. http://dx.doi.org/10.1093/gbe/evu279. Zhang, J., Nielsen, R., Yang, Z., 2005. Evaluation of an improved branch-site likelihood method for detecting positive selection at the molecular level. Mol. Biol. Evol. 22, 2472e2479. http://dx.doi.org/10.1093/molbev/msi237. Zhang, W., Fan, Z., Han, E., Hou, R., Zhang, L., Galaverni, M., Huang, J., Liu, H., Silva, P., Li, P., Pollinger, J.P., Du, L., Zhang, X.Y., Yue, B., Wayne, R.K., Zhang, Z., 2014. Hypoxia adaptations in the grey wolf (Canis lupus chanco) from qinghai-tibet plateau. PLoS Genet. 10 http://dx.doi.org/10.1371/journal.pgen. 1004466.
Contents lists available at ScienceDirect
Biochemical Systematics and Ecology journal homepage: www.elsevier.com/locate/biochemsyseco
Divergent adaptation to Qinghai-Tibetan Plateau implicated from transciptome study of Gymnocypris dobula and Schizothorax nukiangensis Mengchao Yu a, 1, Dongsheng Zhang a, 1, Peng Hu a, Sihua Peng a, Weiwen Li b, Shunping He c, Wanying Zhai a, Qianghua Xu b, **, Liangbiao Chen a, d, * a
Key Laboratory of Aquaculture Resources and Utilization, Ministry of Education, College of Fisheries and Life Science, Shanghai Ocean University, Shanghai, China Key Laboratory of Sustainable Exploitation of Oceanic Fisheries Resources, Ministry of Education, College of Marine Sciences, Shanghai Ocean University, Shanghai, China c The Key Laboratory of Aquatic Biodiversity and Conservation of Chinese Academy of Sciences, Institute of Hydrobiology, Chinese Academy of Sciences, Wuhan, China d Marine Biology and Biotechnology Laboratory, Qingdao National Laboratory for Marine Science and Technology, Qingdao, China b
a r t i c l e i n f o
a b s t r a c t
Article history: Received 29 December 2016 Received in revised form 2 February 2017 Accepted 4 February 2017
The Schizothoracine fishes are widely distributed in the Qinghai-Tibetan Plateau (QTP) area and its peripheral regions, which provide a prime example of adaptation in highland aquatic environments. Recent progresses have revealed various genetic adaptations of these fishes by comparing to distantly related lowerland species, however, comparative studies on closely-related species of different altitudes are still lacking. In this study, we sequenced and annotated a primitive Schizothoracine fish Schizothorax nukiangensis Tsao and a highly specialized one Gymnocypris dobula. We performed evolutionary analyses to investigate the candidate genes and signaling pathways involved QTP highland adaptation in both Schizothoracine fishes. Analysis of the 11,007 one-copy orthologs to the primitive cyprinid species, Danio rerio, revealed that both G. dobula and S. nukiangensis showed elevated evolutionary rates. A large number of genes related to hypoxia, including genes involved metabolic processes and cardiovascular system development, exhibited signatures of positive selection in both Schizothoracine fishes, but very few positively selected genes were found overlapping among these Schizothoracines. Our results indicated divergent genetic adaptation to highland environment for aquatic species living in QTP. © 2017 Elsevier Ltd. All rights reserved.
Keywords: Tibetan plateau Adaptive evolution Schizothoracine fish Hypoxia
1. Introduction The Qinghai-Tibet Plateau (QTP) is the highest and one of the biggest plateaus on earth, covering 2.5 106 square kilometers with an elevation of 3000e5000 m for most parts of the area. The QTP has been uplifting since approximately 45
* Corresponding author. ** Corresponding author. E-mail addresses: [email protected] (Q. Xu), [email protected] (L. Chen). 1 These authors contributed equally. http://dx.doi.org/10.1016/j.bse.2017.02.003 0305-1978/© 2017 Elsevier Ltd. All rights reserved.
98
M. Yu et al. / Biochemical Systematics and Ecology 71 (2017) 97e105
million years ago, resulted from collision of the India plate and the Eurasia plate (Li and Fang, 1999; Favre et al., 2014). The uplifting dramatically changed the environment conditions from an originally humid and warm climate to currently a dry and cold one (Wu et al., 2008). Characterized by factors including hypoxia, cold and strong ultraviolet radiation, the QTP environment posed harsh challenges to the endemic animals (Scheinfeldt and Tishkoff, 2010). Using comparative genomics and transcriptomics analysis, recent studies have identified various genes and signaling pathways that may be responsible for highland adaptation in both terrestrial and aquatic vertebrates. For example, genetic modifications to metabolism and cardiovascular system development are frequently found in the highland adapted animals (Qiu et al., 2012; Qu et al., 2013; Gou et al., 2014). The family Cyprinidae is the most diverse clade in freshwater fishes (Chen and Mayden, 2009). The Schizothoracine fishes (Teleostei: Cyprinidae) dominate the lakes and rivers throughout the Tibetan Plateau and its peripheral regions. More than 70 species in twelve genera in the subfamily Schizothoracinae are endemic to the Tibetan Plateau (He and Chen, 2006). According to the degree of specialization of their morphological traits, the Schizothoracine fishes are divided into three grades: Primitive, specialized, and highly specialized, with Primitive ones live in the peripheral regions at around 1000 m above sea level (a.s.l) and the other two groups are intermingled with each other in the central part of QTP at above 3000 m a.s.l (He and Chen, 2006; Qi et al., 2012). Recent transcriptome studies on Schizothoracines and other endemic Tibetan fishes have identified multiple biological processes and genes involved in highland adaptation, including genes that participate in metabolic processes and responses to hypoxia (Yang et al., 2014; Ma et al., 2015). Despite great progresses made in these investigations, most of current studies are limited to single highland species with comparisons made to distantly related teleosts. Such comparisons are limited in power to decipher the detailed molecular adaptations associated with Schizothoracine expansion in QTP. Therefore, genomewide investigations of adaptive signals among Schizothoracine species of differently graded specialization are highly demanded to further reveal the genetic mechanisms of adaptation in these fishes. The objective of this study was to carry out a comparative genome-wide screen for genes that might be involved in highland adaptation from species living at different altitudes. To achieve these aims, we generated transcriptomes of a primitive Schizothoracine fish Schizothorax nukiangensis Tsao, which lives at about 1,000 m a.s.l, and a highly specialized one Gymnocypris dobula which lives at 3000e4000 m a.s.l, and we performed evolutionary analyses on these data. 2. Materials and methods 2.1. Ethics statement This study was approved by the Ethics Committee for the Use of Animal Subjects of Shanghai Ocean University. 2.2. Sample collection, cDNA library construction and illumina sequencing G. dobula was collected from Lake Duoqing, situated at eastern QTP (28 03.370, 8917.830 ) with an altitude of 4509 m a.s.l, and a temperature at 11.0 C, the concentration of dissolved oxygen in the water was 1.9 mg/L; while S. nukiangensis was collected from Nujiang, Yunnan, situated on the eastern margin of QTP (25 41.24’, 98 53.22), with an altitude of 1201 m a.s.l and at a temperature of 16.0 C, the concentration of dissolved oxygen in the water was 7.8 mg/L. The collection loci were labelled in Fig. 1. For both species, five individuals were captured from the same sampling locations, and fishes with similar size were selected for tissue dissection. The fish were live trapped, anesthetized. Tissues including brain, gill, head kidney, kidney, liver, muscle, skin, spleen were quickly biopsied and placed in RNAlater (QIAGEN), and stored at 80 C on arrival at the laboratory till RNA extraction. Total RNAs from each tissue sample were extracted using TRIzol reagent (Invitrogen Corp., Carlsbad, CA). Library constructions from G. dobula and S. nukiangensis were made using Illumina Hiseq1500 RNA sample preparation kits (Illumina, San Diego, CA) following manufacturer's instructions. RNAseq libraries were quantified using an Agilent 2100 Bioanalyzer (Agilent Technologies). The libraries were then sequenced on an Illumina HiSeq 1500 platform (Illumina, Inc.) with 100 bps paired-end reads. 2.3. De novo assembly and annotation of transcripts Raw reads were cleaned by removing adapter sequences using sequence pre-processing tool Trimmomatic. Reads with Phred score 25 were kept for sequence assembly. The transcriptomes of the tissues were pooled to generate an improved de novo assembly using the Trinity package, and open reading frames (ORFs) were predicted using the transdecoder with a minimum length of 50 amino acids (Haas et al., 2013). The assembled unigenes were first annotated by searching for homologous sequences against National Center for Biotechnology Information (NCBI) nonredundant protein (Nr) database and then were searched against the KOG database using BLASTP (evalue<1e-5). We also used the KOBAS software to annotate the unigenes with Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways (Xie et al., 2011). As zebrafish is currently the most closely related species to the Schizothoracines which has a completed and well annotated genome, we further annotated the unigenes of the two Schizothoracine fishes by mapping them to zebrafish proteins using BLASTP bidirectional best hit (BBH) method. The protein dataset of zebrafish (Danio rerio) was downloaded from
M. Yu et al. / Biochemical Systematics and Ecology 71 (2017) 97e105
99
Fig. 1. Geographical locations of sample collection.
Ensembl online database (ftp://ftp.ensembl.org/pub/release-84/fasta/danio_rerio/pep/) and GO annotation data for zebrafish genome was downloaded from Ensembl BioMart.
2.4. Ortholog determination and sequence alignment Orthologs among the two Schizothoracine fishes and zebrafish were determined using OrthoMCL (version_2.0.9) with default settings (Li et al., 2003). Ortholog groups of single copy genes (one-copy orthologs) that were shared by all three species were retained for further analysis. These one-copy orthologs were aligned at the protein level using the Prank pro€ytynoja and Goldman, 2008) and the corresponding codon alignments were produced using PAL2NAL (Suyama et al., gram (Lo 2006) with the unreliable alignments removed using the Gblocks program (Castresana, 2000). The trimmed alignments with lengths longer than 150 bps were applied for further evolutionary analysis.
2.5. Calculation of evolutionary rates To estimate the evolutionary rates for the two Schizothoracine fishes and zebrafish, we calculated ratios of nonsynonymous (dN) to synonymous (dS) substitution rate (u ¼ dN/dS) using the free-ratio model in codeml program of the PAML package (Yang, 1997). The dN/dS ratios were estimated with both individual one-copy orthologs and concatenated alignments from all one-copy orthologs. The values of dN, dS, and dN/dS ratios were retrieved from output files of the program. Genes were discarded if N*dN or S*dS < 1 or dS >1 according to the method described in a previous study (Goodman et al., 2010). We drew the boxplot for the distributions of dN/dS ratios. We also estimated lineage-specific mean values of dN/ dS using concatenated alignments from all one-copy orthologs. To investigate what molecular processes and functions are under rapid evolution, we estimated evolutionary rates of genes based on their Gene Ontology (GO) terms. Genes were categorized by the GO terms they were annotated, and only GO categories with more than 30 genes were included in the analyses. Lineage-specific mean dN/dS values were estimated by the concatenated alignments of all one-copy orthologs in each GO category, and the average dN/dS values for each GO category were calculated with free-ratio model in codeml program, GO terms experiencing relatively accelerated evolution were identified using a binomial test following the method applied in a previous study (Sequencing and Consortium, 2005).
2.6. Identification of positively selected genes (PSGs) Since the Branch-site model in codeml program has the advantage of detecting positive selection that affects only a few sites on a pre-specified (foreground) branch of the species tree (Zhang et al., 2005), this model was used to detect positive selection in the single-copy orthologous genes. We used likelihood ratio test (LTR) to compare the alternative hypothesis of positive selection on the foreground branch to a null hypothesis with no positive selection on the foreground branches for each orthologous gene. Positively selected genes were inferred only if their P-values were less than 0.05 after FDR adjusting (Benjamini and Hochberg, 1995). GO and KEGG enrichment analyses of PSGs were performed using the Fisher's exact test (Mehta and Patel, 1983).
100
M. Yu et al. / Biochemical Systematics and Ecology 71 (2017) 97e105
3. Results 3.1. De novo assembly and annotation of transcriptomes RNA-sequencing generated 188.70 million raw reads in about 94 Gb of data from S. nukiangensis, and 198.82 million raw reads in about 98 Gb of data of G. dobula. The original sequencing data were deposited in the NCBI and can be accessed in the Short Read Archive (SRA) under the accession numbers (SRP094994 and SRP095655). De novo assembly of these reads generated 91,854 unigenes for G. dobula and 171,555 unigenes for S. nukiangensis. The N50s of the two datasets are 1665 bp and 1661 bp, respectively. A total of 97,648 ORFs in S. nukiangesis and 51,782 ORFs in G. dobula were predicted (Supplementary Table 1). About 94.5% (92,307) and 94.7% (49,302) ORFs in S. nukiangensis and G. dobula had significant matches to protein sequences in the NR database, respectively. There are 87% and 58.2% of the BLASTP top-hits matched zebrafish protein sequences in S. nukiangensis and G. dobula, respectively (Fig. 2). As the Schizothoracine fishes and zebrafish belong to the same family Cyprinidae, the transcripts were further annotated against the protein database of zebrafish (Danio.rerio) using BBH method. In total, 19,520 and 21,043 protein-coding genes of G. dobula and S. nukiangensis mapped to D. rerio, respectively. We also did functional annotation using KOG. A total of 26,706 ORFs (27.35% of 97,648 ORFs) in S. nukiangesis and 35,632 ORFs (68.81% of 51,782 ORFs) in G. dobula were assigned to 26 eukaryotic orthologous groups. The top two categories were “Signal transduction” and “General function prediction only”. 3.2. Identification of orthologous genes and dN/dS analyses A total of 15,623 orthologous groups (also known as orthogroups) were detected for G. dobula, S. nukiangensis and D. rerio. Among these, 13,330 orthogroups were shared by all species, including 11,141 one-copy ortholog groups. After sequence alignment and trimming, 11,007 one-copy orthologous genes were retained for estimating the evolutionary constraints acting on G. dobula, S. nukiangensis and D. rerio. The lineage-specific evolutionary rates were estimated by measuring dN/dS ratios using the free-ratio model in PAML, which allows a separate dN/dS ratio for each branch. The dN/dS ratios for each branch were calculated using a concatenated alignment of 11,007 one-copy orthologs, and the evolutionary rates for the species were G. dobula (0.3101) > S. nukiangensis (0.2167) > D. rerio (0.1456) (Fig. 3A). Furthermore, dN/dS ratios for each ortholog were evaluated at gene level using the free ratio model. There were 7980 genes included in the analysis as they met the criteria N*dN or S*dS > 1 and dS < 1. Both G. dobula and S. nukiangensis had significantly higher dN/dS ratios than D. rerio (Wilcoxon rank sum test, P-values <2.2e-16). There was no significant difference between dN/dS ratios of G. dobula and S. nukiangensis (Wilcoxon rank sum test, P > 0.05) (Fig. 3B). 3.3. Accelerated evolution of the G. dobula lineages and S. nukiangensis lineages To identify lineage-specific accelerated GO categories in G. dobula and S. nukiangensis, the average dN/dS ratios for each GO category were calculated for each branch using concatenated alignment of genes belong to the GO category. We analyzed 1276 GO categories with more than 30 one-copy orthologs (Supplementary Table 2). When comparing S. nukiangensis to D. rerio, there were much more GO categories with significantly higher average dN/dS ratios in S. nukiangensis than in zebrafish
Fig. 2. Top five hit species distributions for S. nukiangensis and G. dobula based on BLASTP against Nr databases.
M. Yu et al. / Biochemical Systematics and Ecology 71 (2017) 97e105
101
Fig. 3. The average dN/dS ratios for each species estimated using concatenated alignments for all one-copy orthologs (A) and each ortholog (B).
lineage (1077 vs. 14); and there were 1008 GO categories with significantly higher average dN/dS ratios in G. dobula while no GO categories showed significantly higher dN/dS ratios in zebrafish. Furthermore, when comparing G. dobula to S. nukiangensis, we found that 936 GO categories had significantly higher dN/dS ratios in G. dobula, including “immune response”, “metabolic process” and “DNA repair”. While only 16 GO categories in S. nukiangensis: including “hemostasis” and “blood coagulation” (Fig. 4). 3.4. Positively selected genes of G. dobula and S. nukiangensis The improved branch-site model in PAML was used to detect positively selected genes (PSGs) among the 11,007 one-copy orthologous genes (Zhang et al., 2005). We identified 315 and 239 PSGs in G. dobula and S. nukiangensis by setting respective branches as the foreground branch. There were 18 PSGs overlapping between the two species, including nup214, which is involved in regulation of glucokinase, and bcam, a member of the immunoglobulin superfamily (Supplementary Table 3). When comparing to a previous study on a specialized Schizothoracine fish Gymnodiptychus pachycheilus (Yang et al., 2014), we found that nup214 was the only one in common for the three species and that there were two PSGs overlapping between G. dobula and G. pachycheilus. We also conducted GO and KEGG functional classification for these PSGs in these two species (Fig. 5). PSGs in G. dobula and S. nukiangensis covered similar GO categories, including “Metabolic process” (84 and 67 PSGs, respectively), “response to stress” (18 and 12), “immune system process” (12 and 10), “cardiovascular system development” (9 and 7) and “cell adhesion” (9 and 6). Some GO categories were specifically enriched in G. dobula or S. nukiangensis. For example, “receptor activity” and “cell adhesion” were enriched in G. dobula, while “ion binding” and “erythrocyte differentiation” were specifically enriched in S. nukiangensis. Among KEGG classification, the class ‘Environmental Information Processing’ contained the largest number of PSGs, for example: there were 26 and 10 PSGs assigned to its subclass “signal transduction” for G. dobula and S. nukiangensis respectively. Another predominant KEGG class was “immune system”, which included 11 and 9 PSGs, respectively. Some genes related to DNA damage, for example, ddb2 in S. nukiangensis, and msh5 in G. dobula, were identified in our study, which may relate to high UV exposure in Schizothoracines. Since hypoxia is one of the most challenging environmental factors for QTP endemic species, we further examined hypoxia-related PSGs by comparing to hypoxiaDB, which is a database of hypoxia-regulated proteins extracted from previous publications (Khurana et al., 2013; Zhang et al., 2014). We identified 23 and 25 PSGs from G. dobula and S. nukiangensis in the hypoxiaDB, respectively. There were no common hypoxia-related PSGs in these two species (Supplementary Table 4). As key regulator in hypoxia response, HIF-1 signaling pathway is believed to play pivotal roles in highland adaptation (Qiu et al., 2012; Qu et al., 2013; Gou et al., 2014). In this study, we identified four PSGs among the 20 genes annotated in the KEGG HIF-1 signaling pathway, including three PSGs in G. dobula and one in S. nukiangesis (Fig. 6). 4. Discussion In this study, we performed multi-tissue transcriptome sequencing of two Schizothoracine fishes and found many genes and functions might be involved in highland adaptation for these fishes. As increased evolutionary rates of genes are regarded as signals of adaptive evolution at molecular level, we investigated the evolutionary rates of Schizothoracines by estimating dN/dS ratios. We found that dN/dS ratios of both Schizothoracines were significantly higher than that of zebrafish. Our results are consistent with some previous studies which found elevated dN/dS ratios in species endemic to Tibet comparing to those living at lower altitudes (Yang et al., 2014; Ma et al., 2015). These results indicated that Schizothoracines have undergone rapid evolution to adapt to the extreme environment of the Tibetan Plateau. On the other hand, as the chromosome numbers in the Schizothoracine fishes are frequently duplicated (Wang et al.,
102
M. Yu et al. / Biochemical Systematics and Ecology 71 (2017) 97e105
Fig. 4. Comparisons of average dN/dS ratios for each GO category with more than 30 orthologs. (A: G. dobula vs. D. rerio; B: S. nukiangensis vs. D. rerio, C: G. dobula vs. S. nukiangensis). GO categories with significantly higher dN/dS ratios for species on Y-axis (red) and X-axis (blue) are highlighted. Grey dots stand for GO categories without significant differences in dN/dS ratios between the two species. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
2016), the increased dN/dS ratios in the Schizothoracines may thus associate with the specific genome duplication events in the subfamily. We also compared dN/dS ratios between the two Schizothoracines. There is no significant difference between them when the comparison was made at individual gene level, but dN/dS ratios estimated by the concatenated alignments of orthologs indicated that G. dobula has increased average dN/dS ratio than S. nukiangensis. The incongruence of the results above may come from the datasets used for the estimation. As many as 3027 genes were not included in the gene level analysis, and thus the mean evolutionary rates estimated from individual genes may be biased comparing with the ones from concatenated alignments. Investigation on evolutionary rates of GO categories also suggested accelerated evolution in Schizothoracines. Most GO categories have higher dN/dS ratios in both Schizothoracines than in zebrafish. In addition, G. dobula have much more GO categories with elevated dN/dS ratios than S. nukiangensis, including many GO categories that were related to the adaptation of high altitude. These results may indicate that extensive genetic adaptation to highland occurred in these Schizothoracines, and that species living at higher altitude in QTP are under stronger selection pressure. The most striking feature of the QTP highland is hypoxia. Adaptation to hypoxia requires decreased oxygen consumption and/or increased oxygen intake. The former is mostly related to metabolism, while the latter involves cardiovascular system development and hemopoiesis. In this study, we identified as many as 84 and 67 PSGs that belong to the GO category
M. Yu et al. / Biochemical Systematics and Ecology 71 (2017) 97e105
103
Fig. 5. A: Distribution of GO classification of PSGs. B: Distribution of KEGG classification of PSGs. (A) Cellular Processes; (B) Environmental Information Processing; (C) Genetic Information Processing; (D) Metabolism; (E) Organismal Systems.
Fig. 6. HIF-1 signaling pathway in KEGG populated with one-copy ortholog genes annotated in this study. The annotated genes are shown in colored boxes: red ones are PSGs identified in G. dobula, blue ones in S. nukiangensis. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
“metabolic process” in G. dobula and S. nukiangensis, respectively. KEGG pathway inference of the PSGs again suggested that various metabolic processes in both Schizothoracines might have been evolved under selection pressure, such as lipid
104
M. Yu et al. / Biochemical Systematics and Ecology 71 (2017) 97e105
metabolism, energy metabolism, carbohydrate metabolism, amino acid metabolism and nucleotide metabolism. Nine and seven PSGs that associated with the GO category “cardiovascular system development” were identified in G. dobula and S. nukiangensis respectively. Three genes involved in “erythrocyte differentiation” was found under positive selection in S. nukiangensis, including erythropoietin, a principal regulator in erythropoiesis, suggesting possible developmental adaptation of this species in hemopoiesis. Furthermore, there were 23 and 25 PSGs listed in hypoxiaDB for G. dobula and S. nukiangensis, respectively, suggesting the two scizothoracines in this study are possibly under selection pressure caused by hypoxia. Hypoxia-inducible factor (HIF) is a master switch that mediates major changes in gene expression under low oxygen. Previous studies have shown that some genes involved in the HIF-1 signaling pathway are under positive selection in scizothoracine fishes living in QTP (Guan et al., 2014; Yang et al., 2014; Shao et al., 2015). In this study, three PSGs were identified to be in the HIF-1 signaling pathway in G. dobula, suggesting hypoxia adaptation in this species. Very few PSGs were found overlapping among different Schizothoracinae species. nup214, which is involved in regulation of glucokinase, was the only common PSG identified in these two Schizothoracines investigated in this study and another specialized schizothoracine fish, G. przewalskii (Yang et al., 2014), indicating divergent evolution during highland adaptation for Schizothoracines. Acknowledgement This work was supported in part by the National Natural Science Foundation of China (31130049), the Research Project of the Chinese Ministry of Education (no. 213013A), the National Natural Science Foundation of China (grant no. 31572598), the “Shuguang Program” supported by the Shanghai Education Development Foundation and the Shanghai Municipal Education Commission (grant no. 13SG51), and the Shanghai Municipal Project for First-Class Discipline of Fishery to Shanghai Ocean University. Appendix A. Supplementary data Supplementary data related to this article can be found at http://dx.doi.org/10.1016/j.bse.2017.02.003. References Benjamini, Y., Hochberg, Y., 1995. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J. R. Stat. Soc. B 57, 289e300. http://dx.doi.org/10.2307/2346101. Castresana, J., 2000. Selection of conserved blocks from multiple alignments for their use in phylogenetic analysis. Mol. Biol. Evol. 17, 540e552. http://dx. doi.org/10.1093/oxfordjournals.molbev.a026334. Chen, W.J., Mayden, R.L., 2009. Molecular systematics of the Cyprinoidea (Teleostei: cypriniformes), the world's largest clade of freshwater fishes: further evidence from six nuclear genes. Mol. Phylogenet. Evol. 52, 544e549. http://dx.doi.org/10.1016/j.ympev.2009.01.006. €ckert, M., Pauls, S.U., Ja €hnig, S.C., Uhl, D., Michalak, I., Muellner-Riehl, A.N., 2014. The role of the uplift of the Qinghai-Tibetan Plateau for the Favre, A., Pa evolution of Tibetan biotas. Biol. Rev. 90 http://dx.doi.org/10.1111/brv.12107 na/na. Goodman, M., Sterner, K.N., Islam, M., Sherwood, C.C., Hof, P.R., Hou, Z., Lipovich, L., Jia, H., Grossman, L.I., Derek, E., Pease, L.R., Goodman, M., Sterner, K.N., Islam, M., Uddin, M., Sherwood, C.C., Hof, P.R., 2010. Phylogenomic analyses reveal convergent patterns of adaptive evolution in elephant and human ancestries. Correction for Goodman et al. Proc. Natl. Acad. Sci. 107, 8498. http://dx.doi.org/10.1073/pnas.1003435107. Gou, X., Wang, Z., Li, N., Qiu, F., Xu, Z., Yan, D., Yang, S., Jia Konga, X., Wei, Z., Lu, S., Lian, L., Wu, C., Wang, X., Li, G., Teng, M., Jiang, Q., Zhao, X., Yang, J., Liu, B., Wei, D., Li, H., Yang, J., Yan, Y., Zhao, G., Dong, X., Li, M., Deng, W., Leng, J., Wei, C., Wang, C., Mao, H., Zhang, H., Ding, G., Li, Y., 2014. Whole-genome sequencing of six dog breeds from continuous altitudes reveals adaptation to high-altitude hypoxia. Genome Res. 24, 1308e1315. http://dx.doi.org/10. 1101/gr.171876.113. Guan, L., Chi, W., Xiao, W., Chen, L., He, S., 2014. Analysis of hypoxia-inducible factor alpha polyploidization reveals adaptation to Tibetan plateau in the evolution of schizothoracine fish. BMC Evol. Biol. 14, 192. http://dx.doi.org/10.1186/s12862-014-0192-1. Haas, B.J., Papanicolaou, A., Yassour, M., Grabherr, M., Blood, P.D., Bowden, J., Couger, M.B., Eccles, D., Li, B., Lieber, M., Macmanes, M.D., Ott, M., Orvis, J., Pochet, N., Strozzi, F., Weeks, N., Westerman, R., William, T., Dewey, C.N., Henschel, R., Leduc, R.D., Friedman, N., Regev, A., 2013. De novo transcript sequence reconstruction from RNA-seq using the Trinity platform for reference generation and analysis. http://dx.doi.org/10.1038/nprot.2013.084. He, D., Chen, Y., 2006. Biogeography and molecular phylogeny of the genus Schizothorax (Teleostei: Cyprinidae) in China inferred from cytochrome b sequences. J. Biogeogr. 33, 1448e1460. http://dx.doi.org/10.1111/j.1365-2699.2006.01510.x. Khurana, P., Sugadev, R., Jain, J., Singh, S.B., 2013. HypoxiaDB: a database of hypoxia-regulated proteins. Database 2013, 1e12. http://dx.doi.org/10.1093/ database/bat074. Li, J., Fang, X., 1999. Uplift of the Tibetan Plateau and environmental changes. Chin. Sci. Bull. 44, 2117e2124. http://dx.doi.org/10.1007/BF03182692. Li, L., Stoeckert Jr., C.J., Roos, D.S., 2003. OrthoMCL: Identification of Ortholog Groups for Eukaryotic Genomes, pp. 2178e2189. http://dx.doi.org/10.1101/gr. 1224503.candidates. €ytynoja, A., Goldman, N., 2008. Phylogeny-aware gap placement prevents errors in sequence alignment and evolutionary analysis. Science 320, Lo 1632e1635. http://dx.doi.org/10.1126/science.1158395. Ma, X., Dai, W., Kang, J., Yang, L., He, S., 2015. Comprehensive Transcriptome Analysis of Six Catfish Species from an Altitude Gradient Reveals Adaptive Evolution in Tibetan Fishes, pp. 1e33. http://dx.doi.org/10.1534/g3.115.024448. Mehta, C.R., Patel, N.R., 1983. A network algorithm for performing Fisher's exact test in r c contingency tables. J. Am. Stat. Assoc. 78, 427e434. http://dx. doi.org/10.2307/2288652. Qi, D., Chao, Y., Guo, S., Zhao, L., Li, T., Wei, F., Zhao, X., 2012. Convergent, parallel and correlated evolution of trophic morphologies in the subfamily schizothoracinae from the qinghai-tibetan plateau. PLoS One 7, 1e10. http://dx.doi.org/10.1371/journal.pone.0034070. Qiu, Q., Zhang, G., Ma, T., Qian, W., Wang, J., Ye, Z., Cao, C., Hu, Q., Kim, J., Larkin, D.M., Auvil, L., Capitanu, B., Ma, J., Lewin, H.A., Qian, X., Lang, Y., Zhou, R., Wang, L., Wang, K., Xia, J., Liao, S., Pan, S., Lu, X., Hou, H., Wang, Y., Zang, X., Yin, Y., Ma, H., Zhang, J., Wang, Z., Zhang, Y., Zhang, D., Yonezawa, T., Hasegawa, M., Zhong, Y., Liu, W., Zhang, Y., Huang, Z., Zhang, S., Long, R., Yang, H., Wang, J., Lenstra, J.A., Cooper, D.N., Wu, Y., Wang, J., Shi, P., Liu, J., 2012. The yak genome and adaptation to life at high altitude. Nat. Genet. 44, 946e949. http://dx.doi.org/10.1038/ng.2343. Qu, Y., Zhao, H., Han, N., Zhou, G., Song, G., Gao, B., Tian, S., Zhang, J., Zhang, R., Meng, X., Zhang, Y., Zhang, Y., Zhu, X., Wang, W., Lambert, D., Ericson, P.G.P., Subramanian, S., Yeung, C., Zhu, H., Jiang, Z., Li, R., Lei, F., 2013. Ground tit genome reveals avian adaptation to living at high altitudes in the Tibetan plateau. Nat. Commun. 4, 2071. http://dx.doi.org/10.1038/ncomms3071.
M. Yu et al. / Biochemical Systematics and Ecology 71 (2017) 97e105
105
Scheinfeldt, L.B., Tishkoff, S.A., 2010. Living the high life: high-altitude adaptation. Genome Biol. 11, 133. http://dx.doi.org/10.1186/gb-2010-11-9-133yrgb2010-11-9-133 [pii]. Sequencing, T.C., Consortium, A., 2005. Initial sequence of the chimpanzee genome and comparison with the human genome. Nature 437, 69e87. http://dx. doi.org/10.1038/nature04072. Shao, Y., Li, J.-X., Ge, R.-L., Zhong, L., Irwin, D.M., Murphy, R.W., Zhang, Y.-P., 2015. Genetic adaptations of the plateau zokor in high-elevation burrows. Sci. Rep. 5, 17262. http://dx.doi.org/10.1038/srep17262. Suyama, M., Torrents, D., Bork, P., 2006. PAL2NAL: robust conversion of protein sequence alignments into the corresponding codon alignments. Nucleic Acids Res. 34 http://dx.doi.org/10.1093/nar/gkl315. Wang, X., Gan, X., Li, J., Chen, Y., He, S., 2016. Cyprininae phylogeny revealed independent origins of the Tibetan Plateau endemic polyploid cyprinids and their diversifications related to the neogene uplift of the plateau. Sci. China Life Sci. 59, 1149e1165. http://dx.doi.org/10.1007/s11427-016-0007-7. Wu, Zhenhan, Barosh, P.J., Zhonghai, W., Daogong, H., Xun, Z., Peisheng, Y., 2008. Vast early Miocene lakes of the central Tibetan plateau. Bull. Geol. Soc. Am. 120, 1326e1337. http://dx.doi.org/10.1130/B26043.1. Xie, C., Mao, X., Huang, J., Ding, Y., Wu, J., Dong, S., Kong, L., Gao, G., Li, C.Y., Wei, L., 2011. KOBAS 2.0: a web server for annotation and identification of enriched pathways and diseases. Nucleic Acids Res. 39 http://dx.doi.org/10.1093/nar/gkr483. Yang, Z., 1997. User Guide PAML: Phylogenetic Analysis by Maximum Likelihood, vol. 6, pp. 1e70. http://dx.doi.org/10.1093/molbev/msm088. Yang, L., Wang, Y., Zhang, Z., He, S., 2014. Comprehensive transcriptome analysis reveals accelerated genic evolution in a Tibet fish, Gymnodiptychus pachycheilus. Genome Biol. Evol. 7, 251e261. http://dx.doi.org/10.1093/gbe/evu279. Zhang, J., Nielsen, R., Yang, Z., 2005. Evaluation of an improved branch-site likelihood method for detecting positive selection at the molecular level. Mol. Biol. Evol. 22, 2472e2479. http://dx.doi.org/10.1093/molbev/msi237. Zhang, W., Fan, Z., Han, E., Hou, R., Zhang, L., Galaverni, M., Huang, J., Liu, H., Silva, P., Li, P., Pollinger, J.P., Du, L., Zhang, X.Y., Yue, B., Wayne, R.K., Zhang, Z., 2014. Hypoxia adaptations in the grey wolf (Canis lupus chanco) from qinghai-tibet plateau. PLoS Genet. 10 http://dx.doi.org/10.1371/journal.pgen. 1004466.