Genetic Structure and Evolutionary History of Rhinopithecus roxellana in Qinling Mountains, Central 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.

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.
Keywords: Rhinopithecus roxellana, population structure, gene flow, ecological niche models, evolutionary history 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., 2002Li et al., , 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).
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.

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(Li et al., , 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. 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 tetranucleotide 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 (A O ), observed heterozygosity (H O ), expected heterozygosity (H E ), polymorphic information content (PIC), allelic richness (A R ), 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 θ (F ST ) 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).

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.

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 (H O ) ranged from 0.563 to 0.624, averaging 0.595; expected heterozygosity (H E ) 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 (F IS ) 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).

Population Differentiation and Structure
The results of F ST 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.

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;

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).

Ecological Niche Modeling
We obtained four distributions tendencies of the Rhinopithecus roxellana, including the LIG, LGM, HOL, and current models by MAXENT (Figure 4)

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 F ST 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 . 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 . 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.