Phylotaxonomic position of Tianzhu white yak (Poephagus grunniens) based on nucleotide sequences of multiple subunits of cytochrome c oxidase

Abstract

To estimate the phylotaxonomic position of Tianzhu white yak (Poephagus grunniens), multiple subunit genes of cytochrome c oxidase (cox) were sequenced: 1545 bp of cox1 (JF946751), 684 bp of cox2 (JN008944) and 781 bp of cox3 (JF946752). Sequence divergences of cox genes between yaks and cattle/zebus (5.7% in cox1, 7.7% in cox2 and 6.45% in cox3) were higher than those between yaks and American bisons (2.4% in cox1, 3.0% in cox2 and 2.1% in cox3). Molecular phylogenetic analysis with these genes also found that yaks and American bisons firstly clustered in one clade, indicating there was higher genetic comparability than that of cattle/zebus. The findings sustained the idea that in the choice of nomenclature yaks belong to the subgenus of Poephagus. According to the nucleotide substitution rate of cox genes among species of Bovinae, the speculated divergence time of yaks from cattle/zebus, American bisons, European bisons, Asian buffaloes/African buffaloes was 1.05–1.50 million years ago (MYA), 2.85–3.89 MYA, 2.90–3.70 MYA, 6.85–7.50 MYA, respectively. The sequential evolution of Bovinae members could be predicted that buffaloes were first to be domesticated during the end of Miocene and the early of Pliocene. In the end of Pliocene, the Bovinae genera were evolved to Bos, Bison and Poephagus. Poephagus was the latest evolved genera among the species of Bovinae.

Keywords:

• cytochrome c oxidase (cox) genes
• taxonomy
• phylogenetic relationship
• divergence time
• Bovinae species
• Tianzhu white yak

1. Introduction

Yaks (Poephagus grunniens) are distributed across the Qinghai-Tibetan Plateau (QTP) and adjacent highlands (Zhang 1989; Wiener et al. 2003). QTP is the largest continuous high-elevation ecosystem in the world, occupying nearly 2.5 million km² of the Asian continent and reaching an average elevation of more than 3000 m as sea level (Li et al. 1995). More than 14 million yaks are currently herded there and they provide many of life's necessities (food, hids, dung fuel and transport power) for the nomadic Tibetan pastoralists living in this extremely harsh region (Wiener et al. 2003).

Bos grunniens, Bos mutus and Poephagus mutus are all species names, which have been assigned to yaks. This variety of nomenclature reflects the uncertain relationship of the yak to other members of the Bovinae. It has been ascertained that yaks are originated from China. Based on the petrifactions found in the northern China, Siberia and Alaska, researchers suggest that the forefather of the today's domestic yak and wild yak (P. mutus) was a primitive wild yak that lived about 2.50 million years ago (MYA; Hou 1991; Cai 1992). Investigation of 20 microsatellite loci indicates that the divergence time between yaks and cattle (Bos taurus) was 0.57–1.53 MYA, between yaks and the American bisons (Bison bison) 0.39–1.04 MYA (Ritz et al. 2000), showing the genetic relationship between the yaks and American bisons is closer than the yaks and cattle. More molecular evidences, either from mitochondrial DNA (mtDNA) or from nuclear DNA, have been obtained, indicating that the yak and American bison form the same clade, separated from other members of Bovinae (Tu et al. 1998; Tu et al. 2002; Hassanin & Ropiquet 2004; Li et al. 2006a; Xie et al. 2010; Zhao et al. 2011; Qiu et al. 2012; Ma et al. 2013). However, so far some envidences derived from archaeology and historical records still have not been reached a consensus about the origns of the yak and its taxonomy. The ancient Tibetans in the northern Tibet had captured and domesticated wild yaks for agricultural use during the Palaeolithic era 5000–10,000 years ago (Olsen 1990; Wiener et al. 2003; Li et al. 2008). Chinese archaeologists discovered fossils in Tibet, Sichan and Qinghai province in the 1950s indicating that the domestic yak appeared 4500 years ago (Zhang 1989; Cai 1992). These findings were consistent with the report by Li (2004) who used microsatellite markers to determine the period of yak origin.

