Skip to main content
Download PDF

Genetic diversity, phylogeography, and maternal origin of yak (Bos grunniens)

BMC Genomics volume 25, Article number: 481 (2024) Cite this article

Abstract

Background

There is no consensus as to the origin of the domestic yak (Bos grunniens). Previous studies on yak mitochondria mainly focused on mitochondrial displacement loop (D-loop), a region with low phylogenetic resolution. Here, we analyzed the entire mitochondrial genomes of 509 yaks to obtain greater phylogenetic resolution and a comprehensive picture of geographical diversity.

Results

A total of 278 haplotypes were defined in 509 yaks from 21 yak breeds. Among them, 28 haplotypes were shared by different varieties, and 250 haplotypes were unique to specific varieties. The overall haplotype diversity and nucleotide diversity of yak were 0.979 ± 0.0039 and 0.00237 ± 0.00076, respectively. Phylogenetic tree and network analysis showed that yak had three highly differentiated genetic branches with high support rate. The differentiation time of clades I and II were about 0.4328 Ma, and the differentiation time of clades (I and II) and III were 0.5654 Ma. Yushu yak is shared by all haplogroups. Most (94.70%) of the genetic variation occurred within populations, and only 5.30% of the genetic variation occurred between populations. The classification showed that yaks and wild yaks were first clustered together, and yaks were clustered with American bison as a whole. Altitude had the highest impact on the distribution of yaks.

Conclusions

Yaks have high genetic diversity and yak populations have experienced population expansion and lack obvious phylogeographic structure. During the glacial period, yaks had at least three or more glacial refugia.

Peer Review reports

Introduction

The Qinghai-Tibet Plateau (QTP), one of the largest and youngest plateaus in the world, was formed around 40 million years ago (Ma) following the collision of the Indian tectonic plate with the Asian plate through several uplift events [1]. A large number of endemic species have appeared in the QTP and adjacent areas [2], due to its unique ecological environment. Yak is one of the representative species of QTP. At present, there are about 17.5 million yaks in the world, of which 94.4% are distributed in China [3]. Contemporary highland pastoralists rely on the strength and hardiness of domestic yak for transportation across vast mountainous terrain and for supplies of milk, meat, fiber, and dung for fuel [4]. Yaks play a vital role in socioeconomic development, pasture ecosystem maintenance, and agricultural biodiversity conservation in the QTP region [5]. Hence, yaks are referred to as “all-round animals” [6]. Morphological data suggest that yaks outside China originated from the Chinese yak [7]. According to previous studies, yaks were first domesticated in Tibet [8]. The combined archaeological and mitochondrial DNA (mtDNA) evidence suggests that Qinghai is one of the places where the yak either originated or was domesticated [7, 9, 10]. Qiu et al. re-sequenced the whole yak genome to find that the domestication of yak occurred 7 300 years ago [11]. However, the available genetic data do not provide a definitive conclusion and it is not known whether yak domestication occurred as a single event or multiple events in a single wild gene pool [7].

Domestication of animals is one of the major achievements of human civilization [12], although there had been many studies on yak ancestry, origin, and domestication, the answers to these questions are not clear. Studies have suggested an association between yak fossils and early human activity in Tibet [13], suggesting that yaks were first domesticated in Tibet [8]. The combined archaeological and mtDNA evidence suggests that Qinghai is one of the places where the yak either originated or was domesticated [7, 9, 10]. Qiu et al. re-sequenced the whole yak genome to find that the domestication of yak occurred 7 300 years ago [11]. However, the available genetic data do not provide a definitive conclusion and it is not known whether yak domestication occurred as a single event or multiple events in a single wild gene pool [7]. mtDNA is a good molecular marker for studying animal origins, evolution, classification, and population genetic diversity [14]. In recent years, the mtDNA D-Loop region sequence has been widely used to evaluate the origin, domestication, and genetic diversity of yaks [15, 16], cytochrome b (Cytb) is also a commonly used marker gene for studying the molecular genetic diversity of populations. Qi et al. [17] conducted a cluster analysis of the mtDNA D-Loop region and Cytb gene of 428 yak individuals from 29 yak populations in China and its surrounding countries. Phylogenetic analysis showed that 29 yak populations were clustered into three categories. In addition to mtDNA D-Loop region and Cytb, other regions can also be used for genetic diversity and phylogenetic analysis. Zhao et al. [18] determined the mtDNA Cytochrome c oxidase polypeptide III (COIII) sequence of 111 yak individuals from 11 yak populations in Tibet, indicating that Tibetan yaks have rich genetic diversity. Recently, researchers have discussed the importance of conducting complete mtDNA sequencing, because such analysis can produce detailed genetic maps when the sample size is large enough [19]. Wang et al. [10] sequenced the whole mitochondrial genome of yak, and used the 10 710 bp protein coding sequence except NADH dehydrogenase6 (ND6) for phylogenetic analysis. It was found that wild yaks were divided into three categories and domestic yaks were divided into two categories. It is speculated that there may be three maternal origins of wild yaks and two maternal origins of domestic yaks. However, most studies to date have been limited to the mitochondrial D-loop region and Cytb gene and have failed to clearly distinguish some important ancient clades in the domestic yak [15, 20, 21], the D-loop region is highly variable and information-rich in determining intraspecific diversity but often has parallel mutations [22, 23]. Recent studies have emphasized the importance of complete mtDNA sequencing [19], as this allows the investigation of 18 times as many sites as the D-loop region [10], and more detailed information can be obtained.

At present, studies using the complete mitochondrial genome of the yak have mainly focused on a single yak genetic resource [24,25,26]. Also, there are few studies on the genetic evolution of the yak that have used the complete mitochondrial genome. Wild yaks that share a common ancestor with domestic yaks are still in existence, making the yak an excellent model species for studying the domestication of large animals. In our previous studies [27, 28] and here, we collected samples of genetic resources from every yak breed (Fig. 1) in China to ensure a large sample size for a detailed genetic map. The genetic diversity and phylogenetic structure of the yak were analyzed by the sequencing of the complete mitochondrial genome sequencing of these samples, which also provides a foundation for the effective protection and utilization of yak genetic resources.

Fig. 1
figure 1

Collection locations of experimental samples. The details of the yak population represented at each sampling site are provided in Supplementary Table S1

Full size image

Results

Identification and analysis of haplotypes

