Skip to main content

Genomic variation in the American pika: signatures of geographic isolation and implications for conservation



Distributional responses by alpine taxa to repeated, glacial-interglacial cycles throughout the last two million years have significantly influenced the spatial genetic structure of populations. These effects have been exacerbated for the American pika (Ochotona princeps), a small alpine lagomorph constrained by thermal sensitivity and a limited dispersal capacity. As a species of conservation concern, long-term lack of gene flow has important consequences for landscape genetic structure and levels of diversity within populations. Here, we use reduced representation sequencing (ddRADseq) to provide a genome-wide perspective on patterns of genetic variation across pika populations representing distinct subspecies. To investigate how landscape and environmental features shape genetic variation, we collected genetic samples from distinct geographic regions as well as across finer spatial scales in two geographically proximate mountain ranges of eastern Nevada.


Our genome-wide analyses corroborate range-wide, mitochondrial subspecific designations and reveal pronounced fine-scale population structure between the Ruby Mountains and East Humboldt Range of eastern Nevada. Populations in Nevada were characterized by low genetic diversity (π = 0.0006–0.0009; θW = 0.0005–0.0007) relative to populations in California (π = 0.0014–0.0019; θW = 0.0011–0.0017) and the Rocky Mountains (π = 0.0025–0.0027; θW = 0.0021–0.0024), indicating substantial genetic drift in these isolated populations. Tajima’s D was positive for all sites (D = 0.240–0.811), consistent with recent contraction in population sizes range-wide.


Substantial influences of geography, elevation and climate variables on genetic differentiation were also detected and may interact with the regional effects of anthropogenic climate change to force the loss of unique genetic lineages through continued population extirpations in the Great Basin and Sierra Nevada.


The impact of global climate change on species with wide geographic distributions is of particular interest given the greater likelihood that refugia from climate-mediated extirpations may exist for these taxa [1]. For widely dispersed alpine mammals that are already constrained by warming ambient air temperatures, the habitat characteristics that correlate with population persistence across mountain ranges will help define those areas where refugia may be found [2,3,4,5]. Although the complex topography and temporal variability of mountain ecosystems may offer refugia across relatively small spatial scales (kilometers), anthropogenic climate change threatens to test the limits of such microhabitat buffering [3, 6,7,8,9]. Therefore, in addition to habitat predictors of persistence, connectivity among habitat patches via gene flow may be critical to maintain genetic diversity and evolutionary potential within species [10,11,12]. Here we examine extant genetic variation within and among American pika (Ochotona princeps) populations sampled from diverse habitats across the western United States. We consider variation across geographic scales shaped by both historical (Pleistocene-era) and contemporary levels of connectivity.

The American pika has become a canary-in-the-coal-mine for the effects of anthropogenic climate change in montane habitats, as a result of extensive extirpations from warmer, lower elevation sites over the past two decades [13,14,15,16,17]. This small lagomorph is found in the mountainous regions of western North America where it lives primarily on talus slopes above timberline [18]. Talus (rocky slopes formed by freeze–thaw processes) provides thermal refuge for the pika, which has a relatively narrow thermal tolerance, and thus relies heavily on the stable thermal microclimate provided by the interstitial spaces of this rocky habitat [19, 20]. The American pika is one of only 30 Ochotonidae species worldwide. Most species are found in Asia and Eastern Europe with one North American congener in Alaska and Canada, O. collaris [18, 21].

Based upon analyses of Sanger sequenced mitochondrial (mtDNA) and nuclear DNA loci [22, 23], as well as allozyme [24], morphological [23] and behavioral data [23, 25, 26], Hafner and Smith [21] suggested collapsing the previously recognized 36 subspecies of O. princeps [18] into five subspecies. These designations are congruent with mitochondrial lineage designations [22] and correspond to distinct mountain ranges and provinces across western North America (i.e., Cascade and coastal ranges; Sierra Nevada and Great Basin; central Utah; Northern Rocky Mountains; and Southern Rocky Mountains). The origin of these lineages dates to 1.3 mya for the basal split between Cascade and Sierra Nevada lineages, and 0.8 mya for the subsequent divergence of Central Utah, Northern Rocky Mountain, and Southern Rocky Mountain lineages [23]. Repeated warming and cooling periods throughout the Pleistocene drove alternating range expansions and contractions. As the climate began warming at the end of the Pleistocene, O. princeps retreated upslope to higher elevations in the Sierra Nevada, Cascade and Rocky Mountains as well as in numerous smaller mountain ranges in the Great Basin physiographic province [23, 27, 28]. Over the ensuing 8000–10000 years, pika populations were lost from many mountain ranges in the Great Basin, particularly from low elevation ranges that lack sufficient talus habitat to provide refuge in a warming climate [29, 30]. Climate change in the twentieth and twenty-first centuries has further accelerated population declines and losses primarily from the Sierra Nevada lineage, O. p. schisticeps [17, 31,32,33]. This is especially true for the mountain ranges of the Great Basin, where resurveys of sites occupied in the mid-twentieth century revealed extensive extirpations by the early to mid-2000s [13, 15, 31].

In the few range-wide genetic studies that have been conducted for pikas, data suggest low levels of diversity within lineages [23, 24, 28]. However, within lineage population genetic data come from a few and mostly geographically distant populations characterized using a small number of traditional molecular markers (allozymes, mtDNA and nuclear introns and microsatellite loci) [22,23,24, 34]. As a result, detailed information on the amount and spatial structuring of variation within lineages is limited. A few studies have examined gene flow among populations separated by shorter distances of 0.5–10 km in both high elevation continuous or semi-continuous habitats and in marginal habitat or at the range periphery [35,36,37,38]. While these studies reveal significant population genetic structure at small spatial scales, patterns are not always consistent with an isolation-by-distance model, thereby strongly suggesting an interaction between habitat features and gene flow [39]. Because of naturally fragmented habitat and low dispersal capabilities, local extinction and colonization suggestive of a metapopulation dynamic has been proposed for this species at multiple spatial and temporal scales [28, 35, 40, 41].

Several recent studies have used more thorough genomic sampling to evaluate the influence of environmental variation and gene flow on fine scale patterns of genetic variation across populations at small geographic scales. Russello et al. [42] generated genotyping-by-sequencing (GBS) data to identify genomic regions responding to selection across four pika populations found along an elevational gradient within the coastal pika lineage in British Columbia. Genome-wide estimates of genetic diversity showed that the warmer, lower elevation populations had significant heterozygote deficiency suggestive of inbreeding and representing potential sink habitat [42]. In an additional study based on similar data, Waterhouse et al. [43] revealed directional gene flow from high to low elevation sites, which suggests that low elevation pika populations and any unique genetic variation they harbor may be lost in local extirpation events. In general, however, we largely lack an understanding of how genetic structure and diversity in pikas varies across ecoregions and habitats. Additional high throughput sequencing data for pika populations sampled across a continuum of spatial scales is likely to improve our understanding of patterns of population genetic structure, standing genetic variation, and their environmental correlates.

The populations sampled for this study represent three of the five O. princeps mitochondrial lineages as well as sampling locations that span the diversity of talus habitat found across the species’ range (Table 1; Fig. 1). We sampled two populations from the Sierra Nevada lineage, which comprises populations in the large Sierra Nevada range as well as the majority of smaller Great Basin mountain ranges. These populations are found in contrasting habitats separated by 38 km: one in high elevation continuous or semi-continuous talus habitat in the Sierra Nevada (3170 m; “Pipet Tarn”), and the other in lower elevation, highly fragmented habitat in the Bodie Hills (2500 m; “Bodie”), a lower elevation spur to the main Sierra Nevada range. Both of these populations have been the subject of earlier genetic studies assessing movement dynamics and gene flow under widely differing habitat spatial structures [35, 38, 44]. The lower elevation site has been the focus of extensive research on metapopulation dynamics in highly fragmented habitat [35, 40, 41, 45].

Table 1 Sample information for 11 pika populations used in this study which represent three subspecies (O. p. princeps, O. p. schisticeps, and O. p. saxitilis) and three of the five mtDNA lineages (Sierra Nevada, Northern and Southern Rocky Mountains)
Fig. 1

Photographs provide by authors KBK (Bodie and Pipet Tarn) and CR (West Knoll)

Photographs for three of the study areas from California and Colorado showing habitat diversity (Bodie Hills, Pipet Tarn, West Knoll).

We also include populations from the Ruby Mountains and the East Humboldt Range in eastern Nevada, which are geographically adjacent and separated by a low elevation pass (~ 2800 m; Secret Pass; Fig. 2). Despite being located in the Great Basin, previous genetic work has shown that these populations represent the western most peripheral extent of the Northern Rocky Mountain mtDNA lineage [22]. The Great Basin contains several hundred, relatively small mountain ranges separated by wide low elevation desert basins, which serve as effective barriers to gene flow for small mammals [29, 46]. These mountain ranges are also strikingly linear and narrow, with a much smaller spatial extent than the Rocky Mountains and Sierra Nevada [29]. However, the Ruby Mountains and the East Humboldt Range constitute one of the largest contiguous montane habitats in the Great Basin (141 km2 above 2280 m) [47], with extensive talus that supports one of the largest pika populations remaining in this physiographic province. The remaining populations sampled for this study are from the Northern and Southern Rocky Mountains and occupy large contiguous talus slopes where pika ecology has been under study for decades [48,49,50].

Fig. 2

Figure created by author KBK

Map illustrating the 11 sampling locations of American pika (Ochotona princeps) populations distributed across the western United States. Shape and color correspond to region and sampling location and match precisely with the color and shapes used in subsequent figures. Populations representing the three sampled subspecies (O. p. princeps, O. p. schisticeps, and O. p. saxatilis) are labeled with text on the map.

