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 . 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–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–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–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 in structuring range-wide genetic resources.
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–17]. This small lagomorph is found in the mountainous regions of western North America where it lives primarily on talus slopes above timberline . 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 , morphological  and behavioral data [23, 25, 26], Hafner and Smith  suggested collapsing the previously recognized 36 subspecies of O. princeps  into five subspecies. These designations are congruent with mitochondrial lineage designations  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 . 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 8,000–10,000 years, pika populations were lost from the relatively low elevation ranges of the Great Basin, which also lack sufficient talus habitat to provide refuge in a warming climate [29, 30]. Climate change in the 20th and 21st centuries has further accelerated population declines and losses primarily from the Sierra Nevada lineage, O. p. schisticeps [17, 31–33]. This is especially true for the mountain ranges of the Great Basin, where resurveys of sites occupied in the mid-20th 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–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–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 . 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.  generated genotyping-by-sequencing (GBS) data to identify genomic regions responding to selection across four pika populations found along an elevational gradient for 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 . In an additional study based on similar data, Waterhouse et al.  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. However, we largely lack an understanding of how genetic structure and diversity vary 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 (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 a fragmented habitat [35, 40, 41, 45].
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. 1). 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 . 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 . 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) , 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–50].
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) across 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 expect lower levels of genetic variation in the peripheral and isolated populations of the Bodie Hills and the Ruby-East Humboldt Mountains. Finally, we expected 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.