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

Abstract

Background

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.

Results

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.

Conclusions

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.

Background

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

Population	N	Elevation (m)	Latitude	Longitude	Sampling years	Subspecies
Bodie, CA	15	2500	38.211861	− 119.005037	2014–2015	O. p. schisticeps
Pipet Tarn, CA	20	3170	37.945285	− 119.28413	2014–2015	O. p. schisticeps
West Knoll, CO	21	3612	40.056753	− 105.59571	2013–2015	O. p. saxatilis
Emerald Lake, MT	21	2896	45.407061	− 110.93998	2011–2015	O. p. princeps
Swan Creek, MT	7	1829	45.374392	− 111.14408	2013	O. p. princeps
Overland Lake, NV (RM)	17	2281	40.458067	− 115.45524	1999–2000	O. p. princeps
Island Lake, NV (RM)	16	3245	40.616662	− 115.37916	1999–2000	O. p. princeps
Hidden Lake, NV (RM)	15	3051	40.746363	− 115.2845	1999–2000	O. p. princeps
Week’s Creek, NV (EH)	10	2825	40.928332	− 115.11273	1999–2000	O. p. princeps
Lizzie’s Basin, NV (EH)	10	2803	40.943973	− 115.11243	1999–2000	O. p. princeps
Smith Lake, NV (EH)	19	2882	41.034497	− 115.094786	1999–2000	O. p. princeps

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

