- Split View
-
Views
-
Cite
Cite
Jan E Janecka, Yuguang Zhang, Diqiang Li, Bariushaa Munkhtsog, Munkhtsog Bayaraa, Naranbaatar Galsandorj, Tshewang R Wangchuk, Dibesh Karmacharya, Juan Li, Zhi Lu, Kubanychbek Zhumabai Uulu, Ajay Gaur, Satish Kumar, Kesav Kumar, Shafqat Hussain, Ghulam Muhammad, Matthew Jevit, Charlotte Hacker, Pamela Burger, Claudia Wultsch, Mary J Janecka, Kristofer Helgen, William J Murphy, Rodney Jackson, Range-Wide Snow Leopard Phylogeography Supports Three Subspecies, Journal of Heredity, Volume 108, Issue 6, September 2017, Pages 597–607, https://doi.org/10.1093/jhered/esx044
- Share Icon Share
Abstract
The snow leopard, Panthera uncia, is an elusive high-altitude specialist that inhabits vast, inaccessible habitat across Asia. We conducted the first range-wide genetic assessment of snow leopards based on noninvasive scat surveys. Thirty-three microsatellites were genotyped and a total of 683 bp of mitochondrial DNA sequenced in 70 individuals. Snow leopards exhibited low genetic diversity at microsatellites (AN = 5.8, HO = 0.433, HE = 0.568), virtually no mtDNA variation, and underwent a bottleneck in the Holocene (∼8000 years ago) coinciding with increased temperatures, precipitation, and upward treeline shift in the Tibetan Plateau. Multiple analyses supported 3 primary genetic clusters: (1) Northern (the Altai region), (2) Central (core Himalaya and Tibetan Plateau), and (3) Western (Tian Shan, Pamir, trans-Himalaya regions). Accordingly, we recognize 3 subspecies, Panthera uncia irbis (Northern group), Panthera uncia uncia (Western group), and Panthera uncia uncioides (Central group) based upon genetic distinctness, low levels of admixture, unambiguous population assignment, and geographic separation. The patterns of variation were consistent with desert-basin “barrier effects” of the Gobi isolating the northern subspecies (Mongolia), and the trans-Himalaya dividing the central (Qinghai, Tibet, Bhutan, and Nepal) and western subspecies (India, Pakistan, Tajikistan, and Kyrgyzstan). Hierarchical Bayesian clustering analysis revealed additional subdivision into a minimum of 6 proposed management units: western Mongolia, southern Mongolia, Tian Shan, Pamir-Himalaya, Tibet-Himalaya, and Qinghai, with spatial autocorrelation suggesting potential connectivity by dispersing individuals up to ∼400 km. We provide a foundation for global conservation of snow leopard subspecies, and set the stage for in-depth landscape genetics and genomic studies.
The snow leopard (Panthera uncia), considered the world’s most elusive large felid, inhabits a vast area (∼1.6 million km2) across 12 countries in Asia (Jackson et al. 2008; McCarthy et al. 2016) (Figure 1). It is a high-altitude specialist that primarily occupies mountains above 3000 m in elevation (Hemmer 1972). This region is characterized by low oxygen levels, temperature extremes, aridity, low productivity, and harsh climatic condition, yet harbors many distinctive taxa, including the Tibetan fox (Vulpes ferrilata, Harris 2014), Chinese desert cat (Felis bieti, Riordan et al. 2015), argali (Ovis ammon, Harris and Reading 2008), markhor (Capra falconeri, Michel and Rosen 2015), urial (Ovis orientalis, Valdez 2008), and Tibetan antelope (Pantholops hodgsonii, IUCN SSC Antelope Specialist Group 2016). The endangered snow leopard is a flagship species for Asia, the largest carnivore in its high-altitude communities, and yet is under substantial threat throughout its range (Jackson et al. 2010; Rosen and Zeller 2016). Research on distribution (Jackson et al. 2006, Janecka et al. 2008, 2011a; McCarthy et al. 2008; Lovari et al. 2009; Karmacharya et al. 2011; Alexander et al. 2016), ecology (Jackson and Ahlborn 1989; McCarthy et al. 2005; Anwar et al. 2011; Lovari et al. 2013; Johansson et al. 2016; Chetri et al. 2017), adaptation to high-altitude (Cho et al. 2013; Janecka et al. 2015), and conservation (Hussain 2000; Mishra et al. 2003; Jackson and Wangchuk 2004; Rosen et al. 2012; Li et al. 2014; Kachel et al. 2016) has provided many insights into snow leopard abundance, habitat use, behavior, movement patterns, and feeding ecology that are important for guiding conservation and management actions needed to ensure its persistence.
Although snow leopards prefer high-altitude mountainous habitat (e.g., Himalaya, Pamir, Alay, Kunlun, Tian Shan) (Hemmer 1972), they also occur in lower, isolated massifs (e.g., Tost and Noyon Uul in the southern Mongolia, Janecka et al. 2011a; Johansson et al. 2016) and have been observed moving through flat or rolling terrain (Schaller 1998; McCarthy et al. 2005; Johansson et al. 2016). However, limited information is available on the level of connectivity among snow leopard populations. Two recent studies modeled snow leopard habitat and connectivity primarily based on topography and climate (Riordan et al. 2015; Li et al. 2016b). Genetic analyses are needed to provide more direct information on permeability of the landscape, dispersal, and demographic fluctuations, and to identify barriers to movement (Avise 1994,, 2000).
Taxonomic classification, phylogeography, and population structure serve as the basis for conservation, management, and research (Wilson and Brown 1953; Avise 1990; Schwartz et al. 2007). Both ecological and molecular data are needed to understand the species, prioritize populations for conservation, and develop recovery or management plans (O’Brien 1991; Avise 1994; Moritz 1994; Haig et al. 2006; Rodgers and Janecka 2013). The snow leopard remains the last large felid to be the subject of a comprehensive subspecies assessment, phylogeographic analysis, and population structure study. Previous range-wide phylogeography studies of felids have primarily relied on samples from captive animals, telemetry studies, hunter-harvested individuals, or museum specimens (Culver et al. 2000; Eizirik et al. 2001; Uphyrkina et al. 2001; Luo et al. 2004). The gap in knowledge for snow leopards is a direct result of the following challenges: 1) they inhabit remote, inaccessible regions that are often politically unstable, 2) opportunities for radio/GPS telemetry are limited because they are difficult to observe and trap in the wild, and 3) most founders of the captive population have an unknown origin.
Noninvasive genetic sampling via collection of scat along wildlife trails and marking sites are an effective and efficient way to survey snow leopard populations and have become an important approach for studying snow leopards and numerous other felids (e.g., Janecka et al. 2008; Gilad et al. 2011; Rodgers and Janecka 2013; Rodgers et al. 2015; Wultsch et al. 2016a, 2016b). Our collaborative efforts to conduct noninvasive surveys of this species have yielded snow leopard DNA samples from all major parts of the range. Here, we present the results of the first, to our knowledge, range-wide snow leopard genetic study to establish a framework for understanding its taxonomy, population history, and phylogeography using noninvasively collected scat samples.
Materials and Methods
Sample Collection, DNA Extraction, and Genotyping
We collected snow leopard scat in 21 localities distributed throughout the range (Figure 1; Supplementary Table S1) using noninvasive genetic surveys following Janecka et al. (2008, 2011a). DNA was extracted with the Qiagen DNA Stool Kit (Qiagen, Inc., Valencia, CA). Snow leopard scats were identified by amplifying and aligning a 96 bp fragment of the mitochondrial cytochrome b gene with reference sequences following Janecka et al. (2008), or via a species-specific PCR assay from Janecka et al. (2014). We initially identified individuals by genotyping 8 microsatellite loci in triplicate using fluorescently labeled primers (Janecka et al. 2011a; 2014). Sex was determined by amplification of the Y-linked AMELY marker (Murphy et al. 1999) following methods in Janecka et al. (2008). Only those scat samples with no detectable errors in the initial microsatellite panel were used in the phylogeographic study and genotyped at 25 additional microsatellite loci (Supplementary Table S2).
Genetic Analysis
The majority of felid phylogeographic studies have relied on microsatellites developed for the domestic cat (Felis catus) linkage and radiation hybrid maps (Menotti-Raymond et al. 1999,, 2003). To improve PCR-based genotyping success in snow leopard scat we sequenced 49 microsatellites in 2 captive-bred snow leopard samples using Sanger sequencing, following the methods in Janecka et al. (2008). We designed new snow leopard-specific primers from repeat-masked microsatellite motif flanking sequences (i.e., no SINES, LINES, or LTRs), which were positioned closer to the repeats so the amplicon size was between 100 and 150 bp. These new primers were designated by the prefix “PUN” and the respective locus number used by Menotti-Raymond et al. (1999,, 2003). We selected a final set of 33 microsatellites based on amplification intensity, unambiguous allele peaks, and chromosomal location, and genotyped these loci in 70 individuals (primers provided in Supplementary Table S2 and PCR conditions in Supplemental Material). Three individuals were genotyped in replicate at all 33 loci, and the genotypes for 2 microsatellite loci were replicated in all 70 individuals to ensure consistency. In addition, we sequenced 3 mtDNA segments in 70 individuals including 96 bp of cytochrome b, 244 bp of the hyper variable region II, and 323 bp of the central conserved region and aligned the sequences with SEQUENCHER 5.5.5 (Gene Codes Corporation, Ann Arbor, MI). The only variable site was found in the central conserved region in one individual, therefore the mtDNA was not informative for population structure.
Genetic Diversity, Population Structure, and Coalescent Simulations
Standard estimates of genetic diversity were derived in GENALEX 6.502 (Peakall and Smouse 2006) including the number of alleles (AN), number of private alleles (AP), effective number of alleles (AE), observed heterozygosity (HO), expected heterozygosity (HE), and fixation index (FI, 1 − (HO/HE). Loci were tested for linkage disequilibrium (LD) and deviations from Hardy–Weinberg equilibrium (HWE) with significant P-values adjusted for multiple comparisons.
Recent changes in the effective population (NE) size in snow leopards were investigated using coalescent simulations implemented in MSVAR 1.3 (Beaumont 1999; Storz 2002; Girod et al. 2011) using 32 polymorphic microsatellite loci (i.e., one monomorphic locus was removed). The model assumes a single stable ancestral population NE1 in the past that experienced a demographic alteration (bottleneck or expansion) starting at time t and subsequently changed exponentially in size to the current population NE0. We simulated 2 different demographic scenarios: 1) larger prior distribution values for the contemporary population size NE0 than the ancestral NE1 (expansion) and 2) larger priors for NE1 than NE0 (bottleneck). We tested various prior distributions for each scenario to assess the independency of the posterior estimates for the parameters NE0, NE1, and t. In lieu of a snow leopard-specific microsatellite mutation rate we choose an average mammalian mutation rate (Brinkmann 1998; Rooney 1999) of 1 × 10−4 sub/site/year allowing rate variation between 10−3 and 10−5. We ran four coalescent simulations for each population with 2.5 × 109 iterations of the Markov chain Monte Carlo (MCMC) algorithm, discarding the first 20% as burn-in. Convergence of the chains from each population simulated with 4 different priors, respectively, were assessed with Gelman–Rubin’s diagnostic (Brooks 1998) implemented in the R package boa (Smith 2007). Gelman–Rubin’s convergence tests of the MCMC algorithm for the independent runs and each variable resulted in values below the threshold of 1.1 (Gelman et al. 2004).
We assessed population structure using several methods. First, we estimated the pairwise fixation index (FST, Weir and Cockerham 1984) to examine gene flow among predefined population groups (NQ, northern Qinghai, this included Aksei in northern Gansu; SQ, southern Qinghai; HIM, Tibet-Nepal-Bhutan; IP, India-Pakistan; TK, Tajikistan-Kyrgyzstan; WM, western Mongolia; SM, southern Mongolia) (Figure 2a). FST values and their significance were estimated using the analysis of molecular variance (AMOVA) framework in GENALEX. FST values were also calculated using ARLEQUIN 3.5 (Excoffier et al. 2007), but no substantial difference was found so we present the results from GENALEX. We also performed population assignment tests with frequency methods to evaluate the level of differentiation in GENALEX. These tests yield probabilities that an individual came from each population, based on its genotype and allele frequencies. If assignment probability was highest for a population in which that individual was not observed, it was considered genetically misassigned. There is a direct correlation between the misassignment rate and dispersal (Rannala and Mountain 1997; Paetkau et al. 2004; Janecka et al. 2011b). In addition, to assess genetic similarity among individual snow leopards without making assumptions regarding HWE and LD, we conducted a principal component analysis (PCA) in adegenet 1.4.2 (Jombart 2008) using R 3.2.4.
We also used individual-based Bayesian clustering approaches that explored the number of genetic clusters (K) within our samples (Pritchard et al. 2000; Guillot et al. 2005). The first approach used STRUCTURE 2.2.4 to identify genetic clusters and estimate ancestry for each individual (Pritchard et al. 2000). This method estimates the likelihoods of different numbers of genetic clusters as well as cluster membership (Q) for each individual. The analysis was done using the following model: admixture, alpha inferred from the data (initial value 1.0), correlated allele frequencies, sampling locations as prior (LOCPRIOR), and 2 × 106 MCMC iterations after a burn-in of 2 × 105 replicates. We varied the number of potential genetic clusters from 1 to 10. The most likely value of K was determined using 10 independent runs for each value of K. We analyzed our results by applying the posterior probability (Pritchard et al. 2000) and the ΔK method (Evanno et al. 2005), as implemented by pophelper (Francis 2017) in R 3.2.4. The Q for each individual was averaged across all ten STRUCTURE runs. To examine hierarchical genetic structure, we conducted additional Bayesian analysis within identified genetic clusters until no further genetic subdivision was detected, or inference was impossible due to low sample sizes (Balkenhol et al. 2014; Wultsch et al. 2016a). We also used a second Bayesian clustering approach that incorporates a spatially explicit model to generate priors as implemented in GENELAND 4.0.6 (Guillot et al. 2005). We applied the spatial model with uncorrelated allele frequencies and simulated the number of K from 1 to 10 using 10 independent runs. Each run consisted of 1 × 106 MCMC iterations with a thinning of 100. The level of spatial uncertainty was set to 50 m.
We also explored whether a smaller number of microsatellites typically used in noninvasive genetic surveys of snow leopards (i.e., 6–8 loci) would be sufficient to assign individuals to the major genetic groups identified in this study. We therefore created a reduced matrix with only six loci (PUN82, PUN100, PUN124, PUN225, PUN229, and PUN327) and included 26 additional samples from Ladakh, India. Probability of identity for unrelated (PID-unr) and related (PID-sib) individuals was estimated in GENALEX. We used the STRUCTURE test for migrants with K = 3 to assign these 26 samples to 1 of the 3 clusters. The following parameters were used: admixture, alpha inferred from the data (initial value 1.0), correlated allele frequencies, sampling locations for reference samples as spatial prior (LOCPRIOR), population information used to test for migrants, and 8 × 105 MCMC iterations after a burn-in of 4 × 105 replicates.
We examined isolation-by-distance (IBD) patterns by correlating genetic and geographic distances via the Mantel test using ecodist 1.2.9 (Goslee and Urban 2007) in R 3.2.4. We also conducted spatial autocorrelation analysis in GENALEX to assess the spatial extent of genetic structure by examining genetic similarity between pairs of snow leopard individuals at several spatial distance classes (25, 50, 100, 250, 500, 750, 1000, 1500, 2000, 2500, and 3000 km). The spatial autocorrelation coefficient (r) was evaluated at each distance class against the null hypothesis of no genetic structure (r = 0) via permutation (10000 simulations) and bootstrapping (1000 repeats).
Results
We obtained scat samples from 21 localities in 7 geographic regions (northern Qinghai [this region included Aksei in northern Gansu], southern Qinghai, Tibet-Nepal-Bhutan, India-Pakistan, Tajikistan-Kyrgyzstan, western Mongolia, and southern Mongolia) (Figures 1 and 2a; Supplementary Table S1). A total of 70 representative individuals were genotyped at 33 microsatellite loci. Measures of genetic diversity were consistently low across the entire snow leopard range (AN of 2.6–3.3; HO of 0.399–0.508; HE of 0.434–0.485; Table 1; Supplementary Table S3). Southern Mongolia and southern Qinghai had the lowest measures of diversity and the highest frequency of private alleles. When all samples were pooled, 16 loci were out of HWE, in contrast to only 3 loci within each of the 7 locations, indicative of the Wahlund effect. The coalescent simulations supported a demographic contraction scenario with the time estimate for the last bottleneck of t = 7782 years ago (ya) (range of 4574–11893; highest probability density, HPD, of 1084–107484; Supplementary Table S4 and Supplementary Figure S1) with an ancestral effective population size NE1 = 8235 (range of 7287–9852; HPD of 1951–47374) and current effective population size NE0 = 1279 (range of 1091–1504; HPD of 249–4902). We performed the analysis for all samples pooled and for the 3 main genetic clusters individually (Supplementary Figure S2). There was only a single variable site in the 683 bp concatenated mtDNA alignment, with 69 individuals having one haplotype, and the second haplotype observed in only a single individual.
Geographic region . | n . | A . | AP . | AE . | HO . | HE . | FI . |
---|---|---|---|---|---|---|---|
Snow leopard | 70 | 5.8 | n. a. | 2.8 | 0.433 | 0.568 | 0.246 |
Central group | 24 | 4.2 | 28 | 2.5 | 0.446 | 0.521 | 0.138 |
Northern Qinghai | 6 | 3.1 | 10 | 2.3 | 0.457 | 0.484 | 0.053 |
Southern Qinghai | 5 | 2.6 | 1 | 2.1 | 0.450 | 0.434 | −0.019 |
Tibet/Nepal/Bhutan | 13 | 3.3 | 12 | 2.2 | 0.441 | 0.452 | 0.012 |
Western group | 16 | 4.0 | 19 | 2.6 | 0.461 | 0.522 | 0.126 |
India/Pakistan | 8 | 3.2 | 10 | 2.3 | 0.508 | 0.473 | −0.082 |
Tajikistan/Kyrgyzstan | 8 | 3.3 | 6 | 2.4 | 0.415 | 0.485 | 0.174 |
Northern group | 30 | 3.9 | 22 | 2.3 | 0.408 | 0.481 | 0.152 |
Western Mongolia | 15 | 3.1 | 5 | 2.2 | 0.416 | 0.442 | 0.058 |
Southern Mongolia | 15 | 3.3 | 14 | 2.1 | 0.399 | 0.450 | 0.110 |
Geographic region . | n . | A . | AP . | AE . | HO . | HE . | FI . |
---|---|---|---|---|---|---|---|
Snow leopard | 70 | 5.8 | n. a. | 2.8 | 0.433 | 0.568 | 0.246 |
Central group | 24 | 4.2 | 28 | 2.5 | 0.446 | 0.521 | 0.138 |
Northern Qinghai | 6 | 3.1 | 10 | 2.3 | 0.457 | 0.484 | 0.053 |
Southern Qinghai | 5 | 2.6 | 1 | 2.1 | 0.450 | 0.434 | −0.019 |
Tibet/Nepal/Bhutan | 13 | 3.3 | 12 | 2.2 | 0.441 | 0.452 | 0.012 |
Western group | 16 | 4.0 | 19 | 2.6 | 0.461 | 0.522 | 0.126 |
India/Pakistan | 8 | 3.2 | 10 | 2.3 | 0.508 | 0.473 | −0.082 |
Tajikistan/Kyrgyzstan | 8 | 3.3 | 6 | 2.4 | 0.415 | 0.485 | 0.174 |
Northern group | 30 | 3.9 | 22 | 2.3 | 0.408 | 0.481 | 0.152 |
Western Mongolia | 15 | 3.1 | 5 | 2.2 | 0.416 | 0.442 | 0.058 |
Southern Mongolia | 15 | 3.3 | 14 | 2.1 | 0.399 | 0.450 | 0.110 |
n, samples size; AN, number of alleles; AP, private alleles; AE, effective number of alleles (1/(∑p2); HO, observed heterozygosity; HE, expected heterozygosity; FI, fixation index (1 − HO/HE).
Geographic region . | n . | A . | AP . | AE . | HO . | HE . | FI . |
---|---|---|---|---|---|---|---|
Snow leopard | 70 | 5.8 | n. a. | 2.8 | 0.433 | 0.568 | 0.246 |
Central group | 24 | 4.2 | 28 | 2.5 | 0.446 | 0.521 | 0.138 |
Northern Qinghai | 6 | 3.1 | 10 | 2.3 | 0.457 | 0.484 | 0.053 |
Southern Qinghai | 5 | 2.6 | 1 | 2.1 | 0.450 | 0.434 | −0.019 |
Tibet/Nepal/Bhutan | 13 | 3.3 | 12 | 2.2 | 0.441 | 0.452 | 0.012 |
Western group | 16 | 4.0 | 19 | 2.6 | 0.461 | 0.522 | 0.126 |
India/Pakistan | 8 | 3.2 | 10 | 2.3 | 0.508 | 0.473 | −0.082 |
Tajikistan/Kyrgyzstan | 8 | 3.3 | 6 | 2.4 | 0.415 | 0.485 | 0.174 |
Northern group | 30 | 3.9 | 22 | 2.3 | 0.408 | 0.481 | 0.152 |
Western Mongolia | 15 | 3.1 | 5 | 2.2 | 0.416 | 0.442 | 0.058 |
Southern Mongolia | 15 | 3.3 | 14 | 2.1 | 0.399 | 0.450 | 0.110 |
Geographic region . | n . | A . | AP . | AE . | HO . | HE . | FI . |
---|---|---|---|---|---|---|---|
Snow leopard | 70 | 5.8 | n. a. | 2.8 | 0.433 | 0.568 | 0.246 |
Central group | 24 | 4.2 | 28 | 2.5 | 0.446 | 0.521 | 0.138 |
Northern Qinghai | 6 | 3.1 | 10 | 2.3 | 0.457 | 0.484 | 0.053 |
Southern Qinghai | 5 | 2.6 | 1 | 2.1 | 0.450 | 0.434 | −0.019 |
Tibet/Nepal/Bhutan | 13 | 3.3 | 12 | 2.2 | 0.441 | 0.452 | 0.012 |
Western group | 16 | 4.0 | 19 | 2.6 | 0.461 | 0.522 | 0.126 |
India/Pakistan | 8 | 3.2 | 10 | 2.3 | 0.508 | 0.473 | −0.082 |
Tajikistan/Kyrgyzstan | 8 | 3.3 | 6 | 2.4 | 0.415 | 0.485 | 0.174 |
Northern group | 30 | 3.9 | 22 | 2.3 | 0.408 | 0.481 | 0.152 |
Western Mongolia | 15 | 3.1 | 5 | 2.2 | 0.416 | 0.442 | 0.058 |
Southern Mongolia | 15 | 3.3 | 14 | 2.1 | 0.399 | 0.450 | 0.110 |
n, samples size; AN, number of alleles; AP, private alleles; AE, effective number of alleles (1/(∑p2); HO, observed heterozygosity; HE, expected heterozygosity; FI, fixation index (1 − HO/HE).
Among these 7 geographic regions, the greatest genetic similarity based on the pairwise FST was between northern Qinghai and southern Qinghai (FST = 0.039, P = 0.174), between India-Pakistan and Tajikistan-Kyrgyzstan (FST = 0.007, P = 0.412), and between western Mongolia and southern Mongolia (FST = 0.057, P < 0.000; Table 2). The most divergent FST values (>0.25) indicated high levels of differentiation for southern Mongolia versus southern Qinghai (FST = 0.308, P < 0.000), western Mongolia versus southern Qinghai (FST = 0.287, P < 0.000), and Tibet-Bhutan-Nepal versus southern Mongolia (FST = 0.258, P < 0.000). The 2 Mongolian regions were the most differentiated with respect to the other regions.
. | Northern Qinghai . | Southern Qinghai . | Tibet/Nepal/Bhutan . | India/Pakistan . | Tajikistan/Kyrgyzstan . | Western Mongolia . | Southern Mongolia . |
---|---|---|---|---|---|---|---|
Northern Qinghai | — | 0.174 | 0.001 | 0.000 | 0.001 | 0.000 | 0.000 |
Southern Qinghai | 0.039 | — | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 |
Tibet/Nepal/Bhutan | 0.121 | 0.177 | — | 0.000 | 0.000 | 0.000 | 0.000 |
India/Pakistan | 0.138 | 0.211 | 0.166 | — | 0.412 | 0.000 | 0.000 |
Tajikistan/Kyrgyzstan | 0.133 | 0.204 | 0.143 | 0.007 | — | 0.000 | 0.000 |
Western Mongolia | 0.218 | 0.287 | 0.227 | 0.142 | 0.092 | — | 0.000 |
Southern Mongolia | 0.220 | 0.308 | 0.258 | 0.161 | 0.122 | 0.057 | — |
. | Northern Qinghai . | Southern Qinghai . | Tibet/Nepal/Bhutan . | India/Pakistan . | Tajikistan/Kyrgyzstan . | Western Mongolia . | Southern Mongolia . |
---|---|---|---|---|---|---|---|
Northern Qinghai | — | 0.174 | 0.001 | 0.000 | 0.001 | 0.000 | 0.000 |
Southern Qinghai | 0.039 | — | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 |
Tibet/Nepal/Bhutan | 0.121 | 0.177 | — | 0.000 | 0.000 | 0.000 | 0.000 |
India/Pakistan | 0.138 | 0.211 | 0.166 | — | 0.412 | 0.000 | 0.000 |
Tajikistan/Kyrgyzstan | 0.133 | 0.204 | 0.143 | 0.007 | — | 0.000 | 0.000 |
Western Mongolia | 0.218 | 0.287 | 0.227 | 0.142 | 0.092 | — | 0.000 |
Southern Mongolia | 0.220 | 0.308 | 0.258 | 0.161 | 0.122 | 0.057 | — |
FST (below diagonal) and associated P-values (above diagonal) were based on 33 microsatellites.
. | Northern Qinghai . | Southern Qinghai . | Tibet/Nepal/Bhutan . | India/Pakistan . | Tajikistan/Kyrgyzstan . | Western Mongolia . | Southern Mongolia . |
---|---|---|---|---|---|---|---|
Northern Qinghai | — | 0.174 | 0.001 | 0.000 | 0.001 | 0.000 | 0.000 |
Southern Qinghai | 0.039 | — | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 |
Tibet/Nepal/Bhutan | 0.121 | 0.177 | — | 0.000 | 0.000 | 0.000 | 0.000 |
India/Pakistan | 0.138 | 0.211 | 0.166 | — | 0.412 | 0.000 | 0.000 |
Tajikistan/Kyrgyzstan | 0.133 | 0.204 | 0.143 | 0.007 | — | 0.000 | 0.000 |
Western Mongolia | 0.218 | 0.287 | 0.227 | 0.142 | 0.092 | — | 0.000 |
Southern Mongolia | 0.220 | 0.308 | 0.258 | 0.161 | 0.122 | 0.057 | — |
. | Northern Qinghai . | Southern Qinghai . | Tibet/Nepal/Bhutan . | India/Pakistan . | Tajikistan/Kyrgyzstan . | Western Mongolia . | Southern Mongolia . |
---|---|---|---|---|---|---|---|
Northern Qinghai | — | 0.174 | 0.001 | 0.000 | 0.001 | 0.000 | 0.000 |
Southern Qinghai | 0.039 | — | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 |
Tibet/Nepal/Bhutan | 0.121 | 0.177 | — | 0.000 | 0.000 | 0.000 | 0.000 |
India/Pakistan | 0.138 | 0.211 | 0.166 | — | 0.412 | 0.000 | 0.000 |
Tajikistan/Kyrgyzstan | 0.133 | 0.204 | 0.143 | 0.007 | — | 0.000 | 0.000 |
Western Mongolia | 0.218 | 0.287 | 0.227 | 0.142 | 0.092 | — | 0.000 |
Southern Mongolia | 0.220 | 0.308 | 0.258 | 0.161 | 0.122 | 0.057 | — |
FST (below diagonal) and associated P-values (above diagonal) were based on 33 microsatellites.
The PCA revealed 3 major groups consistent with the geographic distribution of sampled localities (Figure 2b). Specifically, the snow leopards from the Tibetan Plateau (northern Qinghai, southern Qinghai, and Tibet) and the principal portion of the Himalaya (Bhutan and Nepal) clustered together into a “Central group,” the snow leopards from Western Himalaya (India), Karakorum, Pamir, Alay, and Tian Shan (Tajikistan, Kyrgyzstan) formed a “Western group,” and those from Altai (western Mongolia) and Southern Gobi (southern Mongolia) formed a “Northern group” (Figure 2b). Population assignment tests also supported strong genetic differentiation of the 3 groups. There was an 8.5% misassignment rate when samples were divided into 7 populations. When samples were divided into the Northern, Central, and Western populations there was no misassignment (Supplementary Tables S5 and S6).
Bayesian clustering in STRUCTURE tested for the presence of K of 1–10 genetic clusters using the admixture model with correlated allele frequencies and sampling locations as priors (Supplementary Figures S3 and S4). The Evanno et al. (2005) ad hoc method supported division into 2 groups (K = 2) with the first cluster comprised exclusively of samples from Mongolia (i.e., the Northern group) and the second cluster composed of the individuals from all remaining sites (Figure 2c). The majority of individuals were assigned (Q > 90%) to 1 of these 2 clusters. Snow leopards in Tajikistan and Kyrgyzstan showed evidence of genetic admixture with Mongolia. At K = 3, the Central, Western, and Northern groups supported by PCA analysis and the assignment tests were recovered (Figure 2d). The spatially explicit Bayesian model implemented in GENELAND detected the same 3 clusters, although with more admixture (Figure 3). When a reduced matrix was created with only 6 microsatellite loci (PID-unr = 0.000033, PID-sib = 0.0097), a number commonly used for individual identification in noninvasive surveys, the STRUCTURE test to detect migrants assigned 96% of 26 additional samples from Ladakh (India) to the appropriate genetic cluster (i.e., the Western group) in which they were sampled (Supplementary Table S7).
Bayesian clustering analysis in STRUCTURE was performed at 2 additional hierarchical levels to further explore differentiation within the Central, Western, and Northern groups. The first level analyzed the Central and Western groups together (samples from northern Qinghai, southern Qinghai, Tibet, Bhutan, Nepal, India, Pakistan, Tajikistan, and Kyrgyzstan; Figure 2e) and in a separate analysis the Western and Northern groups were analyzed together (samples from India, Pakistan, Tajikistan, Kyrgyzstan, western Mongolia, and southern Mongolia; Figure 2f). Both analyses resulted in K = 2, with the clusters corresponding to those identified with the full data set for K = 3, and the PCA and population assignment analyses. There was no evidence for genetic admixture between Central and Western groups at this spatial scale (Figure 2e). In contrast, when Western and Northern groups were analyzed together, although the snow leopards sampled from Tajikistan and Kyrgyzstan clustered with India and Pakistan, they did exhibit low levels of admixture with Mongolia (Figure 2f). The second level of hierarchy included separate analyses of the Central (northern Qinghai, southern Qinghai, Tibet, and Himalaya; Figure 2g), Western (India, Pakistan, Tajikistan, and Kyrgyzstan; Figure 2h), and Northern (western Mongolia and southern Mongolia; Figure 2i) groups. Together, these analyses revealed additional structure within these regions. In the Central group, northern and southern Qinghai clustered separately from Tibet, Bhutan, and Nepal with admixture in the central part of the Tibetan Plateau (Figure 2g). The Western group had 2 clusters, one composed exclusively of Kyrgyzstan and the other of India, Pakistan, and Tajikistan (Figure 2h). Finally, in the Northern group, western and southern Mongolia formed 2 distinct clusters, with a few individuals showing mixed ancestry (Figure 2i). The genetic clusters observed in the hierarchical analysis largely corresponded to those recovered with the full dataset for K = 6 (Supplementary Figure S3), which had the highest posterior probability in the STRUCTURE analysis with 70 samples (Supplementary Figure S4a).
The simple Mantel test showed weak IBD (r = 0.183, P = 0.020; Supplementary Figure S5). Spatial autocorrelation analysis detected IBD in the first four distance classes (25, 50, 100, and 250 km; Figure 4). The x-intercept of r was between 250 and 500 km, indicating that snow leopard populations located within these distances are potentially connected by dispersing individuals.
Discussion
The snow leopard, first described by Schreber (1775), is a high-altitude specialist that occupies mountains primarily between 1500 and 4500 m, with confirmed sightings to 6000 m in the Himalaya (Hemmer 1972). Historically, several subspecies have been proposed including the nominate Panthera uncia uncia (Stroganov 1962), Panthera uncia uncioides in Nepal (Horsfield 1855), Panthera uncia schneideri in Sikkim (India, Zukowsky 1950), and Panthera uncia baikalensis-romanii in the Transbaikal (Russia, Medvedev 2000) based on coat color differences. The lack of collection information for many museum specimens and the difficulty of observing and trapping snow leopards in the wild have to date prevented comprehensive taxonomic assessments and therefore this species is considered monotypic. Here, we address the taxonomic question regarding subspecies designation of the snow leopard using genetic data from noninvasive scat samples. Subspecies are generally considered distinct populations that correspond to geographic boundaries and meet discreteness and significance criteria (Wilson and Brown 1953; Haig et al. 2006). Based on the differentiation of nuclear loci that separated samples into 3 discrete and significant genetic clusters (Western, Central, and Northern groups) occurring in nonoverlapping geographic regions we propose the snow leopard be classified into 3 subspecies; P. u. uncia (type locality restricted to Central Asia including Tian Shan, Alay, Pamir, Karakoram, and trans-Himalaya), P. u. uncioides with schneideri as a junior synonym (core Himalaya and the Tibetan Plateau), and P. u. irbis with baikalensis-romanii as a junior synonym (the Altai and Southern Gobi of Mongolia) (Figures 1a and 2b).
Although historically recommendations for subspecies delineations have also included mitochondrial divergence and monophyly (Moritz 1994), phylogeographic studies over the last decade indicate that this strict criterion may be unreasonable because of mitonuclear discordance (Toews and Brelsford 2012). First, mitogenomes are frequently paraphyletic, notably in Felidae, due to introgression from past hybridization events (Roca et al. 2005; Li et al. 2016a). Second, the mtDNA represents a single linked locus (i.e., no recombination) with a smaller NE (i.e., haploid and only passed through females) therefore it is more sensitive to incomplete lineage sorting than nuclear loci (Avise 2000). Although our mtDNA data revealed no haplotype differences across the snow leopard range, the consistent recovery using nuclear markers of 3 discrete allopatric genetic clusters with significant differentiation, each occurring in unique geographic regions, warrants subspecies delineation.
There are several possible reasons why we did not observe different mtDNA haplotypes in the 3 subspecies. First, more extensive mitogenome sequencing may be required to detect polymorphism. Second, the snow leopard mitogenome may have undergone a selective sweep. Mitochondria are responsible for oxidative respiration and therefore may be under selective pressures in hypoxic environments (da Fonseca 2008; Hassanin et al. 2009). In addition, previous studies have shown that the snow leopard lineage underwent mitochondrial replacement after hybridization with the African lion lineage (Li et al. 2016a), and therefore may have been subject to adaptive introgression that resulted in low mtDNA variation (Toews and Brelsford 2012). Finally, mtDNA has a 4-fold smaller NE compared to nuclear DNA (Hudson and Turelli 2003), and therefore more ancestral polymorphism would have been lost during the bottleneck ∼8000 ya detected with microsatellites. The lack of distinct mtDNA lineages is consistent with previous studies in the Tibetan region showing weaker Pleistocene refugia effects (Qu and Lei 2009; Yang et al. 2009; Zhan et al. 2011) compared to those observed in European and North American taxa (Taberlet et al. 1998; Petit et al. 2003; Shafer et al. 2010). Sequencing of mitochondrial and nuclear genomes in the 3 subspecies would shed light on the mechanisms that contributed to the observed mitonuclear discordance.
Snow leopards face threats including low prey densities, retaliatory killing by farmers and herdsmen in response to livestock depredation, illegal wildlife trade, climate change, and development of roads, rails, mining, and hydropower facilities (Jackson et al. 2010; McCarthy et al. 2016). Traditional pastoralism and agro-pastoralism dominate local economies within the snow leopard range often leading to human-wildlife conflict (Mishra et al. 2003,, 2016). Successful community-based conservation initiatives have been implemented in several areas, including Mongolia, Nepal, and Pakistan (Jackson et al. 2010). Recently, there has been an effort to coordinate conservation range-wide with a comprehensive Global Snow Leopard & Ecosystem Protection Program that seeks to secure populations in 20 different landscapes by 2020 (Snow Leopard Secretariat 2013). Our phylogeographic assessment strongly highlights the importance of large-scale international efforts (Rosen and Zahler 2016). Snow leopard populations exhibit cross-boundary connectivity in several important parts of their range, such as between Pakistan and India on the western portion of the Himalaya, and between Nepal, Bhutan, and southern Tibet. It is critical that international corridors between these populations are maintained, and that synchronized conservation actions are realized so that no single area becomes isolated, or a population sink, contributing to decline in neighboring countries.
Presently, the International Union for Conservation of Nature (IUCN) considers the snow leopard a monotypic species and applies criteria for “Endangered (EN) Category 1 (C1)” status range-wide. These include <2500 mature individuals and an estimated 20% decline in 2 generations, corresponding to ∼16 years in snow leopards (IUCN Standards and Petitions Subcommittee 2016). Applying the results from our phylogeographic analysis, we generated a preliminary population size estimate for each subspecies using population size estimates in McCarthy et al (2016b) by summing those within the approximate range of each respective subspecies (Supplementary Table S8 and Supplementary Figure S6). The estimate for P. u. uncia was 2124–3356 individuals, for P. u. uncioides it was 1402–3083, and for P. u. irbis it was 741–1646. This suggests the latter 2 subspecies may meet IUCN EN C1 criteria. These population estimates are preliminary and additional research is needed to both definitively assign populations to subspecies, and obtain abundance information in areas where quantitative data is not available, particularly on the Tibetan Plateau.
In order to determine if assignment of individuals to subspecies could be made with fewer microsatellites we created a reduced matrix of 6 loci and tested 26 additional samples from Ladakh, India. Even with this reduced matrix, we were able to correctly assign 96% of the individuals to P. u uncia, thus illustrating the utility of this reference data set and the substantial level of differentiation between the subspecies. We have performed whole genome amplification of a subset of the samples yielding synthetically derived amplicons not subject to CITES restrictions (Janecka et al. 2006,, 2007). We will make these available to other laboratories to ensure uniform allele designations, thus facilitating direct comparisons with our data set in future studies.
The low microsatellite diversity and lack of mtDNA variation within snow leopards is typical of felids or subspecies that have either been historically isolated, or have undergone recent population bottlenecks, such as the Far Eastern leopard (Panthera pardus orientalis, Uphyrkina et al. 2001), Sumatran tiger (Panthera tigris sumatrae, Luo et al. 2004), North American puma (Puma concolor cougar, Culver et al. 2000), North American ocelot (Leopardus pardalis albescens, Janecka et al. 2011b), and the Asiatic cheetah (Acinonyx jubatus venaticus, Charruau et al. 2011). Although we sequenced a limited amount of mtDNA, Luo et al. (2004) detected 4–6 haplotypes for the same segments in tigers despite analyzing fewer individuals. The microsatellite and mtDNA diversity in snow leopards is consistent with the low genome-wide polymorphism previously reported for a single snow leopard from Mongolia (Cho et al. 2013). We estimated the most recent bottleneck occurred during the middle Holocene, potentially as a consequence of changing climatic and habitat conditions (Zhang et al. 2006; Yang et al. 2009). This finding is reinforced by the uniformly reduced variation across the snow leopard range. Alternatively, if anthropogenic factors were the primary cause of lower variation in extant snow leopards, the diversity would likely vary across populations.
The estimated time of the bottleneck in snow leopards coincides with the start of the Holocene Thermal Maximum (approximately 6000–8000 years ago), a period of warming and increased precipitation in the Tibetan Plateau and a synchronous treeline shift to higher elevations (Zhao et al. 2011). The correspondence of these events with the snow leopard bottleneck has important implications for understanding the potential impact of global climate change. Similar climatological trends are occurring throughout the world, with particularly elevated warming trends in the Tibetan Plateau and the Himalaya (Liu and Chen 2000; Walther et al. 2002; Farrington and Li 2016). Recent studies examining the potential impact of climate change have predicted a substantial reduction in snow leopard habitat and increased fragmentation (Forrest et al. 2012; Li et al. 2016b). Our inferences on the demographic contraction in the Holocene lend support to models that indicate snow leopards are susceptible to global warming.
Major landscape features in Asia correlate with the observed phylogeographic patterns. In the north, the P. u. irbis populations that occupy low-altitude mountains of the Gobi in southern Mongolia are separated from the Qilian Shan in northern Qinghai by the Alashan Plateau, with >400 km of unsuitable habitat. This potential movement barrier corresponds to the greatest observed genetic differentiation within the species, and is consistent with the recent habitat and connectivity models (Riordan et al. 2015; Li et al. 2016b). The admixture observed in Kyrgyzstan, on the other side of China, indicates the presence of more recent introgression between the subspecies or unsampled populations with intermediate allele frequencies. The Dzungarian Basin (∼500 km wide) is likely an impenetrable barrier for snow leopards. However, there are isolated mountains to the west along the boundary of Kazakhstan and China, and to the east in Xinjiang between Tian Shan and the Gobi in southern Mongolia that may act as stepping stones between P. u. uncia and P. u. irbis (Li et al. 2016b). However, both of these possible routes appear too far north and west to connect Qinghai/Gansu (P. u. uncioides) with the northern subspecies. The barrier separating the Central P. u. uncioides and the Western P. u. uncia, between Nepal and India (Ladakh), is not as obvious because the trans-Himalaya form a nearly continuous chain. One possibility is that the combination of their height (i.e., many peaks >6000 m above sea level), linear distance (∼1000 km), and the presence of several major rivers in the region may limit connectivity. Another possibility is that in the past this area may have been covered by extensive uninhabitable glacial fields. These explanations are mutually compatible and both may contribute to the observed differentiation in snow leopards between the Central and Western groups, which is consistent with the Li et al. (2016b) model of snow leopard habitat that indicates more fragmentation in this area than predicted by the Riordan et al. (2015) model.
The observed population genetic structure also provides a coherent and objective basis with which to define management units (MUs). In the north, western Mongolia (Altai MU) and southern Mongolia (Gobi MU) formed separate genetic clusters in the hierarchical analysis. In this area, snow leopards primarily occupy the Altai mountain range, which runs east into the Gobi Altai, with the mountains becoming lower and more isolated, often separated by >100 km of flat desert. Nonetheless, there was evidence for transient dispersal between each of these areas suggesting that the smaller massifs act as stepping stones for migration between larger habitat patches. Long distance movements have been observed among radio and GPS-collared snow leopards in Mongolia (McCarthy et al. 2005; Johansson et al. 2016). In the southern portion of Tibetan Plateau there was connectivity with the Himalaya, consistent with snow leopard habitat models (Riordan et al. 2015; Li et al. 2016b). Northern Qinghai/Gansu and southern Qinghai (Qinghai MU) were genetically divergent from Tibet, Bhutan, and Nepal (Tibet-Himalaya MU) similar to phylogeographic patterns observed in other species (Qu and Lei 2009; Yang et al. 2009; Qu et al. 2010; Zhan et al. 2011). Within Central Asia, the Tajikistan, Pakistan, and India (Pamir-Himalaya MU) region also appeared connected. The Pamirs in Tajikistan are separated from the Tian Shan in Kyrgyzstan (Tian Shan MU) by ∼600 km, some of which includes river valleys (e.g., Vakhsh, Kyzyl-Suu, and Naryn Rivers) which may explain the differentiation in this area. Our subspecies designations, and the delineation of the Qinghai MU and Tibet-Himalaya MU, also correspond with 4 physiography and prey types zones described from a recent meta-analysis of feeding ecology studies (Lyngdoh et al. 2014). Snow leopard conservation efforts need to focus on maintaining natural connectivity within these MUs and to develop context-specific conservation programs. These efforts should take priority over attempts to establish corridors that would cross natural phylogeographic boundaries. Additional sampling is urgently needed to better understand population structure within the MUs and landscape factors that affect connectivity.
Conclusions
We conducted the first range-wide genetic analysis of wild snow leopard populations and delineate 3 subspecies. The criteria for IUCN’s Red List should be applied to each of these individually. The snow leopard underwent a bottleneck ∼8000 ya in the Holocene coinciding with global warming that occurred during the epoch. Population structure within the subspecies indicates a minimum of 6 MUs, 3 of which span multiple countries. Each of these may require different conservation initiatives. Our results serve as a foundation for understanding landscape connectivity of snow leopards and set the stage for more in-depth genomic studies.
Supplementary Material
Supplementary data are available at Journal of Heredity online.
Funding
The work was supported by the National Geographic Society (grant 8369-07 to J.E.J., W.J.M., L.Q., and Z.Y.), the Snow Leopard Conservancy (grant G1400042 to J.E.J., R.J., B.M., and W.J.M.), Britten Foundation (grants 12121 and 11104 to R.J.), Larry Bowman Foundation (grants 16092 and 14031 to R.J.), The College of Veterinary Medicine of Texas A&M University to J.E.J. and W.J.M., The Snow Leopard Conservation Grants Program (grants 120808 and G1500042 to J.E.J.), B.M., N.G., and M.J.), Duquesne University to J.E.J., NSF (grant EF0629849 to W.J.M.), the University of Montana to T.R.W., Council of Scientific and Industrial Research, India (grant BSC0207 to A.G. and S.K.), and Department of Biotechnology, India (grant GAP0374 to A.G. and S.K.)
Data Availability
Microsatellite data was archived in DRYAD (doi:10.5061/dryad.sb20c) and sequences for microsatellite flanking regions and mtDNA segments were deposited in GenBank under accessions KY967522–KY967572.
Acknowledgments
We would like to thank the rangers, biologists, herders, and local communities that helped in collection of samples. All samples were obtained in accordance with the appropriate agencies and all international shipments were in compliance with CITES requirements.