The sequence and de novo assembly of the wild yak genome (2024)

The yak, which can survive extremely cold, harsh and oxygen-poor conditions, is endemic to the Qinghai-Tibetan Plateau (QTP) and adjacent high-altitude regions¹. Yak were domesticated by nomadic people from wild yak at least 7300 years ago². Nowadays, more than 22 million domestic yak (Bos grunniens) provide necessities, such as food, transport, shelter and fuel, for Tibetans and other humans in high-altitude areas¹. In addition, there are still 15–20 thousand wild yak (Bos mutus) surviving in northwestern parts of the QTP^3,4. Due to long-term over-breeding and inbreeding caused by traditional yak breeding practices, the reproductive capacity, growth rate, adult size and milk production of domestic yak have declined and mortality has increased, especially among the newborn and young¹. However, Datong yak, the only artificially cultivated yak breed that is a cross between wild yak and domestic yak, shows excellent growth characteristics and production performance. Datong yak are generally 30 and 50% heavier than domestic yak at birth and six months of age, respectively, produce more than 15% milk, and have 25 and 31% higher carcass weights at ages of 6 and 18 months, respectively^5,6. The high growth, development and production rates of the Datong yak show the feasibility of improving traits of domestic yak with wild yak resources, and potential importance of exploiting wild yak genetic resources for yak breeding in the future. However, the wild yak genome has not been previously sequenced, which has impeded both research and breeding efforts.

Thus, to elucidate genomic features of this vulnerable species, we have constructed a draft genome for wild yak. We extracted genomic DNA from blood tissues, constructed 3 Paired-End (PE) and 3 Mate-Pair (MP) libraries, which were sequenced using the Illumina HiSeq. 2000 platform. After quality filtering and trimming of raw data, Genome Characteristics Estimation (GCE, v1.0.0)⁷ software was employed to evaluate the genome size using PE reads, and Platanus v1.2.4⁸ to assemble the genome using all clean data. In addition, GapCloser v1.12⁹ was used to perform another round of gap closure based on the assembly results. The final genome assembly size was 2.83 Gb, containing 808,541 contigs (N50 = 63.2 kb) and 734,073 scaffolds (N50 = 16.3 Mb), representing 91.5% of the estimated genome. Structural annotation of the genome yielded 22,910 genes, 90.18% of which could be functionally annotated with at least one of the five protein databases (TrEMBL, SwissProt, InterPro, GO and KEGG). The wild yak genome assembled in this study provides a valuable genetic resource for future efforts to protect the vulnerable wild yak and further comparative analysis of genome biology among bovine species to promote breeding research.

Sample collection, library construction and sequencing

Genomic DNA was extracted from a blood sample collected from a female wild yak originally captured from the wild and reared at the Datong Yak Farm of Qinghai Province (37°15′35.6″N, 101°22′24.0″E, altitude around 3200 m) using a standard phenol/chloroform method. The quality and integrity of the extracted DNA were checked by measuring its A260/A280 ratio and agarose gel electrophoresis. For paired-end libraries with insert sizes of 280, 500 and 800 bp, 6 μg portions of genomic DNA were used to generate the corresponding libraries using Illumina TruSeq DNA Nano Preparation Kit (Illumina, San Diego, CA, USA). For mate-pair libraries with insert sizes of 2, 5 and 10 kb, 60 μg portions of DNA were used for circularization and further library construction using Nextera Mate Pair Library Preparation Kit (Illumina, San Diego, CA, USA). Both the sample collection and experimental library construction protocol were approved by the Ethical Committees of Lanzhou University. All libraries were sequenced on an Illumina HiSeq. 2000 platform with 150 bp read length, following the manufacturer’s instructions. Finally, 760.85 Gb of raw data were generated in total (Table1).

Preprocessing and genome size estimation

All the sequencing reads were preprocessed for quality control and filtered with stringent criteria using Lighter v1.1.1¹⁰ software. Firstly, raw data were filtered by removing reads with >10% unknown bases. Then, paired reads with low-quality bases (quality scores ≤7) covering more than 65% of the read length were filtered out. Reads with PCR duplicates or adapter contamination were also removed. Finally, both read 1 and read 2 files were filtered out if they had >10 bp overlap, allowing 10% mismatch. In total, 662.3 Gb of clean reads were obtained after filtering (Table1).

Prior to genome assembly, all the preprocessed sequences from the short-insert library were subjected to genome size estimation using Genome Characteristics Estimation (GCE) with a k value of 21. The genome size of wild yak was estimated to be around 3.09 Gb, using the following formula: genome size = k-mer number/k-mer depth, where the k-mer number refers to the total number of k-mers, and k-mer depth is the depth of the main peak in the k-mer frequency distribution (Fig.1).

Genome assembly

