Main

The yak (B. grunniens) is an iconic symbol of Tibet and of high altitude. More than 14 million domestic yaks provide the basic resources (such as meat, milk, transportation, dung for fuel and hides for tented accommodation) that are necessary for Tibetans and other nomadic pastoralists in high-altitude environments1. In contrast, closely related and cross-fertile taurine cattle (B. taurus) suffer from severe pulmonary hypertension when reared in the yak habitat2,3,4. Yaks have numerous anatomical and physiological traits that equip them for life at high altitude, including large lungs and hearts1, lack of hypoxic pulmonary vasoconstriction5, increased foraging ability6, strong environmental sense1 and high energy metabolism1,7. Thus, comparing yak and cattle contributes to the understanding of evolutionary adaptation to high altitude5,6,7.

Genomic comparisons between closely related species provide insights into the genetic basis of mammalian divergence and adaptation8,9. In this study, we sequenced the genome of a female domestic yak using a whole-genome shotgun strategy and the Illumina HiSeq 2000 platform. De novo assembly of 4.4 billion reads from paired-end libraries (Supplementary Table 1) yielded a draft assembly (65-fold coverage) with a total length of 2,657 Mb, close to the 2,649 Mb of sequence obtained for the cattle genome10 (UMD 3.1), and contig and scaffold N50 sizes of 20.4 kb and 1.4 Mb, respectively (Supplementary Table 2). Approximately 90% of the total sequence was covered by 2,083 scaffolds of >307 kb, with the largest scaffold spanning 8.8 Mb. The assembly metrics of the yak genome were comparable to those of other animal genome assemblies generated by next-generation sequencing technology (Supplementary Table 3). The sequencing depth of 98% of the assembly was more than 20-fold (Supplementary Fig. 1), ensuring high accuracy at the nucleotide level11. The genome assembly covered, without any obvious errors in assembly, 97% of 6 fosmids sequenced by Sanger sequencing and 98% of 81,020 Unigenes assembled from Illumina RNA sequencing (RNA-seq) data for 5 tissues (Supplementary Fig. 2 and Supplementary Tables 4 and 5), indicating coverage of most of the euchromatic regions. GC content distributions were similar to those of the cattle genome (Supplementary Fig. 3 and Supplementary Note). We predicted 22,282 protein-coding genes in yak on the basis of RNA-seq, homology and ab initio gene prediction (Supplementary Fig. 4 and Supplementary Table 6).

High coverage11,12 permitted the identification of 2.2 million heterozygous single-nucleotide variants (SNVs) within the sequenced individual (Supplementary Figs. 5 and 6 and Supplementary Tables 7 and 8). The heterozygosity rate (0.89 × 10−3) was approximately 1.5 times higher than that estimated for cattle (0.59 × 10−3) (refs. 10,13). This may be due to a longer and more systematic selection in cattle and/or to introgression from wild yaks still living on the Qinghai-Tibetan Plateau14,15. Yak and cattle both have 30 chromosomes and have similar karyotypes1, which allows chromosomal assignments on the basis of homology, despite a lack of physical maps for yak chromosomes. Using the human genome as an outgroup, we reconstructed 207 conserved ancestral homologous synteny blocks (aHSBs)16 covering 2.51 Gb (94%) of the yak genome. The existence of these blocks confirms the extensive synteny of the cattle and yak genomes and also allows the detection of breakpoint regions (Supplementary Table 9).

