Climate change refugia
Monthly climate data, specifically temperature, precipitation, and snow water equivalent (SWE), were derived using monthly PRISM data statistically downscaled to 800 m [17] and then, based on a hydrologic process model, downscaled further to 270 m [18]. The following thresholds were used to define climate change refugia, based on comparing the historical period (1910–1939) to the recent period (1970–1999) [4 for more details]: (1) ≤1 °C changes in temperature (mean annual temperature, annual minimum temperature, annual maximum temperature, and annual mean temperature of the coldest quarter); (2) ≤10% change in relative annual mean precipitation and spring snowpack (yearly estimates of SWE); and (3) ≤2 months per year on average exceeding the extreme maximum (95th percentile) historical temperature and precipitation and the extreme minimum (95th percentile) historical precipitation variation.
Habitat connectivity
Maher et al. [4] generated a set of friction surfaces to represent how landscape features may act as barriers or facilitators to species’ dispersal among Sierra Nevada meadows based upon four hypotheses: (1) isolation by distance; (2) isolation by topography, tested with a raster layer of distance weighted by elevation (upslope or downslope) between a layer of delineated meadows [19] with no dispersal permitted if the slope was greater than 45°; (3) isolation by watercourses in two forms: a) rivers and streams weighted with a high resistance value (100) compared to the surrounding matrix (“river_pres”), and b) distance from rivers and streams weighted as an increasing barrier to dispersal (“river_dist”); and (4) isolation by roads, considering distance from primary or secondary roads, which also acted as a proxy for human activity and presences. From these friction surfaces, maps of potential connectivity between meadow polygons buffered by 150 m were created using Circuitscape (ver 3.5.8; [20]). The resulting maps were continuous values of flow through a given pixel and represent the expected frequency of movement across a pixel. The approach was agnostic to species, such that a breadth of surfaces were examined, some of which are not applicable to the natural history of U. beldingi. In addition to these surfaces, we use two friction surfaces representing spatial variation in climatic water deficit (CWD) in the historical (“histCWD”) and modern (“modCWD”) periods, which allows a test of whether less variation in moisture availability facilitated dispersal.
Estimating site occupancy and persistence
A survey of U. beldingi occurrence in 38 meadows in Yosemite National Park was conducted in the summer of 2011 (Fig. 1, inset). Observers walked meadow transects while visually scanning for U. beldingi and listening for alarm calls; a site was considered occupied if a U. beldingi was seen or heard anywhere in the meadow. Meadows were buffered by 2 km using the R package rgeos [21] and mean climate change values were extracted from rasters using the package raster [22]. We used a series of two-sample Wilcoxon tests to determine whether climate change refugia (as defined by the thresholds described above) and connectivity predicted increased probability of occupancy.
Next, as a measure of persistence, we examined results from resurveys (conducted 2003–2011) in or near 74 short-grass sites (either human-modified (Fig. 2a-b) or natural meadows (Fig. 2c)) in the Sierra Nevada, southern Cascades, and Modoc Plateau of California, USA (Fig. 1), where U. beldingi were historically (1902 to 1966) observed and/or trapped (based on detailed field notes from the University of California Museum of Vertebrate Zoology (MVZ)). As reported elsewhere [16], U. beldingi detectability per site averaged 0.96 and they were found to have been extirpated from 42% of the resurvey sites. To test whether climate change refugia were related to site persistence over the twentieth century, we used binomial tests based on the mean value of select climatic variables for each site with a 2 km radius around its centerpoint. Campgrounds, municipal parks, and agricultural fields were grouped as human-modified sites (the effect of connectivity on persistence of populations in them was analyzed separately).
There is a temporal mismatch between the modern climate window that Maher et al. 2017 used to define meadow refugia and the Belding’s ground squirrel resurvey data we used to test them. We believe this is justified because the measures of change in conditions need to precede the survey window in order to measure a potential impact of those changes, especially when trying to measure long-term persistence or extirpation that may show lag effects. If our climate-related variables shifted in magnitude and/or direction after 1999 and prior to the resurvey, then we would expect to find a poor relationship between the environmental data and resurvey data, assuming that climate-related effects are the predominant factor in persistence. Thus at the very least, datasets from Maher et al. (2017) represent a conservative estimate of climate changes relative to the resurvey.
Analysis of genetic diversity and gene flow
We sampled 244 U. beldingi trapped at 15 sites spread across its California range (Fig. 2) from 2003 to 2011 using Sherman or Tomahawk Live Traps, the majority in 2010 and 2011. Genetic samples were collected from 187 adults from an ear snip prior to release or from a liver biopsy in the case of MVZ voucher specimens (archived samples described in Additional file 1: Table S1), extracted using a salt extraction protocol [23], then stored at 4 °C until amplification.
In brief (see Additional file 2: Supplemental Methods for more details), 14 polymorphic microsatellite loci were amplified by PCR using fluorescently labelled primers and conditions optimized from protocols developed for other sciurids (see Additional file 3: Table S2). After adjusting for null alleles and linkage disequilibrium, we calculated the fixation index FST [24, 25], where lower values indicate higher gene flow, to test predictions of meadow connectivity and sample size-corrected allelic richness (Ar) to test the effects of both connectivity and refugia. We focused on allelic richness as a measure of genetic diversity; in both theory and practice, allelic richness is more sensitive to population reduction than is heterozygosity [26, 27].
structure 2.3.4 [28] was used to examine genetic structuring and identify distinct genetic populations based on the 15 sampled sites (Fig. 2; Additional file 4: Table S3). Each run was replicated ten times with K set to run from 1 to 8 populations with an initial burn-in of 100,000 iterations followed by 106 Markov chain Monte Carlo repetitions. We summarized across runs using Structure Harvester [29], the statistic ΔK to choose the best K [30], Clumpp 1.1.2 [31] to cluster results from structure run repetitions, and Distruct 1.1 [32] to visualize population structure.
We examined the pattern of genetic diversity for 11 populations (n = 124) of U. beldingi in the central Sierra Nevada including Yosemite National Park (hereafter “Central Sierra”, Fig. 2) and the western Great Basin using linear regressions to test for correlations between allelic richness and minimal change in climate factors, as well as mean minimum temperature in the recent period.
To test correlations with the connectivity predictions, we calculated the minimum Euclidean distance between populations using the package sp. [33] in R (for isolation by distance) and determined the effective resistances among populations using Circuitscape (for the remaining connectivity hypotheses). We used the Central Sierra populations as focal nodes and the friction surfaces as resistance (topography, presence of watercourses, environmental variation) and conductance (distance from watercourses, distance from roads). The Multiple Matrix Regression with Randomization method [34] and backward selection were used to test whether a combination of distance matrices could explain the FST patterns. The advantage of this approach, as opposed to Mantel tests, is that it allowed the comparison of several distance matrices simultaneously, as well as the amount of variation explained by the given model. We compared connectivity models to patterns in FST among the Central Sierra sites. We also used linear regressions to test for correlations between mean connectivity and allelic richness.
Finally, BayesAss [35] was used to test for directional gene flow between central Sierra populations. It estimates the recent proportion of migrants among proximal populations and thus gives a different, and more precise, measure of gene flow than FST [36]. The number of iterations was 108, with 106 as burn-in period, and a sampling frequency of 2000. Three randomly seeded replicates were run with Δp (allele freq) of 0.80, Δm (mig rate) of 0.80, and ΔF (inbreeding) of 0.80.
Principal component analysis
To summarize overall trends in climate change across the suite of environmental variables and since univariate refugia may be unlikely, we also conducted a principal components analysis (PCA) with data across all of our sites (n = 102). Because variance was unequal among variables, we used a correlation matrix, centering and normalizing all variables before calculations. Eigenvalues of the first three components of the PCA were greater than 1 and together explained 85.4% of the variation in climate change (Additional file 5: Table S4 and Additional file 6: Table S5), with sites separated out in general by increases in precipitation (PC1), increases in SWE vs. warming conditions (PC2), and wetter conditions with more months of extreme maximum temperatures (positive) vs. increases in minimum temperature (PC3).
For these first three components, we compared how well the summarized change explained patterns in occupancy, persistence, and genetic diversity. For occupancy and persistence, we used a series of logistic regressions to test the ability of each set of component scores to explain our observations. For genetic diversity, we compared allelic richness to each component score using linear models.