Noninvasive sampling reveals population genetic structure in the Royle’s pika, Ochotona roylei, in the western Himalaya

Abstract Understanding population genetic structure of climate‐sensitive herbivore species is important as it provides useful insights on how shifts in environmental conditions can alter their distribution and abundance. Herbivore responses to the environment can have a strong indirect cascading effect on community structure. This is particularly important for Royle's pika (Lagomorpha: Ochotona roylei), a herbivorous talus‐dwelling species in alpine ecosystem, which forms a major prey base for many carnivores in the Himalayan arc. In this study, we used seven polymorphic microsatellite loci to detect evidence for recent changes in genetic diversity and population structure in Royle's pika across five locations sampled between 8 and 160 km apart in the western Himalaya. Using four clustering approaches, we found the presence of significant contemporary genetic structure in Royle's pika populations. The detected genetic structure could be primarily attributed to the landscape features in alpine habitat (e.g., wide lowland valleys, rivers) that may act as semipermeable barriers to gene flow and distribution of food plants, which are key determinants in spatial distribution of herbivores. Pika showed low inbreeding coefficients (F IS) and a high level of pairwise relatedness for individuals within 1 km suggesting low dispersal abilities of talus‐dwelling pikas. We have found evidence of a recent population bottleneck, possibly due to effects of environmental disturbances (e.g., snow melting patterns or thermal stress). Our results reveal significant evidence of isolation by distance in genetic differentiation (F ST range = 0.04–0.19). This is the first population genetics study on Royle's pika, which helps to address evolutionary consequences of climate change which are expected to significantly affect the distribution and population dynamics in this talus‐dwelling species.


| INTRODUC TI ON
Over the past century, climate has been changing at a faster pace in alpine habitats than in other ecosystems, for which many species must either adapt or migrate to areas with optimal conditions (Naftz et al., 2002;Shrestha, Gautam, & Bawa, 2012). Among mountaindwelling species, pikas (Ochotona species), in particular, have experienced a significant contraction in their lower elevation range (Moritz et al. 2008). For example, climate-induced heat (Beever, Brussard, & Berger, 2003), cold, or nutritional stress (Wilkening, Ray, Beever, &  Brussard, 2011) led to historical extinctions of local populations and recent range contraction. While dispersal to higher elevation may allow pikas to persist in suitable microclimate, habitat fragmentation or non-availability of preferred food plants can impair their ability to cope with the changing environment and pose a threat to their survival and fitness (Bhattacharyya, Dawson, Hipperson, & Ishtiaq, 2018;Ray, Beever, & Loarie, 2012;Schloss, Nuñez, & Lawler, 2012;Walther et al., 2002).
The Royle's pika (Ochotona roylei) is a widespread small herbivore, occupying much of the Himalayan arc-from northwestern Pakistan to India (Jammu and Kashmir, Himachal Pradesh, and Uttarakhand), Nepal, and adjacent Tibet (Bhattacharyya & Smith, 2018) and inhabits rocky boulder or talus habitat between 2,400 and 5,200 m elevation. It is an asocial mammal (Bhattacharyya & Smith, 2018, Figure 1) and, unlike other talus-dwelling pika species (e.g., American pika Ochotona princeps, Collared pika Ochotona collaris), does not exhibit haying activity to store winter food resources (Bhattacharyya, Adhikari, & Rawat, 2013). The Royle's pika usually occurs at low population densities from 12.5 per ha in Nepal to 16.2 per ha in the western Himalaya in India (Bhattacharyya, Adhikari, & Rawat, 2014a).
The population genetic structure of small mammals is often influenced by landscape features and habitat fragmentation (Gerlach & Musolf, 2000;Peacock & Smith, 1997a;Peakall, Ruibal, & Lindenmayer, 2003). For talus-dwelling species, topographical features (e.g., aspect) and alpine habitat connectivity play a crucial role in determining increased gene flow (e.g., American pikas; Castillo, Epps, Davis, & Cushman, 2014). Similarly, microhabitat characteristics such as rock cover area, availability of rock talus with small crevice size (<15 cm), govern the habitat occupancy in Royle's pika (Bhattacharyya et al., 2015). A large area of talus supports greater plant diversity and results in higher species richness in the pika's diet . Therefore, we expect that Royle's pika populations would be influenced by topographical barriers (e.g., river low land valleys) as well as forage availability (e.g., C 3 plants).
Philopatric settlements (dispersal very close to their natal territory) in pikas often lead to incestuous matings and low intrapopulation genetic variability (Henry, Sim, & Russello, 2012;Peacock & Smith, 1997a, 1997bSmith & Weston, 1990). While we lack information on sex-biased dispersal or philopatry in Royle's pika, being a rocky talus obligate species, Royle's pika could have high inbreeding coefficients and relatedness across spatially close talus habitats or signatures of historical connectivity in habitat.
Population dynamics of Royle's pika are often influenced by snow melting pattern as snow acts as a thermal insulator to pika and their food plants (Bhattacharyya et al., 2014a). With recent  changes in snow patterns and rising temperatures in winters in the Himalaya region (Shrestha et al., 2012), these could result in patchy distribution of pika, loss of genetic diversity, and a demographic bottleneck in the species (Bhattacharyya et al., 2014a).
In the Himalayan arc, seven Ochotona species have been reported (Smith, Formozov, Hoffmann, Changlin, & Erbajeva, 1990); however, detailed long-term studies on pika ecology (see Bhattacharyya et al., 2013;Bhattacharyya et al., 2014a;Bhattacharyya et al., 2014b;Bhattacharyya et al., 2015) and population genetics are yet to be considered. A fine-scale pattern of distribution and habitat selection of the Royle's pika has been studied in the western Himalaya for the past nine years (Bhattacharyya et al., 2013(Bhattacharyya et al., , 2014a(Bhattacharyya et al., , 2014b(Bhattacharyya et al., , 2015. Using genetic markers in combination with fine-scale spatial ecological data provides a unique opportunity to explore the genetic structure, patterns of relatedness, and demographic responses to changing environment. We employed a panel of seven polymorphic markers (see Table 1 for details) to understand the population structure and gene flow in the Royle's pika in the western Himalaya. We estimated genetic diversity between and within populations and explored evidence of inbreeding, isolation by distance, population genetic bottleneck, all of which will provide insights on how shifting F I G U R E 1 Royle's pika (a) and its talus habitat (b) in Kedarnath Wildlife Sanctuary, India. Photograph credit: Sabuj Bhattacharyya