The yak has similar morphological characteristics to American bison and cattle, so researchers focused on the phylogenetic relationship between yaks and other species of Bovinae (Li et al. 2006b; Gu et al. 2007; Xie et al. 2010). Gray (1843) reported that yaks were identified as Poephagus according to their morphological differences between genus Bos and B. bison. Olsen (1990, 1991) found that yaks have the same arrangement of premaxillaries, maxillaries and nasals as in Bison and are different from the structure in Bos. Geraads (1992), using a matrix of 57 cranial characters and 32 taxa of fossil of different Bovinae species, concluded that the yak is close to the Bison. These results were consistent with the studies of Li et al. (2005) and Hassanin and Ropiquet (2004). In contrast, according to different morphological characteristics of Bovinae species, Bohlken (1958, 1961) first grouped the yak with cattle, followed by gaur (Bos gaurus) and banteng (Bos javanicus) in one clade, which then clustered with the American bison and European bison (Bison bonasus). The same conclusion was made by Hartl et al. (1988) and Ritz et al. (2000) from biochemical genetics and microsatellite loci, respectively.

Tianzhu white yak, with 0.8 million, is one of the unique populations of Chinese yaks that centralized distribution in Tianzhu Tibetan Autonomous County, Gansu Province in China, which is located in the eastern end of the Qilian mountain and the northern edge of the QTP. The aim of the present study was to research the taxonomy of the Chinese yak and the phylogenetic relationship between the yaks and other species of Bovinae using nucleotide sequences of mitochondrial cytochrome c oxidase (cox) genes in Tianzhu white yak.

2. Materials and methods

2.1. Animals and sampling

Tianzhu white yaks were selected from Tianzhu White Yak Breeding Farm in Xidatan Village, Xidatan Town, Tianzhu Tibetan Autonomous County, Wuwei Prefecture, Gansu Province, which is located in the eastern end of the Qilian mountain and the northern edge of the QTP (102°02′–103°29′E; 36°29′–37°41′N). Any treatments on animals followed the guidelines stated in the Guide for the Care and Use of Agricultural Animals in Research and Teaching in Gansu Province, China. Ten mL of anticogulated carotid arterial blood was sampled from each adult female yak (n = 4).

2.2. Amplification of cox genes

Genomic DNA extraction was performed using animal blood genomic DNA extraction kit (Tiangen, China) according to the manufacturer's instruction. Primers were designed according to the nucleotide sequence of mitochondrial complete genome of P. grunniens (NC_006380) by Primer Premier 5.0 software and confirmed by primer-BLAST in http://www.ncbi.nlm.nih.gov/tools/primer-blast/ (Table 1). Gene amplification of cox1, cox2 and cox3 was performed with 5 min at 94°C for one cycle, followed by 35 cycles of denaturation for 45 sec at 94°C, 45 sec of annealing at 55°C and 1 min of extension at 72°C, with the final extension step for 10 min at 72°C for one cycle. Analyse 5 μL of each sample on a 0.9% agarose/EtBr gel. Each 20 μL reaction system was comprised of 1 μL genomic DNA (200 ng/μL), 1 μL forward primer (10 μM/μL) and 1 μL reverse primer (10 μM/μL), 1 μL dNTP mix (10 μM/μL), 2μL 10 × buffer, 0.5 μL Taq DNA polymerase (5 U/μL), 13.5 μL ddH₂O. Expected target gene fragment was then subcloned to T-easy vector. Positive plasmid was sequenced in Sangon Biotech (Shanghai, China) Co., Ltd. At least three positive colonies of each fragment from each sample were sequenced.

Table 1. Nucleotide sequence of primers used in the present study (forward primers are listed first).

Name of gene	Symbol	Primers (from 5′ to 3′)	Product size (bp)
Cytochrome c oxidase subunit I	cox1	GGGGTCAAACCTCTATCTTTATCTATGATCGGTGATGTAGCA	1798
Cytochrome c oxidase subunit II	cox2	TCAAGCCAACACCATAACCTCTAGTTGTGGCATGTCACCAAG	876
Cytochrome c oxidase subunit III	cox3	ATCCAAGCCTATGTATTCACCGAGCCTCATCAATAGATAGAAAC	838

2.3. Data analyses

