Diversity Analysis of Chinese Tibetan Naqu Yak (Bos grunniens) Populations Using mtDNA

The domestic yak (Bos grunniens), which was domesticated from the wild yak (Bos mutus), is a long-haired, domesticated bovid that is found throughout the Himalayan region of the Indian subcontinent, the Tibetan Plateau and as far north as Mongolia and Russia. The yak plays important food, economic, and cultural roles in the Asian range and is an important domestic animal for the local people.

Mitochondrial genome sequence (mtDNA) analysis is a popular tool to estimate the phylogenic evolution and migration of humans (e.g., Schaan et al., 2017; Hernández et al., 2017), wild (e.g., Xie et al., 2017; Ming et al., 2017) and domestic animals (Zhang et al., 2016; Kim et al., 2016; Eusebi et al., 2017). In particular, previous studies have investigated yak diversity and population structures in different habitat locations in China, including Gansu (Cheng et al., 2014) and Qing-Hai (e.g., Qian et al., 2013), Xin-Jiang (e.g., Wang et al., 2013) and Tibet (e.g., Zhang et al., 2012; Song et al., 2014), using mtDNA.

Naqu, which is located in the northern part of the Tibet Autonomous Region in the hinterland of the Qinghai-Tibet Plateau, is one of the main farm sources of yaks in Tibet. This region has an average altitude of more than 4500 meters.

Studies of the domestication history of the Naqu yak and estimations of the gene flow among different ecotype yak populations in this area are important and will help improve the genetics of the local yak.

Venous blood samples were obtained from 133 individual yaks from 4 ecotype populations in the Tibetan Naqu; their geographic information is presented in Table I. The blood samples were collected in EDTA tubes and frozen at -20°C prior to extraction. Genomic DNA was isolated using standard procedures (Sambrook and Russell, 2001); the DNA quality was verified on a 1% agarose gel and quantified using a DTX microplate reader (Beckman Coulter, USA).

The high-variability region of the mitochondrial DNA control region (D-loop) was amplified using the mDNA-F (5’- GTA AAG AGC CTC ACC AGT AT -3’) and mDNA-R (5’- GTC GGG AGA CTC ATC TAG GC - 3’) from Mipam et al. (2012). PCR amplification was conducted in a PTC-100TM PCR instrument (MJ Research, Inc., MA, USA) with a total reaction volume of 50 μL containing 150 ng of DNA, 5 μL of 10× PCR standard reaction buffer, 4 μL (10 pmol/μL) of dNTPs, 2 μL (50 mmol/μL) of MgCl2, 1 μL (10 pmol/μL) of each forward and reverse primer, and 2.5 U of Taq DNA polymerase from Promega (Beijing, China). The PCR program was described in Mipam et al. (2012). The PCR products were directly sequenced using mDNA-F with the Genetic Analyzer 3130 xl (Applied Biosystems, USA).