Frequent turnover of gene copy and/or protein domain number has been proposed as a major mechanism underlying the adaptive divergence of closely related species9,17,18,19. First, we used TreeFam20 to identify 13,810 homologous gene families shared by 4 species (yak, cattle, human and dog): 362 gene families were specific to yak and cattle, and 100 were found only in yak (Fig. 1). The sequence depth of these multiple-copy genes was in the same range as for single-copy genes (Supplementary Fig. 7). The yak-specific gene families contained 170 genes, 75 of which have known InterPro domains. Compared with the cattle-specific gene families, the yak-specific families were significantly over-represented in two major functional categories: olfactory sensation (14 genes; P < 0.01) and host defense and immunity (11 genes; P < 0.01) (Supplementary Table 10). Next, we identified 596 gene families that were substantially expanded in yak compared to other mammals (Fig. 2a). Functional categories that were enriched for significant gene family expansions mainly included sensory perception (gene ontology (GO) 0004984, olfactory receptor activity, P < 0.01; GO 0050909, sensory perception of taste, P < 0.01) and energy metabolism (GO 0004129, cytochrome-c oxidase activity, P < 0.01; GO 0015986, ATP synthesis coupled proton transport, P < 0.01) (Supplementary Table 11). Third, matching of ORFs to PFAM domain families at the protein level showed expansion in the numbers of specific domains (Supplementary Table 12). Thus, with regard to sensory perception receptors (olfactory and taste) and other G protein–coupled receptor (GPCR) rhodopsin-like receptors, known to be involved in sensing of the extracellular environment21, we found significantly more GPCR transmembrane domains in yak than in cattle (1,558 versus 1,358). Categories related to hypoxic stress seemed also to have enriched expansions of the corresponding domains in yak. For example, genes with the Hig_1_n (PF04588.6 in Pfam) domain were highly expressed under hypoxic stress22,23 and also had expanded copy numbers in yak (13 copies) relative to cattle (9 copies) and other mammals. Phylogenetic analysis of genes that encoded this domain showed an expansion in the numbers of the closely related Hig_1_n domains in both yak and cattle, as well as three additional copies in yak (Fig. 2b).

Figure 1: Venn diagram showing unique and shared gene families between the yak, cattle, dog and human genomes.
figure 1

The number of gene families is listed in each of the diagram components and the total number for each animal is given in parentheses.

Full size image
Figure 2: Gene expansion and contraction in the yak genome.
figure 2

(a) Dynamic evolution of orthologous gene families. The proportions of expanded (red) and contracted (blue) gene families are shown as pie charts at each branch terminus. MRCA, most recent common ancestor. (b) A neighbor-joining tree of mammalian Hig domain sequences.

Full size image

Adaptive divergence at the molecular level may also be expressed by an increased rate of nonsynonymous changes within genes involved in adaptation9,24. We identified 8,923 high-confidence 1:1:1 orthologous genes in the yak, cattle and human genomes, most of which also correspond to genes in the horse, dog, mouse and chimpanzee genomes (Supplementary Fig. 8). Overall, yak and cattle genes were highly similar, with 45% of encoded proteins identical and mean protein similarity approximating 99.5% (Supplementary Figs. 9 and 10). Average synonymous (Ks) and nonsynonymous (Ka) gene divergence values between yak and cattle were 0.0114 and 0.00207, respectively, close to the values between human and chimpanzee genes (Supplementary Fig. 11). Yak and cattle were estimated to have diverged approximately 4.9 million years ago, which is comparable to the time at which humans and chimpanzees diverged9 (Supplementary Fig. 12). Ka/Ks ratios of nonsynonymous-to-synonymous substitutions for different GO categories revealed an enrichment of elevated pairwise Ka/Ks values in the hypoxia response25,26,27 and energy metabolism28,29 categories, including in 'regulation of blood vessel size', 'regulation of angiogenesis', 'heme binding', 'glycerolipid biosynthetic process' and 'electron carrier activity' (Supplementary Table 13). Analysis of Ka/Ks ratios in the cattle and yak lineages verified that genes with elevated Ka/Ks values in yak were significantly enriched for these functions (Fig. 3a).

Figure 3: Adaptive evolution in the yak genome.
figure 3

(a) Data points represent pairs of yak and cattle mean Ka/Ks ratios by GO category. GO categories with putatively accelerated (P < 0.05, binomial test) nonsynonymous divergence in the yak lineage (red) and in the cattle lineage (blue) are highlighted. A complete list of categories is provided in Supplementary Table 17. (b) Comparison of the proportions of genes showing evidence for positive selection in the yak and cattle lineages. The numbers of positively selected genes are given in parentheses. (c) Five genes involved in integrated nutrition pathways (according to KEGG pathway: map04971 and map00020) were found to show evidence of positive selection in the yak lineage. Solid lines indicate direct relationships between enzymes and metabolites. Dashed lines indicate that more than one step is involved in the process.

Full size image

