Climate change refugia and habitat connectivity promote species persistence

Abstract

Background

Climate change refugia, areas buffered from climate change relative to their surroundings, are of increasing interest as natural resource managers seek to prioritize climate adaptation actions. However, evidence that refugia buffer the effects of anthropogenic climate change is largely missing.

Methods

Focusing on the climate-sensitive Belding’s ground squirrel (Urocitellus beldingi), we predicted that highly connected Sierra Nevada meadows that had warmed less or shown less precipitation change over the last century would have greater population persistence, as measured by short-term occupancy, fewer extirpations over the twentieth century, and long-term persistence measured through genetic diversity.

Results

Across California, U. beldingi were more likely to persist over the last century in meadows with high connectivity that were defined as refugial based on a suite of temperature and precipitation factors. In Yosemite National Park, highly connected refugial meadows were more likely to be occupied by U. beldingi. More broadly, populations inhabiting Sierra Nevada meadows with colder mean winter temperatures had higher values of allelic richness at microsatellite loci, consistent with higher population persistence in temperature-buffered sites. Furthermore, both allelic richness and gene flow were higher in meadows that had higher landscape connectivity, indicating the importance of metapopulation processes. Conversely, anthropogenic refugia, sites where populations appeared to persist due to food or water supplementation, had lower connectivity, genetic diversity, and gene flow, and thus might act as ecological traps. This study provides evidence that validates the climate change refugia concept in a contemporary context and illustrates how to integrate field observations and genetic analyses to test the effectiveness of climate change refugia and connectivity.

Conclusions

Climate change refugia will be important for conserving populations as well as genetic diversity and evolutionary potential. Our study shows that in-depth modeling paired with rigorous fieldwork can identify functioning climate change refugia for conservation.

Background

Attention on climate change refugia, areas relatively buffered from contemporary climate change over time that enable persistence of valued physical, ecological, and socio-cultural resources [1], has increased as scientists attempt to explain anomalies in the expected patterns of population responses to climate change [2,3,4,5]. Although managing climate change refugia is frequently identified as a key climate change adaptation strategy [1, 6], spatially explicit predictions of refugia and their effects on biodiversity are rare [5, 7, 8]. Most studies skip the key step of using independent data to evaluate whether these hypothesized refugia increase species persistence (Although see [9, 10]), which is necessary to confirm that investing limited resources into refugia is warranted [1].

We test whether highly connected montane meadows that have experienced little recent climate change [4] increase persistence of an associated species and thus act as refugia in the face of anthropogenic climate change. Montane meadows are key to alpine and subalpine animal communities, hydrological function, and human recreation, agriculture, and economy [11,12,13]. Recent work [4] used high resolution, dynamically downscaled climate models to map refugial meadow networks to identify where climate-sensitive species could persist, and restoration and conservation actions could be focused. Here we use data on site occupancy, site persistence, and population persistence as indicated by genetic diversity for the Belding’s ground squirrel. Urocitellus beldingi is a 250 g social, semi-fossorial mammal that feeds primarily on grasses and is commonly found in mid to high elevation meadows in western North America [14, 15]. Although ecologically important as key prey items and ecosystem engineers, U. beldingi in California have recently undergone large climate-related declines in occupancy, except where “anthropogenic refugia” were created by supplements with water or food through high human activity (such as campgrounds and municipal parks) [16]. We predicted that populations located in highly connected refugial meadows identified by Maher et al. [4] on the basis of temperature and precipitation changes over the twentieth century would have higher contemporary occupancy, persistence, and, as a result, higher genetic diversity than sites that were not identified as climate change refugia. Our results validate the use of climate refugia as a key aspect of change climate adaptation.

Methods

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.

Fig. 1

Urocitellus beldingi (pictured in inset) study sites mapped across California. Circles in main image indicate resurvey sites, with green indicating persistence and black indicating extirpation since historical surveys. Hatching represents the approximate range of U. beldingi (IUCN 2014). Blue and black circles in inset represent recent one-time surveys in Yosemite National Park (outlined in bold), with blue indicating presence and gray indicating absence

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

Fig. 2