Here, we expand on previous genetic studies of O. princeps by using a reduced representation GBS approach (ddRADseq) [51, 52]. We use these data to quantify patterns of population genetic variation across differing spatial scales: (1) among lineages; (2) among populations within lineages; and (3) at a finer scale within an isolated region of the Great Basin. Based upon previous studies we expected significant genetic divergence among lineages and strong evidence for isolation of the Great Basin populations. We predicted concordance between levels of genetic diversity and overall habitat area within mountain ranges and therefore expected lower levels of genetic variation in the peripheral and isolated populations of the Bodie Hills and the Ruby-East Humboldt Mountains. Finally, we expected the spatial genetic structure of populations to be jointly influenced by geography and environment, both across the range and across finer scales within the isolated Great Basin populations.


Reduced representation GBS libraries sequenced on two lanes of the Illumina HiSeq 4000 platform generated ~ 366 million total reads after initial filtering for potential contaminant reads, including those representing Illumina adaptors and PCR primers. After de-multiplexing by barcodes and matching reads to individual sample IDs, we retained an average of 1,909,370 reads (sd ± 394,102) per individual across 171 individuals for analyses. Assemblies run with bwa aligned 86,751,462 reads across all individuals (mean of 507,318 mapped reads per individual) to the O. princeps reference genome [56]. We identified 144,611 SNPs in the combined alignments of these individuals, but after filtering based on mapping quality, base quality, coverage and minor allele frequency, we retained a subset of 27,670 SNPs, with mean coverage depth per individual per locus of 6.2×, for downstream analyses (Fig. 2).

Spatial genetic structure

We used a hierarchical Bayesian model (entropy) and PCA to characterize regional-scale (CA, MT, NV, and CO) genetic structure using the entire dataset of 171 individuals. entropy analyses supported four differentiated genetic clusters (k = 4) (Additional file 1: Table S1), and all individuals exhibited ancestry coefficients of close to 100% for one of the four (Fig. 3). Three of the clusters represented the subspecies, O. p. princeps (Montana samples; brown), O. p. saxatilis (Colorado samples; green) and O. p. schisticeps (eastern California samples; pink). The fourth cluster was comprised of O. p. princeps individuals from the Ruby-East Humboldt Mountains (Nevada samples; purple). PCA revealed patterns of population genetic structure consistent with the ancestry estimates from entropy (Fig. 3). The first two principal components (PCs) accounted for 59.5 and 19.4% of the variation in the genotype probabilities, respectively, and revealed significant differentiation among all three subspecies (O. p. princeps, O. p. schisticeps, O. p. saxitlis) as well as O. p. princeps from Nevada (Fig. 3). The first PC separated O. p. princeps individuals (Montana and Nevada populations), while PC2 identified a distinct break between individuals sampled from the southwestern region of the species’ range (California and Nevada individuals, O. p. schisticeps and O. p. princeps, respectively) and the northeast region (Colorado and Montana; O. p. saxitilis and O. p. princeps, respectively) (Fig. 3).

Fig. 3

Figure created by author JPJ

Population genetic structure of 11 populations of the American pika (Ochotona princeps) sampled from California (O. p. schisticeps), Colorado (O. p. saxitilis), Nevada (O. p. princeps) and Montana (O. p. princeps) based on 27,670 SNPs. a The first and second principal components (PCs) from the PCA are plotted for each individual (points were jittered to avoid overplotting). b An unrooted neighbor-joining tree analysis produced a tree topology reflecting strong evolutionary differentiation among pika subspecies. c A hierarchical Bayesian model (entropy) shows the assignment probability estimated for each individual for each of four genetic clusters (k = 4). Symbols are only used to delineate different populations within the same region and do not correspond to those found in Fig. 2.

Estimates of Nei’s D indicated the presence of hierarchical genetic differentiation within and among subspecies and geographic regions (Fig. 3b and Additional file 2: Fig. S1). The Sierra Nevada lineage (O. p. schisticeps) was highly differentiated from both the Northern Rockies lineage (O. p. princeps; Nei’s D range: 0.400–0.503) and the Southern Rockies lineage (O. p. saxatilis; Nei’s D range: 0.376–0.390). The Southern Rockies lineage was also well differentiated from the Northern Rockies lineage (Nei’s D range: 0.214–308). Within O. p. princeps, Great Basin populations exhibited modest differentiation from Montana populations (Nei’s D range: 0.111–0.116). Finally, genetic distances between populations from the Ruby Mountains and the East Humboldt Range were relatively small (Nei’s D range: 0.011–0.012).

Despite low overall genetic divergence among populations within each lineage or region, our analyses consistently revealed evidence for fine-scale spatial structure within them. Based on PCA, clear differentiation was observed within lineages where the sampled populations were separated by large geographic distances, including the pair separated by ~ 18 km within the Northern Rockies (Emerald Lake and Swan Creek) and the pair separated by ~ 38 km within the Sierra Nevada (Bodie and Pipet Tarn) (Fig. 4). For the six Great Basin populations sampled within the Ruby-East Humboldt Mountains (separated by a maximum distance of 70 km) fine-scale spatial structure was particularly striking. The first two PCs explained 32.3% of the variance and strongly differentiated between populations separated by the low elevation Secret Pass within this otherwise continuous mountain range (Fig. 5a). The Island Lake population was separated along PC3 axis from the other two populations in the Ruby portion of the range, Overland Lake and Hidden Lake, despite its central location within the range. Smith Lake was separated from Lizzie’s Basin and Week’s Creek in the East Humboldt portion of the range along the PC4 axis (Fig. 5b). The entropy model based on k = 2 (Fig. 5c) was strongly supported, but the k = 6 model had similar support and reflected clear fine scale spatial structure analogous to PC3 and PC4 (Fig. 5b, Additional file 1: Table S2).

Fig. 4

Figure created by author JPJ

Principal components analysis is consistent with population genetic structure for American pika (Ochotona princeps) populations from a eastern California and b Montana. Symbols are only used to delineate different populations within the same region and do not correspond to those found in Fig. 2.

Fig. 5

Figure created by author JPJ

Population genetic structure of six populations of the American pika (Ochotona princeps) from the Ruby Mountains (RM) and the East Humboldt range (EH) in eastern Nevada. The a first, second, b third, and fourth principal components (PCs) from the PCA are plotted for each individual. c A hierarchical Bayesian model (entropy) shows the assignment probability estimated for each individual for each of two genetic clusters (k = 2). Symbols are used only to delineate different populations within the same region and do not correspond to those found in Fig. 2.

Geographic and environmental predictors of spatial genetic structure

Across all study populations, ~ 80% of the climatic variation was explained by the first two PCs. The first climate PC (hereafter Climate 1) had strong loadings for variables related to temperature, whereas precipitation and seasonality variables loaded more heavily on the second climate PC (hereafter Climate 2) (Additional file 1: Table S3). Of the four predictor variables, only geographic distance and Climate 1 distance were strongly associated with one another across all populations (R2 = 0.344; P = 0.010; Additional file 3: Fig. S2). Nei’s D was strongly predicted by Climate 2 distance (R2 = 0.538; P = 0.009) and geographic distance (R2 = 0.342; P = 0.009), but was not influenced by Climate 1 distance (R2 = 0.129; P = 0.079) or elevational distance (R2 = 0.000; P = 0.960) (Table 2; Fig. 6). A model with both Climate 2 distance and geographic distance explained 67.8% of the variation in Nei’s D and had comparable explanatory power relative to models containing more parameters (Table 2).

Table 2 Comparison of MRM tests predicting Nei’s genetic distance (D) among all 11 populations
Fig. 6

Figure created by author JPJ

The relationship between Nei’s D and elevational distance (m), haversine geographic distance (km), climate 1 and climate 2 predictor variables (19 climatic variables downloaded from; see supplementary tables for specific variable identities). Panels a, c, e, and g represent all sampled pika populations, and panels b, d, f, and h represent the subset of six Nevada populations located along an elevational gradient between the Ruby Mountains and East Humboldt Range.

For the Great Basin subset of populations, ~ 91% of the variation in climate among sites was explained by the first two PCs. As in the full dataset above, the first two PCs had strong loadings for variables related to temperature and precipitation, respectively, but seasonality variables loaded more strongly on PC1 for the Nevadan subset (Additional file 1: Table S4). Elevational and geographic distance were correlated with one another for the subset of Nevadan populations (R2 = 0.250; P = 0.046), but neither climate variable was strongly coupled with any other predictor (Additional file 3: Fig. S2). Nei’s D among populations in the Ruby Mountains and East Humboldt Range was strongly affected by geographic distance (R2 = 0.596; P = 0.024), but was not associated with Climate 1 distance (R2 = 0.020; P = 0.399), Climate 2 distance (R2 = 0.094; P = 0.130), or elevational distance (R2 = 0.016; P = 0.284) (Table 3; Fig. 6). Alternative models including geographic distance plus additional predictors did have improved explanatory value, but the gains were fairly modest relative to the more parsimonious model including only geographic distance (Table 3).

Table 3 Comparison of multiple regression on distance matrices (MRM) predicting Nei’s genetic distance (D) among the six Nevadan populations

Levels of genetic diversity

Estimates of nucleotide diversity (π) [77] and Watterson’s theta (θW) [78] were highly variable across the sampled localities (π range: 0.0006–0.0027; θW range: 0.0005–0.0023; Fig. 7; Additional file 1: Table S5). Populations from Montana and Colorado had the highest levels of standing variation (π range: 0.0025–0.0027; θW range: 0.0021–0.0023), while those sampled in the Ruby-Humboldt ranges had by far the lowest (π range: 0.0007–0.0009; θW range: 0.0005–0.0007) (Fig. 7; Additional file 1: Table S5). All populations had positive estimates of Tajima’s D (consistent with recent population contraction), with the largest estimate found in Bodie, CA (Tajima’s D = 0.811) and the lowest estimate in Lizzie’s Basin, NV (Tajima’s D = 0.240) (Additional file 1: Table S5).