In total, 278 haplotypes were defined from 509 yaks. Of these, 28 haplotypes were shared by different breeds, and 250 haplotypes were unique to a specific breed. Among the 278 haplotypes, the H8 haplotype was the most common (found 68 times) and shared by 18 yak breeds except wild yak, Tianzhu white yak, and Pamir yak. In total, 23 haplotypes were identified in 25 wild yaks, of which only the H28 haplotype was shared by the wild yak and Changtai yak; the other 22 haplotypes were unique to wild yaks. In total, 11 haplotypes were defined in Tianzhu white yak; the H20 haplotype was shared by the Tianzhu white, Changtai, Huanhu, and Qinghai Plateau yaks, while the other 10 haplotypes were unique to the Tianzhu white yak. Twenty-two haplotypes were defined in 25 Pamir yaks; only the H150 haplotype was shared by the Pamir yak and Sibu yak, while the other 21 haplotypes were unique to the Pamir yak. Wild yak (96.65%), Pamir yak (95.45%), Tianzhu white yak (90.91%), Jiulong yak (89.47%), Yushu yak (81.25%), and Xueduo yak (80.95%) exhibited higher specific proportions of haplotypes (Table 1); for detailed information please see Supplementary Table S2.

Table 1 Haplotype statistics of yak breeds/populations
Full size table

Phylogenetic analysis of yaks

Phylogenetic analysis was performed using 278 haplotypes (with Bison bison [GU947006.1, GU946996.1] as outgroups). The basic topological structures of the maximum likelihood (ML) and MrBayes phylogenetic trees were the same, showing three highly differentiated genetic clades with high support rates (Fig. 2). Most of the haplotypes were distributed within clades I and II. The yak phylogenetic tree contained six haplogroups (A-F) forming two clades, with clade I containing haplogroups A, C, E, and F and clade II containing haplogroups B and D; A-C were the three main haplogroups. Haplogroup A included all the yak populations, while haplogroup B included all populations except Tianzhu white yak and haplogroup C included all populations except Qinghai plateau yak, Tianzhu white yak, and Pamir yak. Haplogroups D-F contained only a few yak breeds. Among which the haplogroup D only contained the Yushu yak, wild yak, Sunan yak and Pamir yak; Haplogroup E only contained the Yushu yak; Haplogroup F only contained the Yushu yak and wild yak. Among all yak populations, the Yushu yak was common to all haplogroups.

Fig. 2
figure 2

Phylogenetic tree of 278 haplotypes in the yak. Clades I, II, and III represent the three clades of yak. A, B, C, D, E, and F represent the six haplogroups of yak

Full size image

Combined with the geographical distribution of yak populations, a haplotype network map was constructed based on yak mtDNA (Fig. 3). Consistent with the phylogenetic tree, the haplotype network diagram also revealed that yaks were divided into three clades and six haplogroups. The three main haplogroups A-C were distributed in a star-shaped radial pattern; the H8 and H36 haplotypes were shared by multiple individuals and located in the center of the star-shaped radial. The yaks in haplogroup A were distributed in all taxa. The yaks in haplogroup B were distributed in all yak distribution areas except Tianzhu in Gansu. The yaks in haplogroup C were distributed in all yak distribution areas except the Pamir region, Tianzhu in Gansu, and parts of Qinghai. Only the wild and Yushu yaks were widely distributed with the Yushu yak distributed in all haplogroups.

Fig. 3
figure 3

Network diagram of haplotypes. (a) The total network of mtDNA of 509 individuals, with different colors indicating yaks from different provinces. The network diagram of haplogroups (b) A, (c) B, and (d) C

Full size image

Genetic diversity analysis of yak mtDNA

Higher haplotype diversity (Hd) and nucleotide diversity (Pi) values are indicative of greater genetic diversity in a population. Genetic diversity analysis showed that the haplotype and nucleotide diversities of the complete mitochondrial genome sequences of 509 individuals from 21 yak breeds/populations were 0.979 ± 0.0039 and 0.00237 ± 0.00076, respectively, indicating a high overall genetic diversity in the yak. The haplotype diversity of wild yak (0.993 ± 0.013), Sunan yak (0.996 ± 0.015), Pamir yak (0.990 ± 0.014), and Xueduo yak (0.992 ± 0.015) was relatively high, while Muli yak (0.810 ± 0.063) and Tianzhu white yak (0.830 ± 0.068) exhibited the lowest haplotype diversity. The nucleotide diversity of wild yak (0.00352 ± 0.00145), Leiwuqi yak (0.00319 ± 0.00065), Pamir yak (0.00309 ± 0.00018), and Tibet Gaoshan yak (0.00309 ± 0.00067) was relatively high; the lowest nucleotide diversity was of Tianzhu white yak (0.00034 ± 0.00013). The comprehensive analysis of haplotype and nucleotide diversities revealed that the genetic diversity of the wild yak was the highest, and that of the Tianzhu white yak was the lowest. Additional details including variable loci (S), Hd, and Pi are shown in Table 2.

Table 2 Genetic structure and diversity of yaks
Full size table

Population dynamics analysis

By calculating the historical population dynamics of each haplogroup of yak (Table 3), the Tajima’s D value (-1.230) of the Total group was less than 0, and P > 0.05. Fu and Li’s D test (-7.262 (P < 0.05)) and Fu and Li’s F test (-4.543 (P < 0.05)) showed that the neutral test results of the Total group were contradictory. The mismatch analysis showed that the SSD and Hrag values of the yaks in the Total group were 0.0084 (P > 0.05) and 0.0038 (P > 0.05), respectively, this results indicated that the yak population had experienced the expansion. In haplogroups A and C, the results of neutral test showed P < 0.05, and the results of mismatch distribution showed P > 0.05, indicating that haplogroups A and C had experienced population expansion. In the haplogroup B, the results of the neutral test were P > 0.05, and the results of the mismatch distribution were P > 0.05, resulting in a contradiction in the calculation results of the historical population dynamics of the haplogroup B. The historical population dynamics of yaks (Supplementary Fig. S1) were analyzed by Bayesian Skyline Plot (BSP), revealing that each haplogroup of yaks experienced large-scale expansion. Due to the small number of individuals in the haplogroups D, E and F, the neutrality test and mismatch distribution analysis were not performed. We used the PermutCpSSR-2.0 software to analyze the Nst and Gst values at the population level. The results showed that Nst = 0.05225 > Gst = 0.03993 (P > 0.05). Analysis of Molecular Variance (AMOVA) analysis of the yak whole-mtDNA genomes showed that most (94.70%) of the genetic variation in yaks occurred within populations with only 5.30% observed between populations. This indicated the lack of obvious phylogeographical structure in yaks.

Table 3 Analysis of historical population dynamics
Full size table

Estimation of differentiation time

Based on the Bayesian method, the differentiation time of the respective yak clades was calculated. The differentiation time of the two major clades (I and II) was about 0.4328 Ma with a 95% highest posterior density (HPD) of 0.3218–0.5326 Ma. The differentiation time of clades I and II and clade III was about 0.5654 Ma, with a 95% HPD of 0.4283–0.7162 Ma.

Taxonomic status of the yak in the Bovidae