For de novo genome assembly, Platanus software was used for constructing contigs and scaffolds with default parameters, and GapCloser was employed to fill the remaining gaps in the scaffolds with all sequencing reads. These steps finally yielded a wild yak draft genome with a total length of 2.83 Gb, accounting for 91.5% of the estimated genome size (contig and scaffold N50 sizes: 63.2 kb and 16.3 Mb, respectively) (Table2).

To evaluate the completeness of our assembly, we carried out BUSCO¹¹ analyses and the results indicated that 3,974 of the 4,104 conserved single-copy genes in mammals were present in our assembly, of which 3,799 were single, 55 were duplicated, and 120 fragmented matches (Table3). To validate the single-base accuracy of the genome assembly, we aligned the high-quality reads of short-insert libraries to the assembly using Burrows-Wheeler Aligner (BWA, v0.7.15-r1140)¹² software, and the alignment outputs were converted to Binary Alignment Map (BAM) format via SAMtools v1.3¹³. The genome coverage was then calculated by a custom Perl script, which indicated that more than 93.9% of the assembly had >20-fold coverage.

Repeat annotation

Repetitive regions of the wild yak genome were identified using a combination of de novo and hom*ology-based approaches, as applied in a previous analysis of the Ovis ammon polii genome¹⁴. For the de novo prediction, RepeatModeler v1.0.11 was employed first to construct a de novo repeat library, then RepeatMasker v4.0.7¹⁵ was used to identify repeats using both the RepBase¹⁶ library of known transposable elements (TEs) and a self-trained repeat database. Next, we applied RepeatProteinMask (a package in RepeatMasker) to identify repeats at the protein level using the TE protein database. In addition, tandem repeats were further annotated using Tandem Repeat Finder (TRF, v4.0.9)¹⁷. Finally, the non-redundant repeats were checked according to their coordinates in the genome. Overall, we identified 1.41 Gbp of non-redundant repetitive sequences, representing 49.65% of the wild yak genome assembly; of which long interspersed elements (LINE) were the most abundant, accounting for 35.98% of the whole genome (Fig.2; Table4).

Gene prediction and annotation

We employed a combination of hom*ology-based and de novo prediction methods to identify protein-coding genes. For hom*ology-based prediction, protein sequences of seven species (Bos taurus, Equus caballus, hom*o sapiens, Ovis aries, Sus scrofa, Bison bonasus, Bos grunniens) downloaded from Ensembl¹⁸ and GigaDB^19,20 were aligned to the wild yak genome using TBLASTN²¹. Then GeneWise v2.4.1²² software was applied to search for accurately spliced alignments based on the filtered hom*ologous genome sequences. For de novo prediction, we used Augustus²³, Geneid²⁴, GeneMark, GlimmerHMM²⁵ and SNAP²⁶ to predict genes with parameters trained on wild yak and human repeat-masked genomes. EVidenceModeler software (EVM, v1.1.1)²⁷ was employed to generate a consensus gene set by integrating the genes predicted by the hom*ology and de novo approaches. Low-quality genes of short length (proteins with fewer than 30 amino acids) and/or exhibiting premature termination were removed to produce the final gene set, which is composed of 22,910 genes (Fig.3; Table5).

Putative biological functions of these predicted high-quality genes were assigned by searching against five publicly available databases: TrEMBL, Swiss-Prot²⁸, InterPro²⁹, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG)³⁰. Approximately 90.18% of these genes were functionally annotated with at least one of these databases, with 90.05, 88.15, 83.51, 64.20 and 53.01% scoring positive hits in TrEMBL, SwissProt, InterPro, GO and KEGG, respectively (Table6).

Quality assessment of the genome assembly

The assembly presented here is the first wild yak genome version. The contig N50 and scaffold N50 sizes were 63.2 kb and 16.3 Mb respectively, with the longest scaffold 75,900,441 bp. There are 258 scaffolds more than 1 Mb long, with a total length of 2,486,540,864 bp, representing 87.83% of the wild yak genome. By aligning the reads of short insert libraries to the wild yak assembly, we found more than 93.9% of the genome had >20-fold coverage, indicating high accuracy at the nucleotide level. BUSCO analysis carried out to assess the completeness of our assembly resulted in a BUSCO score of 96.8% (complete = 93.8%, single = 92.4%, duplicated = 1.4%, fragmented = 3.0%, missed = 3.2%, genes = 4,104). These results are comparable with those for the published European bison (wisent)³⁴ and domestic yak⁴ genomes, suggesting our assembly has high quality and is quite complete.

Gene prediction and annotation validation

