Population genetics reveals high connectivity of giant panda populations across human disturbance features in key nature reserve

Abstract The giant panda is an example of a species that has faced extensive historical habitat fragmentation, and anthropogenic disturbance and is assumed to be isolated in numerous subpopulations with limited gene flow between them. To investigate the population size, health, and connectivity of pandas in a key habitat area, we noninvasively collected a total of 539 fresh wild giant panda fecal samples for DNA extraction within Wolong Nature Reserve, Sichuan, China. Seven validated tetra‐microsatellite markers were used to analyze each sample, and a total of 142 unique genotypes were identified. Nonspatial and spatial capture–recapture models estimated the population size of the reserve at 164 and 137 individuals (95% confidence intervals 153–175 and 115–163), respectively. Relatively high levels of genetic variation and low levels of inbreeding were estimated, indicating adequate genetic diversity. Surprisingly, no significant genetic boundaries were found within the population despite the national road G350 that bisects the reserve, which is also bordered with patches of development and agricultural land. We attribute this to high rates of migration, with four giant panda road‐crossing events confirmed within a year based on repeated captures of individuals. This likely means that giant panda populations within mountain ranges are better connected than previously thought. Increased development and tourism traffic in the area and throughout the current panda distribution pose a threat of increasing population isolation, however. Maintaining and restoring adequate habitat corridors for dispersal is thus a vital step for preserving the levels of gene flow seen in our analysis and the continued conservation of the giant panda meta‐population in both Wolong and throughout their current range.


| INTRODUC TI ON
Rare and elusive large-bodied mammal populations intrinsically occur at low densities (Mumma, Zieminski, Fuller, Mahoney, & Waits, 2015;Taberlet & Bouvet, 1992) and face increasing threats from climate change and anthropogenic influences Zhu et al., 2013). Managers are frequently tasked with monitoring population sizes, distributions, and connectivity in order to guide management actions. Noninvasive genetic sampling (NGS) is increasingly being used in the conservation and management of threatened animals, as it allows for the estimation of important population parameters such as total size, genetic diversity, and gene flow among populations (Barba et al., 2010;Schregel et al., 2012;Wang et al., 2016;Zhan et al., 2006). Analyses of gene flow grant inference about the functional connectivity of a landscape and have important implications for conservation. Maintaining adequate connectivity both helps to maintain genetic diversity in small subpopulations (Sharma et al., 2012) and allows for recolonization of areas that undergo localized extinctions (Hanski, 1998).
The giant panda (Ailuropoda melanoleuca) (Figure 1) is an example of a species that has faced historical population declines and been the focus of intensive conservation effort through the establishment of protected areas and habitat restoration (Tuanmu et al., 2016).
Although there is evidence of recent population recovery which resulted in the reduction of their extinction risk on the IUCN red list (Swaisgood, Wang, & Wei, 2017), pandas still face ongoing increases in habitat fragmentation and subpopulation isolation (Xu et al., 2017;Yang et al., 2017). Currently, their occupancy has been reduced to the eastern edge of the Tibetan plateau in six separate mountain ranges (Schaller, Hu, Pan, & Zhu, 1985). Within these mountain ranges, major rivers, roads, and habitat loss are estimated to further segregate panda populations into 33 subpopulations (State Council Information Office).
Road development in particular has increased substantially across the giant panda range. While roads cover seemingly small proportions of land surface, they affect the environment in various ways, such as through the loss of suitable habitat, animal mortality, acting as barriers to individual movements, and causing landscape fragmentation (Balkenhol & Waits, 2009;Fahrig, 2004;Zhao et al., 2016). These effects can act to increase genetic structure between populations and decrease genetic diversity within populations, which further reduces population viability (Keyghobadi, 2007). This was found to be the case in the giant panda population occupying the Xiangling Mountains, which exhibited genetic differentiation on either side of a major road (Zhu, Zhan, Meng, Zhang, & Wei, 2010;Zhu, Zhang, Gu, & Wei, 2011). whether the southern and northern subpopulation are connected via effective migration across the road. This would mean that these subpopulations in Wolong act as a single meta-population, which has implications for their persistence and management. Such empirical evidence of population connectivity is also largely lacking across the giant panda range, as analyses within the occupied mountain ranges have rarely employed genetics methods (Shen et al., 2008;Xu et al., 2006). The main objective of this study was to thus evaluate the genetic connectivity, through the presence or absence of migration, of the panda subpopulations on either side of the main road through Wolong. We also endeavored to use NGS methods to determine the population size, distribution, and genetic diversity of the panda population in Wolong to better understand their ecology and inform effective conservation. We hypothesized that the population size would be fairly large due to the widespread availability of understory bamboo habitat in the reserve, and that this would translate to high levels of genetic diversity. That said, we also hypothesized that there would be a detectable effect of the road on gene flow. We expect our results to have implications for the conservation of remaining giant panda populations both within Wolong and throughout their range.