| Field sampling
The Royle's pika is usually found in talus habitats around forest trails with a home range of approximately 50 m 2 (Bhattacharyya et al., 2015;Kawamichi, 1968). We surveyed talus habitats using noninvasive fecal sampling across five locations spanning an eleva-

| Molecular methods
The Royle's pika often feeds on plants with high levels of secondary metabolites (Bhattacharyya et al., 2013 (Alves et al., 2015;Castillo et al., 2014;Li et al., 2009;Mougel et al., 1997;Peacock et al., 2002;Zgurski et al., 2009; see Table 1). The forward sequences of each primer were labeled with fluorescent dyes (FAM, NED, VIC, PET). The markers were multiplexed after the PCR amplification stage, as annealing temperature varied across loci. Each sample for each locus was amplified a TA B L E 1 Genetic diversity in Royle's pika populations in the western Himalaya minimum of four times to control for allelic dropout. Amplified products were diluted and run on an ABI 3730 automated sequencer and analyzed in GENEIOUS version 9.0.4 (https://www.geneious.com/).

| Estimation of genetic diversity
The presence and frequency of null alleles and polymorphism final genotypic dataset, we assigned all samples to unique individuals using probability of identity statistics (PID; Paetkau et al., 1998) in CERVUS version 2.0. Any sample, which showed similarity at five loci, was considered a duplicate sample and was removed from the final dataset (see Cullingham et al., 2016). All loci were tested for departures from Hardy-Weinberg equilibrium (HWE) using CERVUS version 2.0 and used Holm's Bonferroni sequential corrections to adjust p values in HWE estimations (Gaetano, 2013). We tested for linkage disequilibrium using GENEPOP version 4.6 (Raymond, 1995;Rousset, 2008). We estimated likelihood of allelic dropout rate and false allele rate within each population with 10,000 search steps (Johnson & Haydon, 2007), using PEDANT version 1.0. We used ESTIMATES version 9.0 with 10,000 iterations (Colwell et al., 2012) and the nonparametric Chao2 estimator (mean ± SD; Colwell & Coddington, 1994) to estimate our success in sampling alleles in populations (with more than 13 samples). We constructed rarefaction curves plotting the cumulative number of alleles found with increasing sample size.

| Spatial autocorrelation
We examined the relationship between genetic relatedness and geographic distance using Wang's estimator "r" (Wang, 2002) in SPAGeDi version 1.4 (Hardy & Vekemans, 2002). This estimator is preferable among all relatedness indices due to its low sensitivity to sampling error in allele frequency calculation and has low sampling variance due to change in number of loci or alleles (Blouin 2003;Robinson, Simmons, & Kennington, 2013, Wang, 2017 (Hardy & Vekemans, 2002). The standard error for each distance class was estimated using a jackknife procedure over loci (Hardy & Vekemans, 2002), and deviation of Wang's estimate (r) from 0 suggests that Royle's pika individuals within a specific distance class are significantly related (positive values = high relatedness, negative values = low relatedness) and not random.

| Genetic differentiation in pika populations
We used five different methods to assess the population structure and differentiation in Royle's pika. The first two methods used a clustering approach: STRUCTURE and the discriminant analysis of principal components (DAPC) and three methods were based on distances: F ST estimations, the analysis of the molecular variance (AMOVA), and isolation by distance (IBD). We included one locus that showed HWE deviation as the deviation was not observed in all populations.
We first investigated the likelihood of each individual sample belonging to one of several clusters (K) based on allele frequencies, using Bayesian clustering methodology in STRUCTURE version 2.3.4 (Pritchard, Stephens, & Donnelly, 2000). We implemented the admixture ancestry model with correlated allele frequencies. The putative numbers of population clusters (K) were allowed to vary from one to 10 with 50,000 burn-in iterations and 500,000 Markov chain Monte Carlo (MCMC) iterations (Ishtiaq, Prakash, Green, & Johnson, 2015;Singh et al., 2017). A total of 20 independent runs were performed for each K (1-10) to achieve consistency across the runs.
We used Structure Harvester (Earl & VonHoldt, 2012) to calculate Delta K (Evanno, Regnaut, & Goudet, 2005)  We performed DAPC with the R package "Adegenet" (Jombart, 2008). DAPC define clusters using clustering algorithms K-means on transformed data with principal component analysis (PCA). The clustering identifies groups of populations with similar genotypes by specifying the actual number of clusters (K = 5) as a priori information (Putman & Carbone, 2014;Zachos et al., 2016). We retained 32 principal components of PCA, which explained approximately 90% of the total variation of our dataset. We used q > 0.7 to assign individuals to each location.
We used AMOVA to determine hierarchical distribution of genetic variation in ARLEQUIN. To run this analysis, we a priori defined five groups according to sampling sites and talus sampled within these sites. The significance of AMOVA was tested using 1,000 replicate bootstrap of all data.
To determine spatial patterns driving genetic differentiation among locations, linearized pairwise F ST among talus (F ST /(1−F ST ) obtained for seven microsatellite loci were correlated against log-transformed geographic distances between talus using a Mantel test with 1,000 permutations in ARLEQUIN.

| Demographic effects of environmental conditions and habitat (bottleneck)
The detection of evidence for any population bottleneck was primarily based on Wilcoxon's signed rank test and shift in allelic distribution (Luikart, Allendorf, Cornuet, & Sherwin, 1998). We selected locations (TUN, MAD, and NAN) with sample size >10 individuals. We used BOTTLENECK version 1.2.02 to assess the genetic signature of demographic contractions using the heterozygote excess test (Cornuet & Luikart, 1996;Piry, Luikart, & Cornuet, 1999). The heterozygote excess test determines whether sudden decline in population size has led to the loss of rare alleles which can cause an expected heterozygosity excess ; Luikart, Sherwin, Steele, of a recent bottleneck in a population could be indicated by an increase in the number of heterozygotes relative to the number of alleles, whereas heterozygote deficiency indicates population growth (Cornuet & Luikart, 1996). Six out of seven microsatellite primers used in the study were dinucleotide. Hence, a two-phase mutational model (TPM) with 95% single step mutations (SMM) and 12% variance among mutational steps, and 1,000 simulations for each location were used, which is generally recommended practice for dinucleotide microsatellite repeat loci (Di Rienzo et al., 1994;Piry et al., 1999).

| Genetic diversity
A total of 203 fecal samples were collected from five locations. Of these samples, 68.96% were successfully amplified and produced consistent results across all replicates (Supporting Information Table   S1). The average missing data across all loci were 2.95%, ranging from 0.71% (SAT3) to 7.14% (Ocp6). The rarefaction curve for allelic richness for each microsatellite locus indicated that we were able to detect all possible alleles across sampled populations (Supporting Information Table S2). No significant linkage disequilibrium was detected for any locus, whereas deviation from HW equilibrium was detected only in one out of seven loci in RUD and HAK populations after Bonferroni correction (Table 1). We retained all loci for further analysis as no systematic deviation of HWE was observed across all sites.

| Spatial autocorrelation
The pairwise relatedness analysis suggested local-scale genetic structure in pika populations. The negative regression slope (b = −0.22 ± 0.005, p < 0.05) between Wang's estimator (r) and geographic distance suggested neighboring individuals within 1 km distance are more genetically related than any random pair of individuals (Figure 3).

| Genetic differentiation in pika populations
The STRUCTURE results suggested a best fit of K = 2 clusters. The maximum Delta was found at K = 2, followed by K = 3 (Figure 4a,b).
A large proportion of individuals (80%) were successfully assigned (q > 0.7) at K = 2 (Supporting Information Table S3), whereas at K = 3.5%, a large proportion of individuals were assigned (q > 0.7) to specific clusters. In K = 2, individuals in NAN (97%) and HAK (83.34%) were assigned in genetic cluster "A" and MAD (38.46%), TUN (69%), and RUD (57.14%) were assigned in cluster "B" (Figure 4c, Supporting Information Table S3). The DAPC suggested K = 3 as the most probable number of clusters which indicated genetic differentiation in HAK as a distinct cluster, whereas NAN, MAD, TUN, and RUD either partially or completely overlapped (Figure 2b).
The AMOVA suggested that a majority of the genetic variance was attributed to variation within individuals (84.45%), followed by 6.05% of the variation among locations, and 6.84% variation among talus within each location (Table 3). TUN and MAD were the only two locations with significantly low genetic differentiation (F ST = 0.02, p < 0.01, Table 4), whereas other paired locations showed significantly high genetic differentiation (F ST = 0.04−0.16, p < 0.01, Table 4).
The Mantel test showed significant evidence for a decrease in genetic differentiation across sampled individuals with an increase in geographic distance (r = 0.033, p < 0.001; Figure 2c).

| Demographic change (bottleneck) in pika populations
We found significant (p < 0.05) heterozygosity excess across three locations suggesting that demographic change leads to a recent bottleneck (Table 5).

| D ISCUSS I ON
In this first genetic study on the Royle pika in south Asia, we used a noninvasive sampling approach where the quality of samples often relies on environmental conditions (Lucchini et al., 2002;Murphy, Kendall, Robinson, & Waits, 2007), sample preservation method (Nsubuga et al., 2004), and diet of model species (Panasci et al., 2011). Royle's pika occurs in open rocky talus habitat where fecal pellets were often exposed to high temperature fluctuation, high UV radiation, and extreme weather conditions, all of which could potentially degrade the DNA quality. Despite these constraints, we were able to successfully genotype large number (68.96%) of samples which showed high PIC value, low allelic dropout, and false allele rates to establish accuracy and reliability in our dataset.

| Genetic variation and population structure in Royle's pika
Overall, the Royle's pika populations exhibited similar levels of genetic diversity (0.57) to those observed in other studies on talus-dwelling Ochotona species (see Table 2). However, our genetic diversity values are not directly comparable as we used relatively fewer genetic markers than other studies. Nonetheless, the observed patterns in our data reflect pika life-history traits rather than data quality. Natural populations with small numbers of individuals often lose genetic variability due to genetic drift and inbreeding (Lacy, 1997). The Royle's pika has a small litter size (Bhattacharyya et al., 2014a)  coefficients (Glover, Smith, Ames, Joule, & Dubach, 1977;Henry et al., 2012;Moilanen, Smith, & Hanski, 1998;Peacock, 1997;Peacock & Smith, 1997a;Smith & Ivins, 1983). The Royle's pika populations showed high genetic variation within individuals. Distribution of genetic diversity within and among populations is often found to be influenced by migration of individuals (Berthier, Charbonnel, Galan, Chaval, & Cosson, 2006;Hartl & Clark, 1997). Hence, populations with restricted movement and gene flow usually exhibit high genetic differentiation (Peery et al., 2008;Slarkin, 1985). (iii) limited distribution of rock cover to protect from heat stress and predation risk (Bhattacharyya et al., 2014b); (iv) high interannual variation in the snow melting pattern, resulting in low survival and pika abundance (Bhattacharyya et al., 2014a).
Across two cluster approaches, DAPC results were in congruence with F ST estimates due to limitation of STRUCTURE for detecting fine-scale genetic structure (Janes et al, 2017). Landscape features and topography (e.g., elevation, aspect), habitat availability, and temperature gradient are factors which influence pika dispersal and connectivity between metapopulations over 10 km apart (Castillo et al., 2014(Castillo et al., , 2016Henry et al., 2012;Robson, Lamb, & Russello, 2015;Smith, 1974

| Effects of habitat and snow patterns on demography (population bottleneck)
We found evidence of a recent bottleneck, which is possibly a result of significant recent changes in habitat structure and rising winter temperature. The distribution and depth of snow in mountain regions are often influenced by wind and topography (e.g., aspect, slope), and microclimate, causing significant variation in snow cover, even across small geographical distances (Deems, Birkeland, & Hansen, 2002;Gurung et al., 2017). Royle's pika does not hibernate during winter and snow cover acts as a thermal insulator against extreme cold and fluctuating temperatures (Bhattacharyya et al., 2014a). While warm winters with thin snow cover can cause severe mortality in adult pikas, F I G U R E 4 Nonspatial Bayesian individual-based clustering STRUCTURE results from seven polymorphic microsatellite loci and admixture model. Plots showing the modal value of best K using maximum likelihood method (mean L (K) ± SD; a) and using ΔK method (b), the maximum Delta K = 13.70 was found at K = 2, followed by K = 3 with a Delta K = 5.12; visualization of the Bayesian clustering analysis implemented in STRUCTURE for K = 2-K = 6 (c) where each K is shown as a different color. Individuals are shown as single vertical bars, grouped by geographic location; location codes (below the plot) are the same as Figure 1. Data shown are aggregated data from 20 independent runs thin snow cover during the spring reduces reproductive fitness and thereby delays breeding dates (Bhattacharyya et al., 2014a;Morrison & Hik, 2007). In Royle's pika, talus size (area) and connectivity between talus are known to influence the habitat occupancy, whereas predation risk significantly impacts foraging ecology (Bhattacharyya et al., 2013(Bhattacharyya et al., , 2015 and potentially also access to nutritive plants, thereby affecting individual fitness. Furthermore, the decreased precipitation (3-4 mm per year from 1982 to 2006; Shrestha et al., 2012) during spring affects the growth rates of nutritive C 3 plants in summer, which form the main diet of young pika .
Climate-induced habitat loss and range contractions might be the main drivers for the apparent population bottleneck in pika.
Contemporary climate change is expected to significantly affect the distribution and population dynamics in this talus-dwelling species.

| CON CLUS IONS
This is the first genetic study to establish the population structure of an important lagomorph species in the western Himalaya. Using clustering approaches, we found evidence for well-defined population structure. We detected moderate levels of inbreeding and evidence for a recent population bottleneck and genetic isolation by geographic distance. We found little evidence of gene flow was observed between individuals >1 km apart. Therefore, being a climatesensitive small mammal which lives in isolated patches of mountain habitat and having limited dispersal ability, Royle's pika may depend on local adaptation in order to survive changing environmental conditions in the future. Our study highlights the need for an in-depth study of this high-altitude-restricted model species and presents new evidence that it is imminently threatened by climate change.

E TH I C A L A PPROVA L
The field experiments comply with the current laws of India where the study was performed. We thank Uttarakhand Forest Department for ethical approval and permission for collection of fecal samples. Note. Population codes are the same as Figure 1a.

Sum of squares
TA B L E 5 Wilcoxon sign-rank test and change in allelic distribution to detect heterozygosity excess due to significant reduction in Royle's pika population in the western Himalaya TA B L E 3 Analysis of molecular variance (AMOVA) for five sampled locations of Royle's pika

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

AUTH O R CO NTR I B UTI O N
FI and SB designed the experiment. SB conducted the field and laboratory experiments. SB and FI analyzed the data and wrote the manuscript. Both authors approved the final version of the manuscript.

DATA ACCE SS I B I LIT Y
The GPS locations of all sample collection points have already been added in Supporting Information. The microsatellite genotype data(allele frequency