To test the hypothesis that these rapidly evolving genes in yak have been under positive selection, we used the branch-site likelihood ratio test to identify positively selected genes (PSGs) in both the yak and cattle lineages. We identified 85 PSGs (in yak) and 95 PSGs (in cattle) (Fig. 3b and Supplementary Table 14). The PSGs detected in yak were enriched for genes involved in the hypoxia response and energy metabolism (Supplementary Tables 15 and 16). Of 81 genes examined in the response to hypoxia functional category (GO 0001666), 3 (3.7%) showed evidence of positive selection in yak (compared to none in cattle), which is significantly higher than the background level of positive selection across the genome (P < 0.05) (Fig. 3b). The three yak PSGs comprise two important regulators (Adam17 and Arg2) and one target gene (Mmp3) of hypoxia-inducible factor–1α (Hif-1α). As a master regulator of the cellular response to hypoxia, Hif-1α triggers wide transcription of genes involved in angiogenesis, vasodilatation and energy metabolism30,31,32,33,34,35. Notably, alleles of human ADAM17 were previously shown to be present at significant different frequencies in Tibetans and low-altitude dwellers36, indicating a possible role for this gene in altitude adaptation. The Adam17 and Arg2 proteins affect Hif-1α stability and activity by regulating production of tumor necrosis factor α (TNF-α)31,32 and nitric oxide, respectively33,34, whereas Mmp3 has key roles in numerous physiological processes35. In their high-altitude environments, yaks must not only maintain normal energy production under hypoxic pressure7,37 but must also optimize nutritional assimilation, as a consequence of the limited herbal resources available1. Indeed, we found five key genes that show signs of positive selection in yak nutrition pathways (Fig. 3c): Camk2b regulates the secretion of gastric acid in the rumen, which contributes to the assimilation of volatile fatty acids produced by ruminal fermentation38,39,40,41,42, and Gcnt3, Hsd17b12, Whsc1 and Glul have important roles in polysaccharide, fatty acid and amino-acid metabolism, respectively43,44,45,46. In addition, the positively selected changes in Glul may be important for the high level of nitrogen utilization in yak7.

Our evolutionary analyses based on genomic data have provided important insights into adaptation to high altitude in yak. Further understanding may be gained by functional analysis of the identified genes with signs of adaptive evolution in comparative stress studies of yak and other animals living at high altitude. The identification of genes required for natural high-altitude adaptation may help to improve current understanding, treatment and prevention of altitude sickness and other hypoxia-related diseases in humans. In addition, this report of the yak genome sequence, together with the many SNVs identified, will facilitate genetic dissection of agronomically important traits in the species and will accelerate the genetic improvement of milk and meat production in this animal that is essential to the lifestyle and economy of the Tibetan people.

URLs.

SOAP, http://soap.genomics.org.cn/; Ensembl, http://www.ensembl.org/; TimeTree, http://www.timetree.org/.

Methods

Genome sequencing and assembly.

Genomic DNA was extracted from the liver of a female yak with an estimated inbreeding coefficient of 0.094 (ref. 1) that lived above 3,700 m in Huangyuan County of Qinghai Province, China. Sequencing libraries were constructed with multiple insert sizes (200 bp to 20 kb) according to the Illumina protocol. For short insert libraries (200 to 800 bp), 6 μg of DNA was fragmented to the desired insert size, end-repaired and ligated to Illumina paired-end adaptors. Ligated fragments were size selected at 200, 500 and 800 bp on agarose gels and were purified by PCR amplification to yield the corresponding libraries. For long insert sizes (2, 5, 10 and 20 kb) mate-pair library construction, 60 μg of genomic DNA was used; we circularized DNA, digested linear DNA, fragmented circularized DNA and purified biotinylated DNA and then performed adaptor ligation. All libraries were sequenced on an Illumina HiSeq 2000 platform.