Fig. 7

Photo provided by author KBK. Figure created by author JPJ

Estimates of Watterson’s theta (θW) and nucleotide diversity (π) are displayed for all eleven sites. Symbols are used only to delineate different populations within the same region and do not correspond to those found in Fig. 2. Insert: a pika at Bodie Hills, CA.


For many North American taxa, including the pika, persistence in Pleistocene-era refugia enabled survival across repeated glacial cycles and led to strongly divergent populations or lineages [22, 23, 87]. Our results support earlier range wide analyses conducted by Galbreath et al. [22, 23], which show clear genetic differentiation among the three pika mtDNA lineages sampled here—Northern Rockies, Southern Rockies and Sierra Nevada—consistent with glacial cycles and periods of isolation during the Pleistocene. The Galbreath et al. [22] phylogeny based upon mtDNA cytochrome b and D-loop sequence data placed the populations in the Ruby-Humboldt Mountains of Nevada into a clade with the other southern most populations of the Northern Rocky Mountain lineage (O. p. princeps). However, additional sequence data from nuclear introns did not further clarify relationships among the Ruby-East Humboldt populations and other populations in the Northern Rocky Mountain clade despite the geographic isolation of this mountain range [23]. In contrast, increased genomic sampling in our data provides clear evidence of differentiation of the Ruby-East Humboldt populations from other populations in the Northern Rocky Mountain clade. Bayesian clustering analyses, PCA, and pairwise estimates of genetic distance from this study reveal the distinctiveness of the pikas in the Ruby-East Humboldt Mountains.

Spatial genetic structure and its predictors

While geographic distance predicted genetic distance among all populations, a composite climate distance (Climate 2; precipitation and seasonality variables) was more strongly associated with genetic distance (Table 2, Additional file 1: Table S3; Fig. 6). Importantly, geographic and Climate 2 distance were not correlated with one another across all populations (Additional file 3: Fig. S2), suggesting that our results are consistent with both isolation-by-distance and isolation-by-environment [88] and that the evolutionary history of American pikas has been shaped in part by adaptation to climatic variability across the species range. Indeed, the relative importance of the interaction between temperature stress and the seasonality [89], location [90], and form of precipitation (i.e. rain, snowpack, moisture) [33] are increasingly being recognized as important determinants of pika distribution and persistence under a changing climate. Unfortunately, our sequencing dataset—and RADseq/GBS datasets in general—has limited utility for identifying genomic regions potentially associated with local adaptation because the decay of linkage disequilibrium typically occurs at a scale much smaller than typical marker densities ([91,92,93]; but see [94, 95]). Marker density for our dataset was roughly one locus per 80,000 bp, suggesting that more thorough sequencing will be needed to effectively survey the genome.

In contrast, genetic distance at smaller spatial scales within the Ruby Mountains and East Humboldt Range was overwhelmingly associated with geographic distance (Table 3; Fig. 6). A lack of evidence supporting isolation-by-environment for the Great Basin populations in the Ruby-East Humboldt range is not surprising, as they experience similar climatic conditions, diverged relatively recently (Fig. 3b and Additional file 2: Fig. S1), and were potentially connected during the Last Glacial Maximum when the upper reaches of the Ruby Mountains and the East Humboldt Range were glaciated [96]. Isolation-by-distance has been previously reported for Great Basin pikas [97] and has also been documented in other Great Basin alpine organisms, including marmots [98] and forbs [99].

Variation in genetic diversity across lineages and populations

Genetic diversity was an order of magnitude lower in the Ruby-East Humboldt range populations compared to all other populations sampled (Table 3), including the Bodie Hills population, which is found in highly fragmented habitat surrounded by desert-shrub vegetation and isolated from the main Sierra Nevada range. Interestingly, the nucleotide diversity estimates in Ruby-East Humboldt pika populations were similar to those generally found for arctic and subarctic alpine mammals (ddRADseq data: collared pika, O. collaris, π = 0.00011–0.00034; hoary marmot, Marmota caligata, π = 0.00009–0.00063; singing vole, Microtus miurus, π = 0.00032–0.00089; brown lemming, Lemmus trimucronatus, π = 0.00079–0.00108; arctic ground squirrel, Urocitellus parryii, π = 0.00043–0.00069; [100]). The low levels of genetic variation in these arctic species are likely due to the influences of limited glacial refugia during the Pleistocene, as dispersal ability did not explain the observed genetic patterns [100].

In a general sense, the levels of genetic diversity within the pika populations sampled for this study were related to habitat size, insofar as populations from the larger mountain ranges had higher diversity levels (Fig. 7; Additional file 1: Table S5). Although the majority of our sampled populations are located in high elevation habitat with continuous to semi-continuous talus, such habitat characteristics in and of themselves did not ameliorate the loss of genetic variation in the Ruby-East Humboldt populations. In contrast, although talus habitat is limited to widely spaced rocky outcrops across the Bodie Hills plateau with very few sites currently occupied [101], the Bodie Hills population sampled for this study inhabits ~ 100 talus patches formed from ore dumps associated with hard rock gold mining in the later nineteenth century [35, 102, 103]. Genetic diversity at this site was similar to the levels observed in the Rocky Mountains and Sierra Nevada. The difference in levels of genetic diversity found between the Ruby-East Humboldt range and the Bodie Hills, especially given the differences in habitat, suggests potential roles for both the spatial extent of habitat and configuration of local habitat patches in influencing rates of gene flow and thus the maintenance of genetic variation [35].

Metapopulation dynamics, local habitat spatial structure and maintenance of genetic diversity

An extinction and colonization dynamic among semi-independent subpopulations underpins metapopulation theory [104] and predicts reductions in overall genetic diversity in metapopulations over time for scenarios of limited source populations, low recolonization rates, and genetic drift [105,106,107,108]. There are numerous empirical examples from a diverse set of taxa that associate lower genetic diversity with a metapopulation dynamic (e.g., hyrax, Heterohyrax brucei and Procavia johnstoni [109]; moths, Aglaope infausta [110]; trout, Oncorhynchus clarkii henshawi [111]; caddisflys, Allogamus uncatus [112]; treefrogs, Hyla wrightorum [113]).

The fragmented nature of pika habitat, whether at small spatial scales as in the Bodie Hills or larger habitat patches in mountain ranges such as in the Ruby-East Humboldt, Sierra Nevada or Rocky Mountains, suggests pikas across their range may experience metapopulation dynamics albeit at very different temporal and spatial scales, which in turn may result in differing genetic signals. Metapopulation dynamics may therefore partially explain the low genetic diversity observed in pika populations found in the Ruby-East Humboldt range.

The basin and range topography of the Great Basin dates to approximately 30 million years ago and pikas are thought to have colonized these interior basin mountain ranges more recently during the Pleistocene [29]. As the climate warmed, pikas disappeared from most of the Great Basin mountains ranges and are currently found only in the largest ranges with the highest elevational extent [29]. The genetic data from this study suggest that the pika populations in the Ruby-East Humboldt range have long been isolated from other populations in their own lineage. The narrow linear configuration of the Ruby-East Humboldt mountain range limits the amount and spatial extent of talus and may thus constrain both dispersal distance and direction. Such limitations may increase population genetic structure and reduce colonization potential of extirpated talus patches. Although the Bodie Hills plateau was likely colonized by pikas from the main Sierra Nevada also during the Pleistocene, this high elevation spur of the main range has remained connected to the Sierra Nevada and thus ongoing gene flow was at least possible. However, the population sampled for this study is one of the very few extant populations remaining in the Bodie Hills. As a result, the genetic diversity found within the Bodie Hills population might seem anomalous given that it is an isolated metapopulation in highly fragmented habitat [35, 41, 45]. Earlier genetic work using DNA minisatellite markers showed that average heterozygosity for the Bodie Hills population was similar to the Pipet Tarn population inhabiting the larger and contiguous habitat of the main Sierra Nevada [35]. One explanation for the relatively high genetic diversity at Bodie (Table 3), despite the effects of genetic drift, involves the large number of small habitat patches at this site as well as a relatively recent founding event, which is estimated to have occurred sometime during the late nineteenth century [103, 114, 115]. If extinction probabilities are uncorrelated among these patches, then Bodie can act as a highly persistent metapopulation of small local populations that drift toward fixation for different alleles, thereby maintaining genetic diversity at the metapopulation level similar to that of a large un-subdivided population [115]. Furthermore, dispersal among extant patches has been very fluid in this population, driven by juveniles leaving small patches (when territories are unavailable) and having multiple habitat patches within dispersal distance [35, 38]. Therefore, ongoing within-population dispersal may have acted to spread and maintain genetic variation among patches in between extinction events [35].

The role of spatial configuration of talus patches, in contrast to overall habitat size within mountain ranges, has not received as much attention as other potential correlates of extirpation or persistence probabilities for pika populations (but see [116]). Examples from the literature, for metapopulations that have retained genetic variation, suggest that ongoing dispersal and gene flow among extant subpopulations, in between extinction events, can in some cases be critical for maintenance of genetic diversity (e.g., natterjack toad, Bufo calamita [117]; black-tailed prairie dogs, Cynomys ludovicianus [118]; the treacle-mustard, wormseed wallflower Erysimum cheiranthoides [119]; Galapagos warbler finches (Certhidea) [120]: cichlid fishes, Eretmodus cyanostictus, Variabilichromis moorii and Tropheus moorii [121]). Levels of genetic diversity within metapopulations will ultimately depend upon the number of local populations, extirpation-colonization rates and gene flow among local extant populations [105, 113, 114, 122].

Implications for pika persistence