Using complete mtDNA of sheep, Ovis aries (AF010406) as an out-group taxon, phylogenetic analysis of Tianzhu white yak was performed based on the cox1 (JF946751), cox2 (JN008944) and cox3 (JF946752) comparison with representative Bovinae species, such as P. grunniens (Qinghai black breed), B. taurus (cattle), Bos indicus (zebu), B. bison (American bison), B. bonasus (European bison), Bubalus bubalis (Asian buffalo) and Synerus caffer (African buffalo; Table 2). The analysis of sequence alignment, arrangement and the sequence divergence were carried out using DNAStar 5.02 software. The transition: transversion (Ts/Tv) bias was calculated as equation , where A, T, C, G was represented nucleotide substitution frequence, K₁ and K₂ were indicated the Ts/Tv rate ratio of purines and pyrimidines using the maximum parsimony method of MEGA5.0 software (Tamura et al. 2004, 2011) by selecting the Kimura 2-parameter model. Phylogenetic tree was constructed by Neighbour-joining method with bootstrap replication test (1000 replications). The divergence time of different species of Bovinae was estimated by the 2% molecular clock of the cox genes per million years (Birungi & Arctander 2001).

Table 2. Source of sequence for cox subunit genes in species of Bovinae and in O. aries.

Species	Cytochrome c oxidase subunit I (cox1)	Cytochrome c oxidase subunit II (cox2)	Cytochrome c oxidase subunit III (cox3)
Poephagus grunniens (Tianzhu white)	JF946751	JN008944	JF946752
Poephagus grunniens (Qinghai black)	AY684273	AY684273	AY684273
Bos taurus	NC_006853	NC_006853	NC_006853
Bos indicus	AF492350	AF492350	AF492350
Bison bison	NC_012346	BBU62568	NC_012346
Bison bonasus	NC_014044	NC_014044	NC_014044
Bubalus bubalis	NC_006295	NC_006295	NC_006295
Syncerus caffer	NC_020617	NC_020617	NC_020617
Ovis aries	AF010406	AF010406	AF010406

3. Results and discussion

3.1. Cloning of cox1, cox2 and cox3 in Tianzhu white yak

We amplified cox genes of Tianzhu white yak using genomic DNA template and primers listed in Table 1. The expected genes were subcloned and three positive colonies of each fragment from each sample were sequenced. As shown in Figure 1, about 1798 bp, 876 bp and 838 bp fragments containing cox1, cox2 and cox3 were obtained from each sample (Figure 1). Sequencing results indicated that the full length of cox1, cox2 and cox3 in Tianzhu white yak was 1545 bp, 684 bp and 781 bp, which has been conferred GenBank Accession Number with JF946751, JN008944 and JF946752, encoding 515, 228 and 260 amino acids, respectively. Nucleotide information showed that content of bases A and T, bases G and C is 58.06% and 41.94% in cox1, 62.57% and 37.43% in cox2, 55.31% and 44.69% in cox3, showing obviously higher contents of bases A and T than that of bases G and C.

Figure 1. Electrophoresis of multiple subunits of cox genes amplified by PCR. M, 1 kb DNA molecular marker; 1, cytochrome c oxidase subunit I gene (cox1); 2, cytochrome c oxidase subunit II gene (cox2); 3, cytochrome c oxidase subunit III gene (cox3). About 1798 bp, 876 bp and 838 bp fragments containing cox1, cox2 and cox3 in Tianzhu white yak were amplified by PCR using genomic DNA as templates.

3.2. Sequence divergence percentage among Bovinae species based on cox1, cox2 and cox3