The D-loop sequence alignments were constructed using the ClustalX software v2.0 (Larkin et al., 2007). DnaSP 5.10 (Rozas and Rozas, 1995) was used to screen haplotypes and to estimate polymorphisms and the average numbers of nucleotide differences between populations (Kxy). The best fitting model of DNA substitution for BI was obtained using jModelTest (0.1.1) (Posada, 2008). The maximum likelihood phylogenetic network of D-loops among all individuals was constructed with the MEGA (5.0) software (Tamura et al., 2011), and the bootstrap values to support the nodes of the tree were based on 1000 iterations of the heuristic search. Pairwise differences in populations (FST; Slatkin, 1995) were displayed using the Arlequin software version 3.5.1.3 (Excoffier et al., 2010). In addition, a visual haplotype phylogenetic network and frequency distribution were conducted and subjected to median-joining network analysis according to the methods of Bandelt et al. (1999) and Lyimo et al. (2014) using Network 4.1 (http://www.fluxus-engineering.com/sharenet.htm).

All of the Naqu yak sequences were aligned to the complete yak mitochondrial genome (GenBank No. KM223416) (Guangxin et al., 2016). A total of 133 sequences covering bp 1 to 359 and 15872 to 16319 (811 bp total) were obtained as reference sequences (M223416) of the yak mitochondrial D-loop region and submitted to GenBank (MG213580 to MG213712). In addition, 57 polymorphic sites, including 54 single-nucleotide polymorphisms and 3 single-nucleotide copy number variants, were identified in the D-loops of these 133 individuals.

The nucleotide polymorphisms ranged from 0.00894 (SN) to 0.01479 (JL). Tajima’s D ranged from -1.31787 (SN) to 0.26375 (NR), and the P-values of Tajima’s D in all the populations were not significant based on the chi-square test (P > 0.10; Table II).

A total of fifty-nine haplotypes were identified in the 133 individuals. The DX population carried the largest number of haplotypes (28), whereas the SN and JL populations had the smallest number of haplotypes (17). The haplotype diversity of the four populations ranged from 0.9420 (SN) to 0.9770 (JL; Table II). The phylogenetic relationship of the 59 haplotypes was constructed using the maximum likelihood method, and two haplogroups were identified from the 59 haplotypes (Fig. 1). The highest frequency haplotypes were Hapotype_2, Haplotype_4, and Haplotype_8, which were shared by the four populations. Additionally, 13 of the 59 haplotypes were shared by two or three populations, and 43 haplotypes were unique (Fig. 2).

Above diagonal, Kxy; below diagonal, FST.

In the FST analysis of the Naqu yak populations, the largest difference was found between JL and SN (FST=0.01724, P=0.25225±0.0445), and the smallest difference was found between JL and NR (FST=-0.01721, P=0.71171±0.0497 (Table III; Fig. 3). The FST distribution indicated that these four populations were separated into two groups; the first group contained JL, and the second group included SN, DX, and NR, which was consistent with the phylogenetic patterns of these four populations constructed using Kxy (Fig. 4; Table III). However, no significant divergence was observed between the populations according to the chi-square test of FST.

Recently, mitochondrial DNA polymorphisms have been widely used to estimate the gene flow and phylogenetic relationships of maternal lineages in domestic animals (e.g., Nguluma et al., 2017; Jia et al., 2017; Deng et al., 2017; Almarzook et al., 2017) and wild animals (e.g., Elsner et al., 2017; Khaire et al., 2017). The domestic yak (Bos grunniens) is a large and commercially important animal living in the Qinghai-Tibetan Plateau of China with a high-altitude climate (Guangxin et al., 2016). Various research efforts have addressed the diversity and population structure analysis of domestic yaks from different geographic distributions (e.g., Guo et al., 2006; Lai et al., 2007; Mipam et al., 2012; Huang et al., 2012).

Naqu is an important habitat of the Tibetan domestic yak in the northern part of the Tibetan Autonomous Region, as well as a junction between several Chinese yak breeding areas. Therefore, studying the diversity of different Naqu yak populations and assessing the gene flow with other yak populations from neighboring areas will contribute to our understanding of their genetic diversity status and help inform conservation policies.

First, the nucleotide polymorphisms of the D-loop region in this study ranged from 0.00894 to 0.01479; this range was larger than the number of polymorphisms detected in the Zhongdian yak (0.00534) (Tu et al., 2016) but similar to the findings for 8 other Chinese Tibetan yak populations (0.00451 to 0.01438) (Song et al., 2014). Additionally, this value was much larger than the number of polymorphisms found in the mtDNA ND6 gene (Hai et al., 2014). This observation is consistent with the hypothesis that the D-loop region contains the most polymorphisms in the mitochondrial DNA (Tsai and St John, 2016; Gao et al., 2017). Second, in comparison with previous studies, the haplotype diversity (Hi) of the 4 Naqu yak populations (0.9420 to 0.9770) was higher than the diversities among the eight different known Tibetan populations (0.827 to 0.927; Song et al., 2014), the Tianzhu white yak (Gansu, China), and the Jiulong yak (Sichuan, China) but was lower than the diversity of the Maiwa yak (Sichuan, China; Lai et al., 2005). Additionally, the haplotype numbers of the four ecotype populations ranged from 17 to 28, revealing that the Naqu yak populations carried wide and abundant genetic diversity. In particular, the finding that 72.88% of all haplotypes were unique not only indicated high diversity within each ecotype population but also inferred that their genetic characteristics followed their geographical and ecological distributions.

However, no significant divergence in FST was found among the 4 Naqu yak populations, and their genetic distances recapitulated the geographic distances between populations. This finding is consistent with the grazing features on the Tibetan Plateau, including the nomadic process, which results in gene exchange among domestic animal populations accompanied by human migration (China National Commission of Animal Genetic Resources, 2011).

In addition, according to the phylogenetic network constructed using the maximum likelihood method, two D-loop haplotype haplogroups were identified from the 59 Naqu yak haplotypes; this finding was consistent with previous studies of the two known domestic sites in the Chinese Tibetan yak (e.g., Song et al., 2014; Lai et al., 2005). Finally, no significant differences were found based on Tajima’s D test (P > 0.10) of D-loop sequences in those populations, indicating that no historical population expansion occurred in the Naqu yak population.

The Naqu yak has multiple maternal origins and high diversity within different geographic ecotype yak populations. Additionally, the yaks are highly diverse within populations, although they are not heavily divergent, because all populations share the highest frequency haplotypes. Our work underlines the importance of Naqu yak genetic diversity studies not only to obtain a better understanding of the current domestic Naqu yak status but also to enhance and provide data for the development of biological genetic conservation strategies in the Qinghai-Tibet Plateau.

This work was supported by National Technical System for Beef and Yak Industry (CARS-37), and the National Natural Science Foundation of China (No. 31172195).

Authors have declared no conflict of interest.