Gene models in the wild yak assembly were predicted using a combination of hom*ology-based and ab initio gene approaches. Then EVM software was employed to integrate the gene prediction results to produce a consensus gene set. To enhance the quality of the gene prediction, we removed low-quality genes of short length (proteins with fewer than 30 amino acids) and/or exhibiting premature termination. The final gene set consisted of 22,910 genes, and the distributions of gene length, CDS length, exon length, intron length and exon number were similar to those of other mammals (Fig.3). BUSCO analysis was also performed to assess the completeness of these predicted genes, resulting in a BUSCO value of 97.7% (complete = 94.8%, single = 93.1%, duplicated = 1.7%, fragmented = 2.9%, missed = 2.3%, genes = 4,104) (Table6). In addition, functional annotation of these predicted genes indicated that 90.18% of them could be assigned to at least one functional term (Table5). These results clearly indicate that the annotated gene set of the wild yak genome is quite complete.

The software versions, settings and parameters used are described below.

(1) GCE, version 1.0.0, parameters used: kmer_freq_hash -k 21 -l reads.list -t 24 -i 5000000 -o 0 -p wild_yak & > kmer_freq.log; gce -f wild_yak.freq.stat -c 46 -g 155866493014 -m 1 -D 8 -b 1 > wild_yak.Table2> wild_yak.log.

(2) Lighter, version 1.1.1, parameters used: -k 17 3000000000 -trim -t 20.

(3) Platanus, version 1.2.4, parameters used: platanus assemble -o wild_yak -f <insert size 280 bp pair-end reads> <insert size 500 bp pair-end reads> <insert size 800 bp pair-end reads> -t 30 -m 500 -tmp temp; platanus scaffold -o wild_yak -c wild_yak_contig.fa -b wild_yak_contigBubble.fa -IP1 < insert size 280 bp pair-end reads> -a1 280 -IP2 <insert size 500 bp pair-end reads> -a2 500 -IP3 < insert size 800 bp pair-end reads > -a3 800 -OP4 <insert size 2 k pair-end reads> -a4 2000 -OP5 <insert size 5 k pair-end reads> -a5 5000 -OP6 <insert size 10 k pair-end reads> -a6 10000; platanus gap_close -o wild_yak -c wild_yak_scaffold.fa -IP1 <insert size 280 bp pair-end reads> -IP2 <insert size 500 bp pair-end reads> -IP3 <insert size 800 bp pair-end reads> -OP4 < insert size 2 k pair-end reads > -OP5 <insert size 5 k pair-end reads> -OP6 <insert size 10 k pair-end reads>.

(4) Gap Closer, version 1.12, parameters used: -l 150 -t 30, in configFile: max_rd_len = 100; Paired-end libs: reverse_seq = 0, asm_flags = 3; Mate-pair libs: reverse_seq = 11, asm_flags = 2.

(5) BUSCO, version 3: mammal default parameters, mammalia_odb9.

(6) BWA, version 0.7.15-r1140: default parameters.

(7) SAMtools, version 1.3; default parameters.

(8) RepeatMasker, version 4.0.7 (with RepBase library release-20170127).

(9) RepeatModeler, RepeatModeler-open-1.0.11.

(10) TRF, version 4.09, parameters used: trf wild_yak.gapclose.fa 2 7 7 80 10 50 500 -d -h.

(11) TBLASTN, version 2.5.0, parameter used: –e 1E-5.

(12) GeneWise, version 2.4.1, parameters used: -tfor/-trev (-rfor for genes on forward strand and -trev for reverse strand) -gff.

(13) Augustus, version 3.2.3, parameter used: -species = human.

(14) Geneid, version 1.0, parameters used: -3 -P.

(15) Genemark, version 3.9, parameter used: -f gff3.

(16) Snap, version 2006-07-28, parameter used: –gff.

(17) GlimmerHMM, version 3.0.4, default parameters.

(18) EVM, version 1.1.1, default parameters.

(19) InterProScan, version 5.25-64.0, parameters: -f tsv -iprlookup -goterms -pa -t p.

Authors and Affiliations

1. State Key Laboratory of Grassland Agro-Ecosystems, School of Life Sciences, Lanzhou University, Lanzhou, China

   Yanbin Liu,Jiayu Luo,Jiajia Dou,Biyao Yan,Qingmiao Ren,Bolin Tang&Qiang Qiu

2. Research Center for Ecology and Environmental Sciences, Northwestern Polytechnical University, Xi’an, China

   Kun Wang

Contributions

Q. Qiu and K. Wang designed and supervised the project; B. Yan and Q. Ren prepared the samples; Y. Liu, J. Luo, J. Dou and B. Tang analyzed the data; Y. Liu wrote the manuscript with other authors’ help and Q. Qiu and K. Wang revised the manuscript. All authors read and approved the final manuscript.

Corresponding authors

Correspondence to Kun Wang or Qiang Qiu.

Competing interests

The authors declare no competing interests.

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as 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 images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

The Creative Commons Public Domain Dedication waiver http://creativecommons.org/publicdomain/zero/1.0/ applies to the metadata files associated with this article.

Reprints and permissions