The percentage divergences of cox1, cox2 and cox3 in Bovinae species were summarized in Table 3, Table 4 and Table 5. Tianzhu white yak had the highest identity to Qinghai black breed with the minimum sequence divergence of 0.5% in cox1. The following was yak and American bison (2.4%), yak and cattle/zebu (5.7%), yak and European bison (6.2–6.3%), yak and Asian buffalo (13.8–14.1%), yak and African buffalo (14.4–14.7%). Sequence divergence between cattle and zebu, cattle/zebu and American bison/European bison, cattle/zebu and Asian buffalo/African buffalo was 1.2%, 5.3–5.7% and 14.5–15.1%, respectively (Table 3). The percentage divergence of cox2 sequence of Tianzhu white yak and Qinghai black breed was 0.6%, with the minimum value, the following order of divergence between yak and American bison, yak and European bison, yak and cattle/zebu, yak and Asian buffalo, yak and African buffalo was 3.0%, 7.4%, 7.5–7.9%, 14.4–14.8%, 15.4%, respectively (Table 4). In Table 5, the nucleotide divergence of cox3 between yak and American bison, yak and European bison, yak and cattle/zebu, yak and African buffalo, yak and Asian buffalo was 2.1%, 5.8%, 6.4–6.5%, 11.9%, 15.5%, respectively. In general, the sequence divergences between yaks and cattle/zebus (5.7% in cox1, 7.7% in cox2 and 6.5% in cox3) were higher than those between yaks and American bisons (2.4% in cox1, 3.0% in cox2 and 2.1% in cox3), indicating the genetic relationship between yaks and American bisons is closer than yaks and cattle/zebus.

Table 3. Nucleotide divergence percentage of cytochrome c oxidase subunit I (cox1) gene in species of Bovinae.

	Poephagus grunniens (Tianzhu white)	Poephagus grunniens (Qinghai black)	Bos taurus	Bos indicus	Bison bison	Bison bonasus	Bubalus bubalis	Syncerus caffer
Poephagus grunniens (Tianzhu white)
Poephagus grunniens (Qinghai black)	0.5
Bos Taurus	5.7	5.7
Bos indicus	5.7	5.7	1.2
Bison bison	2.4	2.4	5.4	5.5
Bison bonasus	6.3	6.2	5.7	5.3	7.0
Bubalus bubalis	13.8	14.1	15.1	14.9	14.0	14.4
Syncerus caffer	14.4	14.7	14.5	14.5	14.5	13.9	11.5
Ovis aries	15.6	15.6	15.3	15.7	14.9	14.9	15.7	15.2

Table 4. Nucleotide divergence percentage of cytochrome c oxidase subunit II (cox2) in species of Bovinae.

	Poephagus grunniens (Tianzhu white)	Poephagus grunniens (Qinghai black)	Bos taurus	Bos indicus	Bison bison	Bison bonasus	Bubalus bubalis	Syncerus caffer
Poephagus grunniens (Tianzhu white)
Poephagus grunniens (Qinghai black)	0.6
Bos Taurus	7.9	7.8
Bos indicus	7.5	7.5	0.9
Bison bison	3.0	3.0	7.2	7.2
Bison bonasus	7.4	7.4	5.1	4.7	6.4
Bubalus bubalis	14.8	14.4	14.8	14.4	14.0	12.8
Syncerus caffer	15.4	15.4	15.0	15.4	16.2	13.9	10.3
Ovis aries	18.0	18.2	16.0	15.6	18.4	15.9	18.0	14.5

Table 5. Nucleotide divergence percentage of cytochrome c oxidase subunit III (cox3) in species of Bovinae.

	Poephagus grunniens (Tianzhu white)	Poephagus grunniens (Qinghai black)	Bos taurus	Bos indicus	Bison bison	Bison bonasus	Bubalus bubalis	Syncerus caffer
Poephagus grunniens (Tianzhu white)
Poephagus grunniens (Qinghai black)	0
Bos Taurus	6.4	6.4
Bos indicus	6.5	6.5	0.9
Bison bison	2.1	2.1	6.4	6.5
Bison bonasus	5.8	5.8	5.1	5.5	6.4
Bubalus bubalis	15.5	15.5	14.6	15.1	15.0	15.0
Syncerus caffer	11.9	11.9	13.9	13.6	12.7	13.0	10.1
Ovis aries	15.6	15.5	15.3	15.6	16.0	15.7	18.7	16.7

3.3. Phylogenetic analysis among Bovinae species based on cox1, cox2 and cox3

Based on the nucleotide sequence of 1545 bp of cox1, 684 bp of cox2 and 838 bp of cox3 in Tianzhou white yak, a phylogenetic tree of Bovinae was constructed using the Neighbour-joining tree method and listed as Figure 2A, Figure 2B and Figure 2C. Bovinae was clustered in one group, while the out-group of sheep was in another group. Among eight species of Bovinae, yaks were clustered with American bisons, but not with cattle/zebus.