Whole-genome shotgun assembly of the yak was performed using short oligonucleotide analysis package (SOAP)denovo47 (Supplementary Note). First, a series of strict filtering steps were performed before assembly to avoid artifactual duplication, adaptor contamination and low-quality reads. Second, reads from the short-insert (≤800-bp) libraries were assembled into distinct contigs on the basis of k-mer overlap information. Third, reads from the long-insert (≥800-bp) libraries were aligned to the contig sequence, and the paired-end relationships between reads were used to construct scaffolds. We used a hierarchical assembly strategy in which we added data step by step from short paired-end reads to long paired-end reads. Finally, in order to fill gaps between scaffolds, we used the paired-end information to retrieve read pairs that had one read well aligned on the contigs and another read located within the gap region. We then performed a local assembly of the collected reads.

Transcriptome and fosmid sequencing and assembly.

From the same yak animal, RNA was extracted from fresh heart, liver, brain, stomach and lung tissues for the generation of transcriptome data. The quality and integrity of the RNA samples were examined using the Agilent 2100 Bioanalyzer, their RNA integrity number (RIN) values ranged from 8.6–10.0, with no sign of degradation. Approximately 20 μg of total RNA (at a concentration of ≥400 ng/μl) from each tissue was used to construct cDNA libraries. Poly(A) mRNA was isolated using beads conjugated to oligo(dT), mRNA was fragmented, and cDNA was synthesized using random hexamer primers and reverse transcriptase (Invitrogen). After end repair, adaptor ligation and PCR amplification, the libraries were sequenced using the Illumina HiSeq 2000 platform. Libraries that gave reads that were unevenly distributed among the gene regions (for example, showing a strong bias toward 5′ or 3′ regions) were discarded and replaced. Transcripts were assembled using SOAPdenovo. Reads for the assembly were filtered as for the genome assembly, and duplicate reads were removed. For further study, we used only those transcripts that were longer than 150 bp and that were covered at least twice.

We also constructed a fosmid (inset size of ∼40 kb) library from the same DNA resource and randomly selected six clones for Sanger sequencing. These fosmids were assembled using the Celera Assembler48. We then evaluated the completeness and accuracy of the genome assembly by comparing the assembled scaffolds with 6 complete fosmid clones and 81,020 Unigene sequences using BLAST searches.

Heterozygous SNV detection.

To evaluate the heterozygosity rate and its distribution, high-quality reads (average quality score of >30) from short-insert libraries were realigned to the assembly with SOAP (see URLs). The probabilities of each possible genotype at every position on the reference genome were calculated, and a statistical model based on Bayesian theory and the Illumina quality system was used to call SNVs. The allelic sequence with the highest probability was used as the reference sequence, and heterozygous SNVs were called if other alleles also had high probability11. To estimate the accuracy of the identified heterozygous SNVs, we randomly selected 150 heterozygous SNVs and validated them by PCR amplification and Sanger sequencing. In 146 of the 150 sequences, SNVs were validated by double sequence peaks (Supplementary Fig. 6 and Supplementary Table 8).

Annotation.

Transposable elements in the yak genome assembly were first identified using a combination of homology-based and de novo approaches at both the DNA and protein levels47. We then used homology and ab initio prediction, as well as RNA-seq to identify protein-coding genes, building a consensus gene set by merging all predicted genes. For homology-based gene prediction, we aligned human and cattle protein sequences to the repeat-masked yak genome using TBLASTN and Genewise49 for fast alignment and accurate spliced alignment, respectively. Next, we used the ab initio gene prediction methods Genscan50 and Augustus51 to predict protein-coding genes, using parameters trained from a set of high-quality homolog prediction proteins. RNA-seq–based gene prediction was performed by aligning all RNA-seq data against the assembled genome using TopHat52, and Cufflinks53 was used to predict cDNAs from the resultant data. The final gene set was generated by merging these three gene prediction resources by GLEAN54. Gene functions were assigned according to the best match of the alignment to the SwissProt and Translated EMBL Nucleotide Sequence Data Library (TEMBL) databases, using BLASTP. For yak reference genes, motifs and domains were determined by searches in InterProScan55 of the sequences against publicly available databases, including Pfam, PRINTS, PROSITE, ProDom and SMART. Gene Ontology56 IDs for each gene were obtained from the corresponding InterPro entry. We also mapped yak reference genes to KEGG57 pathway databases and identified the best match for each gene.

Gene families.