Our results agree with previous analyses suggesting vulnerability of Great Basin pika populations [13,14,15, 31, 32]. Pika populations in the Great Basin ecoregion represent two mtDNA lineages including the Sierra Nevada and Northern Rocky Mountains. The topography of the Great Basin effectively isolates populations within these ranges and while extirpations have been occurring since the end of the Pleistocene, they have accelerated over the late twentieth and early twenty-first centuries [13,14,15, 29]. Few of the extant pika populations in the Great Basin have been characterized genetically, but the distinctness of the populations in the Ruby-East Humboldt Mountains suggests these ranges may harbor unique genetic variation. The accelerated rate of pika population losses in recent decades in this ecoregion could certainly be associated with range isolation, mountain range size, as well as mean elevation and a warming climate [13, 15, 46]. This is especially concerning considering that the most evolutionarily distinct lineage sampled in this study (O. p. schisticeps; Fig. 2b) is also the lineage that occurs in both eastern California and the Great Basin [22, 23]. These two regions, which have produced most of the evidence for climate change-induced pika population extirpations [17, 33, 123, 124], are also predicted to lose more of the currently suitable pika habitat during this century [16, 22, 125].

Habitat fragmentation together with increasing temperatures and expansion of invasive species are reducing the habitable landscape for diverse alpine taxa including alpine vascular flora [126, 127], mammals [12, 128,129,130], and aquatic invertebrates [131, 132]. The Grinnell resurvey project in Yosemite National Park showed an upslope contraction of lower elevational limits for half of 28 small mammal species surveyed, consistent with the ~ 3 °C increase in minimum temperatures observed over the past century [128]. Despite the maintenance of genetic variation in some of the populations sampled for this study, Tajima’s D estimates indicate that all have suffered recent population contractions, with the Bodie Hills population showing the greatest contraction (Additional file 1: Table S5). Earlier work across California, which included the pika population at the Bodie Hills site, revealed that both temperature and habitat area influenced pika occupancy [17, 38]. There has also been a concomitant reduction in allelic diversity and effective population size in the Bodie Hills population over the past 25 years [35, 38].

For the American pika, climate change is likely to reduce the already low rates of gene flow among local populations [44, 133] through increased mortality and reduced dispersal [116], thereby intensifying the impacts of genetic drift on standing genetic variation. The majority of pika populations sampled had levels of nucleotide diversity (Fig. 7; Additional file 1: Table S5) below the average estimate for a wide range of taxa (green plants, animals, fungi and Chromalveolates; mean = 0.0065, π = 0.0005–0.05; [134]). However, mammalian species in general tend to have lower levels of diversity, compared to the other taxa, where diversity is likely to co-vary with geographic range and life history [100, 134]. The very low genetic diversity of the Ruby-East Humboldt populations raises concerns about the effects of genetic drift on genetic variation in isolated populations and highlights the need for more extensive sampling of this and other insular ranges.


Expanding the genetic characterization of population level genetic diversity will be necessary to define fully the genetic variation within pika lineages and how this diversity is partitioned across the landscape. Our results suggest that the natural fragmentation inherent to rugged, mountainous terrain may have disproportionate effects [135] on pikas and other small alpine mammals by constraining connectivity critical for maintaining local demographic stability [41]. The patterns of climate-associated, range-wide differentiation observed in this study, along with the significant role of geography at finer scales, reinforces the need to tease apart the specific climatic or non-climatic factors that can shape genetic variation across small to moderate spatial scales within each mountain system. Indeed, while these findings provide important context regarding the effects of geography and climate on the amount and distribution of extant genetic diversity, our understanding of the potential adaptive diversity that these genetic data represent in pikas is unknown. Furthermore, while adaptive genetic diversity is likely to play a role in persistence, recent work has emphasized the relevance of landscape context at the ecoregional scale, not genetic lineage, in structuring intraspecific variation in response to climate change [136]. Given these insights, genomic analyses, like those conducted here, at within-mountain range scales, may be most useful for identifying those portions of the pika’s range that contain talus of sufficient quality, quantity and spatial configuration to support maintenance of species-wide genetic diversity.


DNA resources

All individual pikas included in this study were live-trapped and released at point of capture after processing. No animals were euthanized for this study. The majority of tissue samples were collected over the last ten years as part of long-term monitoring efforts or recent live-trapping studies (Table 1). Additional DNA samples curated from earlier studies were obtained for a subset of sites (Ruby Mountains [Overland Lake, Island Lake, Hidden Lake] and East Humboldt Range [Week’s Creek, Lizzie’s Basin and Smith Lake]; Montana sites, Emerald Lake and Swan Creek). Samples were collected from 11 distinct populations, which represent three subspecies (O. p. princeps, O. p. schisticeps, and O. p. saxitilis) and three of the five mtDNA lineages (Sierra Nevada, Northern and Southern Rocky Mountains) (Table 1; Fig. 2).

Live-trapping for pikas from the Bodie and Pipet Tarn sites was conducted as follows. Pikas were trapped using Tomahawk live traps (16″ × 6″ × 6″) baited with native alpine vegetation and covered with rocks to prevent exposure to weather and predators [35, 37]. Traps were set predawn and checked within 2 h. Trapped individuals were treated with an inhalant anesthetic (isoflurane) before a small sample of ear tissue (< 50 mg) was collected and placed in cryovials and stored in liquid nitrogen. Animals were released at the point of capture after 15 min of recovery. We followed the guidelines of the American Society of Mammalogists for live animal research [53]. Permission for live capture and release of pikas was approved by Institutional Animal Care and Use Committees (IACUC) at the University of Nevada Reno and University of Colorado, Boulder. Trapping permits were obtained from California Department of Parks and Recreation; California Department of Fish and Wildlife; Colorado Parks and Wildlife; Montana Fish, Wildlife and Parks; Nevada Department of Wildlife, and US Department of Agriculture Forest Service.

Illumina sequencing of restriction fragment libraries

DNA was extracted from tissue samples representing 35 O. p. schisticeps, 21 O. p. saxitilis, and 115 O. p. princeps individuals (Fig. 2; Table 1) using a Qiagen DNeasy Blood and Tissue Kit and quantified on a QIAxpert UV/VIS spectrophotometer (Qiagen, Inc., Valencia, CA, USA). Reduced representation libraries were prepared for Illumina sequencing following a genotyping-by-sequencing (GBS) protocol [51] analogous to ddRADseq [52]. Genomic DNA was digested with two restriction endonucleases (EcoRI and MseI) and double-stranded adaptor oligonucleotides were ligated to the digested fragments. Adaptor ligation occurred simultaneously with restriction digestion using adaptor sequences that incorporated unique 8–10 base pair (bp) barcodes for each individual, as well as the priming sites for Illumina sequencing. Ligated fragments were PCR amplified using a high-fidelity proofreading polymerase, uniquely barcoded products were pooled into a single library and then size selected for fragments ranging from 350 to 425 bp in length using a Blue Pippin unit (Sage Science, Beverly, MA, USA). Single end 150 bp reads were generated with two lanes of sequencing on an Illumina HiSeq 4000 at the University of Texas Genomic Sequencing and Analysis Facility (Austin, TX, USA).

DNA sequence assembly, and variant calling

We first filtered out low-quality reads, those representing common bacterial contaminants, and those associated with Illumina oligos using bowtie_db2 [54] and the tapioca pipeline ( Individual bar codes and the adjacent six bases corresponding to the EcoRI cut sites were then trimmed from each read, and individual identifiers were incorporated into separate fastq files for each individual (available at Dryad; [55]). We aligned all reads to the draft O. princeps genome (NCBI identifier: GCF_000292845.1; OchPri3.) [56] using the aln and samse algorithms of bwa v0.7.8 [57] with an edit distance of four and otherwise default parameter settings.

We used the mpileup command in samtools v0.1.19 [58] to merge the .bam alignment files of all individuals into .bcf formatted files. We subsequently used bcftools v0.1.19 [58]) to identify bi-allelic single nucleotide polymorphisms (SNPs) and to calculate genotype likelihoods (-e option). We used the default scaled substitution mutation rate, a full prior setting (-P option) and called SNPs if P(ref|D) < 0.01 (vcftools v0.1.14) [59] was used to perform additional quality filtering. We retained SNPs in our final data set when 92% of the individuals had at least one read at the locus. We excluded variable sites with more than one alternate allele and loci with minor allele frequencies less than 5%. Finally, we randomly selected a single SNP per every 1000 bp to increase the independence of loci used in subsequent analyses. Additional parameter settings for these analyses are available at Dryad [55].

Spatial genetic structure

Analyses were performed on the full set of all eleven populations as well as a subset of six populations of O. p. princeps sampled along semi-continuous talus habitat across the Ruby–East Humboldt mountain range in eastern Nevada. To infer individual ancestry proportions (q) and the number of ancestral populations represented (k) in the dataset, we used a hierarchical Bayesian model (entropy v1.2) [60] that is based upon the commonly implemented model of structure [61, 62]. Entropy simultaneously estimates genotype probabilities and individual ancestry coefficients while incorporating genotype uncertainty inherent to low-to-medium coverage depth sequencing designs, as well as error associated with sequencing or subsequent mapping [60, 63, 64]. We ran five replicate MCMC chains for models ranging from k = 2 to k = 11 to determine the most probable k for the entire set of individuals, and for k = 2 to k = 8 for the Ruby-East Humboldt dataset. To speed the stabilization of MCMC chains, we initialized ancestry coefficients using cluster assignment probabilities generated from k-means clustering and linear discriminant analysis of principal component scores using the mass package in R v3.6.3 [65, 66]. All MCMC chains were run for 40,000 steps following 10,000-step burn-ins and thinning for one out of ten steps. Deviance Information Criterion (DIC) was used for model comparison. We further summarized patterns and levels of population genetic variation across all samples using the estimated genotype probabilities produced by entropy as the input variables for principal component analysis (PCA) performed using the prcomp function in R. The matrix of genotype probabilities is available at Dryad [55].