Figure 2. Phylogenetic trees of species of Bovinae based on cox genes of Tianzhu white yak by Neighbour-joining method. Figure 2A, Figure 2B and Figure 2C is derived from cox1, cox2 and cox3, respectively. The bootstrap consensus trees inferred from 1000 replicates are taken to represent the evolutionary history of the taxa analysed. Branches corresponding to partitions reproduced in less than 50% bootstrap replicates are collapsed. The percentage of replicate trees in which the associated taxa clustered together in the bootstrap test (1000 replicates) are shown next to the branches. The evolutionary distances were computed using the p-distance method and are in the units of the number of base differences per site. The analysis involved nine nucleotide sequences. Codon positions included were 1st + 2nd + 3rd + Noncoding. All positions containing gaps and missing data were eliminated.

3.4. Estimation of divergence time among different Bovinae species

The Ts/Tv rate of cox1, cox2 and cox3 among the different species of Bovinae was 9.042, 7.646 and 6.052, higher than the critical value 2.0 (Knight & Mindell 1993) and showed higher transition bias. Therefore, the molecular clock of the cox1, cox2 and cox3 gene sequence was 2% per million years in Bovinae and chosen to estimate the divergence time between different species of Bovinae. Based on cox1, the divergence time between Tianzhu white yak and Qinghai black yak, yak and American bison, yak and cattle/zebu, yak and European bison, yak and buffalo was 0.25 MYA, 1.20 MYA, 2.85 MYA, 3.15 MYA and 7.13 MYA, respectively. Based on cox2, the divergence time between Tianzhu white yak and Qinghai black yak, yak and American bison, yak and European bison, yak and cattle/zebu, yak and buffalos was 0.30 MYA, 1.50 MYA, 3.70 MYA, 3.84 MYA and 7.50 MYA, respectively. Based on cox3, the divergence time between yak and American bison, yak and European bison, yak and cattle/zebu, yak and African buffalo, yak and Asian buffalo was 1.05 MYA, 2.90 MYA, 3.25 MYA, 5.95 MYA and 7.75 MYA, respectively. These results indicated that the divergence time between yak and American bison was much later than other species of Bovinae.

The earliest report of Linnaeus (1766) listed yak as Bos grunniens of Bos, which was consistent with Bohlken's results from morphological characteristics that yak and cattle first grouped, followed by gaur and banteng in one clade, then clustered with the American bison and European bison (Bohlken 1958, 1961). Gray (1843) identified yak as Poephagus according to the morphological differences between genus Bos and B. bison, which was consistent with the opinions of Groves (1981), Olsen (1990, 1991), Geraads (1992) from the head bone features and archaeological study, and Li et al. (2006a) from the molecular research. Since the late 1980s, the sequences of mtDNA of yak have been used for the phylogenetic studies of Bovinae. Miyamoto et al. (1989) sequenced the 12S rRNA and three tRNA genes, and a 247 bp partial hypervariant D-loop fragment of four taxa in the Bovinae and the results showed a similar topology of phylogeny with Groves (1981), the yak grouped with the B. bison first, with an average divergence of 2.6% for the conservative rRNA/tRNA genes and 9.1% for the D-loop fragment and then followed by B. taurus (Kraus et al. 1992). Ward et al. (1999) using a partial mtDNA control region of 667 base pairs, found that the percentage nucleotide divergence of B. indicus and B. taurus from the yak was 24.23% and 29.53%, respectively, but that of American bison and European bison from the yak was only 12.28% and 16.59%, respectively. Similar results have been obtained by Hassanin and Douzery (1999), Li et al. (2006a), Xie et al. (2010) from complete sequence of cytochrome b (cytb), Zhao et al. (2011) from cox3, Zhang et al. (2012) from mtDNA D-loop, Hassanin and Douzery (1999) from mtDNA 12S rRNA, aromatase cytochrome P-45. Within the nuclear genome alone, Wall et al. (1992) confirmed the position of yak as more closely related to the bison and wisent than to B. taurus and B. indicus as earlier suggested by Groves (1981). Hassanin and Douzery (1999) sequenced the promoter segment of the lactoferrin-encoding genes of nuclear genomic DNA in the yak alongside with other Bovidae species, showing the yak clustered with the American bison first and then grouped with B. taurus in the phylogenetic tree. More similar molecular evidences have been made by Ritz et al. (2000) from 20 bovine microsatellite markers, Buntjer et al. (2002) from 361 fingerprinting Bovini markers of AFLP, Fan et al. (2000) from exon 4 of kappa casein, Li et al. (2005) from exon 2 of MHC-DRB3.