The Caprinae, including Ovis ammon (NC_047196.1) and Ovis aries (NC_001941.1), were selected as an outgroup for phylogenetic analysis. The results showed that Bos grunniens and Bos mutus were located together, and yaks as a whole clustered with Bison bison. The overall clustering relationship is shown in Fig. 4: (((((((Bos grunniens + Bos mutus) + Bison bison) + (Bos gaurus + Bos javanicus)) + ((Bos taurus + Bos primigenius + Bos indicus) + Bison bonasus)) + (((Bubalus bubalis + Bubalus arnee) + (Bubalus depressicornis + Bubalus quarlesi)) + (Syncerus caffer)) + Pseudoryx nghetinhensis) + ((Boselaphus tragocamelus + Tetracerus quadricornis) + Tragelaphus spekii)) + (Ovis ammon + Ovis aries)). Based on the differentiation time of the respective Bovidae species, the differentiation time of Bos grunniens and Bison bison was 2.2011 Ma, the differentiation time of Bos grunniens, Bos gaurus, and Bos javanicus was 3.9735 Ma, and the differentiation time of Bos grunniens and Bos taurus was 4.7392 Ma.

Fig. 4
figure 4

Taxonomic status of the yak in the Bovidae family

Full size image

Population distribution dynamics of the yak

The dynamic distribution of the yak in different periods was simulated by the MaxEnt model based on 11 selected environmental factors. Similar average Area under the curves (AUCs) of training (0.918) and test (0.873) data in different periods, indicated high accuracy of MaxEnt model simulation (Supplementary Table S3). The contribution rate of each environmental factor was tested by the knife-cut method. Elevation (74.6%) and annual temperature range (0.0009%) contributed the highest and lowest, respectively (Supplementary Table S4). A dynamic distribution map of yaks from the last interglacial (LIG) to 2100 was constructed (Fig. 5), showing that yaks were mainly located at the edge of the QTP during the last glacial maximum (LGM).

Fig. 5
figure 5

Dynamic distribution of yaks from the last interglacial period to 2100. LIG: last interglacial. LGM: last glacial maximum. MH: Mid-Holocene

Full size image

Discussion

In the neutral test, values of Tajima’s D > 0 and P < 0.05 indicated that the population underwent a bottleneck effect and equilibrium selection, while Tajima’s D < 0 and P < 0.05 indicated that the population experienced the expansion and directional selection [29]. Fu and Li’s D*&F* indices can be used for neutral test [30]. However, the historical dynamics of the population cannot be inferred only by the Tajima’s D value, which is usually analyzed by a combination of neutral test and mismatch distribution. In this experiment, the neutral test results of the Total group were contradictory, and the mismatch analysis showed that the Total group experienced population expansion. In haplogroup B, the neutral test results are contrary to the mismatch distribution results. Only when the neutral test was negative and P < 0.05, it indicated that the population was significantly deviated from the neutral mutation, indicating that it was intervened by artificial or natural selection, and the curve of mismatch distribution analysis was unimodal distribution, indicating that the population was in an expanding state. Due to the neutral test results of Total group and haplogroup B are contradictory and the neutral test results are contrary to the mismatch distribution results, which makes it difficult to judge whether the population has expanded during evolution. The result of mismatch distribution is based on the ideal state, which is limited in practical applications. In fact, the historical population dynamics are often more complex than the parameter models involved in these methods. The BSP method is based on the clustering theory and is used to quantify the relationship between gene sequence lineages and population geographic history. Using molecular clock or fossil correction, with the help of BEAST series software, the Markov Chain Monte Carlo (MCMC) algorithm based on Markov chain is used to calculate the change of effective population size with time. Especially for the analysis of different genes and a small number of individuals, this method can better estimate the effective population size. Therefore, the neutral test results deviate from 0 and P < 0.05 is only a prerequisite. If the neutral test results are inconsistent with the BSP results, the BSP results (Supplementary Fig. S1) should be used.

Genetic variation in yaks originated largely within populations without any obvious geographical structure. A study in goats also showed no obvious phylogeographical structure between the highly differentiated genetic clades and haplogroups due to the large-scale migration of goats around the world [31]. Notably, the domestication of yaks started 7300 years ago [11]. After domestication, yaks migrated with herdsmen on a large scale, enabling integration with yak populations in distant geographical regions which would account for the lack of an obvious pedigree geographical structure in yaks.

Genetic diversity is a central facet of biological diversity [32]. High Hd and Pi values indicated high genetic diversity in wild, Pamir, Xueduo, Qinghai Plateau, Yushu, Tibet Gaoshan, Sunan, and Jiulong yaks. The Gannan, Niangya, Bazhou, and Sibu yaks had high Hd but low Pi values. While a single base mutation can generate a new haplotype, this has little impact on nucleotide diversity. Compared with haplotype diversity, increases in nucleotide diversity take longer. Therefore, the bottleneck effect caused by repeated changes in the environment and the rapid population expansion and variation accumulation after the bottleneck effect can lead to high Hd and low Pi values [33]. The low Hd and high Pi values of Leiwuqi and Muli yaks can be attributed to selective pressure from a new environment after migration or the coming into contact of two relatively independent populations or simply to an insufficient number of samples [33]. This requires further analysis. Zhongdian and Tianzhu white yaks exhibited low Hd and Pi values, indicating the influence of the founder effect, i.e., the re-establishment of a new population from a few individuals. The numbers in this population would increase but without an increase in genetic diversity [33]. Low levels of genetic diversity lead to inbreeding and decreased population fitness [34]. Which is consistent with the artificial selection process seen in the Tianzhu white yak. For 130 years, the Tianzhu white yak has been bred by strict selection of fur color [35], and this strict artificial selection and control have led to low gene flow between the Tianzhu white yak and other yak populations, resulting in a high degree of genetic differentiation, consistent with the earlier findings that the Tianzhu white yak differed from other domestic yak populations [11].

Previous studies based on mtDNA D-loop fragments showed that yaks were divided into two clades (I and II), with a differentiation time of 100 000–130 000 years ago [7]. The differentiation time of the main clades in the yak phylogenetic tree was also calculated based on the third codon of the protein-coding gene of mtDNA. This new method and new fossil marker (2.5 Ma) revealed the differentiation time of three yak clades to be between 420 000 and 580 000 years ago [10]. Analysis of the whole mtDNA also revealed three main clades in the yak evolutionary tree. Differentiation between clades I and II occurred 0.4328 Ma with a 95% HPD of 0.3218–0.5326 Ma. Differentiation between clade III and clades (I, II) occurred before 0.5654 Ma with a 95% HPD of 0.4283–0.7162 Ma. This is in agreement with previous findings [10]. In addition, our estimated differentiation time is also consistent with the records of glacial events in the middle and late Pleistocene in the Tibetan Plateau [36], suggesting that glacial activity in the middle and late Pleistocene may have triggered yak migration and therefore genetic differentiation.