| DNA extraction and amplification
Total genomic DNA was extracted from fecal samples using QIAamp DNA stool mini kits (Qiagen, Germany), according to the manufacturer's instructions. We used seven tetra-microsatellite loci to distinguish among individuals. These were as follows: GPL-60, gpz-20, GPL-29, gpz-6, GPL-53, GPL-44, and gpz-47 (Huang et al., 2015). The probability of identity across these loci in the target population was estimated using GIMLET 1.3.3 (Valière, 2002). PCR amplifications were carried out in 25 μl reaction mixtures comprising approximately 50 ng of template DNA, 2 mm MgCl 2 , 200 μmol of dNTP each, 15 pmol of each primer, 1.0 μg of bovine serum albumin (BSA), and 0.3 units of Hotstart DNA polymerase (Takara). Amplifications were performed using the following PCR procedure: an initial denaturation step for 5 min at 95°C, followed by 35 cycles of 95°C for 45 s, 30 s at locus-specific annealing temperature and 50 s at 72°C, and a final elongation for 10 min at 72°C. For genotyping, the PCR amplification products were separated by capillary electrophoresis using a denaturing acrylamide gel matrix on an ABI 3730xl Genetic Analyzer. Alleles were detected using Genemapper 3.2 software.

| Quality control
Genotyping errors caused by amplification of poor quality DNA from fecal samples such as allelic dropout and false alleles can severely bias estimates of population parameters (Broquet & Petit, 2004;Lampa, Henle, Klenke, Hoehn, & Gruber, 2013). Therefore, we performed control measures to ensure the quality of our genetic data. All fecal samples were amplified at least three times for each marker. A single-locus genotype was not accepted until our replicates resulted in at least three identical homozygote profiles or two identical heterozygote profiles. These criteria were based on a pilot study, where genotypes obtained from feces versus blood samples were compared (Huang et al., 2015). Huang et al. (2015) concluded that the seven loci used in this study always featured exact matches between blood and fecal samples of n = 15 captive pandas, and that results from feces exposed to the natural environment for up to 5 weeks (longer than the estimated 2-week cutoff of our study) were consistent. As an additional quality control, we used MICRO-CHECKER to search for loci with large allele dropout and scoring errors caused by stutter peaks (Oosterhout, Hutchinson, Wills, & Shipley, 2004). No evidence of allelic dropout or scoring error due to stuttering was found for any locus. Finally, we used FreeNA to estimate null allele frequencies for each locus (Chapuis & Estoup, 2007).
There was an average null allele frequency of <0.04 across the 7 loci.

| Estimation of population size
Individual genotypes were identified with the MStools plugin for Microsoft Excel using the following rules: (a) Genotypes from different samples were believed to represent the same individual if all alleles in all loci were identical. (b) If only one allele was found to differ between individuals, DNA was re-extracted and three more PCR replication was performed. If the allele was still different, we judged the samples as belonging to different individuals. (c) If there were differences of two or more alleles, the samples were accepted as belonging to different individuals.
The noninvasive records of individual genotypes throughout an area can be used to estimate the total population size via capturemark-recapture methods (Gervasi et al., 2008;Lukacs & Burnham, 2005;Mumma et al., 2015). We used the identification of different F I G U R E 3 Sampling locations of giant panda feces in Wolong nature reserve, China individuals through the fecal genetics data to build a CMR model and estimate the giant panda population size in Wolong using the "CAPWIRE" package (Pennell, Stansbury, Waits, & Miller, 2013) in the R programming environment. CAPWIRE performs population size estimation as well as or better than other abundance estimators when the data contain multiple observations of an individual within a session and there are <200 individuals (Miller, Joyce, & Waits, 2005;Mumma et al., 2015). Because our fecal collection efforts focused on all the giant panda's suitable habitat, we inferred that recapture probability was even among all individuals. We thus ran models under the assumption of equal capture (ECM) probabilities in CAPWIRE to estimate the population size. Because our study population was not closed during the study period and there was likely migration across the Northern and Southern borders, we also used the R package "secr" to employ spatially explicit capture-recapture (SECR) methods to estimate a density of pandas per square km across our study area (Efford, 2013). We used the polygon trap format corresponding to the 520 survey cells and grouped the data into 1 session of 30 sampling occasions based on proximity of collection time. We then multiplied the estimated density by the sample area to get an estimate of the number of pandas in Wolong.