To evaluate the potential for population genetic structure across a range of geographic scales, we conducted separate PCA analyses on: (1) the full set of 171 individuals; (2) the subset of individuals sampled from the Ruby-East Humboldt Mountains in Nevada; (3) the subset of individuals from California (Bodie and Pipet Tarn); and (4) the subset of individuals sampled from Montana (Emerald Lake and Swan Creek). PCA is an effective model-free tool for illustrating spatial genetic structure across a continuum of differentiation, including among weakly differentiated populations [67,68,69,70], where even axes explaining lower levels of variance can reveal spatial genetic structure and population identifiability [71]. We further summarized patterns of genetic differentiation based on allele frequency variation by calculating Nei’s genetic distance (D) [72] for each pairwise comparison among populations.

Geographic and environmental predictors of spatial genetic structure

To quantify the relative importance of climatic and geographic factors as predictors of spatial genetic structure, we implemented a series of multiple regressions on distance matrices (MRMs) within a model comparison framework. Separate analyses were conducted for both the entire eleven population dataset and the subset of six populations from the Ruby Mountains and East Humboldt Range. Haversine geographic distances between sampling locations were calculated based on the midpoints of latitude and longitude for each sampling region using the fossil package in R [73], whereas elevational distance was defined as the absolute difference in elevation between sites. Climatic variation was characterized for sampling localities using PCA based on 19 climatic variables downloaded for each site (1 km2 scale; [74]. For each of the first two principal components (PCs) from the climate PCA (proportion of variance explained: full model = 46.9% and 33.8%; Ruby-East Humboldt subset = 60.0% and 31.0%) climatic distance among sites was calculated using the dist function in R (hereafter referred to as the Climate 1 and Climate 2 predictor variables). MRMs were performed using the ecodist package in R [75, 76], with Nei’s D as the response variable. All possible combinations of MRMs including 1, 2, 3, and 4 predictor variables (geographic distance, elevational distance, Climate 1 distance, Climate 2 distance) were analyzed.

Quantifying levels of genetic diversity

We estimated nucleotide diversity (π) [77] and Watterson’s theta (θW) [78] with ANGSD v0.921 [79] following the empirical Bayes method described by Korneliussen et al. [80]. Diversity metrics were calculated for each population independently using individual.bam files. First, site allele frequency likelihoods were estimated using the "-doSAF 1" command and specifying the samtools genotype likelihood model. Next, the folded site frequency spectrum (SFS) was estimated based on allele frequencies using the "realSFS" function. π and θW were calculated with the "doThetas 1" command using the SFS likelihoods as priors for each site. The "thetaStat" command was used to summarize values for each scaffold, and mean π and θW across scaffolds were calculated in R.

We also estimated Tajima’s D [81] using ANGSD. Tajima’s D is based on the difference between π and θW, which are expected to be equivalent for populations at mutation-drift equilibrium [82]. However, these estimates of θ inherently differ in their response to selection or demographic events: π > θW when the SFS is skewed towards common variants (positive Tajima’s D), whereas π < θW when the SFS is skewed towards rare variants (negative Tajima’s D) [81]. During a bottleneck, low-frequency variants are disproportionately lost from a population [83], resulting in a SFS with more common variants than expected (positive Tajima’s D). Following a bottleneck, a population can rapidly expand, resulting in an increased number of mutations per generation and thus a SFS containing more rare variants than expected (negative Tajima’s D). It is worth noting that this is a broad oversimplification of the effects of a bottleneck on Tajima’s D, and that the resulting estimates can be strongly influenced by population size, the strength of the bottleneck, the time since the bottleneck, and the rate of subsequent population growth [82, 84,85,86]. For each population, Tajima’s D was calculated for 50 kbp windows across the genome (50 kbp step size) using the estimates of π and θW generated above. The mean value of Tajima’s D across all windows was calculated in R.

Availability of data and materials

A directory containing the compressed individual.fastq files and a matrix of genotype probabilities used for downstream analyses can be found at Dryad Digital Repository



Double digest restriction-site associated DNA sequencing


Deviance information criterion




Multiple regression on distance matrices


Principal component


Principal component analysis


Single nucleotide polymorphism


  1. 1.

    Keppel G, Van Niel KP, Wardell Johnson GW, Yates CJ, Byrne M, Mucina L, et al. Refugia: identifying and understanding safe havens for biodiversity under climate change. Glob Ecol Biogeogr. 2012;21:393–404.

    Article  Google Scholar 

  2. 2.

    Epps CW, Wehausen JD, Bleich VC, Torres SG, Brashares JS. Optimizing dispersal and corridor models using landscape genetics. J Appl Ecol. 2007;44:714–24.

    Article  Google Scholar 

  3. 3.

    Ashcroft MB. Identifying refugia from climate change. J Biogeogr. 2010;37:1407–13.

    Article  Google Scholar 

  4. 4.

    Varner J, Dearing MD. The importance of biologically relevant microclimates in habitat suitability assessments. PLoS ONE. 2014;9:e104648.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  5. 5.

    Li J, McCarthy TM, Wang H, Weckworth BV, Schaller GB, Mishra C, et al. Climate refugia of snow leopards in High Asia. Biol Conserv. 2016;203:188–96.

    Article  Google Scholar 

  6. 6.

    Hampe A, Jump AS. Climate relicts: past, present, future. Annu Rev Ecol Evol Syst. 2011;42:313–33.

    Article  Google Scholar 

  7. 7.

    Scherrer D, Körner C. Topographically controlled thermal-habitat differentiation buffers alpine plant diversity against climate warming. J Biogeogr. 2011;38:406–16.

    Article  Google Scholar 

  8. 8.

    Williams CM, Henry HAL, Sinclair BJ. Cold truths: how winter drives responses of terrestrial organisms to climate change. Biol Rev. 2015;90:214–35.

    Article  PubMed  Google Scholar 

  9. 9.

    Rodhouse TJ, Hovland M, Jeffress MR. Variation in subsurface thermal characteristics of microrefuges used by range core and peripheral populations of the American pika (Ochotona princeps). Ecol Evol. 2017;7:1514–26.

    Article  PubMed  PubMed Central  Google Scholar 

  10. 10.

    Graham CH, VanDerWal J, Phillips SJ, Moritz C, Williams SE. Dynamic refugia and species persistence: tracking spatial shifts in habitat through time. Ecography. 2010;33:1062–9.

    Article  Google Scholar 

  11. 11.

    Gunderson AM, Lanier HC, Olson LE. Limited phylogeographic structure and genetic variation in Alaska’s arctic and alpine endemic, the Alaska marmot. J Mammal. 2012;93:66–75.

    Article  Google Scholar 

  12. 12.

    Rubidge EM, Patton JL, Lim ML, Burton AC, Brashares JS, Moritz C. Climate-induced range contraction drives genetic erosion in an alpine mammal. Nat Climate Change. 2012;2:285–8.

    Article  Google Scholar 

  13. 13.

    Beever EA, Brussard PF, Berger J. Patterns of apparent extirpation among isolated populations of pikas (Ochotona princeps) in the Great Basin. J Mammal. 2003;84:37–54.;2.

    Article  Google Scholar 

  14. 14.

    Beever EA, Ray C, Wilkening JL, Brussard PF, Mote PW. Contemporary climate change alters the pace and drivers of extinction. Glob Change Biol. 2011;17:2054–70.

    Article  Google Scholar 

  15. 15.

    Wilkening JL, Ray C, Beever EA, Brussard PF. Modeling contemporary range retraction in Great Basin pikas (Ochotona princeps) using data on microclimate and microhabitat. Quat Int. 2011;235:77e88.

    Article  Google Scholar 

  16. 16.

    Calkins MT, Beever EA, Boykin KG, Frey JK, Andersen MC. Not-so-splendid isolation: modeling climate-mediated range collapse of a montane mammal Ochotona princeps across numerous ecoregions. Ecography. 2012;35:780–91.

    Article  Google Scholar 

  17. 17.

    Stewart JAE, Perrine JD, Nichols LB, Thorne JH, Millar CI, Goehring KE, et al. Revisiting the past to foretell the future: summer temperature and habitat area predict pika extirpations in California. J Biogeogr. 2015;42:880–90.

    Article  Google Scholar 

  18. 18.

    Smith AT, Weston ML. Ochotona princeps. Mammal Sp Acc. 1990;352:1–8.

    Google Scholar 

  19. 19.

    MacArthur RA, Wang LCH. Physiology of thermoregulation in the pika Ochotona princeps. Can J Zool. 1973;51:11–6.

    CAS  Article  PubMed  Google Scholar 

  20. 20.

    Wilkening JL, Ray C, Varner J. Relating sub-surface ice features to physiological stress in a climate sensitive mammal, the American Pika (Ochotona princeps). PLoS ONE. 2015;10:e0119327.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  21. 21.

    Hafner DJ, Smith AT. Revision of the subspecies of the American pika, Ochotona princeps (Lagomorpha: Ochotonidae). J Mammal. 2010;91:401–17.

    Article  Google Scholar 

  22. 22.

    Galbreath KE, Hafner DJ, Zamudio KR. When cold is better: climate-driven elevation shifts yield complex patterns of diversification and demography in an alpine specialist (American pika, Ochotona princeps). Evolution. 2009;63:2848–63.

    CAS  Article  PubMed  Google Scholar 

  23. 23.

    Galbreath KE, Hafner DJ, Zamudio KR, Agnew K. Isolation and introgression in the Intermountain West: contrasting gene genealogies reveal the complex biogeographic history of the American pika (Ochotona princeps). J Biogeogr. 2010;37:344–62.

    Article  Google Scholar 

  24. 24.

    Hafner DJ, Sullivan RM. Historical and ecological biogeography of neararctic pikas (Lagomorpha: Ochotonidae). J Mammal. 1995;76:302–21.

    Article  Google Scholar 

  25. 25.

    Somers P. Dialects in southern Rocky-Mountain pikas, Ochotona-princeps (Lagomorpha). Anim Behav. 1973;21:124–37.

    Article  Google Scholar 

  26. 26.

    Conner DA. Seasonal changes in activity patterns and the adaptive value of haying in pikas (Ochotona princeps). Can J Zool. 1983;61:411–6.

    Article  Google Scholar 

  27. 27.

    Mead JI. Quaternary records of pika, Ochotona, in North America. Boreas. 1987;16:165–71.

    Article  Google Scholar 

  28. 28.

    Hafner DJ. Pikas and permafrost: post-Wisconsin historical zoogeography of Ochotona in the southern Rocky Mountains, U.S.A. Arctic Alpine Res. 1994;26:375–82.

    Article  Google Scholar 

  29. 29.

    Grayson DK. A brief history of Great Basin pikas. J Biogeogr. 2005;32:2103–11.

    Article  Google Scholar 

  30. 30.

    Grayson DK. The Late Quaternary biogeographic histories of some Great Basin mammals (western USA). Quaternary Sci Rev. 2006;25:2964–91.

    Article  Google Scholar 

  31. 31.

    Beever EA, O’Leary J, Mengelt C, West JM, Julius S, Green N, et al. Improving conservation outcomes with a new paradigm for understanding species’ fundamental and realized adaptive capacity. Conserv Lett. 2016;9:131–7.

    Article  Google Scholar 

  32. 32.

    Nichols LB, Klingler KB, Peacock MM. American pikas (Ochotona princeps) extirpated from the historic masonic mining district of eastern California. West N Am Nat. 2016;76:163–71.

    Article  Google Scholar 

  33. 33.

    Stewart JAE, Wright DH, Heckman KA. Apparent climate mediated loss and fragmentation of core habitat of the American pika in the Northern Sierra Nevada, California, USA. PLoS ONE. 2017;12:e0181834.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  34. 34.

    Castillo Vardaro JA, Epps CW, Frable BW, Ray C. Identification of a contact zone and hybridization for two subspecies of the American pika (Ochotona princeps) within a single protected area. PLoS ONE. 2018;3:e0199032.

    CAS  Article  Google Scholar 

  35. 35.

    Peacock MM, Smith AT. The effect of habitat fragmentation on dispersal patterns, mating behavior, and genetic variation in a pika (Ochotona princeps) metapopulation. Oecologia. 1997;112:524–33.

    Article  PubMed  Google Scholar 

  36. 36.

    Henry P, Henry A, Russello MA. Variation in habitat characteristics of American pikas along an elevation gradient at their northern range margin. Northwest Sci. 2012;86:346–50.

    Article  Google Scholar 

  37. 37.

    Robson KM, Lamb CT, Russello MA. Low genetic diversity, restricted dispersal, and elevation-specific patterns of population decline in American pikas in an atypical environment. J Mammal. 2016;97:464–72.

    Article  Google Scholar 

  38. 38.

    Klingler KB. An integrated investigation of the population genetics, physiological stress, and movement patterns in the American pika (Ochotona princeps). Doctoral dissertation. University of Nevada, Reno. 2017.

  39. 39.

    Castillo JA, Epps CW, Jeffress MR, Ray C, Rodhouse TJ, Schwalm D. Replicated landscape genetic and network analyses reveal wide variation in functional connectivity for American pikas. Ecol Appl. 2016;26:1660–76.

    Article  PubMed  Google Scholar 

  40. 40.

    Smith AT, Gilpin ME. Spatially correlated dynamics in a pika metapopulation. In: Hanksi I, Gilpin M, editors. Metapopulation biology: ecology, genetics and evolution. San Diego: Academic; 1997. p. 407–28.

  41. 41.

    Moilanen A, Smith AT, Hanski IA. Long-term dynamics in a metapopulation of the American pika. Am Nat. 1998;152:530–42.[2503:MDEOHQ]2.0.CO;2.

    CAS  Article  PubMed  Google Scholar 

  42. 42.

    Russello MA, Waterhouse MD, Etter PD, Johnson EA. From promise to practice: pairing non-invasive sampling with genomics in conservation. PeerJ. 2015;3:e1106.

    Article  PubMed  PubMed Central  Google Scholar 

  43. 43.

    Waterhouse MD, Erb LP, Beever EA, Russello MA. Adaptive population divergence and directional gene flow across steep elevational gradients in a climate-sensitive mammal. Mol Ecol. 2018;27:2512–3228.

    Article  PubMed  Google Scholar 

  44. 44.

    Peacock MM. Determining natal dispersal patterns in a population of North American pikas (Ochotona princeps) using direct mark-resight and indirect genetic methods. Behav Ecol. 1997;8:340–50.

    Article  Google Scholar 

  45. 45.

    Clinchy M, Haydon DT, Smith AT. Pattern does not equal process: what does patch occupancy really tell us about metapopulation dynamics? Am Nat. 2002;159:351–62.

    Article  PubMed  Google Scholar 

  46. 46.

    Beever EA, Ray C, Mote PW, Wilkening JL. Testing alternative models of climate-mediated extirpations. Eco Appl. 2010;20:164–78.

    Article  Google Scholar 

  47. 47.

    McDonald KA, Brown JH. Using montane mammals to model extinctions due to global change. Conserv Biol. 1992;6:409–15.

    Article  Google Scholar 

  48. 48.

    Southwick CH, Golian SC, Whitworth MR, Halfpenny JC, Brown R. Population density and fluctuations of pikas (Ochotona princeps) in Colorado. J Mammal. 1986;67:149–53.

    Article  Google Scholar 

  49. 49.

    Dearing MD. Disparate determinants of summer and winter diet selection of a generalist herbivore Ochotona princeps. Oecologia. 1996;108:467–78.

    Article  PubMed  Google Scholar 

  50. 50.

    Wilkening JL, Ray C. Characterizing predictors of survival in the American pika (Ochotona princeps). J Mammal. 2016;97:1366–75.

    Article  Google Scholar 

  51. 51.

    Parchman TL, Gompert Z, Mudge J, Schilkey F, Benkman CW, Buerkle CA. Genome wide association genetics of an adaptive trait in lodgepole pine. Mol Ecol. 2012;21:2991–3005.

    CAS  Article  PubMed  Google Scholar 

  52. 52.

    Peterson BK, Weber JN, Kay EH, Fisher HS, Hoekstra HE. Double digest RADseq: an inexpensive method for de novo SNP discovery and genotyping in model and non-model species. PLoS ONE. 2012;7:e37135.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  53. 53.

    Sikes RS, Gannon WL. Guidelines of the American Society of Mammalogists for the use of wild mammals in research. J Mammal. 2011;92:235–53.

    Article  Google Scholar 

  54. 54.

    Langmead B, Salzberg SL. Fast gapped-read alignment with Bowtie 2. Nat Methods. 2012;9:357–9.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  55. 55.

    Klingler KB, Jahner JP, Parchman TL, Ray C, Peacock MM. Data from: “Genomic variation in the American pika: signatures of geographic isolation and implications for conservation.” Dryad Dig Repos. 2020.

    Article  Google Scholar 

  56. 56.

    Fontanesi L, Di Palma F, Flicek P, Smith AT, Thulin C-G, Alves PC, Lagomorph Genomics Consortium. LaGomiCs—lagomorph genomics consortium: an international collaborative effort for sequencing the genomes of an entire mammalian order. J Hered. 2016;107:295–308.

    Article  PubMed  PubMed Central  Google Scholar 

  57. 57.

    Li H, Durbin R. Fast and accurate short read alignment with Burrows–Wheeler transform. Bioinformatics. 2009;25:1754–60.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  58. 58.

    Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. The sequence alignment/map format and SAMtools, 1000 genome project data processing subgroup author notes. Bioinformatics. 2009;25:2078–9.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  59. 59.

    Danecek P, Auton A, Abecasis G, Albers CA, Banks E, DePristo MA, et al. The variant call format and VCFtools. Bioinformatics. 2011;27:2156–8.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  60. 60.

    Gompert Z, Lucas LK, Buerkle CA, Forister ML, Fordyce JA, Nice CC. Admixture and the organization of genetic diversity in a butterfly species complex revealed through common and rare genetic variants. Mol Ecol. 2014;23:4555–73.

    Article  PubMed  Google Scholar 

  61. 61.

    Pritchard JK, Stephens M, Donnelly P. Inference of population structure using multilocus genotype data. Genetics. 2000;155:945–59.

    CAS  PubMed  PubMed Central  Google Scholar 

  62. 62.

    Falush D, Stephens M, Pritchard JK. Inference of population structure using multilocus genotype data: linked loci and correlated allele frequencies. Genetics. 2003;164:1567–87.

    CAS  PubMed  PubMed Central  Google Scholar 

  63. 63.

    Nielsen R, Paul JS, Albrechtsen A, Song YS. Genotype and SNP calling from next-generation sequencing data. Nat Rev Genet. 2011;12:443–51.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  64. 64.

    Buerkle A, Gompert Z. Population genomics based on low coverage sequencing: how low should we go? Mol Ecol. 2013;22:3028–35.

    CAS  Article  Google Scholar 

  65. 65.

    Venables WN, Ripley BD. Random and mixed effects. In: Modern applied statistics with S. Statistics and computing. Springer; 2002. p. 271–300.

  66. 66.

    R Core Team. R: a language and environment for statistical computing. Vienna: R Foundation for Statistical Computing; 2020.

  67. 67.

    Leslie S, Winney B, Hellenthal G, Davison D, Boumertit A, Day T, et al. The fine-scale genetic structure of the British population. Nature. 2015;519:309–14.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  68. 68.

    Jahner JP, Gibson D, Weitzman CL, Blomberg EJ, Sedinger JS, Parchman TL. Fine-scale genetic structure among greater sage-grouse leks in central Nevada. BMC Evol Biol. 2016;16:127.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  69. 69.

    Szulkin M, Gagnaire PA, Bierne N, Charmantier A. Population genomic footprints of fine-scale differentiation between habitats in Mediterranean blue tits. Mol Ecol. 2016;25:542–58.

    CAS  Article  PubMed  Google Scholar 

  70. 70.

    Vendrami DLJ, De Noia M, Telesca L, Handal W, Charrier G, Boudry P, et al. RAD sequencing sheds new light on the genetic structure and local adaptation of European scallops and resolves their demographic histories. Sci Rep-UK. 2019;9:7455.

    CAS  Article  Google Scholar 

  71. 71.

    McVean G. A genealogical interpretation of principal components analysis. PLoS Genet. 2009;5(10):e1000686.

    Article  PubMed  PubMed Central  Google Scholar 

  72. 72.

    Nei M. Genetic distance between populations. Am Nat. 1972;106:283–92.

    Article  Google Scholar 

  73. 73.

    Vavrek MJ. fossil: Palaeoecological and palaeogeographical analysis tools. Palaeontol Electron. 2011;14:16p.

    Google Scholar 

  74. 74.

    Fick SE, Hijmans RJ. WorldClim 2: new 1-km spatial resolution climate surfaces for global land areas. Int J Climatol. 2017;37:4302–15.

    Article  Google Scholar 

  75. 75.

    Goslee SC, Urban DL. The ecodist package for dissimilarity-based analysis of ecological data. J Stat Softw. 2007;22:1–19.

    Article  Google Scholar 

  76. 76.

    Lichstein JW. Multiple regression on distance matrices: a multivariate spatial analysis tool. Plant Ecol. 2007;188:117–31.

    Article  Google Scholar 

  77. 77.

    Tajima F. Evolutionary relationship of DNA sequences in finite populations. Genetics. 1983;105:437–60.

    CAS  Article  Google Scholar 

  78. 78.

    Watterson GA. On the number of segregating sites in genetic models without recombination. Theor Popul Biol. 1975;7:256–76.

    CAS  Article  PubMed  Google Scholar 

  79. 79.

    Korneliussen TS, Albrechtsen A, Korneliussen RN. ANGSD: analysis of next generation sequencing data. BMC Bioinform. 2014;15:2–13.

    Article  Google Scholar 

  80. 80.

    Korneliussen TS, Moltke I, Albrechtsen A, Nielsen R. Calculation of Tajima’s D and other neutrality test statistics from low depth next-generation sequencing data. BMC Bioinform. 2013;14:2–14.

    CAS  Article  Google Scholar 

  81. 81.

    Tajima F. Statistical method for testing the neutral mutation hypothesis by DNA polymorphism. Genetics. 1989a;123:585–95.

    CAS  Article  Google Scholar 

  82. 82.

    Hahn MW. Molecular population genetics. 1st ed. New York: Sinauer Associates; 2019.

    Google Scholar 

  83. 83.

    Nei M, Maruyama T, Chakraborty R. The bottleneck effect and genetic variability in populations. Evolution. 1975;29:1–10.

    Article  Google Scholar 

  84. 84.

    Tajima F. The effect of change in population size on DNA polymorphism. Genetics. 1989b;123:597–601.

    CAS  Article  Google Scholar 

  85. 85.

    Fay JC, Wu C-I. A human population bottleneck can account for the discordance between patterns of mitochondrial versus nuclear DNA variation. Mol Biol Evol. 1999;16:1003–5.

    CAS  Article  PubMed  Google Scholar 

  86. 86.

    Gattepaille L, Jakobsson M, Blum M. Inferring population size changes with sequence and SNP data: lessons from human bottlenecks. Genetics. 2016;110:409–19.

    Article  Google Scholar 

  87. 87.

    Galbreath KE, Hoberg EP. Return to Beringia: parasites reveal cryptic biogeographic history of North American pikas. Proc R Soc B Biol Sci. 2012;279:371–8.

    Article  Google Scholar 

  88. 88.

    Wang IJ, Bradburd GS. Isolation by environment. Mol Ecol. 2014;23:5649–62.

    Article  PubMed  Google Scholar 

  89. 89.

    Morrison SF, Hik DS. Demographic analysis of a declining pika Ochotona collaris population: linking survival to broad-scale climate patterns via spring snowmelt patterns. J Anim Ecol. 2007;76:899–907.

    Article  PubMed  Google Scholar 

  90. 90.

    Jeffress MR, Rodhouse TJ, Ray C, Wolff S, Epps CW. The idiosyncrasies of place: geographic variation in the climate distribution relationships of the American pika. Ecol Appl. 2013;23:864–78.

    Article  PubMed  Google Scholar 

  91. 91.

    Hoban S, Kelley JL, Lotterhos KE, Antolin MF, Bradburd G, Lowry DB, et al. Finding the genomic basis of local adaptation: pitfalls, practical solutions, and future directions. Am Nat. 2016;188:379–97.

    Article  PubMed  PubMed Central  Google Scholar 

  92. 92.

    Lowry DB, Hoban S, Kelley JL, Lotterhos KE, Reed LK, Antolin MF, et al. Breaking RAD: an evaluation of the utility of restriction site-associated DNA sequencing for genome scans of adaptation. Mol Ecol Resour. 2017a;17:142–52.

    CAS  Article  Google Scholar 

  93. 93.

    Lowry DB, Hoban S, Kelley JL, Lotterhos KE, Reed LK, Antolin MF, et al. Responsible RAD: striving for best practices in population genomic studies of adaptation. Mol Ecol Resour. 2017b;17:366–9.

    Article  Google Scholar 

  94. 94.

    Catchen JM, Hohenlohe PA, Bernatchez L, Funk WC, Andrews KR, Allendorf FW. Unbroken: RADseq remains a powerful tool for understanding the genetics of adaptation in natural populations. Mol Ecol Resour. 2017;17:362–5.

    CAS  Article  Google Scholar 

  95. 95.

    McKinney GJ, Larson WA, Seeb LW, Seeb JE. RADseq provides unprecedented insights into molecular ecology and evolutionary genetics: comment on breaking RAD by Lowry et al. (2016). Mol Ecol Resour. 2017;17:356–61.

    CAS  Article  Google Scholar 

  96. 96.

    Laabs BJC, Munroe JS, Best LC, Caffee MW. Timing of the last glaciation and subsequent deglaciation in the Ruby Mountains, Great Basin, USA. Earth Planet Sc Lett. 2013;361:16–25.

    CAS  Article  Google Scholar 

  97. 97.

    Merideth SJ. The impact of habitat spatial structure on pika (Ochotona princeps) dispersal dynamics. Master of Science thesis. University of Nevada, Reno. 2002.

  98. 98.

    Floyd CH, Van Vuren DH, May B. Marmots on Great Basin mountaintops: using genetics to test a biogeographic paradigm. Ecology. 2005;86:2145–53.

    Article  Google Scholar 

  99. 99.

    Kramer AT, Fant JB, Ashley MV. Influences of landscape and pollinators on population genetic structure examples from three Penstemon (Plantaginaceae) species in the Great Basin. Am J Bot. 2011;98:109–21.

    Article  PubMed  Google Scholar 

  100. 100.

    Knowles LL, Massatti R, He Q, Olson LE, Lanier HC. Quantifying the similarity between genes and geography across Alaska’s alpine small mammals. J Biogeogr. 2016;43:1464–76.

    Article  Google Scholar 

  101. 101.

    Nichols LB. Fecal Pellets of American Pikas (Ochotona princeps) provide a crude chronometer for dating patch occupancy. West N Am Nat. 2011;70:500–7.

    Article  Google Scholar 

  102. 102.

    Severaid JH. Natural history of the pikas (Mammalian genus Ochotona). Doctoral dissertation. University of California, Berkeley. 1955.

  103. 103.

    Smith AT. The distribution and dispersal of pikas: consequences of insular population structure. Ecology. 1974;55:1112–9.

    Article  Google Scholar 

  104. 104.

    Hanski I. Habitat connectivity, habitat continuity, and metapopulations in dynamic landscapes. Oikos. 1999;87:209–19.

    Article  Google Scholar 

  105. 105.

    Gilpin ME. The genetic effective population size of a metapopulation. In: Gilpin M, Hanski I, editors. Metapopulation dynamics: empirical and theoretical investigations. London: Academic Press; 1991. p. 165–75.

  106. 106.

    Whitlock MC, Barton NH. The effective size of a subdivided population. Genetics. 1997;146:427–41.

    CAS  Article  Google Scholar 

  107. 107.

    Hedrick PW, Gilpin ME. Genetic effective size of a metapopulation, In: Hanksi I, Gilpin M, editors. Metapopulation biology: ecology, genetics, evolution. San Diego: Academic Press; 1997. p. 165–81.

  108. 108.

    Whitock ME. Selection and drift in metapopulations. In: Hanski I, Gaggiotti OE, editors. Ecology, genetics and evolution of metapopulations. Academic Press; 2004. p. 153–173.

  109. 109.

    Gerlach G, Hoeck HN. Islands on the plains: metapopulation dynamics and female biased dispersal in hyraxes (Hyracoidea) in the Serengeti National Park. Mol Ecol. 2001;10:2307–17.

    CAS  Article  PubMed  Google Scholar 

  110. 110.

    Schmitt T, Seitz A. Low diversity but high differentiation: the population genetics of Aglaope infausta (Zygaenidae: Lepidoptera). J Biogeogr. 2004;31:137–44.

    Article  Google Scholar 

  111. 111.

    Neville HM, Dunham JB, Peacock MM. Landscape attributes and life history variability shape genetic structure of trout populations in a stream network. Landsc Ecol. 2006;21:901–16.

    Article  Google Scholar 

  112. 112.

    Shama LNS, Kubow KB, Jokela J, Robinson CT. Bottlenecks drive temporal and spatial genetic changes in alpine caddisfly metapopulations. BMC Evol Biol. 2011;11:278.

    Article  PubMed  PubMed Central  Google Scholar 

  113. 113.

    Mims MC, Hauser L, Goldberg CS, Olden JD. Genetic differentiation, isolation-by- distance, and metapopulation dynamics of the Arizona treefrog (Hyla wrightorum) in an isolated portion of its range. PLoS ONE. 2016;11:e0160655.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  114. 114.

    Ray C. Maintaining genetic diversity despite local extinctions: effects of population scale. Biol Conserv. 2001;2001(100):3–14.

    Article  Google Scholar 

  115. 115.

    Peacock MM, Ray C. Dispersal in pikas (Ochotona princeps): combining genetic and demographic approaches to reveal spatial and temporal patterns. In: Clobert J, Danchin E, Dhondt AA, Nichols JD, editors. Dispersal. Oxford: Oxford University Press; 2001. p. 43–56.

    Google Scholar 

  116. 116.

    Schwalm D, Epps CW, Rodhouse TJ, Monahan WB, Castillo JA, Ray C, et al. Habitat availability and gene flow influence diverging local population trajectories under scenarios of climate change: a place-based approach. Glob Change Biol. 2016;22:1572–84.

    Article  Google Scholar 

  117. 117.

    Rowe G, Beebee TJC, Burke T. A microsatellite analysis of natterjack toad, Bufo calamita, metapopulations. Oikos. 2000;88:641–51.

    Article  Google Scholar 

  118. 118.

    Antolin MF, Savage LT, Eisen RJ. Landscape features influence genetic structure of black-tailed prairie dogs (Cynomys ludovicianus). Landsc Ecol. 2006;21:867–75.

    Article  Google Scholar 

  119. 119.

    Honnay O, Jacquemyn H, Van Looy K, Vandepitte K, Breyne P. Temporal and spatial genetic variation in a metapopulation of the annual Erysimum cheiranthoides on stony river banks. J Ecol. 2009;97:131–41.

    Article  Google Scholar 

  120. 120.

    Farrington HL, Petren K. A century of genetic change and metapopulation dynamics in the galapagos warbler finches (Certhidea). Evolution. 2011;65:3148–61.

    Article  PubMed  Google Scholar 

  121. 121.

    Nevado B, Mautner S, Sturmbauer C, Verheyen E. Water-level fluctuations and metapopulation dynamics as drivers of genetic diversity in populations of three Tanganyikan cichlid fish species. Mol Ecol. 2013;22:3933–48.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  122. 122.

    Harrison S, Hastings A. Genetic and evolutionary consequences of metapopulation structure. TREE. 1996;11:180–3.

    CAS  Article  PubMed  Google Scholar 

  123. 123.

    Stewart JAE, Wright DH. Assessing persistence of the American pika at historic localities in California’s Northern Sierra Nevada. Wildl Soc B. 2012;36:759–64.

    Article  Google Scholar 

  124. 124.

    Wright DH, Stewart JAE. Within-talus temperatures are not limiting for pikas in the northern Sierra Nevada, California, USA. Calif Fish Game. 2018;104:180–95.

    Google Scholar 

  125. 125.

    Ray C, Beever E, Loarie S. Retreat of the American pika: up the mountain or into the void. In: Brodie JF, Post ES, Doak DF, editors. Wildlife conservation in a changing climate. Chicago: University of Chicago Press; 2012. p. 245–68.

    Google Scholar 

  126. 126.

    Kobiv Y. Response of rare alpine plant species to climate change in the Ukrainian Carpathians. Folia Geobot. 2017;52:217–26.

    Article  Google Scholar 

  127. 127.

    Halloy SRP, Mark AF. Climate-change effects on alpine plant biodiversity: a New Zealand perspective on quantifying the threat. Arct Antarct Alp Res. 2018;35:248–54.[0248:CEOAPB]2.0.CO;2.

    Article  Google Scholar 

  128. 128.

    Moritz C, Patton JL, Conroy CJ, Parra JL, White GC, Beissinger SR. Impact of a century of climate change on small-mammal communities in Yosemite National Park USA. Science. 2008;322:261–4.

    CAS  Article  PubMed  Google Scholar 

  129. 129.

    Rowe RJ, Finarelli JA, Rickart EA. Range dynamics of small mammals along an elevational gradient over an 80-year interval. Glob Change Biol. 2010;16:2930–43.

    Article  Google Scholar 

  130. 130.

    Bhattacharyya S, Dawson DA, Hipperson H, Ishtiaq F. A diet rich in C3 plants reveals the sensitivity of an alpine mammal to climate change. Mol Ecol. 2019;28:250–65.

    CAS  Article  PubMed  Google Scholar 

  131. 131.

    Holzapfel AM, Vinebrooke RD. Environmental warming increases invasion potential of alpine lake communities by imported species. Glob Change Biol. 2005;11:2009–15.

    Article  Google Scholar 

  132. 132.

    Giersch J, Jordan S, Luikart G, Jones LA, Hauer FR, Muhlfeld CC. Climate-induced range contraction of a rare alpine aquatic invertebrate. Freshw Sci. 2015;2015(34):53–65.

    Article  Google Scholar 

  133. 133.

    Smith AT, Ivins BT. Colonization in a pika population: dispersal versus philopatry. Behav Ecol Sociobiol. 1983;13:37–47.

    Article  Google Scholar 

  134. 134.

    Leffler EM, Bullaughey K, Matute DR, Meyer WK, Se’gurel L, Venkat A, et al. Revisiting an old riddle: what determines genetic diversity levels within species? PLoS Biol. 2012;10:e1001388.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  135. 135.

    Sauvajot RM, Buechner M, Kamradt DA, Schonewald CM. Patterns of human disturbance and response by small mammals and birds in chaparral near urban development. Urban Ecosyst. 1998;2:279–97.

    Article  Google Scholar 

  136. 136.

    Smith AB, Beever EA, Kessler AE, Johnston AN, Ray C, Epps CW, et al. Alternatives to genetic affinity as a context for within-species response to climate. Nat Climate Change. 2019;9:787–94.

    Article  Google Scholar 

Download references


We thank two anonymous reviewers who provided valuable comments that we used to improve the manuscript. KBK and MMP give special thanks to BSHP staff for logistical support especially CSP Supervising Rangers J. Heitzmann, T. Peters, and R. Peek. Field work and sample collection at BSHP would not have been possible without K. King and L. Nichols. CR also gives special thanks A. Craighead (Craighead Institute) for help with Montana sample collection. Thanks to J. Kastner for field assistance and help with generating a sample map and V. Kirchoff for laboratory expertise.


This work was supported by a grant from The Eppley Foundation for Research awarded to MMP, which funded the laboratory work and generation of DNA sequence data. Colorado and Montana sampling was funded by awards to CR from the Niwot Ridge Long Term Ecological Research project through NSF Cooperative Agreements DEB-634 1027341 and DEB-1637686. Bodie Foundation Grant, UNR-EECB summer stipend, UNR GSA Travel Grant and a UNR Graduate Research Grant to KBK funded fieldwork in California. TLP was supported by startup funds from the University of Nevada, Reno. The funding bodies played no role in the design of the study and collection, analysis, and interpretation of data and in writing the manuscript.

Author information




KBK and MMP designed the study. KBK, CR, and MMP collected pika samples from the field. KBK and TLP generated GBS libraries for sequencing. JPJ and TLP analyzed the data. All authors contributed to writing, and read and approved the final manuscript.

Corresponding author

Correspondence to Mary M. Peacock.

Ethics declarations

Ethics approval and consent to participate

Institutional Animal Care and Use Committees (IACUC) at the University of Nevada, Reno (protocol No. 00557) and University of Colorado, Boulder (protocol No. 2681), approved research protocols. Scientific collecting permits were obtained from the California Department of Parks and Recreation (No. 2014-03); California Department of Fish and Wildlife (SC-012184); US Department of Agriculture (USDA) Forest Service, Pacific Southwest Region, Harvey Monroe Hall RNA (No. PSW-4000-12, by letter); and USDA Forest Service, Inyo National Forest (No. LVD14021 and No. RMT143); Montana Department of Fish, Wildlife and Parks (2020-038-W); Colorado Division of Wildlife (20TR2014); USDA Forest Service permit (BOU621).

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1: Table S1.

Five replicate entropy (Gompert et al. 2014) chains were run for k = 2–11 for the full dataset of all pika populations. Table S2. Five replicate entropy (Gompert et al. 2014) chains were run for k = 2–8 for the Nevada subset of pika populations. Table S3. The loadings of bioclimatic variables onto the first two PCs for all 11 populations. Table S4. The loadings of bioclimatic variables onto the first two PCs for the six Nevadan populations. Table S5. Sample sizes (N) and estimates of genetic diversity for each population.

Additional file 2: Fig. S1.

Pairwise comparisons of Nei’s D (Nei 1972) based on allele frequencies for each sampled pika population. The distribution of Nei’s D for all pairwise comparisons is represented by a heat map with warmer colors indicating greater genome-wide genetic differentiation. Figure created by authors KBK and JPJ.

Additional file 3: Fig. S2.

The relationships among the four variables (elevational distance, geographic distance, (Climate 1 Distance, Climate 2 distance) used to predict genetic distance (see Fig. 6; Tables 2 and 3 in the main text) are depicted for all populations (All; upper triangle) and the subset of Nevadan populations (NV; lower triangle). Figure created by author JPJ.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Klingler, K.B., Jahner, J.P., Parchman, T.L. et al. Genomic variation in the American pika: signatures of geographic isolation and implications for conservation. BMC Ecol Evo 21, 2 (2021).

Download citation


  • Alpine
  • Climate
  • Conservation
  • Genetic diversity
  • Genotyping-by-sequencing
  • Great basin
  • Metapopulation
  • Ochotona princeps
  • Rocky Mountains
  • Sierra Nevada