There have been three great glacial epochs in the evolutionary history of the earth. The Quaternary glacial epoch was the most recent in geological history and had a major impact on modern biology [37]. Due to specific buffering environmental characteristics, the glacial refugia contained unique genetic lineages during a series of climatic fluctuations occurring during the Tertiary and Quaternary epochs [38], allowing animals and plants to escape from the harsh climatic conditions of the glacial epoch [39]. Glacial refuges are the starting point for the post-glacial redistribution of species after deglaciation [40]. Large glacial refugia may last for hundreds of thousands of years or even longer, and, therefore, the isolation of glacial refugia may have accelerated the differentiation of a species’ populations [41]. The Hengduan Mountains in the southeastern part of the QTP rose rapidly between the Late Miocene and Late Pliocene [42]. The QTP was never completely covered by glaciers during the Quaternary [43] and is one of the most important biodiversity research hotspots in the world [44]. The QTP also formed an important refuge and place of origin for species during the glacial epoch, resulting in rich species diversity and unique geological characteristics [45, 46]. Studies on Juniperus przewalskii [47], Metagentiana striata [48], and Pedicularis longiflora [49] indicated that some species may have retreated to refuges on the edge of the QTP during the glacial epoch and recolonized the plateau and adjacent highlands at the end of epoch. In addition, some glacial refugia on the QTP supported the survival of plant species, such as Hippophae rhamnoides [50] and Spiraea alpina [51], during climate change [1, 52]. Research on Aconitum gymnandrum [53], Hippophae tibetana [54], Rhodiola alsia [55], and Rhodiola chrysanthemifolia [56] has shown that there were several miniature refugia on the QTP. The multiple glacier refuge scenario fits well with the hypothesis of multiregional and multiscale glaciation during the Pleistocene [49, 57]. The mtDNA phylogeny analysis of yak revealed three clades. Wang et al. [53] found four Aconitum gymnandrum populations that possibly evolved from four independent glacial refugia during the LGM. The single refuge hypothesis generally advocates a network of stellate haplotypes in the entire population [58]. The network diagram of the yak mitochondrial genome (Fig. 3) indicated at least three stellate haplotype network structures in the yak, suggesting at least three glacial refuges in the yak evolutionary history. The MaxEnt model was used to predict the dynamic distribution of yaks from the LIG to the present time (Fig. 5). Yaks were mainly but not entirely located on the edge of the QTP during the LGM, suggesting that the yak refuges during the glacial epoch were associated mainly with the marginal areas of the QTP. In Alopex lagopus, the population distribution did not change to follow available habitats during the post-glacial contraction phase but instead went extinct [59]. This suggests that arctic species did not respond to climate change in the same way as temperate species, i.e., their distribution ranges contracted rather than expanded during interglacial warming. It is still unclear whether species migrated to the refuge to cope with climate change or populations outside the refuges were directly made extinct [58]. With the lack of detailed information about some sampling points, yak fossil information, and other data, it is impossible to assess the location of glacial refugia during the glacial epoch, or whether yaks outside the glacial refugia underwent direct extinction.

The yak classification in the Bovidae is very different. Linnaeus placed the yak together with Bos taurus and Bos indicus in the Bos genus. Based on its morphological and skeletal differences from other cattle species, Gray, Olsen, and Geraads classified the yak into an independent yak genus [60]. This study explored the taxonomic status of yak from the direction of maternal inheritance. The representative mtDNA of yaks from three clades and six haplogroups were selected and those of other cattle were downloaded from the GenBank database to determine the phylogenetic relationships among bovids. Each species was represented by at least two individuals. The clustering results showed that the individuals of each species were first clustered into one category; the yaks as a whole were clustered together with Bison bison, and then with other cattle. This is consistent with a previous report [60] on the Cytb gene of domestic yak mtDNA and another report [61] on the mtDNA of wild yak. The time of separation between the yak and other bovid species was calculated by the Bayesian method, finding a value of 2.2011 Ma for the separation between yak and Bison bison, which was consistent with the phylogenetic tree structure. This study supports the classification of Li [60] and Zhong [61] who classified the yak into a separate genus that included two species, namely, domestic and wild yak.

Conclusions

In conclusion, yaks generally showed high genetic diversity, among which wild yaks had the highest genetic diversity and Tianzhu white yaks had the lowest genetic diversity. The mtDNA genome-wide genetic variation of yak mainly occurs within the population and lacks obvious phylogeographic structure. Phylogenetic analysis revealed that yak had three highly differentiated genetic branches with high support rate. The differentiation time of clades I and II were about 0.4328 Ma, and the differentiation time of clades (I and II) and III were 0.5654 Ma. Yak contains 6 haplogroups. In all yak populations, Yushu yak is distributed in all haplogroups, and the three major haplogroups A-C are distributed in a star-shaped radial pattern. Yaks have experienced population expansion. At the last LGM, the yak’s glacial refugia was mainly located at the edge of QTP. The effects of altitude and annual temperature range on the dynamic distribution of yaks were the highest and lowest, respectively. The classification showed that yaks and wild yaks were first clustered together, and yaks were clustered with American bison as a whole.

Methods

Animals and sample collection

The mitochondrial genomes of 372 yaks were sequenced. These yaks were from 15 breeds identified by the National Livestock and Poultry Genetic Resources Commission. Supplementary Table S1 lists information on the breed, numbers, mitochondrial genome accession numbers and geographical coordinates of the sampling points of the yaks. Blood samples were collected from 20 to 26 yaks from each breed/population. There was no sex restriction, and the yaks were three to eight years old without any genetic relationships. Blood samples were collected from the jugular veins of the yaks into EDTA anticoagulant tubes, immediately stored in the vehicle refrigerator, and transported to the laboratory within 24 h. The samples were then stored at -80 °C in the Key Laboratory of Yak Breeding Engineering of Gansu Province. The blood storage numbers were R-5-1-001, R-5-1-002, and R-5-1-003, and the DNA was extracted from the samples within one month. No yaks were sacrificed in this experiment, and all blood samples were taken from live yaks. The jugular vein area of the yaks was locally disinfected with alcohol before blood sample collection. After collected blood samples, wiped with iodine to prevent any possible wounds being infected.

Extraction, amplification, and sequencing of the mitochondrial genome

The primers designed by Wang [10] and synthesized by Xi’an Qingke Biotechnology Co., Ltd. (Xi’an, China) were used to amplify the whole mitochondrial genome of the yaks. The reagents, methods, and reaction conditions used for the extraction and amplification of the mitochondrial genome were as previously described [27]. After amplification, the quality of the products was assessed using 1% agarose gel electrophoresis before sequencing at Xi’an Qingke Biotechnology Co., Ltd. (Xi’an, China).

