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 (Tables 2, S3; Fig. 6). Importantly, geographic and Climate 2 distance were not correlated with one another across all populations (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–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 (Figs. 3B and 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; 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 19th 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-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 19th 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–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 20th and early 21st centuries [13–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-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 (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; 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.