Skip to main content

Forest fragmentation causes an isolated population of the golden takin (Budorcas taxicolor bedfordi Thomas, 1911) (Artiodactyla: Bovidae) in the Qinling Mountains (China)


Budorcas taxicolor bedfordi is a rare animal uniquely distributed in the Qinling Mountains (China). Human disturbance and habitat fragmentation have directly affected the survival of B. t. bedfordi. It is urgent to clarify the genetic diversity and genetic structure of the B. t. bedfordi population and implement effective conservation measures. In this study, 20 new polymorphic microsatellite loci were isolated by Illumina sequencing. The genetic diversity and population structure of 124 B. t. bedfordi individuals from three populations (Niubeliang population, Zhouzhi population, and Foping population) were analysed according to these 20 microsatellite loci. Our results indicated that B. t. bedfordi had a low level of genetic variability and that there was inbreeding in the three populations. The population genetic structure analyses showed that the Niubeliang population had a trend of differentiation from other populations. National roads can affect population dispersal, while ecological corridors can promote population gene exchange. None of the three B. t. bedfordi populations experienced bottleneck effects. For conservation management plans, the Zhouzhi population and Foping population should be considered one management unit, and the Niubeliang population should be considered another management unit. We suggest building an ecological corridor to keep the habitat connected and formulating tourism management measures to reduce the influence of human disturbance on B. t. bedfordi.

Peer Review reports


Biodiversity is the summation of all variations in life, and it is the material basis for human survival and development. With the growth of the population and the development of technology after the Industrial Revolution, human beings have become increasingly capable of utilizing and transforming nature. Human activities have seriously damaged the ecosystem by changing the climate or the environment or by encroaching on habitats. Forest habitat degradation and fragmentation have made a large number of species extinct or endangered [1,2,3]. The population size of many species is decreasing rapidly, and the genetic diversity within species is also declining [4]. Due to human activities, habitat fragmentation and other factors, many endangered species have low genetic diversity [5,6,7,8,9]. Consequently, there is an urgent need to study the genetic polymorphisms of each endangered species, especially those with very small distribution areas, to prevent the further decline of their populations. Determining how to repair the forest habitat environment and protect wild animals reasonably and effectively has become an urgent task for the survival and development of human beings.

Golden takin, Budorcas taxicolor bedfordi, belongs to the order Artiodactyla and family Bovidae. The forest habitat of this species is uniquely located in the Qinling Mountains, which covers a 24,773 km2 area between 32° N to 34° N and 106° E to 110° E with elevations of 1,550 m to 3,600 m above sea level in China [10,11,12,13]. The population quantity is 5500 ~ 8000 [8], and this species is a national first-class protected wild animal in China and is listed in CITES Appendix II [14]. The International Union for the Conservation of Nature and Natural Resources (IUCN) classified B. t. bedfordi as vulnerable (VU) [15]. The nature reserves in the distribution area of B. t. bedfordi are small and disconnected. Moreover, in the distribution area of B. t. bedfordi, there are 108, 210, 312, and 316 National Roads; 102 and 212 provincial roads; Xi ‘an to Hanzhong, Xi ‘an to Ankang, and Xi ‘an to Wuhan expressways; Baoji to Chengdu and Xi ‘an to Ankang railway special lines; and the recently built Xi ‘an to Chengdu high-speed railway. Although these roads connect the northern and southern regions of the Qinling Mountains and promote economic development, they divide the B. t. bedfordi habitat into multiple patches. Due to large-scale forest fragmentation and habitat fragmentation, the B. t. bedfordi populations in the Qinling Mountains are isolated from each other, and their living conditions are not optimal [16].

To formulate scientific and effective protection measures, the scientific evaluation of the population genetic diversity and genetic structure of B. t. bedfordi is necessary. Since microsatellite loci can provide high-resolution genetic information, they are very suitable for population genetic studies of wild animals such as B. t. bedfordi with a small population distribution area [1718]. Compared with traditional microsatellite loci isolation methods, next-generation sequencing (NGS) techniques are more efficient and less costly [19, 20]. In this research, we isolated 20 microsatellite markers by NGS techniques and analysed the population genetic diversity and genetic structure of B. t. bedfordi in detail. The purposes of this research were as follows: (1) to examine the genetic variability of the species; (2) to appraise the population genetic structure; and (3) to analyse the influence of human disturbance on gene flow among populations to investigate the phylogeographical and landscape genetic patterns of B. t. bedfordi.


Tissue sample, DNA extraction and genome sequencing

Muscle samples (n = 124) were collected from dead B. t. bedfordi within three geographical regions: Zhouzhi National Nature Reserve (ZZ population, 108° 11′ E 33°45′ N; n = 40), Foping National Nature Reserve (FP population, 107° 53′ E 33°32′ N; n = 48) and Niubeiliang National Nature Reserve (NBL population, 109° 02′ E 33°91′ N; n = 36) (Fig. 1). The animals died of natural causes in the wild and were found by the reserve staff while they were patrolling the mountains. They had been dead for five to thirty days when they were found. The All muscle tissue samples collected from dead individuals were fixed in 95% ethanol. A QIAGEN DNeasy Blood & Tissue Kit (Qiagen, Valencia, CA) was used to extract whole genomic DNA [21].