Analysis of genetic structure and genetic diversity

The complete mitochondrial genomes (accession numbers OK375501–OK375872) of the 372 yak individuals were sequenced. After the inclusion of the mitochondrial sequences obtained in an earlier study (accession numbers MW414100–MW414210, MK124955.1), a total of 484 domestic yak individuals were used for experimental analysis. This final sample set covered all morphological groups and distribution ranges of the domestic yaks. In addition, the complete mitochondrial genome sequences (accession numbers GQ464246.1–GQ464266.1, MK033130.1, KR106993.1, KY829451.1, NC_025563.1) of 25 wild yaks were downloaded from National Center for Biotechnology Information (NCBI), increasing the total number of yak mitochondrial complete genome sequences to 509.

MAFFT 7.0 was used for sequence alignment [62]. AMOVA were performed using Arlequin v 3.5 software to determine the level of differentiation among the yak populations; the number of permutations was set to 10 000 [63]. DNAsp 6.0 [64] software was used to calculate the Fst among populations, detect the degree of genetic differentiation among populations, and perform genetic diversity analysis and neutrality tests. PermutCpSSR-2.0 software was used to calculate the Nst and Gst values at the species level (parameters set to 1 000 substitutions) to explore the relationship between genetic and geographical distances and examine whether the species distribution had a genetic geographical structure [65]. The yak network diagram was constructed with Popart [66] software. Default settings were used for all the software not described in detail.

Construction of phylogenetic tree

The yak phylogenetic tree was constructed with PhyloSuite [67]. The IQ-tree module was used to construct the ML phylogenetic tree. The optimal base replacement model was K3Pu + F + R5, “Bootstrap” was Standard, “Num of Bootstrap” was 1 000, and the SH-alRT test was enabled. The default value of repeat sampling times was 1 000. The Mbayes module was used for the construction of 4 MCMC chains of a Bayesian phylogenetic tree, and the optimal base replacement model was HKY + F + I + G4. The number of generations was 100 000 000, the sampling frequency was 1000, the number of runs was 2, and 25% of aging samples were discarded.

Estimation of differentiation time

Beast v1.10 [68] was used to estimate the differentiation time, with a Bayesian Information Criterion optimal model of TIM2 + F + I + G4. The corrected Akaike Information Criterion optimal model was GTR + F + R5. According to the fossil time query website (http://www.timetree.org/), the point of divergence between Bos mutus and Bison bison was 1.5 (0.2–3.9) Ma; the relaxed molecular clock model from Beast v 1.10 was used to calculate the differentiation time. A total of 5 000 000 generations, obtaining a tree every 1 000 generations, and ESS > 200 were used.

Research on population distribution dynamics

The climate data from all periods were downloaded from the WorldClimate database (http://www.worldclim.org/). Chinese soil data were obtained from the World Soil Database (HWSD) and downloaded from the National Ice Frozen Desert Science Data Center and the National Special Environment and Special Function Observation Research Station sharing service platform (http://www.crensed.ac.cn/portal/). The map data were downloaded from the National Basic Geographic Information Center (http://www.ngcc.cn/ngcc/). The mask method in ArcGIS 10.5 (https://developers.arcgis.com/) was used to extract the variable data from all environmental factors based on the map data and environmental data of China. MaxEnt v 3.4.1 (https://biodiversityinformatics.amnh.org/open_source/maxent/) was used to construct the distribution dynamics of species in different periods. AUC < 0.60 indicated the failure of the prediction model; AUC 0.60–0.70 indicated that the model prediction results are poor; AUC 0.70–0.80 indicated that the model prediction results are general; AUC 0.80–0.90 indicated that the model prediction results are accurate; AUC 0.90–1.00 indicated that the model prediction results are very accurate [69].

Data availability

The datasets generated and analyzed during the present study are available in the GenBank repository, under the accession number: OK375501–OK375872 (https://www.ncbi.nlm.nih.gov/genbank/).The data supporting the conclusions of this study are available in the supplementary table.

Abbreviations

D-loop:

Displacement loop

QTP:

Qinghai-Tibet Plateau

Ma:

Million years ago

mtDNA:

Mitochondrial DNA

Cytb :

Cytochrome b

Hd:

haplotype diversity

Pi:

Nucleotide diversity

HPD:

Highest posterior density

LIG:

Last interglacial

LGM:

last glacial maximum

NCBI:

National Center for Biotechnology Information

AUC:

Area under the curves

References

  1. Xia M, Tian Z, Zhang F, Khan G, Gao Q, Xing R, Zhang Y, Yu J, Chen S. Lancea tibeticaDeep Intraspecific divergence in the endemic Herb (Mazaceae) distributed over the Qinghai-Tibetan Plateau. Front Genet. 2018;9:492.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  2. Liu JQ, Sun YS, Xue-Jun GE, Gao LM, Qiu YX. Phylogeographic studies of plants in China: advances in the past and directions in the future. J Syst Evol. 2012;50(4):9.

    Article  CAS  Google Scholar 

  3. Li Y, Zong W, Zhao S, Qie M, Yang X, Zhao Y. Nutrition and edible characteristics, origin traceability and authenticity identification of yak meat and milk: a review. Trends Food Sci Technol. 2023;139:104133.

    Article  CAS  Google Scholar 

  4. Chen N, Zhang Z, Hou J, Chen J, Gao X, Tang L, Wangdue S, Zhang X, Sinding MS, Liu X, et al. Evidence for early domestic yak, taurine cattle, and their hybrids on the Tibetan Plateau. Sci Adv. 2023;9(50):eadi6857.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  5. Shah M, Xu C, Wu S, Zhao W, Luo H, Yi C, Liu W, Cai X. Isolation and characterization of spermatogenic cells from cattle, yak and cattleyak. Anim Reprod Sci. 2018;193:182–90.

    Article  PubMed  Google Scholar 

  6. Wiener G, Han J, Long R. The yak 2nd edn. (Regional Office for Asia and the Pacific Food and Agriculture Organization of the United Nations The yak 2nd edn. (Regional Office for Asia and the Pacific Food and Agriculture Organization of the United Nations, Bangkok.

  7. Guo S, Savolainen P, Su J, Zhang Q, Qi D, Zhou J, Zhong Y, Zhao X, Liu J. Origin of mitochondrial DNA diversity of domestic yaks. BMC Evol Biol. 2006;6(1):73.

    Article  PubMed  PubMed Central  Google Scholar 

  8. Li J. Genetic Diversity and phylogenetic analysis of Yak mtDNA in Karakoram Pamir Region. Master Kashgar Univ 2019, pp:12.

  9. Ma Z, Xia X, Chen S, Zhao X, Zeng L, Xie Y, Chao S, Xu J, Sun Y, Li R, et al. Identification and diversity of Y-chromosome haplotypes in Qinghai yak populations. Anim Genet. 2018;49(6):618–22.

    Article  CAS  PubMed  Google Scholar 

  10. Wang Z, Shen X, Liu B, Su J, Yonezawa T, Yu Y, Guo S, Ho SYW, Vila C, Hasegawa M, et al. Phylogeographical analyses of domestic and wild yaks based on mitochondrial DNA: new data and reappraisal. J Biogeogr. 2010;37(12):2332–44.

    Article  Google Scholar 

  11. Qiu Q, Wang L, Wang K, Yang Y, Ma T, Wang Z, Zhang X, Ni Z, Hou F, Long R, et al. Yak whole-genome resequencing reveals domestication signatures and prehistoric population expansions. Nat Commun. 2015;6:10283.

    Article  CAS  PubMed  Google Scholar 

  12. Zhang S, Liu W, Liu X, Du X, Zhang K, Zhang Y, Song Y, Zi Y, Qiu Q, Lenstra J, et al. Structural variants selected during yak domestication inferred from Long-Read whole-genome sequencing. Mol Biol Evol. 2021;38(9):3676–80.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  13. Brantingham PJ, Olsen JW, Schaller GB. Lithic assemblages from the Chang Tang region, Northern Tibet. Antiquity. 2001;75(288):319–27.

    Article  Google Scholar 

  14. Ma ZJ, Zhong JC, Han JL, Xu JT, Liu ZN, Bai WL. Research progress on molecular genetic diversity of yaks. Genetics. 2013;35(02):151–60.

    CAS  Google Scholar 

  15. Ma ZJ, Zhong JC, Han JL, Xu JT, Dou QL, Chang HP. Genetic diversity of mtDNA D-Loop region in wild yaks (Bos grunniens mutus). J Ecol. 2009;29(09):4798–803.

    CAS  Google Scholar 

  16. Song QQ, Zhong JC, Zhang CF, Xin JW, Ji QM, Cai ZX. Analysis on genetic diversity and phyletic evolution of mitochondrial DNA from tibetan yaks. Acta Theriol Sinica. 2014;34(04):356–65.

    Google Scholar 

  17. Qi XB, Han JL, Blench R, et al. Understanding the yak pastoralism in central Asian highlands: genetic evidence for origin, domestication and dispersion of domestic yak. In: Sanchez-Mazas A, BlenchR. Ross MD, Peiros I, Lin M, editors. Past Human migrations in East Asia: matching Archaeology, Linguistics and Genetics[M]. London and New York: Routledge, Taylor & Francis Group; 2008. pp. 427–42.

    Google Scholar 

  18. Zhao SJ, Chen ZH, Ji QM, Cai ZX, Zhang CF, Xin JW, Zhong JC. Sequence analysis of mtDNA COIII of tibetan yaks. Scientia Agricultura Sinica. 2011;44(23):4902–10.

    CAS  Google Scholar 

  19. Derenko M, Denisova G, Malyarchuk B, Dambueva I, Bazarov B. Mitogenomic diversity and differentiation of the buryats. J Hum Genet. 2018;63(1):71–81.

    Article  CAS  PubMed  Google Scholar 

  20. Chang GB, Chang H, Chen GH, Chen R, Zhao J, Zhuo Y, Guan YP. Analysis of genetic diversity and phylogenetic status of Bazhou yak based on partial sequence of Cytb gene. Chin J Anim Sci. 2010;46(17):19–21.

    CAS  Google Scholar 

  21. Zhang CF, Xu LJ, Ji QM, Xin JW, Zhong JC. Genetic diversity and evolution relationship on mtDNA D-loop in tibetan yaks. Acta Ecol Sin. 2012;32(05):1387–95.

    Article  CAS  Google Scholar 

  22. Excoffier L, Yang Z. Substitution rate variation among sites in mitochondrial hypervariable region I of humans and chimpanzees. Mol Biol Evol. 1999;16(10):1357–68.

    Article  CAS  PubMed  Google Scholar 

  23. Ingman M, Gyllensten U. Analysis of the complete human mtDNA genome: methodology and inferences for human evolution. J Heredity. 2001;92(6):454–61.

    Article  CAS  Google Scholar 

  24. Bao P, Pei J, Ding X, Wu X, Chu M, Xiong L, Liang C, Guo X, Yan P. Characterisation of the complete mitochondrial genome of the Jinchuan Yak (Bos grunniens). Mitochondrial DNA B Resour. 2019;4(2):3856–7.

    Article  PubMed  PubMed Central  Google Scholar 

  25. Liang CN, Wu X, Ding X, Wang H, Guo X, Chu M, Bao P, Yan P. Characterization of the complete mitochondrial genome sequence of wild yak (Bos mutus). Mitochondrial DNA Part DNA Mapp Sequencing Anal. 2016;27(6):4266–7.

    Google Scholar 

  26. Guo X, Wu X, Chu M, Bao P, Xiong L, Liang C, Ding X, Pei J, Yan P. Characterization of the complete mitochondrial genome of the Pamir yak (Bos grunniens). Mitochondrial DNA B Resour. 2019;4(2):3165–6.

    Article  PubMed  PubMed Central  Google Scholar 

  27. Wang X, Pei J, Bao P, Cao M, Guo S, Song R, Song W, Liang C, Yan P, Guo X. Mitogenomic diversity and phylogeny analysis of yak (Bos grunniens). BMC Genomics. 2021;22(1):325.

    Article  PubMed  PubMed Central  Google Scholar 

  28. Wu X, Zhou X, Ding X, Liang C, Guo X, Chu M, Wang H, Pei J, Bao P, Yan P. Characterization of the complete mitochondrial genome of the Huanhu Yak (Bos Grunniens). Mitochondrial DNA Part B. 2019;4(1):1235–6.

    Article  Google Scholar 

  29. Zhong D, Ding L. Rising process of the Qinghai-Xizang (Tibet) Plateau and its mechanism. Sci China (Series D) 1996(04):289–95.

  30. Chassot P, Nemomissa S, Yuan YM, Kupfer P. High paraphyly of Swertia L. (Gentianaceae) in the Gentianella-lineage as revealed by nuclear and chloroplast DNA sequence variation. Plant Syst Evol. 2001;229(1–2):1–21.

    Article  CAS  Google Scholar 

  31. Luikart G, Gielly L, Excoffier L, Vigne J, Bouvet J, Taberlet P. Multiple maternal origins and weak phylogeographic structure in domestic goats. Proc Natl Acad Sci USA. 2001;98(10):5927–32.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  32. Guo SC, Qi DL, Chen GH, Xu SX, Zhao XQ. Genetic diversity and classification of mitochondrial DNA (mtDNA) in yaks. J Ecol 2008(09):4286–94.

  33. Hua Y. Population genetics and phylogeny of scorpionflies based on mitochondrial DNA. Doctor Northwest A&F Univ 2020, pp: 20–1.

  34. Jump A, Marchant R, Peñuelas J. Environmental change and the option value of genetic diversity. Trends Plant Sci. 2009;14(1):51–8.

    Article  CAS  PubMed  Google Scholar 

  35. Zhang MQ, Xu X, Luo SJ. The genetics of brown coat color and white spotting in domestic yaks (Bos grunniens). Anim Genet. 2014;45(5):652–9.

    Article  CAS  PubMed  Google Scholar 

  36. Zheng BX, Xu QQ, Shen YP. The relationship between climate change and quaternary glacial cycles on the Qinghai-Tibetan Plateau: review and speculation. Quatern Int 2002, 97 – 8:93–101.

  37. Bandelt H, Forster P, Röhl A. Median-joining networks for inferring intraspecific phylogenies. Mol Biol Evol. 1999;16(1):37–48.

    Article  CAS  PubMed  Google Scholar 

  38. Tang CQ, Matsui T, Ohashi H, Dong YF, Momohara A, Herrando-Moraira S, Qian S, Yang Y, Ohsawa M, Luu HT, et al. Identifying long-term stable refugia for relict plant species in East Asia. Nat Commun. 2018;9(1):4488.

    Article  PubMed  PubMed Central  Google Scholar 

  39. Beck RA, Burbank DW, Sercombe WJ, Riley GW, Barndt JK, Berry JR, Afzal J, Khan AM, Jurgen H, Metje J, et al. Stratigraphic evidence for an early collision between northwest India and asia. Nature. 1995;373(6509):55–8.

    Article  CAS  Google Scholar 

  40. Beheregaray LB. Twenty years of phylogeography: the state of the field and the challenges for the Southern Hemisphere. Mol Ecol. 2008;17(17):3754–74.

    Article  PubMed  Google Scholar 

  41. Molnar P, Boos WR, Battisti DS. Orographic Controls on Climate and Paleoclimate of Asia: Thermal and Mechanical Roles for the Tibetan Plateau. In: Annual Review of Earth and Planetary Sciences, Vol 38 Edited by Jeanloz R, Freeman KH, vol. 38; 2010: 77–102.

  42. Sun BN, Wu JY, Liu YS, Ding ST, Li XC, Xie SP, Yan DF, Lin ZC. Reconstructing Neogene vegetation and climates to infer tectonic uplift in western Yunnan, China. Palaeogeography Palaeoclimatology Palaeoecology. 2011;304(3–4):328–36.

    Article  Google Scholar 

  43. Shi Y. Characteristics of late quaternary monsoonal glaciation on the Tibetan Plateau and in East Asia. Quatern Int 2002, 97 – 8:79–91.

  44. Myers N, Mittermeier RA, Mittermeier CG, da Fonseca GA, Kent J. Biodiversity hotspots for conservation priorities. Nature. 2000;403(6772):853–8.

    Article  CAS  PubMed  Google Scholar 

  45. He K, Jiang X. Sky islands of southwest China. I: an overview of phylogeographic patterns. Chin Sci Bull. 2014;59(7):585–97.

    Article  Google Scholar 

  46. Marchese C. Biodiversity hotspots: a shortcut for a more complicated concept. Global Ecol Conserv. 2015;3:297–309.

    Article  Google Scholar 

  47. Zhang Q, Chiang TY, George M, Liu JQ, Abbott RJ. Phylogeography of the Qinghai-Tibetan Plateau endemic Juniperus przewalskii (Cupressaceae) inferred from chloroplast DNA sequence variation. Mol Ecol. 2005;14(11):3513–24.

    Article  CAS  PubMed  Google Scholar 

  48. Chen S, Wu G, Zhang D, Gao Q, Duan Y, Zhang F, Chen S. Potential refugium on the Qinghai-Tibet Plateau revealed by the chloroplast DNA phylogeography of the alpine species Metagentiana striata (Gentianaceae). Bot J Linn Soc. 2008;157(1):125–40.

    Article  Google Scholar 

  49. Yang F-S, Li Y-F, Ding X, Wang X-Q. Extensive population expansion of Pedicularis longiflora (Orobanchaceae) on the Qinghai-Tibetan Plateau and its correlation with the quaternary climate change. Mol Ecol. 2008;17(23):5135–45.

    Article  PubMed  Google Scholar 

  50. Jia DR, Abbott RJ, Liu TL, Mao KS, Bartish IV, Liu JQ. Out of the Qinghai-Tibet Plateau: evidence for the origin and dispersal of eurasian temperate plants from a phylogeographic study of Hippophaë rhamnoides (Elaeagnaceae). New Phytol. 2012;194(4):1123–33.

    Article  PubMed  Google Scholar 

  51. Khan G, Zhang F, Gao Q, Fu P, Zhang Y, Chen S. Spiroides shrubs on Qinghai-Tibetan Plateau: Multilocus phylogeography and palaeodistributional reconstruction of Spiraea alpina and S. Mongolica (Rosaceae). Mol Phylogenetics Evol. 2018;123:137–48.

    Article  CAS  Google Scholar 

  52. Wang B, Xie F, Li J, Wang G, Li C, Jiang J. Phylogeographic investigation and ecological niche modelling of the endemic frog species Nanorana pleskei revealed multiple refugia in the eastern tibetan Plateau. PeerJ. 2017;5:e3770.

    Article  PubMed  PubMed Central  Google Scholar 

  53. Wang L, Abbott RJ, Zheng W, Chen P, Wang Y, Liu J. History and evolution of alpine plants endemic to the Qinghai-Tibetan Plateau: Aconitum Gymnandrum (Ranunculaceae). Mol Ecol. 2009;18(4):709–21.

    Article  PubMed  Google Scholar 

  54. Wang H, Qiong L, Sun K, Lu F, Wang Y, Song Z, Wu Q, Chen J, Zhang W. Phylogeographic structure of Hippophae Tibetana (Elaeagnaceae) highlights the highest microrefugia and the rapid uplift of the Qinghai-Tibetan Plateau. Mol Ecol. 2010;19(14):2964–79.

    Article  CAS  PubMed  Google Scholar 

  55. Gao Q, Zhang D, Duan Y, Zhang F, Li Y, Fu P, Chen S. Intraspecific divergences of Rhodiola alsia (Crassulaceae) based on plastid DNA and internal transcribed spacer fragments. Bot J Linn Soc. 2012;168(2):204–15.

    Article  Google Scholar 

  56. Gao QB, Zhang FQ, Xing R, Gornall RJ, Fu PC, Li Y, Gengji ZM, Chen SL. Phylogeographic study revealed microrefugia for an endemic species on the Qinghai-Tibetan Plateau: Rhodiola Chrysanthemifolia (Crassulaceae). Plant Syst Evol. 2016;302(9):1179–93.

    Article  Google Scholar 

  57. Owen LA, Finkel RC, Barnard PL, Ma HZ, Asahi K, Caffee MW, Derbyshire E. Climatic and topographic controls on the style and timing of late quaternary glaciation throughout Tibet and the Himalaya defined by Be-10 cosmogenic radionuclide surface exposure dating. Q Sci Rev. 2005;24(12–13):1391–411.

    Article  Google Scholar 

  58. Provan J, Bennett KD. Phylogeographic insights into cryptic glacial refugia. Trends Ecol Evol. 2008;23(10):564–71.

    Article  PubMed  Google Scholar 

  59. Dalen L, Nystrom V, Valdiosera C, Germonpre M, Sablin M, Turner E, Angerbjorn A, Arsuaga JL, Gotherstrom A. Ancient DNA reveals lack of postglacial habitat tracking in the arctic fox. Proc Natl Acad Sci USA. 2007;104(16):6726–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  60. Li QF, Li YF, Zhao XB, Liu ZS, Zhang QB, Song DW, Qu XG, Li N, Xie Z. Sequencing of mitochondrial DNA cytochrome b gene of yak and its origin, classification and position. J Anim Husb Veterinary Med 2006(11):1118–23.

  61. Zhong JC, Chai ZX, Ma ZJ, Wang Y, Yang WY, La H. Sequencing and phylogenetic analysis of the mitochondrial genome of wild yaks. J Ecol. 2015;35(05):1564–72.

    Google Scholar 

  62. Katoh K, Rozewicki J, Yamada K. MAFFT online service: multiple sequence alignment, interactive sequence choice and visualization. Brief Bioinform. 2019;20(4):1160–6.

    Article  CAS  PubMed  Google Scholar 

  63. Excoffier L, Smouse P, Quattro J. Analysis of molecular variance inferred from metric distances among DNA haplotypes: application to human mitochondrial DNA restriction data. Genetics. 1992;131(2):479–91.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  64. Rozas J, Ferrer-Mata A, Sánchez-DelBarrio J, Guirao-Rico S, Librado P, Ramos-Onsins S, Sánchez-Gracia A. DnaSP 6: DNA sequence polymorphism analysis of large data sets. Mol Biol Evol. 2017;34(12):3299–302.

    Article  CAS  PubMed  Google Scholar 

  65. Hartnup K, Huynen L, Te Kanawa R, Shepherd L, Millar C, Lambert D. Ancient DNA recovers the origins of Māori feather cloaks. Mol Biol Evol. 2011;28(10):2741–50.

    Article  CAS  PubMed  Google Scholar 

  66. Leigh JW, Bryant D. POPART: full-feature software for haplotype network construction. Methods Ecol Evol. 2013;6:1110–6.

    Article  Google Scholar 

  67. Zhang D, Gao F, Jakovlić I, Zou H, Zhang J, Li W, Wang G. PhyloSuite: an integrated and scalable desktop platform for streamlined molecular sequence data management and evolutionary phylogenetics studies. Mol Ecol Resour. 2020;20(1):348–55.

    Article  PubMed  Google Scholar 

  68. Suchard M, Lemey P, Baele G, Ayres D, Drummond A, Rambaut A. Bayesian phylogenetic and phylodynamic data integration using BEAST 1.10. Virus Evol. 2018;4(1):vey016.

    Article  PubMed  PubMed Central  Google Scholar 

  69. Bellard C, Bertelsmeier C, Leadley P, Thuiller W, Courchamp F. Impacts of climate change on the future of biodiversity. Ecol Lett. 2012;15(4):365–77.

    Article  PubMed  PubMed Central  Google Scholar 

Download references

Acknowledgements

The authors would like to thank all the reviewers who participated in the review and the brothers and sisters for their help in the experiment.

Funding

This work was supported by the China Agriculture Research System of MOF and MARA (CARS-37) and the Innovation Project of Chinese Academy of Agricultural Sciences (25-LZIHPS-01).

Author information

Authors and Affiliations

  1. Key Laboratory of Yak Breeding in Gansu Province, Lanzhou Institute of Husbandry and Pharmaceutical Sciences, Chinese Academy of Agricultural Sciences, Lanzhou, 730050, P.R. China

    Xingdong Wang, Jie Pei, Lin Xiong, Pengjia Bao, Min Chu, Xiaoming Ma, Yongfu La, Chunnian Liang, Ping Yan & Xian Guo

  2. Key Laboratory of Animal Genetics and Breeding on Tibetan Plateau, Ministry of Agriculture and Rural Affairs, Lanzhou, 730050, P.R. China

    Xingdong Wang, Jie Pei, Lin Xiong, Pengjia Bao, Min Chu, Xiaoming Ma, Yongfu La, Chunnian Liang, Ping Yan & Xian Guo

Authors
  1. Xingdong Wang
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Jie Pei
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Lin Xiong
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Pengjia Bao
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Min Chu
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Xiaoming Ma
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. Yongfu La
    View author publications

    You can also search for this author in PubMed Google Scholar

  8. Chunnian Liang
    View author publications

    You can also search for this author in PubMed Google Scholar

  9. Ping Yan
    View author publications

    You can also search for this author in PubMed Google Scholar

  10. Xian Guo
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

X.G. and X.W. conceptualized this study. P.B. and P.Y. helped in the investigation. J.P., L.X., and X.W. helped in methodology and software. X.W., Y.L., and X.M. performed data curation. X.W. and M.C. helped in writing the original draft. X.G., C.L., and X.W. helped in writing, reviewing, and editing the manuscript. X.G. helped in funding acquisition. All authors contributed to the interpretation of the results and writing of the manuscript.

Corresponding author

Correspondence to Xian Guo.

Ethics declarations

Ethics approval and consent to participate

All procedures involving animals were performed according to the guidelines of the China Council on Animal Care and the Ministry of Agriculture of the People’s Republic of China and approved by the Animal Ethics Committee of Lanzhou Institute of Husbandry and Pharmaceutical Sciences, Chinese Academy of Agricultural Sciences.

Consent for publication

Not applicable.

Competing interests

The authors declare no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Electronic supplementary material

Below is the link to the electronic supplementary material.

Supplementary Material 1

Supplementary Material 2

Rights and permissions

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 licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence 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 licence, 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 data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Wang, X., Pei, J., Xiong, L. et al. Genetic diversity, phylogeography, and maternal origin of yak (Bos grunniens). BMC Genomics 25, 481 (2024). https://doi.org/10.1186/s12864-024-10378-z

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12864-024-10378-z

Keywords

BMC Genomics

ISSN: 1471-2164

Contact us