Photographs provide by authors KBK (Bodie and Pipet Tarn) and CR (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 km² 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.

Figure created by author KBK

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.

Results

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.

Figure created by author JPJ

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.

Figure created by author JPJ

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.

Figure created by author JPJ

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 (R² = 0.344; P = 0.010; Additional file 3: Fig. S2). Nei’s D was strongly predicted by Climate 2 distance (R² = 0.538; P = 0.009) and geographic distance (R² = 0.342; P = 0.009), but was not influenced by Climate 1 distance (R² = 0.129; P = 0.079) or elevational distance (R² = 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).

Model	R2	P
Elevation + Geography + Climate 1 + Climate 2	0.716	0.005
Geography + Climate 1 + Climate 2	0.708	0.003
Elevation + Geography + Climate 2	0.695	0.002
Geography + Climate 2	0.678	0.002
Elevation + Climate 1 + Climate 2	0.673	0.001
Climate 1 + Climate 2	0.673	0.001
Elevation + Climate 2	0.538	0.005
Climate 2	0.538	0.002
Elevation + Geography + Climate 1	0.385	0.057
Elevation + Geography	0.384	0.038
Geography + Climate 1	0.342	0.026
Geography	0.342	0.009
Elevation + Climate 1	0.130	0.184
Climate 1	0.129	0.079
Elevation	0.000	0.960

1. The two climate variables correspond to the first two principal components from a PCA on 19 climate variables from BioClim, which explained 46.9 and 33.8 percent of the variance in the data, respectively (see Table S4 for loadings and Fig. S2 for relationships among predictor variables)

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 https://www.worldclim.org; 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.

Figure created by author JPJ

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 (R² = 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 (R² = 0.596; P = 0.024), but was not associated with Climate 1 distance (R² = 0.020; P = 0.399), Climate 2 distance (R² = 0.094; P = 0.130), or elevational distance (R² = 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).

Model	R2	P
Elevation + Geography + Climate 1 + Climate 2	0.727	0.024
Elevation + Geography + Climate 1	0.722	0.024
Elevation + Geography + Climate 2	0.695	0.034
Elevation + Geography	0.686	0.025
Geography + Climate 1 + Climate 2	0.685	0.013
Geography + Climate 1	0.670	0.017
Geography + Climate 2	0.628	0.021
Geography	0.596	0.024
Elevation + Climate 1 + Climate 2	0.149	0.243
Elevation + Climate 2	0.126	0.144
Climate 1 + Climate 2	0.102	0.227
Climate 2	0.094	0.130
Elevation + Climate 1	0.055	0.242
Climate 1	0.020	0.399
Elevation	0.016	0.284

1. The two climate variables correspond to the first two principal components from a PCA on 19 climate variables from BioClim, which explained 60.0 and 31.0 percent of the variance in the data, respectively (see Additional file 1: Table S5 for loadings and Additional file 3: Fig. S2 for relationships among predictor variables)

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.

Photo provided by author KBK. Figure created by author JPJ

Discussion

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.

Conclusions

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.

Methods

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 (https://github.com/ncgr/tapioca). 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 km² scale; https://www.worldclim.org) [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 https://doi.org/10.5061/dryad.mcvdncjww.

Abbreviations

ddRADseq: Double digest restriction-site associated DNA sequencing

DIC: Deviance information criterion

GBS: Genotyping-by-sequencing

MRM: Multiple regression on distance matrices

PC: Principal component

PCA: Principal component analysis

SNP: Single nucleotide polymorphism

References

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. https://doi.org/10.1111/j.1365-2486.2012.02729.x.

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. https://doi.org/10.1111/j.1365-2664.2007.01325.x.

3. Ashcroft MB. Identifying refugia from climate change. J Biogeogr. 2010;37:1407–13. https://doi.org/10.1111/j.1365-2699.2010.02300.x.

4. Varner J, Dearing MD. The importance of biologically relevant microclimates in habitat suitability assessments. PLoS ONE. 2014;9:e104648. https://doi.org/10.1371/journal.pone.0104648.

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. https://doi.org/10.1016/j.biocon.2016.09.026.

6. Hampe A, Jump AS. Climate relicts: past, present, future. Annu Rev Ecol Evol Syst. 2011;42:313–33. https://doi.org/10.1146/annurev-ecolsys-102710-145015.

7. Scherrer D, Körner C. Topographically controlled thermal-habitat differentiation buffers alpine plant diversity against climate warming. J Biogeogr. 2011;38:406–16. https://doi.org/10.1111/j.1365-2699.2010.02407.x.

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. https://doi.org/10.1111/brv.12105.

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. https://doi.org/10.1002/ece3.2763.

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. https://doi.org/10.1111/j.1600-0587.2010.06430.x.

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. https://doi.org/10.1644/10-MAMM-A-380.1.

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. https://doi.org/10.1038/NCLIMATE1415.

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. https://doi.org/10.1644/1545-1542(2003)084%3c0037:POAEAI%3e2.0.CO;2.

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. https://doi.org/10.1111/j.1365-2486.2010.02389.x.

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. https://doi.org/10.1016/j.quaint.2010.05.004.

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. https://doi.org/10.1111/j.1600-0587.2011.07227.x.

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. https://doi.org/10.1111/jbi.12466.

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

19. MacArthur RA, Wang LCH. Physiology of thermoregulation in the pika Ochotona princeps. Can J Zool. 1973;51:11–6. https://doi.org/10.1139/z73-002.

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. https://doi.org/10.1371/journal.pone.0119327.

21. Hafner DJ, Smith AT. Revision of the subspecies of the American pika, Ochotona princeps (Lagomorpha: Ochotonidae). J Mammal. 2010;91:401–17. https://doi.org/10.1644/09-MAMM-A-277.1.

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. https://doi.org/10.1111/j.1558-5646.2009.00803.x.

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. https://doi.org/10.1111/j.1365-2699.2009.02201.x.

24. Hafner DJ, Sullivan RM. Historical and ecological biogeography of neararctic pikas (Lagomorpha: Ochotonidae). J Mammal. 1995;76:302–21. https://doi.org/10.2307/1382343.

25. Somers P. Dialects in southern Rocky-Mountain pikas, Ochotona-princeps (Lagomorpha). Anim Behav. 1973;21:124–37. https://doi.org/10.1016/S0003-3472(73)80050-8.

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. https://doi.org/10.1139/z83-054.

27. Mead JI. Quaternary records of pika, Ochotona, in North America. Boreas. 1987;16:165–71. https://doi.org/10.1111/j.1502-3885.1987.tb00768.x.

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. https://doi.org/10.1080/00040851.1994.12003082.

29. Grayson DK. A brief history of Great Basin pikas. J Biogeogr. 2005;32:2103–11. https://doi.org/10.1111/j.1365-2699.2005.01341.x.

30. Grayson DK. The Late Quaternary biogeographic histories of some Great Basin mammals (western USA). Quaternary Sci Rev. 2006;25:2964–91. https://doi.org/10.1016/j.quascirev.2006.03.004.

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. https://doi.org/10.1111/conl.12190.

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. https://doi.org/10.3398/064.076.0203.

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. https://doi.org/10.1371/journal.pone.0181834.

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. https://doi.org/10.1371/journal.pone.0199032.

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. https://doi.org/10.1007/s004420050341.

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. https://doi.org/10.3955/046.086.0410.

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. https://doi.org/10.1093/jmammal/gyv191.

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. 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. https://doi.org/10.1890/15-1452.1.

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. https://doi.org/10.1016/B978-0-12-323445-2.X5000-7.

41. Moilanen A, Smith AT, Hanski IA. Long-term dynamics in a metapopulation of the American pika. Am Nat. 1998;152:530–42. https://doi.org/10.1890/0012-9658(1998)079[2503:MDEOHQ]2.0.CO;2.

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. https://doi.org/10.7717/perej.1106.

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. https://doi.org/10.1111/mec.14701.

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. https://doi.org/10.1093/beheco/8.3.340.

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. https://doi.org/10.1086/338990.

46. Beever EA, Ray C, Mote PW, Wilkening JL. Testing alternative models of climate-mediated extirpations. Eco Appl. 2010;20:164–78. https://doi.org/10.1890/08-1011.1.

47. McDonald KA, Brown JH. Using montane mammals to model extinctions due to global change. Conserv Biol. 1992;6:409–15. https://doi.org/10.1046/j.1523-1739.1992.06030409.x.

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. https://doi.org/10.2307/1381011.

49. Dearing MD. Disparate determinants of summer and winter diet selection of a generalist herbivore Ochotona princeps. Oecologia. 1996;108:467–78. https://doi.org/10.1007/BF00333723.

50. Wilkening JL, Ray C. Characterizing predictors of survival in the American pika (Ochotona princeps). J Mammal. 2016;97:1366–75. https://doi.org/10.1093/jmammal/gyw097.

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. https://doi.org/10.1111/j.1365-294X.2012.05513.x.

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. https://doi.org/10.1371/journal.pone.0037135.

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. https://doi.org/10.1644/10-MAMM-F-355.1.

54. Langmead B, Salzberg SL. Fast gapped-read alignment with Bowtie 2. Nat Methods. 2012;9:357–9. https://doi.org/10.1038/NMETH.1923.

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. https://doi.org/10.5061/dryad.mcvdncjww.

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. https://doi.org/10.1093/jhered/esw010.

57. Li H, Durbin R. Fast and accurate short read alignment with Burrows–Wheeler transform. Bioinformatics. 2009;25:1754–60. https://doi.org/10.1093/bioinformatics/btp324.

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. https://doi.org/10.1093/bioinformatics/btp352.

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. https://doi.org/10.1093/bioinformatics/btr330.

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. https://doi.org/10.1111/mec.12811.

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

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.

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. https://doi.org/10.1038/nrg2986.

64. Buerkle A, Gompert Z. Population genomics based on low coverage sequencing: how low should we go? Mol Ecol. 2013;22:3028–35. https://doi.org/10.1111/mec.12105.

65. Venables WN, Ripley BD. Random and mixed effects. In: Modern applied statistics with S. Statistics and computing. Springer; 2002. p. 271–300. https://doi.org/10.1007/978-0-387-21706-2_10.

66. R Core Team. R: a language and environment for statistical computing. Vienna: R Foundation for Statistical Computing; 2020. https://www.R-project.org/.

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. https://doi.org/10.1038/nature14230.

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. https://doi.org/10.1186/s12862-016-0702-4.

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. https://doi.org/10.1111/mec.13486.

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. https://doi.org/10.1038/s41598-019-43939-4.

71. McVean G. A genealogical interpretation of principal components analysis. PLoS Genet. 2009;5(10):e1000686. https://doi.org/10.1371/journal.pgen.10006869.

72. Nei M. Genetic distance between populations. Am Nat. 1972;106:283–92. https://doi.org/10.1086/282771.

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

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. https://doi.org/10.1002/joc.5086.

75. Goslee SC, Urban DL. The ecodist package for dissimilarity-based analysis of ecological data. J Stat Softw. 2007;22:1–19. https://doi.org/10.18637/jss.v022.i07.

76. Lichstein JW. Multiple regression on distance matrices: a multivariate spatial analysis tool. Plant Ecol. 2007;188:117–31. https://doi.org/10.1007/s11258-006-9126-3.

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

78. Watterson GA. On the number of segregating sites in genetic models without recombination. Theor Popul Biol. 1975;7:256–76. https://doi.org/10.1016/0040-5809(75)90020-9.

79. Korneliussen TS, Albrechtsen A, Korneliussen RN. ANGSD: analysis of next generation sequencing data. BMC Bioinform. 2014;15:2–13. https://doi.org/10.1186/s12859-014-0356-4.

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. https://doi.org/10.1186/1471-2105-14-289.

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

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

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

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

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. https://doi.org/10.1093/oxfordjournals.molbev.a026175.

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. https://doi.org/10.1534/genetics.115.185058.

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. https://doi.org/10.1098/rspb.2011.0482.

88. Wang IJ, Bradburd GS. Isolation by environment. Mol Ecol. 2014;23:5649–62. https://doi.org/10.1111/mec.12938.

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. https://doi.org/10.1111/j.1365-2656.2007.01276.x.

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. https://doi.org/10.1890/12-0979.1.

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. https://doi.org/10.1086/688018.

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.

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.

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.

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.

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. https://doi.org/10.1016/j.epsl.2012.11.018.

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. Floyd CH, Van Vuren DH, May B. Marmots on Great Basin mountaintops: using genetics to test a biogeographic paradigm. Ecology. 2005;86:2145–53. https://doi.org/10.1890/04-1227.

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. https://doi.org/10.3732/ajb.1000229.

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. https://doi.org/10.1111/jbi.12728.

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. https://doi.org/10.3398/064.070.0410.

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

103. Smith AT. The distribution and dispersal of pikas: consequences of insular population structure. Ecology. 1974;55:1112–9. https://doi.org/10.2307/1935464.

104. Hanski I. Habitat connectivity, habitat continuity, and metapopulations in dynamic landscapes. Oikos. 1999;87:209–19. https://doi.org/10.2307/3546736.

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. https://doi.org/10.1016/B978-0-12-284120-0.X5001-3.

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

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. https://doi.org/10.1016/B978-012323445-2/50011-0.

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. https://doi.org/10.1016/B978-012323448-3/50009-X.

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. https://doi.org/10.1046/j.0962-1083.2001.01369.x.

110. Schmitt T, Seitz A. Low diversity but high differentiation: the population genetics of Aglaope infausta (Zygaenidae: Lepidoptera). J Biogeogr. 2004;31:137–44. https://doi.org/10.1046/j.0305-0270.2003.01003.x.

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. https://doi.org/10.1007/s10980-005-5221-4.

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. https://doi.org/10.1186/1471-2148-11-278.

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. https://doi.org/10.1371/journal.pone.0160655.

114. Ray C. Maintaining genetic diversity despite local extinctions: effects of population scale. Biol Conserv. 2001;2001(100):3–14. https://doi.org/10.1016/S0006-3207(00)00202-0.

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.

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. https://doi.org/10.1111/gcb.13189.

117. Rowe G, Beebee TJC, Burke T. A microsatellite analysis of natterjack toad, Bufo calamita, metapopulations. Oikos. 2000;88:641–51. https://doi.org/10.1034/j.1600-0706.2000.880321.x.

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. https://doi.org/10.1007/s10980-005-5220-5.

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. https://doi.org/10.1111/j.1365-2745.2008.01452.x.

120. Farrington HL, Petren K. A century of genetic change and metapopulation dynamics in the galapagos warbler finches (Certhidea). Evolution. 2011;65:3148–61. https://doi.org/10.1111/j.1558-5646.2011.01385.x.

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. https://doi.org/10.1111/mec.12374.

122. Harrison S, Hastings A. Genetic and evolutionary consequences of metapopulation structure. TREE. 1996;11:180–3. https://doi.org/10.1016/0169-5347(96)20008-4.

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. https://doi.org/10.1002/wsb.220.

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.

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.

126. Kobiv Y. Response of rare alpine plant species to climate change in the Ukrainian Carpathians. Folia Geobot. 2017;52:217–26. https://doi.org/10.1007/s12224-016-9270-z.

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. https://doi.org/10.1657/1523-0430(2003)035[0248:CEOAPB]2.0.CO;2.

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. https://doi.org/10.1126/science.1163428.

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. https://doi.org/10.1111/j.1365-2486.2009.02150.x.

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. https://doi.org/10.1111/mec.14842.

131. Holzapfel AM, Vinebrooke RD. Environmental warming increases invasion potential of alpine lake communities by imported species. Glob Change Biol. 2005;11:2009–15. https://doi.org/10.1111/j.1365-2486.2005.001057.x.

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. https://doi.org/10.1086/679490.

133. Smith AT, Ivins BT. Colonization in a pika population: dispersal versus philopatry. Behav Ecol Sociobiol. 1983;13:37–47. https://doi.org/10.1007/BF00299675.

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. https://doi.org/10.1371/journal.pbio.1001388.

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. https://doi.org/10.1023/A:1009588723665.

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. https://doi.org/10.1038/s41558-019-0584-8.