Fig. 1
figure 1

The habitat of B. t. bedfordi and distribution of samples. The map of China was based on the standard map service system GS (2019)1673 of the Ministry of Natural Resources ( and the map of Shaanxi Province was based on Shaanxi Provincial Platform for Common GeoSpatial Information Services ( The base map has not been modified

Whole genomic DNA was enzyme digested into 400–500 bp fragments. The B. t. bedfordi DNA library was prepared by terminal repaired, adding poly A tail, ligating sequencing adapter, purification and PCR amplification. The purified DNA library was quantified using Agilent 2100 Bioanalyzer and ABI StepOnePlus Real-Time PCR System, and sequenced using Illumina NovaSeq 6000 platform.

Microsatellite characterization and development

The mean quality score, Q20 values, Q30 values and the GC content of the sequencing raw data were estimated by seqtk 1.3 [22]. Microsatellite repeat units (di-, tri-, tetra-, penta-, hexanucleotide repeats) were scanned by Galaxy server pipeline [23, 24]. Primers were designed by Primer Premier 5.0 [25], using design parameters as follows: GC content between 50% and 70%, primer length range between 18 and 26 bp, target amplicon length range between 80 and 400 bp, annealing temperature between 55 and 66 °C. We designed three pairs of primer for each microsatellite locus and chose the primer pair with the highest score. From the novel microsatellite markers, we randomly selected 100 pairs of primers for further studies. The primers were synthesized by Beijing Tsingke Biotech Co., Ltd. (Beijing, China).

Primer testing and selection

PCR amplification was carried out with 100 primer pairs using mixed DNA pool of golden takin. After agarose-gel electrophoresis, primer pairs that can amplify target bands were used as primary primers. The PCR amplification system was as follows: 1 μL B. t. bedfordi DNA, 12.5 μL PCR Master mix (TIANGEN, Beijing), 2 μL primer pairs, 9.5 μL ddH2O. The PCR reactions were as follows: initial denaturation at 94 °C for 8 min; 38 cycles of denaturation for 30 s at 94 °C, specified annealing temperature for 30s, and 72 °C for 45s; and 72 °C extension for 10 min. The primary primers were validated by direct sequencing and capillary electrophoresis (Beijing Tsingke Biotech Co., Ltd. Beijing, China).

Table 1 Information of twenty new microsatellite loci for B. t. bedfordi

Data analysis

CONVERT 1.3 [26] and GenALEx 6.5 software [27, 28] were used to convert the format of the genetic data. The null alleles were detected by Microchecker, and the null allele frequency (NO) was calculated by Cervus 3.0 software [29]. We also used Cervus 3.0 to estimate the polymorphic information content (PIC) of each microsatellite locus. Hardy–Weinberg equilibrium (PHWE) was calculated by Microchecker. GenALEx 6.5 [27]was used to analyse microsatellite loci polymorphisms and population genetic diversity, including the number of different alleles (NA), number of effective alleles (NE), Shannon’s information index (I), expected heterozygosity (HE), observed heterozygosity (HO), and Wright’s inbreeding coefficient (FIS). Nei’s genetic distance among individuals was calculated by PHYLIP 3.6 [30]. The principal coordinates of three B. t. bedfordi populations were evaluated by GenALEx 6.5. INEST 2.0 [31] was used to estimate FIS per population and the parameter setting was the default Bayesian approach. We used GENEPOP V4 [32] to estimate the pairwise population genetic differentiation index (FST). The individual and population UPGMA trees were constructed by MEGA 5.2 [33, 34].

In addition, the genetic structure and the relationship among the three populations of B. t. bedfordi were analysed by using Structure 2.3.1 software [35]. The parameters were set as follows: number of populations (K), 1 to 10; length of burn-in period, 20000; and number of MCMC repeats after burn-in, 100000. The most likely value of the genetic Cluster K value is selected according to the principle of maximum likelihood value or calculation of the inflection point of the ΔK peak value.

Furthermore, we used two statistical methods to examine potential bottlenecks, heterozygote excess of populations and the mode-shift indicator test. Heterozygote excess was detected by Wilcoxon tests (one tail) using BOTTLENECK 1.2.02 software [36]. Wilcoxon tests were calculated under infinite allele (IAM), stepwise mutation (SMM), and two-phase (TPM) models, and each model was performed with 10000 simulation repeats. In the mode-shift indicator test, rare allele frequency will decrease if the population experiences a bottleneck effect. Consequently, when the allele frequency distribution follows the normal L-shaped form, it indicates that the population has not experienced a bottleneck effect.

To determine whether the isolation between populations was caused by distance factors or human interference, we calculated the correlation between genetic distance (FST among populations) and geographic distance of samples with the Mantel test of GenALEx [27]. The geographical distance between each sampling point was obtained by ArcGIS 10.5.


Microsatellite detection and test of neutrality

A total of 374,645 microsatellite motifs were identified. Among these motifs, there were 259,891 (69.37%) dinucleotide repeats, 60,468 (16.14%) trinucleotide motifs, 46,681 (12.46%) tetranucleotide motifs, 4,983 (1.33%) pentanucleotide motifs and 2,622 (0.7%) hexanucleotide motifs. As a result, 20 microsatellite loci (GenBank accession numbers MN733993-MN734006 and OP970799-OP970804) with high polymorphism were isolated.

The value of NA ranged from 5 to 14, the value of PIC ranged from 0.421 to 0.851, the value of average HO was 0.3260, and the value of average HE was 0.7360 (Table 1). After Bonferroni correction, there were thirteen loci (65%) that significantly deviated from Hardy-Weinberg equilibrium and seven loci (35%) that conformed to Hardy-Weinberg equilibrium (Table 1).

Genetic diversity

In view of its small geographical distribution area and habitat fragmentation, we estimated the genetic diversity level of the B. t. bedfordi population to be low. The results of B. t. bedfordi genetic diversity from the three populations are shown in Table 2. The values of I, HE and UHE were lowest in the NBL population (I = 1.402 ± 0.083, HE = 0.704 ± 0.029, UHE = 0.714 ± 0.029). The values of the inbreeding coefficient over polymorphic loci among populations ranged from 0.462 ± 0.080 to 0.568 ± 0.091.

Table 2 Genetic diversity of three B. t. bedfordi populations (mean ± SE)

Genetic structure

The results of the FST values among the three populations of B. t. bedfordi show that the FST value was lowest between the ZZ and FP populations (FST =0.026) and highest between the FP and NBL populations (FST =0.051) (Table 3). The results demonstrated that the NBL population had higher genetic differentiation than the ZZ and FP populations.

Table 3 Values of FST among three populations of B. t. bedfordi

The UPGMA phylogenetic tree of 124 individuals is shown in Fig. 2. The 124 individuals formed two large clades, clade I and clade II, of which clade II can be further divided into two small clades, clade II-1 and clade II-2.

Using different colours to label 124 individuals based on the collected regions, it was found that individuals from the NBL population were all clustered in clade (I) The individuals from the ZZ population and FP population were all clustered in clade (II) The ZZ population mostly clustered in clade II-1, and 3 individuals were distributed in clade II-2. Most of the FP population was clustered in clade II-2, and 6 individuals were distributed in clade II-1.

Fig. 2
figure 2

UPGMA phylogenetic tree of 124 individuals of B. t. bedfordi

A UPGMA phylogenetic tree of the three populations of B. t. bedfordi is shown in Fig. 3 (a). The three populations first gathered into two large branches. The ZZ population and FP population had a closer genetic relationship. The NBL population had the farthest genetic relationship with other populations.

The principal coordinate analysis showed that 124 individuals could be divided into three groups (Fig. 3(b)), and the clustering results were basically consistent with the phylogenetic tree. The NBL population was gathered in Group A, the ZZ population was mostly gathered in Group B, and the FP population was mostly gathered in Group C. The ZZ population and FP population had sporadic individuals distributed in Group C and Group B.

Fig. 3
figure 3

(a) UPGMA phylogenetic tree of three B. t. bedfordi populations; (b) principal coordinates of three B. t. bedfordi populations

It can be seen from the population genetic structure analysis that when K = 2, ΔK had the maximum value (Supplementary Figure S1). All individuals from the Niubeiliang region were assigned to the red cluster, and all individuals from the Zhouzhi and Foping regions were assigned to another cluster (green) (Fig. 4 (a)). When K = 3, the green cluster was divided into two clusters (blue and green) (Fig. 4(b)). Except for the red cluster, the blue and green clusters were mixed with other colours, indicating that there was genetic exchange between the FP and ZZ populations.

Fig. 4
figure 4

STRUCTURE bar plots for three populations of B. t. bedfordi. (a) K = 2, (b) K = 3

Test for bottleneck effects

Under the three assumptions of the IAM, TPM and SMM models, no significant heterozygosity excess was detected in any of the three populations (p > 0.05) and displayed a lack of heterozygosity. The results of the mode-shift indicator test are shown in Supplementary Table S1 and Figure S2. The values of allele frequency ranged from 0.1 to 0.7 and were mainly distributed at 0.1. Each population was characterized by an L-shaped distribution of allelic frequencies (i.e., many alleles with low frequencies). The results indicated that the B. t. bedfordi populations did not experience bottleneck effects.

Table 4 Test of bottleneck effects for B. t. bedfordi


Genetic variation of B. t. bedfordi

PIC is a commonly used indicator in genetic diversity research and is used to evaluate the identification ability of microsatellite markers and reflect the reliability of information provided by microsatellite loci [37, 38]. When PIC > 0.50, the polymorphism of microsatellite markers was higher, and the information provided was rich, which could well reflect the genetic diversity. When 0.25 < PIC < 0.50, microsatellite markers had high polymorphism and information content, which can provide more reasonable information. When PIC < 0.25, microsatellite markers had low polymorphism and provided less information [39]. All twenty microsatellite loci in this study had high PIC values, greater than 0.25, and the highest value was 0.851. Therefore, these 20 microsatellite markers were highly informative and were sufficient for discriminating among B. t. bedfordi individuals and populations.

Genetic diversity is essential because it is the raw material for natural selection allowing species to adapt to new environmental conditions. The loss of genetic diversity decreases species adaptive potentials [39]. The allele number and heterozygosity value are important indicators of population genetic diversity. The higher the values are, the higher the genetic polymorphism of the populations. We observed a total of 174 alleles from 20 microsatellite loci in this study, with an average number of 8.7 alleles. The number of effective alleles in the three sampling populations was 4.395, slightly lower than that of other endangered wild animals in the same region, giant panda (Ailuropoda melanoleuca) 4.58 [40], golden monkey (Rhinopithecus roxellana) 5.05 [41], forest musk deer (Moschus berezovskii) 6.37 [42], goat (Capra hircus) 5.23 [43], and cattle (Bos taurus)6.21 [44]. The average observed heterozygosity (HO) of B. t. bedfordi was 0.326, which was lower than that of the giant panda (0.488) [45], golden monkey (0.620) [46], goat (0.66) [47] and cattle (0.74) [44]. Our results confirm that the genetic variety of B. t. bedfordi was quite low. In addition, the genetic diversity of the B. t. bedfordi population was also studied based on the mtDNA D-loop region and single nucleotide polymorphisms (SNPs). The nucleotide diversity of B. t. bedfordi analysed by the mtDNA D-loop region was 0.0021 [48], indicating low genetic diversity [49]. The single nucleotide variants of B. t. bedfordi analysed by SNP were lower than those of the giant panda, golden monkey, goat, and cattle [50,51,52]. These results indicated that the B. t. bedfordi population suffered from poor genetic diversity, which was in accordance with the results of this study. At present, the number of golden takins in the Qinling Mountains is only approximately 5000 [53]. In approximately 1960, the quantity of B. t. bedfordi was approximately 600. Later, the Chinese government implemented relevant protection policies and established protected areas, which gradually increased the quantity of B. t. bedfordi. However, the protected areas were not connected, and human disturbances were severe, which hindered communication between individuals from different regions. Therefore, although the population has increased in recent decades, the genetic diversity is low, and the inbreeding coefficient within the population is high. The low genetic diversity and small population size of B. t. bedfordi could have led to its near extinction. Consequently, there is an urgent need to protect B. t. bedfordi.

In the Hardy-Weinberg equilibrium test, 65% of the loci were out of equilibrium, and the three populations also deviated significantly from Hardy-Weinberg equilibrium. The three assumptions of IAM, TPM and SMM by Bottleneck showed that three populations of B. t. bedfordi presented a deficit of heterozygotes. This result indicated that the deviation of the B. t. bedfordi population from Hardy-Weinberg equilibrium was caused by heterozygote deficiency. Heterozygous deficiency is mainly caused by null alleles, inbreeding, and the Wahlund effect [54, 55]. In population genetics studies, null allele frequencies below 0.2 did not affect the result analysis [56, 57]. In this study, out of 20 loci, only 1 locus had a null allele frequency of 0.31, thus indicating that the null allele did not affect the accuracy of the analysis in this study [58]. In addition, the FIS of the three groups were all greater than 0, and the average FIS value reached 0.53, indicating that there was significant inbreeding among the B. t. bedfordi populations [59]. Furthermore, it can be seen from the results of the FST values, phylogenetic analyses, and population genetic structure analysis that B. t. bedfordi had significant population structure (Table 4; Figs. 3 and 4). The significant population structure can lead to a potential Wahlund effect [60, 61]. Consequently, we suspected that heterozygote deficiency was mainly due to inbreeding or population substructure. Although the B. t. bedfordi population did not experience a bottleneck effect and the population increased year by year, the high inbreeding coefficient and heterozygote deficiency will cause B. t. bedfordi to easily suffer from the bottleneck effect in the case of abrupt habitat change.

Effects of forest change on the population genetic structure of B. t. bedfordi

The genetic differentiation among the B. t. bedfordi populations was low. However, the FST value of the NBL population and the FP population reached 0.051, indicating that the NBL population had a trend of differentiation. At the same time, the Structure results also showed that all B. t. bedfordi individuals were divided into two groups, and the NBL population was a separate group. It can also be seen from the Structure results that when K = 3, there is a certain differentiation between the FP and ZZ populations. We examined the correlation between genetic distance and geographical distance and found that the correlation was not significant (R = 0.358, P = 0.480) (Figure S3), suggesting that the genetic differences between B. t. bedfordi populations were incompletely caused by distance. As seen from the satellite map, the three B. t. bedfordi populations were separated from each other by 108 National Roads and 210 National Roads (Fig. 1). On both sides of the road, the vegetation was seriously damaged, and the crown density decreased, which seriously affected the population dispersal of B. t. bedfordi. 108 National Road, built in 1969, obstructed individual communication between the ZZ population and FP population in the short term. In addition, in 2000, the Qinling Tunnel of 108 National Road was put into use, while the mountaintop section was abandoned. Vegetation restoration based on forest care and renewal and planting of native tree species was carried out at the same time. Therefore, the gene flow with the ZZ population and FP population was facilitated, and the genetic differentiation was small. 210 National Road, built in 1930, blocked individuals of the NBL population for a long time. Moreover, in 2002, the Xihan Expressway was built between the habitats of the NBL and ZZ populations, which increased the difficulty of individual communication between the NBL population and the other two populations. These two roads also produced certain barriers to the gene flow of golden monkey populations [46]. Thus, it can be inferred that the genetic differentiation between B. t. bedfordi populations was caused not only by distance but also by the influence of roads.

Conservation for B. t. bedfordi

The current results indicated that B. t. bedfordi had low genetic diversity and inbreeding within the population, and its distribution range was small. The species is vulnerable to abrupt habitat change or infectious disease and faces an uncertain future. Thus, we propose the following protection recommendations.

  1. (1)

    The gene flow between the ZZ population and FP population was facilitated, so B. t. bedfordi inhabiting these two reserves can be protected as a whole population. The NBL population can be protected as a separate unit. In addition, it was found that national roads and expressways may interfere with the migration of B. t. bedfordi, blocking its gene flow. We suggest building ecological corridors in the meadow area of the NBL reserve to reduce or remove the impact of habitat fragmentation on the migration of B. t. bedfordi and gradually improve the quality of the habitat of B. t. bedfordi. Therefore, the habitat available for B. t. bedfordi can be connected and intact, and genes can be exchanged freely [62].

  2. (2)

    Because there is a small number of B. t. bedfordi distributed in the surrounding areas of each reserve, we propose that the surrounding areas should be designated as the extension range of the reserve and the habitat range should be gradually expanded to satisfy the population increase.

  3. (3)

    We suggest that the tourist administration bureau and relevant departments scientifically plan the norms of forest parks, the length of time scenic spots are open and control the size of tourist flow [63, 64]. These measures can reduce the interference with the daily activities and migration of B. t. bedfordi caused by roads and tourists as much as possible.

  4. (4)

    It is necessary to carry out ex situ conservation. The rescued B. t. bedfordi can be released in different places to support the genetic exchange with different populations.


In conclusion, the results obtained in this study indicate that the genetic diversity of B. t. bedfordi, an endangered wild animal, is quite low. Due to forest fragmentation and road construction, genetic differentiation among populations has occurred. However, ecological corridor construction and forest restoration were beneficial to individual exchange. We have proposed some measures for the conservation of B. t. bedfordi population diversity in the Qinling Mountains. This study provides a solid foundation for further research on the endangerment mechanism of B. t. bedfordi.

Data availability

The datasets generated and analysed during the current study are available in the KNB repository,


  1. Gehrig-Fasel J, Guisan A, Zimmermann NE, Gehrig-Fasel J. Tree line shifts in the Swiss Alps: climate change or land abandonment? J Veg Sci. 2007;18(4):571–82.

    Article  Google Scholar 

  2. Coppes J, Ehrlacher J, Thiel D, Suchant R, Braunisch V. Outdoor recreation causes effective habitat reduction in capercaillie Tetrao urogallus: a major threat for geographically restricted populations. J Avian Biol. 2017;48(12):1583–94.

    Article  Google Scholar 

  3. Immitzer M, Nopp-Mayr U, Zohmann M. Effects of habitat quality and hiking trails on the occurrence of Black Grouse (Tetrao tetrix L) at the northern fringe of alpine distribution in Austria. J Ornithol. 2014;155(1):173–81.

    Article  Google Scholar 

  4. Frankham R, Ballou JD, Briscoe DA. Introduction to conservation genetics. 2nd ed. United Kingdom: Cambridge University Press; 2010.

    Book  Google Scholar 

  5. Kang M, Wang J, Huang H. Demographic bottlenecks and low gene flow in remnant populations of the critically endangered Berchemiella Wilsonii var. Pubipetiolata (Rhamnaceae) inferred from microsatellite markers. Conserv Genet. 2008;9(1):191–9.

    Article  Google Scholar 

  6. Ge S, Hong DY, Wang HQ, Liu ZY, Zhang CM. Population genetic structure and conservation of an endangered conifer, Cathaya argyrophylla (Pinaceae). Int J Plant Sci. 1998;159(2):351–7.

    Article  Google Scholar 

  7. Li YY, Chen XY, Zhang X, Wu TY, Lu HP, Cai YW. Genetic differences between wild and artificial populations of metasequoia glyptostroboides: implications for species recovery. Conserv Biol. 2005;19(1):224–31.

    Article  Google Scholar 

  8. Kaneko S, Isagi Y, Nobushima F. Genetic differentiation among populations of an oceanic island: the case of Metrosideros Boninensis, an endangered endemic tree species in the Bonin Islands. Plant Species Biol. 2008;23(2):119–28.

    Article  Google Scholar 

  9. Moreira RG, McCauley RA, Cortés-Palomec AC, Fernandes GW, Oyama K. Spatial genetic structure of Coccoloba cereifera (Polygonaceae), a critically endangered microendemic species of Brazilian rupestrian fields. Conserv Genet. 2009;11(4):1247–55.

    Article  CAS  Google Scholar 

  10. Yang L, Wei FW, Zhan XJ, Fan HZ, Zhao PP, Huang GP, et al. Evolutionary conservation genomics reveals recent speciation and local adaptation in threatened takins. Mol Biol Evol. 2022;39(6):msac111.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  11. Zeng ZG, Song YL. Ecology and protection of Qinling Takin. Bull Biology. 2008;43(8):1–4. 63. (In Chinese).

    Google Scholar 

  12. Ma YT, Wang XF. Current status and protection measures of golden takin (Budorcas taxicolor bedfordi) in Qinling mountain ranges. Shaanxi for Sci Technol. 2008;47(2):80–3. (In Chinese).

    Google Scholar 

  13. Yuan ZH, Sun JC. Monitoring report of population size of giant panda, golden monkey and takin in changqing nature reserve. J Shaanxi Normal Univ. 2007;35(1):100–3. (In Chinese with English abstract).

    Google Scholar 

  14. Zeng Z, Song Y. Golden Takin (Budorcas taxicolor bedfordi). Chin J Zool. 2002;37(1):5. (In Chinese).

    Google Scholar 

  15. Song YL, Smith AT, MacKinnon J. Budorcas taxicolor. The IUCN red list of threatened species 2008: e.T3160A9643719.

  16. Kunz F, Kohnen A, Nopp-Mayr U, Coppes J. Past, present, future: tracking and simulating genetic differentiation over time in a closed metapopulation system. Conserv Genet. 2021;22(3):355–68.

    Article  Google Scholar 

  17. Park CS, Lee SY, Cho GJ. Evaluation of recent changes in genetic variability in thoroughbred horses based on microsatellite markers parentage panel in Korea. Anim Bioscience. 2022;35(4):527–32.

    Article  CAS  Google Scholar 

  18. Kerry RG, Montalbo FJP, Das R, Patra S, Mahapatra GP, Maurya GK, et al. An overview of remote monitoring methods in biodiversity conservation. Environ Sci Pollut Res. 2022;29(53):8017–221.

    Article  Google Scholar 

  19. Fernandez-Silva I, Whitney J, Wainwright B, Andrews KR, Ylitalo-Ward H, Bowen BW, et al. Microsatellites for next-generation ecologists: a post-sequencing bioinformatics pipeline. PLoS ONE. 2013;8(2):e55990.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  20. Yang LL, Hafiz A. Development and characterization of novel polymorphic microsatellite markers for Tapinoma indicum (Hymenoptera: Formicidae). J Insect Sci. 2021;21(6):1–6.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  21. Iyavoo S, Hadi S, Goodwin W. Evaluation of five DNA extraction systems for recovery of DNA from bone. Forensic Sci Int Genet Supplement. 2013;4(1):e174–5.

    Article  Google Scholar 

  22. Li H. seqtk: Toolkit for processing sequences in FASTA/Q formats. 2018.

  23. Grifths SM, Fox G, Briggs PJ, Donaldson IJ, Hoodet S, Richardson P, et al. A galaxy-based bio-informatics pipeline for optimised streamlined microsatellite development from illumina next-generation sequencing data. Conserv Genet Resour. 2016;8(4):481–6.

    Article  Google Scholar 

  24. Kerkhove T, Hellemans B, Troch MD, Backer AD, Volckaert F. Isolation and characterisation of 14 novel microsatellite markers through next generation sequencing for the commercial atlantic seabob shrimp Xiphopenaeus kroyeri. Mol Biol Rep. 2019;46(6):6565–9.

    Article  PubMed  CAS  Google Scholar 

  25. Lalitha S. Primer Premier 5. Biotech Softw Internet Rep. 2000;1(6):270–2.

    Article  Google Scholar 

  26. Glaubitz JC. Mol Ecol Notes. 2004;4(2):309–10. CONVERT: A user-friendly program to reformat diploid genotypic data for commonly used population genetic software packages.

  27. Peakall R, Smouse PE. GenAlEx 6.5: genetic analysis in Excel. Population genetic software for teaching and research-an update. Bioinformatics. 2012;28(19):2537–9.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  28. Selionova MI, Aibazov MM, Mamontova TV, Petrov SN, Kharzinova VR, Arsen DV, et al. 42 genetic differentiation of Russian goats and wild relatives based on microsatellite loci. J Anim Sci. 2020;98(4):19–20.

    Article  PubMed Central  Google Scholar 

  29. Pagnotta MA. Comparison among methods and statistical software packages to analyze germplasm genetic diversity by means of codominant markers. J. 2018;1(1):197–215.

    Article  Google Scholar 

  30. Challa S, Neelapu NRR. Phylogenetic trees: applications, construction, and assessment. In: Essentials of Bioinformatics.; Hakeem KR, Shaik NA, Banaganapalli B, Elango R. Eds. Switzerland: Springer nature switzerland AG. 2019; Volume III: pp.167–192.

  31. Chybicki IJ, Burczyk J. Simultaneous estimation of null alleles and inbreeding coefficients. J Hered. 2009;100(1):106–13.

    Article  PubMed  CAS  Google Scholar 

  32. Raymond M, Rousset F. GENEPOP (version 1.2): population genetics software for exact tests and ecumenicism. J Hered. 1995;86(3):248–9.

    Article  Google Scholar 

  33. Tamura K, Dudley J, Nei M, Kumar S. Molecular evolutionary genetics analysis (MEGA) software version 4.0. Mol Biol Evol. 2007;24(8):1596–9.

    Article  PubMed  CAS  Google Scholar 

  34. Tamura K, Peterson D, Peterson N, Stecher G, Nei M, Kumar S. MEGA5: molecular evolutionary genetics analysis using maximum likelihood, evolutionary distance, and maximum parsimony methods. Mol Biol Evol. 2011;28(10):2731–9.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  35. Hubisz MJ, Falush D, Stephens M, Pritchard JK. Inferring weak population structure with the assistance of sample group information. Mol Ecol Resour. 2009;9(5):1322–32.

    Article  PubMed  PubMed Central  Google Scholar 

  36. Piry S, Luikart G, Cornuet JM. BOTTLENECK: a computer program for detecting recent reductions in the effective population size using allele frequency data. J Heredity. 1999;90(4):502–3.

    Article  Google Scholar 

  37. Barkley NA, Roose ML, Krueger RR, Federici CT. Assessing genetic diversity and population structure in a citrus germplasm collection utilizing simple sequence repeat markers (SSRs). Theor Appl Genet. 2006;112(8):1519–31.

    Article  PubMed  CAS  Google Scholar 

  38. Uzun A, Yesilogu T, Polat I, Aka-Kacar Y, Gulsen O, Yildirim B et al. Evaluation of genetic diversity in lemons and some of their relatives based on SRAP and SSR markers. Plant Mol Biology Report 201;29(3):693–701.

  39. Serrote CML, Reiniger LRS, Silva KB, Rabaiolli SMDS, Stefanel CM. Determining the polymorphism information content of a molecular marker. Gene. 2020;726:144175.

    Article  PubMed  CAS  Google Scholar 

  40. Wei FW, Hu YB, Yan L, Nie YG, Wu Q, Zhang ZJ. Giant pandas are not an evolutionary cul-de-sac: evidence from multidisciplinary research. Mol Biol Evol. 2015;32(1):4–12.

    Article  PubMed  CAS  Google Scholar 

  41. Li YL, Huang K, Tang SY, Feng L, Yang J, Li ZH, Li BG. Genetic structure and evolutionary history of Rhinopithecus roxellana in qinling mountains, central China. Front Genet. 2020;11:611914.

    Article  PubMed  Google Scholar 

  42. Wang D, Xu G, Wang YH, He S, Bu SH, Zheng XL. Study on polymorphisms of microsatellites DNA of Chinese captive forest musk deer (Moschus berezovskii). Acta Theriol Sinica. 2019;39(9):599–607.

    Article  Google Scholar 

  43. Asroush F, Mirhoseini SZ, Badbarin N, Seidavi A, Tufarelli V, Laudadio V, Dario C, Selvaggi M. Genetic characterization of Markhoz goat breed using microsatellite markers. Archives Anim Breed. 2018;61(4):469–73.

    Article  Google Scholar 

  44. Svishcheva G, Babayan O, Lkhasaranov B, Tsendsuren A, Abdurasulov A, Stolpovsky Y. Microsatellite diversity and phylogenetic relationships among East Eurasian Bos taurus breeds with an emphasis on rare and ancient local cattle. Anim (Basel). 2020;24(9):1493.

    Article  Google Scholar 

  45. National Forestry and Grassland Administration. The 4th national survey report on giant panda in China. China: Science Press; 2021.

    Google Scholar 

  46. Huang K. Genetic structure and phylogeography of Rhinopithecus roxellana qinlingensis. Ph.D. Dissertation, North West University, China; 2014.

  47. Whannou HRV, Spanoghe M, Dayo G-K, Demblon D, Lanterbecq D, Dossa LH. Genetic diversity assessment of the indigenous goat population of Benin using microsatellite markers. Front Genet. 2023;14:1079048.

    Article  PubMed  PubMed Central  Google Scholar 

  48. Feng H, Huang Y, Ren Y, Feng CL, Liu XN. Structure of mitochondrial DNA control region and genetic diversity of Budorcas taxicolor in Shaanxi province. Acta Agriculturae Boreali-occidentalis Sinica. 2014;23(1):1–6. (In Chinese with English abstract).

    Google Scholar 

  49. Lan H, Shi LM. The origin and genetic differentiation of native breeds of pigs in South West China: an approach from mitochondrial DNA polymorphism. Biochem Genet. 1993;31(2):51–60.

    Article  PubMed  CAS  Google Scholar 

  50. Zhao SC, Zheng PP, Dong SS, Zhan XJ, Wu Q, Guo XS, et al. Whole-genome sequencing of giant pandas provides insights into demographic history and local adaptation. Nat Genet. 2013;45(1):67–71.

    Article  PubMed  CAS  Google Scholar 

  51. Zhou XM, Meng XH, Liu ZJ, Chang J, Wang BS, Li MZ, et al. Population genomics reveals low genetic diversity and adaptation to hypoxia in snub-nosed monkeys. Mol Biol Evol. 2016;33(10):2670–81.

    Article  PubMed  CAS  Google Scholar 

  52. Li A, Yang Q, Li R, Dai X, Cai K, Lei Y, Jia K, Jiang Y, Zan L. Chromosome-level genome assembly for takin (Budorcas taxicolor) provides insights into its taxonomic status and genetic diversity. Mol Ecol. 2023;32(6):1323–34.

    Article  PubMed  CAS  Google Scholar 

  53. Zeng GZ, Zhong WQ, Song YL, Li JS, Zhao LG, Gong HS. Present status of studies on eco-biology of Takin. Acta Theriol Sinica. 2003;23(2):161–7. (In Chinese with English abstract).

    Google Scholar 

  54. Paxton RJ, Thorén PA, Gyllenstrand N, Tengö J. Microsatellite DNA analysis reveals low diploid male production in a communal bee with inbreeding. Biol J Linn Soc. 2000;69(4):483–502.

    Article  Google Scholar 

  55. Wang ZF, Hamrick JL, Godt MJW. High genetic diversity in Sarracenia leucophylla (Sarraceniaceae), a carnivorous wet-land herb. J Hered. 2004;95(3):234–43.

    Article  PubMed  CAS  Google Scholar 

  56. Abdelkader AA, Ata N, Benyoucef MT, Djaout A, Azzi N, Yilmaz O, et al. New genetic identification and characterisation of 12 Algerian sheep breeds by microsatellite markers. Italian J Anim Sci. 2018;17(1):38–48.

    Article  CAS  Google Scholar 

  57. Al-Atiyat RM, Aljumaah RS, Alshaikh MA, Abudabos AM. Microsatellite-based genetic structure and diversity of local Arabian Sheep breeds. Front Genet. 2018;9:408.

    Article  PubMed  PubMed Central  Google Scholar 

  58. Séré M, Thévenon S, Belem AMG, De Meeûs T. Comparison of different genetic distances to test isolation by distance between populations. Heredity. 2017;119:55–63.

    Article  PubMed  PubMed Central  Google Scholar 

  59. Weir BS, Cockerham CC. Estimation F-statistics for the analysis of pupulation structure. Evolution. 1984;38(6):1358–70.

    Article  PubMed  CAS  Google Scholar 

  60. Jump AS, Peñuelas J. Genetic effects of chronic habitat fragmentation in a wind-pollinated tree. Proc Natl Acad Sci USA. 2006;103(21):8096–100.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  61. Meeûs TD. J Hered. 2018;109(4):446–56. FST, WahlundEffects, and Null Alleles.

    Article  PubMed  Google Scholar 

  62. Kettunen M, Terry A, Tucker G, Jones A. Guidance on the maintenance of landscape connectivity features of major importance for wild flora and fauna - Guidance on the implementation of article 3 of the birds directive (79/409/EEC) and article 10 of the habitats directive (92/43/EEC). Belgium: Institute for European Environmental Policy (IEEP); 2007. pp. 1–97.

    Google Scholar 

  63. Tang FL, Yan Y, Liu WG. Construction progress of national park system in China. Biodivers Sci. 2019;27(2):123–7.

    Article  Google Scholar 

  64. Rule A, Dill SE, Sun G, Chen A, Khawaja S, Li I, Zhang V, Rozelle S. Challenges and opportunities in aligning conservation with development in China’s national parks: a narrative literature review. Int J Environ Res Public Health. 2022;19(19):12778.

    Article  PubMed  PubMed Central  Google Scholar 

Download references


This research was funded by the National Natural Science Foundation of China (No: 32000353), Natural Science Foundation of Shaanxi Province (2020NY-008, 2022SF-191), Science and Technology Program of Shaanxi Academy of Sciences (2019 K-11, 2021 K-06). The funding body played no role in the design of the study and collection, analysis, interpretation of data, and in writing the manuscript.

Author information

Authors and Affiliations



Conceptualization, F.H.; formal analyses, F.H., C.F.J. and W.L.; resources, W.L.; writing—original draft, F.H.; review and editing, F.H. and J.T.Z.; funding acquisition, F.H. All authors have read and agreed to the published version of the manuscript.

Corresponding author

Correspondence to Hui Feng.

Ethics declarations

Ethical approval and consent to participate

This study and included experimental procedures were approved by Institution of Animal Care and the Ethics Committee of Shaanxi Institute of Zoology. The muscle samples used in this study were collected from dead takin in the wild. The main causes of death were diseases, fighting injuries, predation by black bears and leopards, and other non-human factors. The sample collection was permitted by the Management Division of Zhouzhi National Nature Reserve, Foping National Nature Reserve and Niubeiliang National Nature Reserve and the sampling process was followed by the Regulations of the People’s Republic of China on the Protection of Wild Animals.

Ethical guidelines

Field studies and other non-experimental research on animals in this study comply with the IUCN Policy Statement on Research Involving Species at Risk of Extinction and the Convention on the Trade in Endangered Species of Wild Fauna and Flora.

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: Figure S1

Criteria for selecting the optimal K for assigning individuals during the Bayesian cluster analysis. Table S1. The proportion of alleles in different allelic frequency classes for different populations of B. t. bedfordi. Figure S2. Graphic representation of proportion of alleles in different allelic frequency classes for different populations of B. t. bedfordi. Figure S3. Relationship between pairwise FST/(1-FST) and geographic distance among populations of B. t. bedfordi

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 The Creative Commons Public Domain Dedication waiver ( 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

Feng, H., Cao, F., Jin, T. et al. Forest fragmentation causes an isolated population of the golden takin (Budorcas taxicolor bedfordi Thomas, 1911) (Artiodactyla: Bovidae) in the Qinling Mountains (China). BMC Zool 9, 2 (2024).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: