Yuli Li1†
Jia Yang
Zhonghu Li
Baoguo Li- 1Shaanxi Key Laboratory for Animal Conservation, College of Life Sciences, Northwest University, Xi'an, China
- 2School of Pharmacy, Xi'an Jiaotong University, Xi'an, China
- 3Key Laboratory of Resource Biology and Biotechnology in Western China, Ministry of Education, College of Life Sciences, Northwest University, Xi'an, China
- 4Center for Excellence in Animal Evolution and Genetics, Chinese Academy of Sciences, Kunming, China
The Qinling mountainous region is one of the world's biodiversity hotspots and provides refuges for many endangered endemic animals. The golden snub-nosed monkeys (Rhinopithecus roxellana) are considered as a flagship species in this area. Here, we depicted the genetic structure and evolutionary history via microsatellite markers and combination with the ecological niche models (ENMs) to elucidate the intraspecific divergent and the impacts of the population demography on our focal species. Our results revealed three distinct subpopulations of R. roxellana and also uncovered asymmetric historical and symmetric contemporary gene flow that existed. Our evolutionary dynamics analyses based on diyabc suggested that the intraspecific divergence accompanied with effective population sizes changes. The ENM result implied that the distribution range of this species experienced expansion during the last glacial maximum (LGM). Our results highlighted that geological factors could contribute to the high genetic differentiation within the R. roxellana in the Qinling Mountains. We also provided a new insight into conservation management plans with endangered species in this region.
Background
Historical processes such as geographic changes and human activities have greatly affected patterns of dispersal and genetic migration among populations; these processes could restrict gene flow and accelerate genetic differentiation through inbreeding and genetic drift (Waage and Greathead, 1988; Frankham, 2005) and consequently shape population structure and evolutionary history (Balkenhol et al., 2009). In many species, empirical studies have indicated that mountains (Lait and Burg, 2013), rivers (Chambers and Garant, 2010), and deserts (Astrid et al., 2012) could act as physical or ecological obstacles to prevent animal dispersal to mitigate gene flow. In addition, as humans continue to expand across the world, many once-continuous landscapes become divided into separate fragmented patches (Hanski, 1994; Newbold et al., 2015). Low genetic variations will decrease the adaptive abilities of species to circumstances, and then this species become more vulnerable when environmental conditions changes, which will finally lead to their extinction (Frankham, 2005; Willi et al., 2006; Burkey, 2010).
The golden snub-nosed monkey, Rhinopithecus roxellana, is an Asian colobine endemic to the temperate forests of the mountainous regions in central China (Li et al., 2002; Kirkpatrick and Grueter, 2010), which has a typical multi-level social structure discovered in primate (Kirkpatrick and Grueter, 2010). The society of this monkey consists of four levels: unit, band, herd, and troop (Grueter et al., 2012). Furthermore, the bands can be classified into the breeding band (BB) and the all-male band (AMB) (Kirkpatrick and Grueter, 2010; Grueter et al., 2012). This primate once ranged widely across southern, southwestern, and central China (Li et al., 2003). However, due to the historical changes and continuous expansion of human populations and their corresponding activities, numbers and ranges of R. roxellana have been significantly reduced. In consequence, this primate today only distribute within the three isolated areas among four provinces (Shaanxi, Sichuan, Gansu, and Hubei), and its current estimated population size is no more than 20,000 individuals (Ren et al., 1998; Li et al., 2002, 2007).
In Shaanxi, R. roxellana occurs only in the Qinling Mountains (Li et al., 2000; Wang et al., 2014). This area is a major biodiversity hotspot of China with many unique and rare plant and/or animal species (Norton et al., 2011; Zhang et al., 2016). Within the mountains, 39 R. roxellana troops comprising a total of 4,000 individuals are distributed across the counties of Zhouzhi, Ningshan, Yangxian, Taibai, and Foping (Figure 1; Li et al., 2000). However, during the past few years, land reclamation, tree-felling, illegal hunting, and infrastructural development have all posed critical threats to the habitats of golden snub-nosed monkeys in these regions (Feng et al., 2009; Zhang et al., 2016). Previous research has shown that expanding human impact in nearby the Qinling Mountains has resulted in local extinction and population decline of R. roxellana (Li et al., 2002; Wang et al., 2014).
Figure 1. Sampling sites: the range of coordinates were 107.52-108.38°E, 33.64-33.82°N.
Based on the microsatellite data, Pan et al. (2005) showed the R. roxellana had relatively high intraspecific genetic diversity and formed different genetic structures at fine scale in China, e.g., the Minshan Mountains, the Shennongjia Mountains, and the Qinling Mountains. Also, Huang et al. (2016a) had detected that the genetic diversity and effective population size of R. roxellana in Qinling Mountains were not significantly decreased, which is in conflict with their expectations. In addition, a recent study has that, even under threat of habitat fragmentation, internal adjustment mechanisms in multilevel social systems can effectively avoid inbreeding due to the bridging role of certain social units (non-breeding groups: AMB) at a small scale within the Qinling Mountains (Li et al., 2020). However, the evolutionary history and demography of the R. roxellana inhabiting the Qinling Mountains remain poorly understood (Pan et al., 2005).
In this study, using microsatellite data, we evaluate the genetic diversity and gene flow within the Qinling populations and also investigated its complicated evolutionary history. We attempt to figure out the following: (i) whether the golden snub-nosed monkey was threatened by genetic homogeneity in fragmented habitat and, if so, (ii) whether the historical shifts facilitated the current genetic differentiation within R. roxellana, and (iii) whether and how the gene flow impacted its genetic structure.
Methods
Genetic Sampling
We collected a total of 344 fecal and hair samples of five representative populations scattered among five counties (Zhouzhi, Foping, Taibai, Yangxian, and Ningshan), which spread over five National Nature Reserve [the Zhouzhi National Nature Reserve (ZNNR), the Foping National Nature Reserve (FNNR), the Taibaishan National Nature Reserve (TNNR), the Changqing Natioanl Nature Reserve (CNNR), and Ningshan Nature Reserve (NNR); Figure 1: 107.52-108.38°E, 33.64-33.82°N]. These regions have semi-humid montane climate, and their elevations range from 1,100 to 2,930 m above sea level (Li et al., 2000, 2003). The average annual temperature of these areas is 8.59°C, with a minimum temperature in January (−17.7°C) and a maximum temperature in August (34.6°C). The average rainfall per annum is 577.64 mm. Details of the sample information are presented in Table 1; sampling locations are shown in Figure 1.
Table 1. Genetic diversities of five populations.
Because our focal species inhabit mountainous highlands and shy away from humans, it is quite hard to follow them across steep cliffs and precipitous valleys; it took us more than 1 year to collect the samples from March 2016 to April 2017. The percentage of fecal samples was 90 in our all samples, and they were stored in DETs [20% dimethyl sulfoxide (DMSO), 0.25 M of sodium-EDTA, 100 mM of Tris·HCl, pH 7.5, and NaCl to saturation] solution at −20°C (Allen et al., 1998). The hair samples, collected with a stick with adhesive tape by glue and fruit bait on an 80 × 6 cm wooden board, were mainly obtained from the Zhouzhi BB (because this band was half-habituated for 20 years on field observation) and then stored in silica gel for drying at room temperature.
Molecular Methods
Follicle DNA was extracted with proteinase K digestion in a PCR-compatible buffer, while fecal DNA was extracted using QIAamp DNA Stool Mini Kits (Qiagen, German). All samples were amplified based on the primers of 19 tetra-nucleotide microsatellite loci (see Supplementary Table 1 for locus profile) in an ABI Veriti Thermal Cycler with the following protocol: 95°C for 5 min, followed by 30 cycles (94°C for 30 s, 55–60°C for 45 s, 72°C for 45 s), and 72°C for 10 min. PCR products were segregated with an ABI PRISM 3100 Genetic Analyser, and their sizes relative to the internal size standard (ROX-labeled HD400) were determined using GENEMAPPER V3.7 (Applied Biosystems). To prevent genotyping errors such as false allele and allelic dropout, homozygote genotypes were confirmed by five independent replicates, with all heterozygotes observed and confirmed by at least three separate reactions (Taberlet et al., 1996). Replicates were detected by POLYRELATEDNESS V1.6 (Huang et al., 2016b) and excluded from before subsequent analyses.
Genetic Diversity
We used MICROCHECKER v2.2.3 to test the presence of null alleles for all loci (van Oosterhout et al., 2010). For each population, we calculated the genetic diversity indices that included number of alleles (AO), observed heterozygosity (HO), expected heterozygosity (HE), polymorphic information content (PIC), allelic richness (AR), and Wright's inbreeding coefficient at each locus utilizing GENAIEX V6.5 (Peakall and Smouse, 2012). We performed a Hardy–Weinberg equilibrium (HWE) test for each band at each locus using Fisher's exact tests in GENEPOP V4.3 (Rousset, 2008). Significance thresholds were adjusted for multiple tests by the sequential Bonferroni procedure (Rice, 1989).
We employed two methods to test the presence of bottleneck effects within different populations. The first method was based on deviations of allele frequencies in calculations of heterozygosity, where we used the signed test and two-tailed Wilcoxon test in BOTTLENECK V1.2.02 (Piry et al., 1999). We considered two types of mutation models: (i) a two-phase model (TPM) with 95% stepwise mutations and a variance of 12 and (ii) a stepwise mutation model (SMM) with iteration number set to 1,000. The second method involved calculation of the GW coefficient in ARLEQUIN V3.6 (Excoffier and Lischer, 2010).
Population Differentiation and Structure
Genetic differentiation among populations was evaluated using θ (FST) across loci with the 100,000 permutations with ARLEQUIN V3.6 (Excoffier and Lischer, 2010). The number of steps in the Markov chain was set 100,000, and that of burn-in steps was 10,000. Moreover, we performed AMOVA analysis in ARLEQUIN V3.6 (Excoffier and Lischer, 2010) to estimate genetic variations among the populations. The significance of fixation indices was tested using 10,000 permutations.
Bayesian clustering was performed using STRUCTURE V2.3.4 (Pritchard et al., 2000) to examine the population genetic structure. We first assessed the occurrence of population subdivision under an admixture model and with allele frequencies correlated. The program was run for K from 1 to 10. For each run, we used 1,000,000 MCMC cycles following 500,000 burn-in cycles. Two alternative methods were utilized to estimate the most likely number (K) of genetic clusters with the program STRUCTURE HARVSTER (Earl and Vonholdt, 2012, i.e., by tracing the change in the average of log-likelihood L(K) (Pritchard et al., 2000) and by calculating ΔK (Evanno et al., 2005).
Genetic Migration Analyses
To explore the historical genetic gene flow of Rhinopithecus roxellana in the Qinling Mountains, we used the program MIGRATE-N v3.6 (Beerli, 2006) to assess pairwise estimates of migration rates (Nm) among the five populations based on the microsatellite datasets. We performed maximum-likelihood analyses in MIGRATE-N v 3.6 using three long chains (100,000 runs) and 10 short chains (10,000 runs), and the burn-in was set as 10,000. We repeated this procedure five times to obtain the average maximum-likelihood estimates with 95% confidence intervals (CIs).
Moreover, in order to obtain the contemporary genetic migration among intraspecific lineages, we employed BAYESASS v3.0 (Wilson and Rannala, 2002) to calculate the intra-population migration rates. We conducted the analyses with burn-in of 1,000,000 iterations, setting 1,000 as the sampling frequency. Ten independent runs were executed to seek the model convergence.
Demographic History
In order to explore the plausible scenarios of divergence and population dynamics within R. roxellana, approximate Bayesian computation (ABC) was used in DIYABC v 2.04 (Cornuet et al., 2014). Based on the STRUCTURE results (Figure 2), 10 historical population divergence scenarios for these lineages were simulated by DIYABC analysis (Supplementary Table 2). We assumed uniform priors on all parameters, and then we used a goodness-of-fit test to check the priors of all parameters before implementing the simulations (Figure 3A). To estimate the divergence times among the lineages, the average generation time of R. roxellana was assumed to be 5 years, following Oleksyk et al. (2010).
Figure 2. Structure analyses using the LOCPRIOR model. (A,B) The mean estimated L(K) and ΔK as a function of number of clusters (K). (C) The bar plots for K = 2 and 3.
Figure 3. (A) The 10 scenarios of population history of five populations in Rhinopithecus roxellana in Qinling Mountains with DIYABC V 2.04. N1, N2, and N3 denote the current population lineages, while NA represents the ancestral population lineage. t1 and t2 are divergence times for the depicted event. (B) Schematic representation of changes in population size tested within the three population lineages.
Ecological Niche Modeling
Ecological niche modeling (ENM) analyses were performed with MAXENT v3.3.3k (Phillips et al., 2006) to assess the ecological niche of each lineage and to predict their potential range based on their georeferenced localities and environmental variables thereof. Only wild occurrences have been taken into account for the ENM analyses. The occurrence data of R. roxellana were obtained from our field observations, literature (Wen and Wen, 2006; Wen, 2009), and the records from two sources: the Global Biodiversity Information Facility (GBIF; https://www.gbif.org/) and the National Specimen Information Infrastructure (NSII; www.nsii.org.cn). In total, 131 georeferenced points were obtained.
We obtained 19 bioclimatic variables at 2.5 arc-minute resolutions from WorldClim (www.worldclim.org) (Hijmans et al., 2005) for four periods: the current, the Holocene (HOL; 12 ka–current), the last glacial maximum (LGM; 18–21 ka), and the last interglacial period (LIG; 120–140 ka). To avoid model over-fitting linked to correlated climatic parameters, only those seven variables that had low correlation coefficients with one another (r < 0.8) were retained for subsequent analyses (Supplementary Table 3). All ecological distribution models were visualized in ARCGIS 10.5.
Results
Genetic Diversity
A total of 361 samples were collected. After exclusion of 44 repeated samples as determined by microsatellite profiles, we eventually established a dataset consisting 317 individuals. With the use of MICRO-CHECKER, the frequencies of null alleles at each of the 19 loci were found to be lower than the threshold frequency (v = 0.15) across the five populations. The location, population size, and sample size of each band are presented in Table 1.
The genetic diversity indices of each band are also presented in Table 1. On the whole, genetic variability was moderate at the population level. The number of alleles per locus ranged from 4 to 6, averaging 5.047; observed heterozygosity (HO) ranged from 0.563 to 0.624, averaging 0.595; expected heterozygosity (HE) was from 0.558 to 0.672, averaging 0.623; the polymorphism information content (PIC) was from 0.473 to 0.538, averaging 0.511; and allelic richness ranged from 2.342 to 2.642, averaging 2.511 for the populations. The values of Wright's inbreeding coefficient (FIS) showed that the effect of inbreeding to five populations was weak.
No evidence of bottleneck effects has been found from the microsatellite data for each population. The results of bottleneck effects tests are shown in Table 2; none of the sign or Wilcoxon tests suggest heterozygosity excess or deficiency in either the SMM or TPM mutation models. The lowest P-value was 0.060 (Wilcoxon test under TPM model in DPY). The GW coefficients of all bands were above the empirical value 0.132 (Garza and Williamson, 2001), and the lowest GW coefficient was 0.846 ± 0.177 (CYG).
Table 2. Bottleneck effect tests of five populations.
Population Differentiation and Structure
The results of FST and permutation test revealed a relatively low but significant genetic divergence among those five bands (P < 0.05, Table 3). The AMOVA analyses showed that genetic variation among and within five populations was 15.74 and 84.26%, respectively (Table 4). According to the STRUCTURE analysis of microsatellite data, the most likely number of genetic clusters was at 3, and the scenario of K = 2 was also given for comparison. The L(K) (the probability of the data given K and the model) and ΔK (using the method of Evanno et al., 2005) as a function of selecting most likely K-value, and the bar plot of assignment matrix of LOCPRIOR analysis are shown in Figure 2. Monkeys from three populations (e.g., GNG, DPY, and SZZ) were collectively classified into a single cluster, while the remaining two monkey populations, HTP and CYG, were each assigned to a separate cluster.
Table 3. Pairwise FST (lower triangular) and permutation test.
Table 4. Results of AMOVA.
Gene Flow
Our study revealed asymmetrical historical gene flow among the five populations by MIGRATE, with the greatest gene migration between CYG and DPY (54.1218; Table 5) and the lowest between GNG and CYG (6.8760; Table 5). Moreover, BAYESASS analysis revealed contemporary genetic migration was symmetrical among the 10 related pairs. Estimated population sizes according to Bayesian modes, with 95% CI, were 0.3074 (95% CI: 0.2936–0.3661) for GNG, 0.4262 (95% CI: 0.4038–0.4502) for DPY, 0.2567 (95% CI: 0.2427–0.2719) for SZZ, 0.6696 (95% CI: 0.6341–0.7078) for HTP, and 1.4484 (95% CI: 1.3051–1.5423) for CYG.
Table 5. Rates of contemporary and historical gene flows among five populations using microsatellite data with the programs BAYESASS, MIGRATE, and bidirectional gene flow (Nm = immigrants per generation).
Evolutionary Dynamics History
In the DIYABC analysis, the posterior probability for scenario 8 (with 95% CI) was 0.914 (95% CI: 0.906–0.942), much higher than that of the other nine scenarios (Supplementary Figure 1). Our ABC demographic analysis detected that the first divergence in our focal species formed N1 and N3 lineages and the second divergence formed N2, which originated from the second contact of the N1 and N3 lineages (Figure 3B). The median values of the effective population sizes of N1, N2, and N3 were 2.18 × 105, 1.69 × 105, and 1.44 × 105, respectively, whereas NA was 1.85 × 106 (Table 6). In addition, the mean values of effective populations of N1, N2, and N3 were 2.61 × 105, 2.07 × 105, and 1.75 × 105, respectively, whereas NA was 1.96 × 106 (Table 6). The estimated median time of divergence between N1 (GNG, DPY, and SZZ populations) and N3 (CYG populations) (t1) from the common ancestor was 1.95 × 105 generations ago, N2 (HTP populations) diverged from the N1, and N3 (t2) was 4.51 × 104 generations ago (Table 6).
Table 6. Demographic approximate Bayesian computation (ABC) models for Rhinopithecus roxellana in Qinling Mountains.
Ecological Niche Modeling
We obtained four distributions tendencies of the Rhinopithecus roxellana, including the LIG, LGM, HOL, and current models by MAXENT (Figure 4). All models had high predictive ability [area under the receiver operating characteristic curve (AUC) > 0.9]. The predicted model for the present potential range of R. roxellana was fairly congruent with the current distribution of the species. The ENMs suggested that the R. roxellana had population expansion during the LIG to LGM. Compared with the HOL and the current distributions, the current pattern had little diminished.
Figure 4. Predicted past and current distributions of Rhinopithecus roxellana based on ecological niche modeling using MAXENT: the current, the Holocene (HOL; 12 ka–current), the last glacial maximum (LGM; 18–21 ka), and the last interglacial period (LIG; 120–140 ka).
Discussion
The present distributions of animal populations depend upon the historical processes and human activities. The habitat range of Rhinopithecus roxellana is a consequence of variation in ecological factors such as climatic change and intensity of human disturbance. In this study, we firstly used the interdisciplinary approaches to address the evolutionary dynamics of R. roxellana in Qinling Mountains. We examined genetic diversity, gene flow, genetic structures, and evolutionary history combination with the ENMs to elucidate ecological factors on the population demography.
Genetic Diversity and Genetic Differentiation
We used microsatellite data to analyze the genetic diversity and population structure of R. roxellana. A relatively high level of genetic diversity had been found. Compared with previous studies based on the microsatellite profiles, the expected heterozygosity in all populations from our investigation (0.628) was close to that of the other studies of this species (0.559 in Li et al., 2020; 0.625 in Huang et al., 2016a; 0.591 in Chang et al., 2013; and 0.631 in Pan et al., 2005). In addition, there was no evidence of any past genetic bottlenecks in any of the sampled monkey bands (Table 2), revealing that these populations have not suffered significant population size reduction. This result was consistent with that of Li et al. (2020) and Huang et al. (2016a), which also focused on the monkeys in the Qinling Mountains. The FST and permutation test results revealed a genetic differentiation among different populations (P < 0.05) associated with topographical barriers such as mountain ridges, which might lead to the reduced gene flow between populations (Chang et al., 2013).
Genetic Admixture and Gene Flow
The Bayesian clustering revealed three major subpopulations that strongly coincide with the major topographical ridge features in the study area. Our results indicated that the contemporary gene flow among the five populations was symmetric and low, whereas the historical gene flow for all pairs seemed asymmetric (Table 5). The possible explanations for the low contemporary gene flow were the more and more frequent human activities in this area. Human-constructed barriers such as rivers, villages, logging roads, and farmlands have limited the ability of BBs to freely move across these fragment landscapes for decades (Li et al., 2002; Wang et al., 2014). Moreover, the topography of the Qinling Mountains, which is dominated by high altitude temperate forest habitats, and open areas of fragmented and cleared forests containing villages and agricultural fields, may greatly increase the genetic differentiation of R. roxellana and restrict gene flow by decreasing the opportunities for successful dispersal (Wang et al., 2014). Nevertheless, satellite telemetry data indicated that some neighbor populations had small overlaps of their home ranges with occasional group fission–fusions (Qi et al., 2017). Meanwhile, the AMBs of R. roxellana played critical roles, which connected, gathered, and allocated gene flow among nearby populations (Li et al., 2020). But these only happened at a fine scale, which may explain the genetic migration of the bands between GNG, DPY, and SZZ (Figure 2C).
In general, gene flow is critical for mitigating the impacts of local adaptation by homogenizing populations from differing conditions or by spreading deleterious alleles across populations, and it also could spread favorable alleles to populations and increase genetic diversity (Welt et al., 2015; Epps and Keyghobadi, 2016). Previous research had reported that historical and contemporary gene migration was low among populations that often existed in highly fragmented habitats (Chiucchi and Gibbs, 2010), but if genetic structure was admixture, it might cause intense historical gene flow (Epps et al., 2013). However, whether the genetic structure of the Qinling Mountains was due to historical gene flow rather than contemporary ones needs further investigation.
Demographic History
Our ABC demographic analysis detected three lineages (N1, N2, and N3) existing within the monkeys of the Qinling Mountains (N1: GNG, DPY, and SZZ; N2: CYG; and N3: HTP). The first divergence in our focal species that formed N1 and N3 lineages from ancestral linage was at 0.975 Ma [95% highest posterior density (HPD): 0.33–2.86 Ma]. Based on the genetic analysis and fossil data from existing populations, snub-nosed monkeys are distributed across eastern, central, southern, and southwestern China in the past 2 million years (Liedigk et al., 2003). In addition, our initial divergence time in the Qinling Mountains might occur at the Late Pleistocene. Due to climate change and orogenic movement, lots of numbers of animals have declined population sizes and even became extinct during this period (Buuveibaatar et al., 2016; Olsoy et al., 2016). Some species have been forced to shift their range to resist environmental changes and/or habit conversion (Ceballos, 2002; Cleland et al., 2012). The first divergence time estimation coincided with this historical period.
Moreover, the second divergence was at 0.23 Ma (95% HPD: 0.14–0.25 Ma), and formed N2, which originated from the second contact of the N1 and N3 lineages. Since the habitat fragmentation and ecological barriers could lead to restrict gene migration, some species would resist the threat of genetic homogeneity (Lenormand, 2002). It had been reported the monkeys would disperse long distance and build their families in another population (Huang et al., 2017). But so far, since our understanding of the mechanism against genetic homogeneity is still limited, further research would be needed. In addition, according to previous studies (Yu et al., 2016; Zhou et al., 2016; Kuang et al., 2019), one lineage of the Sichuan (SG) population was closely clustered with the Qinling (QL) population. This result indicated that the SG population may have contributed to the genetic structure of the QL population; for our subsequent research, this part would also be taken into account.
The ENM analyses also suggested that the golden snub-nosed monkeys had expanded their ranges during the LIG (0.14–0.12 Ma) to LGM. At the end of the largest glaciation (ca. 1.2–0.6 Ma), the temperature increased, and also, the relatively cold climate may have continued until the late Ionian stage (0.3–0.126 Ma) (Goldewijk et al., 2001). Thus, because the golden snub-nosed monkey could adapt to cold habitats (Yu et al., 2016), it is feasible that they expanded their ranges and increased their population sizes during this period.
It is necessary for us to raise the issues on a conservation strategy for the R. roxellana, in addition to carrying out further survey and studies on its ecology, behavior, and dietary selection. At the moment, what we know about the most possible threats could include road accidents (traffic is heavy in the areas), habitat loss, human's intentional killing due to economic reasons from the local culture, and natural disasters (Li et al., 2002; Wang et al., 2014). Although a series of natural reserve have been established in the Mt. Qinling, rapid tourism development would inevitably lead to the expanding of roads and other artificial constructions, while increasing human disturbance in this area. A monitoring system on their dynamic population profiles should be established, so that emergency strategies for their conservation could be applied if a significant population decline was detected. On the other hand, more technical programs for captive breeding are also required.
The Qinling Mountains are famous for their global biodiversity hotspot and refugia for many endangered mammals and birds, particularly during the last glaciation (Wang et al., 2014). The Chinese Government has invested a great deal of manpower to the conservation of “four precious species in the Mts. Qinling” (the giant panda, Ailuropoda melanoleuca; the golden snub-nosed monkey, Rhinopithecus roxellana; the crested ibis, Nipponia nippon; and the takin, Budorcas taxicolor). Our new results of the evolutionary history with the R. roxellana could also shed light on species with the same geographical distribution, as well as other threatened or endangered amphibians and reptiles in this area.
Data Availability Statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.
Ethics Statement
The animal study was reviewed and approved by Wildlife Protection Society of China.
Author Contributions
BL and YL conceived the study. YL and KH performed the experiments. YL, LF, JY, and ZL contributed materials and analysis tools. BL, YL, and ST wrote the manuscript. LF and JY revised the manuscript. All authors approved the final version of the manuscript.
Funding
This work was funded by the Strategic Priority Research Program of the Chinese Academy of Sciences (XDB31020302), the National Natural Science Foundation of China (31730104, 31572278, and 31770411), the National Key Program of Research and Development, the Ministry of Science and Technology of China (2016YFC0503200), the Young Elite Scientists Sponsorship Program by CAST (2017QNRC001), and the Natural Science Foundation of Shaanxi Province (2018JM3024).
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Acknowledgments
We thanked four National Nature Reserves of the Qinling Mountains for the permission of this study.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2020.611914/full#supplementary-material
Supplementary File 1. Genepop genotype file.
Supplementary Figure 1. The posterior probabilities of the 10 scenarios.
Supplementary Table 1. Microsatellite marker information.
Supplementary Table 2. Prior distributions for model parameters used in divergence model comparisons (Figure 3A).
Supplementary Table 3. The seven bioclimatic variables that were used in ecological niche modeling.
References
Allen, M., Engström, A. S., Meyers, S., Handt, O., and Gyllensten, U. (1998). Mitochondrial DNA sequencing of shed hairs and saliva on robbery caps: sensitivity and matching probabilities. J. Forensic Sci. 43, 453–464. doi: 10.1520/JFS16169J
Astrid, V. S., Nathan, H. S., Graham, J. F., Paul, C. P., and Ryan, K. B. (2012). Landscape resistance to dispersal: simulating long-term effects of human disturbance on a small and isolated wolf population in southwestern Manitoba, Canada. Environ. Monit. Assess. 184, 6923–6934. doi: 10.1007/s10661-011-2469-9
Balkenhol, N., Gugerli, F., Cushman, S. A., Waits, L. P., Coulon, A., Arntzen, J. W., and Participants of the landscape genetics research agenda workshop 2007. (2009). Identifying future research needs in landscape genetics: where to from here? Landsc. Ecol. 24, 455–463. doi: 10.1007/s10980-009-9334-z
Beerli, P. (2006). Comparison of Bayesian and maximum-likelihood inference of population genetic parameters. Bioinformatics 22, 341–345. doi: 10.1093/bioinformatics/bti803
Burkey, T. V. (2010). Extinction rates in archipelagoes: implications for populations in fragmented habitats. Conserv. Biol. 9, 527–541. doi: 10.1046/j.1523-1739.1995.09030527.x
Buuveibaatar, B., Mueller, T., Strindberg, S., Leimgruber, P., Kaczensky, P., and Fuller, T. K. (2016). Human activities negatively impact distribution of ungulates in the Mongolian Gobi. Biol. Conserv. 203, 168–175. doi: 10.1016/j.biocon.2016.09.013
Ceballos, G. (2002). Mammal population losses and the extinction crisis. Science 296, 904–907. doi: 10.1126/science.1069349
Chambers, J. L., and Garant, D. (2010). Determinants of population genetic structure in eastern Chipmunks (Tamias striatus): the role of landscape barriers and sex-biased dispersal. J. Hered. 101, 413–422. doi: 10.1093/jhered/esq029
Chang, Z. F., Yang, B. H., Vigilant, L., Liu, Z. J., and Li, M. (2013). Evidence of male-biased dispersal in the endangered Sichuan snub-nosed monkey (Rhinopithexus roxellana). Am. J. Primatol. 76, 72–83. doi: 10.1002/ajp.22198
Chiucchi, J. E., and Gibbs, H. L. (2010). Similarity of contemporary and historical gene flow among highly fragmented populations of an endangered rattlesnake. Mol. Ecol. 19, 5345–5358. doi: 10.1111/j.1365-294X.2010.04860.x
Cleland, E. E., Allen, J. M., Crimmins, T. M., Dunne, J. A., Pau, S., Travers, S. E., et al. (2012). Phenological tracking enables positive species responses to climate change. Ecology 93, 1765–1771. doi: 10.1890/11-1912.1
Cornuet, J. M., Pudlo, P., Veyssier, J., Dehne-Garcia, A., Gautier, M., Leblois, R., et al. (2014). DIYABC v2.0: a software to make approximate Bayesian computation inferences about population history using single nucleotide polymorphism, DNA sequence and microsatellite data. Bioinformatics 30, 1187–1189. doi: 10.1093/bioinformatics/btt763
Earl, D. A., and Vonholdt, B. M. (2012). Structure Harvest: a website and program for visualizing structure output and implementing the Evanno method. Conserv. Genet. Resour. 4, 359–361. doi: 10.1007/s12686-011-9548-7
Epps, C. W., Castillo, J. A., Anne, S. K., Pierre, D. P., Greg, S. H., Mark, J., et al. (2013). Contrasting historical and recent gene flow among African buffalo herds in the Caprivi Strip of Namibia. J. Hered. 104, 172–181. doi: 10.1093/jhered/ess142
Epps, C. W., and Keyghobadi, N. (2016). Landscape genetics in a changing world: disentangling historical and contemporary influences and inferring change. Mol. Ecol. 24, 6021–6040. doi: 10.1111/mec.13454
Evanno, G. S., Regnaut, S. J., and Goudet, J. (2005). Detecting the number of clusters of individuals using the software STRUCTURE: a simulation study. Mol. Ecol. 14, 2611–2620. doi: 10.1111/j.1365-294X.2005.02553.x
Excoffier, L., and Lischer, H. E. (2010). Arlequin suite ver 3.5: a new series of programs to perform population genetics analyses under Linux and Windows. Mol. Ecol. Resour. 10, 564–567. doi: 10.1111/j.1755-0998.2010.02847.x
Feng, T. T., Manen, F. T., Zhao, N. X., and Wei, L. F. (2009). Habitat assessment for giant pandas in the qinling mountain region of China. J. Wildl. Manage. 73, 852–858. doi: 10.2193/2008-186
Frankham, R. (2005). Genetics and extinction. Biol. Conserv. 126, 131–140. doi: 10.1016/j.biocon.2005.05.002
Garza, J. C., and Williamson, E. G. (2001). Detection of reduction in population size using data from microsatellite loci. Mol. Ecol. 10, 305–318. doi: 10.1046/j.1365-294x.2001.01190.x
Goldewijk, K. K., Beusen, A., van Drecht, G., and de Vos, M. (2001). The HYDE 3.1 spatially explicit database of human-induced global land-use change over the past 12,000 years. Glob. Ecol. Biogeogr. 20, 73–86. doi: 10.1111/j.1466-8238.2010.00587.x
Grueter, C. C., Chapais, B., and Zinner, D. (2012). Evolution of multilevel social systems in nonhuman primates and humans. Int. J. Primatol. 33, 1002–1037. doi: 10.1007/s10764-012-9618-z
Hanski, I. (1994). Patch-occupancy dynamics in fragmented landscapes. Trends Ecol. Evol. 9, 131–135. doi: 10.1016/0169-5347(94)90177-5
Hijmans, R. J., Cameron, S. E., Parra, J. L., Jones, P. G., and Andy, J. (2005). Very high resolution interpolated climate surfaces for global land areas. Int. J. Biometeorol. 25, 1965–1978. doi: 10.1002/joc.1276
Huang, K., Guo, S. T., Cushman, S. A., Dunn, D. W., Qi, X. G., Hou, R., et al. (2016a). Population structure of the golden snub-nosed monkey Rhinopithecus roxellana in the Qinling Mountains, central China. Integr. Zool. 11, 350–360. doi: 10.1111/1749-4877.12202
Huang, K., Ritland, K., Dunn, D. W., Qi, X. G., Guo, S. T., and Li, B. G. (2016b). Estimating relatedness in the presence of null alleles. Genetics 201, 247–260. doi: 10.1534/genetics.114.163956
Huang, Z. P., Bian, K., Liu, Y., Pan, R. L., Qi, X. G., and Li, B. G. (2017). Male dispersal pattern in golden snub-nosed monkey (Rhinopithecus roxellana) in qinling mountains and its conservation implication. Sci. Rep. 7:46217. doi: 10.1038/srep46217
Kirkpatrick, R. C., and Grueter, C. C. (2010). Snub-nosed monkeys: multilevel societies across varied environments. Evol. Anthropol. Issues News Rev. 19, 98–113. doi: 10.1002/evan.20259
Kuang, W. M., Ming, C., Li, H. P., Wu, H., Frantz, L., Roos, C., et al. (2019). The origin and population history of the endangered golden snub-nosed monkey (Rhinopithecus roxellana). Mol. Biol. Evol. 36, 487–499. doi: 10.1093/molbev/msy220
Lait, L. A., and Burg, T. M. (2013). When east meets west: population structure of a high-latitude resident species, the boreal chickadee (Poecile hudsonicus). Heredity 1, 321–329. doi: 10.1038/hdy.2013.54
Lenormand, T. (2002). Gene flow and the limits to natural selection. Trends Ecol. Evol. 17, 183–189. doi: 10.1016/S0169-5347(02)02497-7
Li, B. G., Chen, C., Ji, W. H., and Ren, B. P. (2000). Seasonal home range changes of the Sichuan snub-nosed monkey (Rhinopithecus roxellana) in the Qinling mountains of China. Folia Primatol. 71, 375–386. doi: 10.1159/000052734
Li, B. G., Jia, Z. Y., Pan, R. L., and Ren, B. P. (2003). “Changes in distribution of the snub-nosed monkey in China,” in Primates in Fragments: Ecology and Conservation. ed L. K. Marsh (New York, NY: Kluwer Academic/Plenum Press), 29–51. doi: 10.1007/978-1-4757-3770-7_4
Li, B. G., Pan, R. L., and Oxnard, C. E. (2002). Extinction of snub-nosed monkeys in China during the Past 400 Years. Int. J. Primatol. 23, 1227–1244. doi: 10.1023/A:1021122819845
Li, M., Liu, Z. J., Gou, J. X., Ren, B. P., Pan, R. L., Su, Y. J., et al. (2007). Phylogeography and population structure of the golden monkeys (Rhinopithecus roxellana): inferred from mitochondrial DNA sequences. Am. J. Primatol. 69, 1195–1209. doi: 10.1002/ajp.20425
Li, Y. L., Wang, L., Wu, J. W., Ye, X. P., Garber, G. A., Yan, Y., et al. (2020). Bachelor groups in primate multilevel society facilitate gene flow across fragmented habitats. Curr. Zool. 66, 113–122. doi: 10.1093/cz/zoaa006
Liedigk, R., Yang, M. Y., Jablonski, N. G., Momberg, F., Geissmann, T., Lwin, N., et al. (2003). Evolutionary history of the odd-nosed monkeys and the phylogenetic position of the newly described Myanmar snub-nosed monkey Rhinopithecus strykeri. PLoS ONE 7:e37418. doi: 10.1371/journal.pone.0037418
Newbold, T., Hudson, L. N., Hill, S. L., Contu, S., Lysenko, I., Senior, R. A., et al. (2015). Global effects of land use on local terrestrial biodiversity. Nature 520, 45–50. doi: 10.1038/nature14324
Norton, C. J., Jin, C. Z., Wang, Y., and Zhang, Y. Q. (2011). “Rethinking the palearctic-oriental biogeographic boundary in Quaternary China,” in Asian Paleoanthropology: From Africa to China and Beyond, eds C. J. Norton, and D. R. Braun (Dordrecht: Springer Press), 81–100. doi: 10.1007/978-90-481-9094-2_7
Oleksyk, T. K., Smith, M. W., and O'Brien, S. J. (2010). Genome-wide scans for footprints of natural selection. Philos. Trans. R. Soc. Lond. Ser. B. 365, 185–205. doi: 10.1098/rstb.2009.0219
Olsoy, P. J., Zeller, K. A., Hicke, J. A., Quigley, H. B., Rabinowitz, A. R., and Thornton, D. H. (2016). Quantifying the effects of deforestation and fragmentation on a range-wide conservation plan for jaguars. Biol. Conserv. 203, 8–16. doi: 10.1016/j.biocon.2016.08.037
Pan, D., Li, Y., Hu, H. X., Meng, S. J., Men, Z. M., Fu, Y. X., et al. (2005). Microsatellite polymorphisms of Sichuan golden monkeys. Chin. Sci. Bull. 50, 2850–2855. doi: 10.1007/BF02899655
Peakall, R., and Smouse, P. E. (2012). GenAlEx 6.5: genetic analysis in Excel. Population genetic software for teaching and research-an update. Bioinformatics 28, 2537–2539. doi: 10.1093/bioinformatics/bts460
Phillips, S. J., Anderson, R. P., and Schapire, R. E. (2006). Maximum entropy modeling of species geographic distributions. Ecol. Modell. 190, 231–259. doi: 10.1016/j.ecolmodel.2005.03.026
Piry, S., Luikart, G., and Cornuet, J. M. (1999). Bottleneck: a computer program for detecting recent reductions in the effective size using allele frequency data. J. Hered. 90, 502–503. doi: 10.1093/jhered/90.4.502
Pritchard, J. K., Stephens, M., and Donnelly, P. (2000). Inference of population structure using multilocus genotype data. Genetics 155, 945–959. Available online at: https://www.genetics.org/content/155/2/945
Qi, X. G., Huang, K., Fang, G., Grueter, C. C., Dunn, D. W., Li, Y. L., et al. (2017). Male cooperation for breeding opportunities contributes to the evolution of multilevel societies. Proc. R. Soc. B 284:20171480. doi: 10.1098/rspb.2017.1480
Ren, R. M., Su, Y. J., Yan, K. H., Li, J. J., and Hu, Y. F. (1998). Preliminary Survey of the Social Organization of Rhinopithecus Roxellana in Shennongjia National Natural Reserve. Hubei: World Scientific Press. doi: 10.1142/9789812817020_0014
Rice, W. R. (1989). Analyzing tables of statistical tests. Evolution 43, 223–225. doi: 10.1111/j.1558-5646.1989.tb04220.x
Rousset, F. (2008). Genepop'007: a complete re-implementation of the genepop software for Windows and Linux. Mol. Ecol. Resour. 8, 103–106. doi: 10.1111/j.1471-8286.2007.01931.x
Taberlet, P., Griffin, S., Goossens, B., Questiau, S., Manceau, V., Escaravage, N., et al. (1996). Reliable genotyping of samples with very low DNA quantities using PCR. Nucleic Acids Res. 24, 3189–3194. doi: 10.1093/nar/24.16.3189
van Oosterhout, C., Hutchinson, W. F., Wills, D. P., and Shipley, P. (2010). Micro-checker: software for identifying and correcting genotyping errors in microsatellite data. Mol. Ecol. Resour. 4, 535–538. doi: 10.1111/j.1471-8286.2004.00684.x
Waage, J. K., and Greathead, D. J. (1988). Biological control: challenges and opportunities. Philos. Trans. R. Soc. Lond. B. Biol. Sci. 318, 111–128. doi: 10.1098/rstb.1988.0001
Wang, C. L., Wang, X. W., Qi, X. G., Guo, S. T., zhao, H. T., Wei, W., et al. (2014). Influence of human activities on the historical and current distribution of Sichuan snub-nosed monkeys in the Qinling Mountains, China. Folia Primatol. 85, 343–357. doi: 10.1159/000368398
Welt, R. S., Amy, L., and Franks, S. J. (2015). Analysis of population genetic structure and gene flow in an annual plant before and after a rapid evolutionary response to drought. AoB Plants 7, 632–641. doi: 10.1093/aobpla/plv026
Wen, H., and Wen, R. (2006). The Change of the Plant and Animal in China During Different Historical Period. Chongqing: Chongqing Press.
Wen, R. (2009). The Distributions and Changes of Rare Wild Animals in China. Jinan: Shandong Science and Technology Press.
Willi, Y., Van Buskirk, J., and Hoffmann, A. A. (2006). Limits to the adaptive potential of small populations. Annu. Rev. Ecol. Evol. Syst. 37, 433–458. doi: 10.1146/annurev.ecolsys.37.091305.110145
Wilson, G., and Rannala, B. (2002). Bayesian inference of recent migra-tion rates using multilocus genotypes. Genetics 163, 1177–1191. Available online at: https://www.genetics.org/content/163/3/1177
Yu, L., Wang, G. D., Ruan, J., Chen, Y. B., Yang, C. P., Cao, X., et al. (2016). Genomic analysis of snub-nosed monkeys (Rhinopithecus) identifies genes and processes related to high-altitude adaptation. Nat. Genet. 48, 947–952. doi: 10.1038/ng.3615
Zhang, Y. B., Wang, Y. Z., Phillips, N., Ma, K. P., Li, J. S., and Wang, W. (2016). Integrated maps of biodiversity in the Qinling Mountains of China for expanding protected areas. Biol. Conserv. 210, 64–71. doi: 10.1016/j.biocon.2016.04.022
Keywords: Rhinopithecus roxellana, population structure, gene flow, ecological niche models, evolutionary history
Citation: Li Y, Huang K, Tang S, Feng L, Yang J, Li Z and Li B (2021) Genetic Structure and Evolutionary History of Rhinopithecus roxellana in Qinling Mountains, Central China. Front. Genet. 11:611914. doi: 10.3389/fgene.2020.611914
Received: 29 September 2020; Accepted: 30 November 2020;
Published: 20 January 2021.
Edited by:
Deyan Ge, Chinese Academy of Sciences (CAS), ChinaReviewed by:
Li Yu, Yunnan University, ChinaWei Wang, Chinese Research Academy of Environmental Sciences, China
Copyright © 2021 Li, Huang, Tang, Feng, Yang, Li and Li. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Baoguo Li, YmFvZ3VvbGkmI3gwMDA0MDtud3UuZWR1LmNu
†These authors have contributed equally to this work