The protein-coding genes from 6 mammalian species (Canis familiaris, Homo sapiens, Mus musculus, B. taurus, Equus caballus and Monodelphis domestica) downloaded from Ensembl (release 56; see URLs) were used in addition to yak genes to define gene families that descended from a single gene in the last common ancestor20. The longest translation form was chosen to represent each gene, and stretches of genes encoding fewer than 30 amino acids were filtered out. The 9,393 single-copy families obtained from this analysis were used to reconstruct phylogenies and estimate the times since divergence. Data from fourfold-degenerate sites were extracted from each family and attributed to one 'super gene' for each species. Modeltest58 was used to select the best substitution model (GTR + gamma + I), and Mrbayes59 was used to reconstruct the phylogenetic tree. The MCMCtree program implemented in the Phylogenetic Analysis by Maximum Likelihood (PAML)60 package was used to estimate the time since cattle-yak divergence. Calibration time was obtained from the TimeTree database (see URLs). To identify gene families that had undergone expansion or contraction, we applied the Café program, which is based on a probabilistic graphical model61, to infer the rate and direction of change in gene family size over a given phylogeny.

Evolutionary analyses.

We used conserved genome synteny methodology62 to establish a high-confidence orthologous gene set that included yak, cattle (UMD 3.1), horse (EquCab2.0), dog (CanFam2.0), mouse (mm9), chimpanzee (panTro2) and human (hg19) genes. Briefly, whole-genome multiple alignments were constructed for the relevant genomes using the MULTIZ63 alignment pipeline, with the human genome serving as the reference genome. To minimize the effect of annotation errors, variations in sequence quality and changes in gene structure on subsequent evolutionary rate analyses, we mapped all the human protein-coding genes from RefSeq64, KnownGene65 and VEGA66 to each of the other species via their syntenic alignments, then passed the resulting blocks through a series of rigorous filters that selected for large-scale synteny, high alignment quality and conservation of exon-intron structure. All orthologs were aligned using the codon option in the Probabilistic Alignment Kit (PRANK)67 program, and alignments shorter than 150 bp were discarded. The values of Ka and Ks and the Ka/Ks ratio were estimated for each gene using the Codeml program with the free-ratio model in the PAML package, and 10,000 concatenated alignments constructed from 150 randomly chosen genes were used to estimate lineage-specific mean values. The human GO annotation download from Ensembl was used to assign GO categories to 8,923 orthologs. The binomial test9 was used to identify GO categories with more than 20 orthologs that had an excess of nonsynonymous changes in either yak or cattle lineages. To detect genes evolving under positive selection in either yak or cattle (Supplementary Table 17), we used the optimized branch-site model68 in which likelihood ratio test (LRT) P values were computed assuming that the null distribution was a 50:50 mixture of a chi-squared distribution with 1 degree of freedom and a point mass at zero. Fisher's exact tests were used to test for over-represented functional categories among positively selected genes8. For each category C and set of PSGs S, a 2 × 2 contingency table was constructed for the numbers of genes assigned or not assigned to C and within or outside S. Then, (one-sided) P values for the independence of rows and columns were computed by Fisher's exact test. In addition, the distributions of LRT P values among the genes assigned and not assigned to C were compared by a (one-sided) Mann-Whitney U test. To test whether the unique mutations in yak, which resulted in the detected signal of positive evolution in the yak lineage, were specific to yak, we amplified and sequenced the DNA fragments encompassing the candidate yak-specific mutations in 15 genes (including the 8 genes shown in Fig. 3b,c) in 5 randomly selected yaks and 5 cattle. All mutations were confirmed to be specific to yak.

Accession codes.

The yak whole-genome shotgun project has been deposited at the DNA Data Bank of Japan (DDBJ), the European Molecular Biology Laboratory (EMBL) nucleotide sequencing database and GenBank under the same accession, AGSK00000000. The version of the genome described in this paper is the first version, AGSK01000000 (available at DDBJ, EMBL and GenBank). The mitochondrial sequence has been deposited at GenBank under accession JQ692071. All short-read data have been deposited at the Short Read Archive (SRA) under accession SRA047288. Raw sequencing data for the transcriptome have been deposited in the Gene Expression Omnibus (GEO) under accession GSE33300.