Sites where Urocitellus beldingi were sampled genetically. Colors in pie charts show results of structure run, separating 15 California sites into four populations adjusted by sample size. Analyses primarily focused on the central and southern populations. Images indicate (a) an example of an agricultural “anthropogenic refugium” in northern California, where Belding’s ground squirrels are persisting and acting as crop pests despite extirpation at surrounding natural sites; (b) Mono Lake County Park, an example of another “anthropogenic refugium” (note the sprinkler in the foreground) where Belding’s ground squirrels are persisting despite extirpation from nearby historical sites, which are currently dominated by sagebrush; and (c) a natural meadow site in Yosemite National Park (bold line)

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 F_ST [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 10⁶ 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 F_ST 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 F_ST 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 F_ST [36]. The number of iterations was 10⁸, with 10⁶ 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.

Results

Occupancy

The 20 meadows occupied by U. beldingi in Yosemite National Park in 2011 experienced less change over the last century than the 18 unoccupied meadows based on the following five climate change indices (Fig. 3, Additional file 7: Table S6): change in annual minimum temperature (occupied meadow average = +1.30 °C, unoccupied = +1.66 °C; P < 0.001), change in mean annual precipitation (occupied = +8.38 mm, unoccupied = +50.06 mm; P < 0.001), number of months warmer than the historical extreme minimum temperature (occupied = 39.87, unoccupied = 62.94; P < 0.001), number of months wetter than the historical precipitation extreme (occupied = 25.45, unoccupied = 27.74; P = 0.007), and change in April 1st SWE (occupied = +11.68, unoccupied = +69.12; P < 0.001). According to these indices, U. beldingi were more likely to be found in climate change refugia. Occupied sites within Yosemite National Park also had higher values of connectivity than unoccupied sites in relation to specific predicted barriers (topography: P < 0.001; rivers: P _{river_pres} = 0.007, W = 250, P _{river_dist} = 0.021; roads: P _road < 0.001; and climate-water deficit: W_histCWD = 310, P < 0.001, W_modCWD = 282, P = 0.001).

Fig. 3

U. beldingi presence or absence in Yosemite National Park compared to changes in climate: annual minimum temperature (°C) and mean annual precipitation (mm), number of months warmer than the historical minimum temperature extremes, number of months wetter than historical precipitation extremes, change in SWE (mm), and PC1 (primarily representing increases in precipitation; in white) and PC3 (wetter conditions with more months of extreme maximum temperatures (positive) vs. increases in minimum temperature; in gray) scores

In addition, logistic regression for individual principal components of the PCA, which accounted for correlation among climate variables, strongly explained patterns in occupancy (all P < 0.05; estimated PC1 = 0.9271, estimated PC2 = −3.6038, estimated PC3 = 2.1670). Considering all components using a multiple logistic regression followed by stepwise removal leads to essentially equal support of patterns representing PC1 and PC3 (estimated PC1 = 0.6676, P = 0.054 and estimated PC3 = 1.5326, P = 0.055 when considered together). Logistic regressions of individual components indicated higher occupancy at cooler, drier sites (P < 0.05 in each case; estimated PC1 = 0.9271, estimated PC2 = −3.6038, estimated PC3 = 2.1670; Fig. 3, Additional File 8: Figure S1).

Persistence

U. beldingi populations in California were also more likely to persist throughout the twentieth century in sites that were defined as climate change refugia from a suite of temperature and precipitation variables. Specifically, binomial tests showed that persistence was significantly correlated with refugial status (P < 0.05) for 7 of 8 climatic thresholds, with the exception of change in SWE (Fig. 4). Based on logistic regression, only PC1 predicted persistence (est. = 0.3598, P < 0.01).

Fig. 4

U. beldingi persistence in climate change refugia as defined by 8 climate thresholds: changes in mean annual temperature and mean temperature of the coldest quarter no greater than 1°C; relative annual mean precipitation changes within 10%; and no more than 2 months per year on average exceeding the extreme historical temperature and precipitation variation; and relative annual mean precipitation and spring snowpack (SWE) changes within 10%. Significance indicated at P < 0.05 (*), P < 0.01 (**), and P < 0.001 (***)

Connectivity did not predict persistence of U. beldingi across all sampled sites over the last century in the Central Sierra, but did affect persistence when the 21 anthropogenic refugia were excluded from the analysis (W_topo = 61, P < 0.0001). As with occupancy, each connectivity layer contributed to this effect (W_{river_pres} = 65, P < 0.0001; W_{river_dist} = 63, P < 0.0001; W_{road_road} = 64, P < 0.0001; W_histCWD = 56, P < 0.0001; W_modCWD = 51, P < 0.0001; Fig. 5). In contrast, populations that persisted in anthropogenic refugia had the lowest connectivity values.

Fig. 5

U. beldingi population extirpation and persistence over the last century in natural vs. human-modified short-grass sites (“Human Mod.”, e.g., municipal parks and agricultural fields) according to a number of connectivity metrics. Connectivity measures were based upon friction surfaces representing the elevationally informed distance between meadows (a), distance from roads (b), presence of watercourses (c), and the distance from watercourses (d)

Genetic diversity

California populations of U. beldingi fell into four distinct genetic clusters, although the distinction between the northern and central sites could be due to lack of samples from intermediate populations. Mean allelic richness when adjusted for sample size was 3.3 ± 0.3 alleles per locus for the Central Sierra sites (range 2.28–3.85; Additional file 7: Table S3). Most sites where there were enough U. beldingi to trap for a robust genetic analysis ranked as refugia; thus a comparison of the effect of refugial status on genetic diversity was not possible for most climate indices. However, Central Sierra sites with lower current mean minimum temperature had higher values of allelic richness (β = −0.24, t₉ = −4.48, R² _adj = 0.66, P = 0.002; Fig. 6) at temperature-buffered sites. As predicted, higher allelic richness of Central Sierra meadows was also correlated with higher meadow connectivity (β = 0.48, t₉ = 2.92, R² _adj = 0.43, P = 0.02; Fig. 6), but a multiple regression suggested that mean minimum temperature contributed more strongly (β_a = −0.28, t₇ = −4.83, P = 0.002, β_Conn = 2.3e-06, t₇ = 1.52, P > 0.1, β_Int = 4.0e-07, t₇ = 1.62; R² _adj = 0.82, F_3,7; = 10.39, P = 0.006). There was no significant effect of either meadow area or change in environmental conditions (i.e., environmental change PC axes) on allelic richness per site. Meadow area also was not predictive of persistence at these sites (Morelli et al. 2012).

Fig. 6

The log₁₀ of connectivity (averaged for topography, historic and modern CWD, presence of a river barrier and distance to a river, and presence of roads) and effect of minimum temperature in the modern time period (1970–1999) among the meadows sampled across the central Sierra on genetic diversity as measured by allelic richness. The best-fit line with confidence intervals are shown in light gray; prediction intervals are shown in dark grey. Each of the sample sites is plotted with error bars representing ±1 SE

Gene flow (INTENDED AS A HEADER) across the entire region, pairwise F_ST values (mean = 0.092) ranged from 0.012 to 0.23. The highest F_ST values were between sites in the north and those in the central and southern region (mean F_ST = 0.18) and the lowest values were among the central sites (mean F_ST = 0.092). Analysis of F_ST across all sites indicated that higher connectivity, based on absence of a river barrier, was correlated with higher gene flow (β_{river_pres} = −0.011, β_{river_dist} = 133.8, F = 35.78, P = 0.001, R² = 0.68). A more parameter-rich model supporting multiple connectivity hypotheses was found (β_pathdist = −1.04e-05, β_{river_pres} = 20.76, β_{river_dist} = −0.027, β_modCWD = −1.20e-06, F = 61.69, P = 0.001, R² = 0.89), but the extra coefficients, while significant, were close to unity. Focusing only on the Central Sierra sites, there was again support for watercourses as barriers to gene flow (β_{river_pres} = 54.99, F = 53.32, P = 0.001, R² = 0.34). See Fig. 7 for an example that considered two equidistant pairs of sites, with higher connectivity predictions correlated with lower values of F_ST and thus higher levels of gene flow. This result was corroborated by the migrant rates estimated by BayesAss, with the lowest rates of migration into and out of anthropogenic refugia (Table 1).

Fig. 7

An example of a connectivity surface map for survey points based on isolation related to river or stream presence. Two pairs of sites that are each approximately 45 km apart have very different gene flow levels correlated with different levels of connectivity

Table 1 Bayesian assessment of migration rates among Central Sierra sites, derived by BayesAss

Discussion

This study reveals clear support that climate change refugia, here identified as areas that have changed little over the twentieth century in terms of temperature and precipitation, can predict persistence, with a focus on the climate-sensitive U. beldingi. Yosemite National Park meadows that were predicted to be refugial based on several climate change indices were more likely to contain U. beldingi. Across California, persistence over the last century was higher in sites that were more stable climatically, with the exception of sites that were artificially supplemented with food or water. Allelic richness, an index of effective metapopulation size and population persistence, was highest for the coolest sites across the central Sierra Nevada. Although sample size precluded a direct test of the correlation between change in temperature and precipitation with genetic diversity across sites, this indicates that refugial sites that warm less in the future will better maintain U. beldingi populations. Thus, we show that evidence from occupancy surveys can be combined with simple genetic data to test and validate the buffering qualities of climate change refugia.

We considered several different delineations of climate change refugia in this analysis in order to explore the idea of using a single species as a test of whether an area is a climate refugium. However, species undoubtedly differ in the climate factors that are most relevant for their persistence [37]. Our previous research [38] indicated that climate change has resulted in increased body size in U. beldingi, perhaps related to longer growing seasons; thus less access to food does not appear to be driving population extirpations. We found that minimum temperature was a good predictor of persistence based on occupancy and genetic measures. Other studies have shown that mammals living in temperate montane systems could be negatively affected by warming in the winter [16, 39]. It is hypothesized that homeothermic species such as American pika (Ochotona princeps) [40] could freeze to death if the insulating snow layer disappears due to higher minimum temperatures in winter. This would be particularly relevant for hibernating species like U. beldingi that can lose necessary fat stores trying to maintain body temperature over winter. Shifts in winter timing also appear to be associated with declines in fitness of Columbian ground squirrels [41].

Conversely, it’s possible that our refugia measures were correlated with other variables and thus are just indirect predictors of persistence in the face of climate change. For example, predators, parasites, or habitat may be directly related to the change in climate, and U. beldingi may actually be responding to these other factors. Further study, and experimental manipulations, would be required to partition this variation.

Our study also provides support for biologically relevant habitat connectivity as a predictor of functioning climate change refugia. Among the connectivity hypotheses examined, the impact of watercourses appeared to have the most explanatory power for U. beldingi. As expected, we found that the small, semi-fossorial U. beldingi finds a river to be a substantial barrier to movement, whereas even the largest roads crossing the mid-elevations of the Sierra Nevada are not complete impediments to their dispersal. Genetic diversity and gene flow were also predicted by the habitat connectivity models.

Previous work showed that anthropogenic refugia had the highest densities of U. beldingi [16]; however, our study found they had the lowest gene flow and habitat connectivity values. For example, the Mono Lake County Park population was genetically distinct from the nearest Central Sierra Nevada populations (Fig. 2b); other known historical populations of U. beldingi had all been extirpated from the area [16]. Thus, this research indicates that Mono Lake County Park may be acting as a refugium for a previously more ubiquitous, genetically distinct U. beldingi population in an area that is particularly low and warm for the species. Though there was almost no immigration into these sites, they may be acting as a source for nearby populations. However, this complex example highlights the need to distinguish anthropogenic refugia from ecological traps [42].

Conclusions

Climate change refugia will be important for conserving populations as well as genetic diversity and evolutionary potential [43]. As during past epochs, they may act as safe havens of genetic diversity over the forthcoming decades and centuries of drastic change across the landscape, which is consistent with our finding across the full set of study sites that recent environmental shifts and isolation were negatively associated with persistence. As such, they are key targets with potential for management action [1, 44, 45], particularly while more long-term, sustainable climate change adaptation strategies are developed. Our study shows that in-depth modeling paired with rigorous fieldwork can identify functioning climate change refugia for conservation.