At present study, the average percent nucleotide divergence between yaks (Tianzhu white, Qinghai black) and cattle/zebus was 5.7% in cox1, 7.7% in cox2 and 6.45% in cox3, which was considerably higher than that between yaks and American bisons (2.4% in cox1, 3.0% in cox2 and 2.1% in cox3), indicating that genetic identity between yaks and American bisons is higher than that between yaks and cattle/zebus. Phylogenetic analysis based on cox1, cox2 and cox3 also showed that yaks were clustered with American bisons, but not with cattle/zebus. According to the 2% per million years mutation rate of cox genes, we speculated the divergence time between yaks and American bisons was 1.20 MYA based on cox1, 1.50 MYA based on cox2 and 1.05 MYA based on cox3, which is later than that between yaks and cattle/zebus (2.85 MYA in cox1, 3.70 MYA in cox2 and 3.25 MYA in cox3, respectively). These findings showed that the yak has the closest genetic relationship related to American bison. Therefore, it seems clear that for the domestic yak, the subgenus, Poephagus, seems more appropriate than Bos. This finding is in agreement with the results made by Wiener et al. (2003), Li et al. (2005), Li et al. (2006a).

Using subunit genes of cytochrome c oxidase, we estimated that the divergence time of buffaloes from other species of Bovinae was 7.13–7.75 MYA, whereas the yaks from cattle/zebus, B. bisons from cattle/zebus was 2.85–3.84 MYA and 1.05–1.50 MYA, respectively, which is in agreement with the report of Xie (1985), Li et al. (2006a). The divergence time of cattle from zebus was 0.45–0.60 MYA, consenting with the result of Ritz et al. (2000). Similar to the result of Li et al. (2006a), the divergence time between yaks and American bisons was 1.7–2.05 MYA. Combination with paleontological evidence (Osborn 1990) and morphological characteristics (Xie 1985), we estimated that buffalos were first to be divided into Asian buffalos and African buffalos among the species of Bovinae during the end of Miocene (23.00–5.3 MYA) and the early of Pliocene (5.2–1.8 MYA). In the end of Pliocene, the genera of Bovinae were evolved to Bos, Bison and Poephagus. Poephagus, which branched off from the middle of Pleiocene, was the latest evolved genera among the species of Bovinae.

4. Conclusion

We cloned and sequenced multiple subunit genes of cox of Tianzhu white yak (P. grunniens): 1545 bp of cox1 (JF946751), 684 bp of cox2 (JN008944) and 781 bp of cox3 (JF946752). Sequence divergences between yaks and cattle/zebus (5.7% in cox1, 7.7% in cox2 and 6.45% in cox3) were higher than those between yaks and American bisons (2.4% in cox1, 3.0% in cox2 and 2.1% in cox3). Molecular phylogenetic analysis with these genes also found that yaks and American bisons firstly clustered in one clade, indicating there was higher genetic comparability than that of cattle/zebus. We concluded that in the choice of nomenclature yaks belong to the subgenus of Poephagus. We also speculated that buffaloes were first to be divided into Asian buffalos and African buffalos among the species of Bovinae during the end of Miocene and the early of Pliocene. In the end of Pliocene, the Bovinae genera were evolved to Bos, Bison and Poephagus. Poephagus, which branched off from the middle of Pleiocene, was the latest evolved genera among the species of Bovinae.

Funding

The work was partially financed by the project of National Natural Science Foundation of China, NSFC [grant number 31260533] and the Programme for Changjiang Scholars and Innovative Research Team in University of Ministry of Education, China [grant number IRT13091].

Additional information

Funding

Funding: The work was partially financed by the project of National Natural Science Foundation of China, NSFC [grant number 31260533] and the Programme for Changjiang Scholars and Innovative Research Team in University of Ministry of Education, China [grant number IRT13091].