| Population genetics analysis
The number of alleles (A), observed heterozygosity (Ho), expected heterozygosity (He) and polymorphic information content (PIC) were calculated at individual loci and across loci using the software CERVUS 3.0 (Marshall, Slate, Kruuk, & Pemberton, 1998).
A Bayesian clustering method implemented in Structure 2.3.1 (Pritchard et al., 2000) was used to determine the most likely number of genetic clusters. The admixture model was chosen, allele frequencies were assumed correlated, and analysis was conducted with a burn-in of 100,000 and followed by 1,000,000 MCMC repetitions.
Ten independent runs were carried out for each cluster set (K), from 1 to 4. The most likely K value was determined by evaluating the log likelihood [In P(X/K)] of the posterior probability of the data for each value of K (Pritchard et al., 2000). In addition, the △K statistic, the second-order rate of change in the log probability of the data between successive values of K, was estimated and used to determine the most likely number of genetic clusters (Evanno, Regnaut, & Goudet, 2005). To cross-validate the results of STRUCTURE, we also conducted a principal coordinates analysis (PCoA) using GenALEx V6.5 (Peakall & Smouse, 2012). In this analysis, multivariate genetic distances between individuals (Smouse & Peakall, 1999) are decomposed through PCoA to find sources of genetic variation across the population.
For quantifying genetic differentiation between populations, we estimated an F ST and its significance value through resampling 10,000 permutations of the genotypes between populations to derive a null distribution using Arlequin 3.5 (Excoffier & Lischer, 2010 (Sun & Chang, 2016). Simulations have shown that F ST performs better than other indices of population differentiation, as it is more sensitive in detecting population genetic processes when the mutation rate is high relative to the migration rate (Whitlock, 2011). GeneClass v.2.0 was used to detect first-generation migration of individuals across the road. Specifically, we assigned the two populations on either side different identities before applying Bayesian likelihood-based test statistics to compute the probability of an individual originating in one of the populations with Monte Carlo resampling of 10,000 simulated individuals at an alpha value of 0.01 (Piry et al., 2004). In this analysis, we estimated the ratio of the likelihood that an individual is of the same population from which it was sampled (L_home) divided by the ratio of the highest likelihood of the individual's assignment to any sampled population (L_max) to detect migrants.
The Triadic maximum likelihood (TrioML) estimator and the QuellerGt moment estimator, implemented in Coancestry 1.0, were used to calculate the inbreeding coefficient (f) for each individual and pairwise relatedness value (r) between two individuals, respectively (Wang, 2011). The individual inbreeding coefficient reflects the extent to which their parents are genetically related: f < 0.125 has been defined as low inbreeding, 0.25 > f > 0.125 as moderate, and f ≥ 0.25 as high inbreeding (Marshall et al., 2002). A smaller negative pairwise relatedness value (r) suggests distant kinship, while a larger positive value suggests closer kinship. Offspring of two individuals with a high pairwise relatedness value (r) have a high risk of inbreeding deficiencies.

| Spatial density pattern
We used the kernel density estimation (KDE) function in ArcGIS 10.0 to quantify the spatial pattern of giant panda density (Bailey & Gatrell, 1995). Previous studies in Wolong have estimated the diameter of giant panda home ranges to fall between 1 and 3 km (Guan et al., 2016;Hu, 2001;Schaller et al., 1985). Supposing that a given giant panda's home range is a circle, the 142 identified panda's GPS locations were used to denote the center of the circle with a radius of 3 km (for pandas with multiple recaptures, we only used the site of the first discovery for this analysis). This circle represented the maximum likely area that a giant panda might utilize. Closer regions to the observed panda locations represent areas of higher probable activity frequencies, which are reflected in the kernel density output. The density map was divided into three tiers (low, medium, and high), indicating different levels of density of inferred space use.

| Sampling and molecular identification
Of the 520 survey cells, we found fresh fecal samples in 140. In total, 539 fresh fecal samples were noninvasively collected for genetic analysis from the entire study area during two sampling sessions in the years 2015-2016 (Figure 3). Successful genotyping of 6 or more microsatellite loci was obtained for 322 samples (with three samples that were successfully genotyped at only 5 loci). The probability of two individuals who were full siblings sharing an identical multi-locus genotype was 0.00808 based on 6 loci, indicating that this subset of the original 7 loci was enough for accurate individual identification (PIDsib < 0.01) (Waits, Luikart, & Taberlet, 2001). Although using 5 loci resulted in a PIDsib of 0.015, the three samples that were only successfully genotyped at 5 loci were included in the analysis because of large spatial distance between them and other samples (>2 km).

| Spatial density pattern
There were four areas with relatively high densities of giant pandas in Wolong: the Tiantaishan, Hetaoping-Niutoushan, Wuyipeng, and Xihe areas, ordered from north to south (Figure 4). The large home range overlap and spatial proximity of separate areas of activity indicate that in the absence of strong resistance or barriers to movement, the giant pandas in Wolong constitute a relatively continuous population. Genetic relatedness analysis revealed that 68.84% of genotyped pairs had an estimated relatedness value of r < 0.125, with an average of -0.00013 for the whole population. These results indicate that no significant genetic differentiation has occurred around the G350 national road.

| D ISCUSS I ON
Previous capture/recapture studies of giant pandas that have used genetic markers have resulted in population size estimates that have exceeded those of other methods, with Zhan et al. (2006) estimating nearly double the population size in Wanglang Nature Reserve compare to the 3rd national survey. These methods have been criticized in the past for potential violations of CMR model assumptions, including population closure and genotyping error (Garshelis et al., 2008). Our CMR model estimate of giant panda population size within Wolong was also larger than results from the latest national survey, though not as drastically so (slightly over 50%). The potential for genotyping errors was explicitly addressed in our analysis (see  Estimating and evaluating genetic variation is critical for the effective evaluation and management of endangered populations and species (Caniglia, Fabbri, Galaverni, Milanesi, & Randi, 2014;Du et al., 2016;Wang et al., 2016). Populations with higher genetic diversity are often inferred to have greater capacity to adapt to environmental change (Frankham, 2005). Our analysis of genetic variation in the giant panda population in Wolong revealed relatively high levels of genetic diversity with large MNA, He, and PIC values. Although the data are not directly comparable because different microsatellite markers were used, comparisons of Wolong populations to other mountain populations suggested that genetic diversity in Wolong ranks relatively high (Supporting information Table S1). On a larger scale, a high level of genetic variation was also confirmed by genomewide SNP analysis of 34 wild pandas' whole genomes, suggesting pandas have large evolutionary potential (Li, Fan, Tian, & Zhu, 2010;Wei et al., 2015;Zhao et al., 2013).
Conservation geneticists emphasize the need to prevent the occurrence of inbreeding in endangered species because it is typically associated with decreased fertility and survival (Deborah & John, 2009;Keller, F., Waller, & Donald, M., 2002;Stevenr et al., 2008).
As reflected by our estimated metrics, inbreeding is at a moderate Our findings that the pandas were able to cross the national road bisecting Wolong, and that this road has not resulted in genetic differentiation between the populations on either side, differ from road effects found in the Xiangling Mountains. Zhu et al. (2010) found that the national road G108 has resulted in a significant degree of genetic differentiation in the giant panda populations there. The smaller effect of road G350 on local panda populations in Wolong (Fst = 0.021) compared to those of G108 (Fst = 0.033) could be due to smaller traffic volumes of G350, which was a provincial road S303 before 2017. Generally, wider roads with greater volumes of highspeed traffic affect wildlife populations more strongly than small, less travelled roads (Clevenger, Chruszcz, & Gunson, 2001;Jaarsma, Langevelde, & Botma, 2006).
Our results are inconsistent with previous studies that have suggested that due to the impact of major roads coupled with the destruction of vegetation nearby, the habitat and panda populations in the Qionglai Mountains have been fragmented into four blocks (Xu et al., 2006). This is directly related to previous assumptions that there has been a lack of gene exchange between the two subpopulations separated by unsuitable habitat and the road G350 in Wolong (Hu, 2001;Loucks et al., 2001;Schaller et al., 1985).   Xu et al., 2006), full subdivision of presently connected populations is likely an ongoing process.
The high levels of genetic diversity frequently seen in giant panda populations (Wei et al., 2015) should thus be seen as a resource to preserve, as well as supplement with reintroduction efforts (Yang et al., 2018). Recent molecular and behavioral investigations suggest that giant pandas rely primarily on adequate dispersal opportunities to avoid inbreeding through sex-biased dispersal (Hu et al., 2017). This study thus provides an effective base from which to continue to monitor the giant panda population in Wolong for small-scale dispersal and population connectivity patterns in the face of environmental changes. Continued study will also allow for the investigation of long-term population dynamics in order to better understand panda ecology and inform management both within the reserve and throughout the current panda distribution area in order to ensure the continued recovery of one of the world's foremost conservation icons.

CO N FLI C T O F I NTE R E S T
None declared.

DATA ACCE SS I B I LIT Y
The microsatellite data available from the Dryad Digital Repository: https://doi.org/10.5061/dryad.hf03sm4.