Spatio-temporal variability in the strength, directionality, and relative importance of climate on occupancy and population densities in a philopatric mammal, the American pika (Ochotona princeps)
non-stationarity Spatio-Temporal Analysis Rocky Mountains Species-climate relationships distribution models (naive occupancy) density American pikas
Peter D Erik A Marie L Dylan K
Species distribution models (SDMs) have been widely employed to evaluate species–environment relationships. However, when extrapolated over broad spatial sca...
https://brand.frontiersin.org/m/ed3f9ce840a03d7/favicon_16-tenantFavicon-Frontiers.png
Frontiers
1 Introduction

Contemporary climate change continues to shape species’ distributions and biological communities globally. Observational and experimental studies have now documented responses spanning shifts in distributions, abundance, body size, phenology, behavior, and biotic interactions, among many others (Eastman et al., 2012; Thurman and Garcia, 2019; Goodman et al., 2022; Calabrese and Pfennig, 2023; Kerner et al., 2023). Shifts in distributions and declines in abundance have caused particular concern, given that both indicate species may not be evolving rapidly enough to cope with the impacts of recent climate change, thereby increasing risks of extirpation and extinction (Radchuk et al., 2019). Accordingly, understanding the underlying biological processes (i.e., mechanisms) that influence species responses remains a priority for management decision-making and species-distribution modeling efforts (Urban et al., 2016).

Pinpointing the factors that limit distributions, and how these factors vary in importance over space and time, has long been of interest in ecology (Darwin, 1859). To understand and model multidimensional niche spaces, species-specific life histories should be well understood; simultaneously, investigators should also consider how factors constraining populations might vary across temporal and spatial scales. For example, Brown (1984), emphasized that the abiotic factors constraining distributions differ between cool- and warm-edge range limits. Although rising temperatures may cause declines in occupancy at warm limits, they may concurrently allow greater survival and expansion beyond current cold-edge range limits (e.g., Sultaire et al., 2022). Given the recent motivation to identify evolutionary “winners’’ and “losers” under climate change based on biological traits, a binary classification of species’ vulnerability would be an oversimplification, considering that species often simultaneously benefit and suffer from climate shifts, depending on populations’ range position and numerous other factors (Somero, 2010). This, along with spatial variability in local adaptation, microrefugia, disturbances, and biotic interactions, may partially explain why species–climate relationships tend to be markedly heterogeneous over space and, thus, difficult to predict.

Given the large variability in the relative importance of environmental constraints across species ranges, species distribution models (SDMs) have recently sought to account for some of this heterogeneity (Peterson et al., 2011; Guisan et al., 2017; Martínez-Minaya et al., 2018). Broadly, correlative SDMs aim to statistically and spatially model distributions using known or suspected relationships between organisms and environmental conditions. One recent meta-analysis concluded that SDMs generally perform best when predicting species occurrences, but sequentially decline in predictive ability for population abundance, fitness, and genetic diversity (Lee-Yaw et al., 2022). Although SDMs have had increasing success (i.e., predictive accuracy) in recent years, one major limitation is they generally cannot be applied elsewhere in a species’ range other than where the training data originated (Charney et al., 2021; Rousseau and Betts, 2022). Mechanistic modeling offers an alternative to correlative approaches by modeling distributions based on individual-level, physiological processes that are hypothesized to constrain species (Riddell et al., 2017; Briscoe et al., 2022), such as water requirements necessary to dissipate heat (Riddell et al., 2019). Nonetheless, mechanistic models require high-quality biological data that explain how species experience their environment and are therefore limited by our understanding of how environmental factors actually constrain species (Kearney et al., 2010; Urban et al., 2016).

Both correlative and mechanistic modeling approaches currently struggle with accommodating nonstationarity, a phenomenon wherein ecological processes and relationships vary in strength, direction, and relative importance over space and time (Rollinson et al., 2021). Nonstationarity can arise from shifts in species’ relationships with environmental variables both spatially and temporally, as well as from shifts in environmental conditions. For example, one recent study found that snowshoe hares (Lepus americanus) are primarily constrained in their distribution by snow-cover duration, but the degree of the snow’s importance interacts strongly with maximum temperature across the species’ range (Sultaire et al., 2022). Mounting evidence suggests that due to widespread nonstationarity in ecological relationships, inferences from one time period and region should not be used to infer or predict range dynamics elsewhere or in other time periods (Pardi et al., 2020); accordingly, heavy reliance on range-wide distribution models could create the possibility of overgeneralizations and inaccurate models for more-localized climate-adaptation management and decision-making. Collectively, these observations raise the following questions: 1) What are the most appropriate scales to investigate to best understand species–environment relationships (e.g., Guisan et al., 2007), especially to be most useful for resource managers (Carroll et al., 2022)?; 2) Should future monitoring efforts ideally focus on local, regional, or continental scales?, and, accordingly, 3) Where should model training and validation data be acquired from (e.g., Scherrer et al., 2021)?; 4) When are SDMs transferable over time in a reliable manner (e.g., Smith et al., 2013)?; and 5) What range of environmental conditions are necessary to survey to ensure models are robust (e.g., Thuiller et al., 2004)? Ecological forecasts remain limited today by our overall nascent understanding of varying degrees of nonstationarity across both space and time, and climate is only one niche component that makes SDM transferability difficult (e.g., these models often ignore biotic interactions and habitat fragmentation and loss, Díaz et al., 2019).

Further complicating advances in forecasting climate-driven responses, rates of warming are remarkably heterogeneous over earth, with high-elevation and high-latitude regions generally experiencing much faster rates than elsewhere (Wang Q et al., 2016). For example, warming across four of earth’s largest mountain chains (the Alps, Andes, Himalayas, and Rockies) has been notably faster than surrounding lowlands, whereas differences in precipitation have been less consistent (Pepin et al., 2022). The variable rates of warming may therefore leave mountain-dwelling populations at greater risk relative to their low-elevation counterparts under shifting conditions. Complex topography in mountains provides a high diversity of local climates for species to exploit, highlighting how mountain systems may now act as exceptional testing grounds for questions related to ecological nonstationarity in species responses to climate over local and regional scales, as well as through time.

Here, we evaluate species–climate relationships using occupancy and population density data to evaluate the degree of spatiotemporal nonstationarity in a small mammal. To do so, we tested relationships at several scales including within ecoregions, across ecoregions, and over time, using the American pika (Ochotona princeps, Richardson, 1828) as our focal species. American pikas (“pikas,” hereafter) are small, climate-sensitive lagomorphs that occupy patchily distributed, broken-rock (e.g., talus) habitats in western North America, predominantly in mountainous areas (Hafner and Sullivan, 1995). We used a newly integrated dataset of occupied and historically occupied patches that were surveyed between 2011 and 2021. These data cover three expansive ecoregions, and all field surveys used identical sampling protocols. Because shifts in climate often cause contrasting demographic responses across seasons (e.g., Cordes et al., 2020), we examined both summer and winter metrics that we hypothesized act as proxies of specific, mechanism-based pathways that shape population dynamics (SI, Table 1). Accordingly, we incorporated seasonal metrics for temperature, precipitation, and moisture. In this study, we specifically predicted the following:

1. Across ecoregions, the strength and direction of climate relationships of both response metrics (occupancy and density) will exhibit high spatial nonstationarity, with coefficients displaying large variability within each climate metric.

2. Within an ecoregion, species–climate relationships will show higher congruence in strength and directionality among mountain ranges that are in close geographic proximity, compared to those further away.

3. Between two closely spaced time periods, the strength and direction of species–climate relationships will not change within each climate metric. In other words, such relationships will exhibit relative temporal stationarity.

2 Methods

2.1 Model organism

We focus on American pikas here given their often-high detectability (p > 0.90, Beever et al., 2010; Ray et al., 2016) as well as our ability to quantify population density from patch-level surveys using visual and auditory detections. Due to their limited physiological thermoregulation, pikas occupy broken-rock features, like talus fields and lavascapes, that buffer against fluctuating ambient temperatures (Benedict et al., 2020). This habitat is also readily delineated using satellite imagery prior to field surveys and persists over centuries to millennia with little structural change. Studies that have incorporated mechanistically-informed variables suggest that heat and cold stress are often the best predictors of pika distribution and density across various ecoregions (e.g., Beever et al., 2003; Wilkening et al., 2011; Johnston et al., 2019; Billman et al., 2021), and abundance has often been linked strongly to precipitation and water-balance (Beever et al., 2013; Johnston et al., 2019), highlighting the species’ overall responsiveness to the direct effects of weather and climate. In addition to direct exposure, the indirect effects of climate have also been found to shape occupancy and/or population density, such as precipitation’s role on the quality, diversity, and availability of forage (e.g., Erb et al., 2014; Ray et al., 2016; Varner et al., 2023). Although temperate mountain regions exhibit high seasonality throughout the year, pikas do not hibernate, which can leave them vulnerable to both winter conditions, such as low snowpack, and summer conditions. Due to numerous life-history characteristics, pikas appear to have a comparatively lower adaptative capacity in the face of climate change than many other montane mammals in western North America, making them an ideal model organism for testing climate-responsiveness theory and questions (Beever et al., 2023).

2.2 Field surveys: occupancy and population density

We surveyed for evidence of current or historical pika occupancy in talus patches across the Rocky Mountains (USA) between 2011 and 2021. Patches – each defined as a contiguous area of broken-rock habitat separated from all other such areas by >25 m – were identified using satellite imagery on the online mapping platform CalTopo.com prior to fieldwork. Field surveys occurred from June to September and were restricted to crepuscular hours, to achieve the highest detectability (few surveys between 12:00 and 16:00; these occurred only when conditions were cooler, with fully overcast skies). To survey, we walked 50-m transects across each talus patch, spaced 15 m apart elevationally, to obtain estimates of population density while avoiding double-counting individuals (# of pikas/50 m, sensu Beever et al., 2011). We considered a patch occupied if we unequivocally detected pikas visually, aurally, or if we found fresh haypiles with at least 10 fresh sprigs (green with chlorophyll) of herbaceous vegetation in them. A patch was considered historically occupied or extirpated if we only found old evidence such as old fecal pellets or degrading haypiles with no fresh vegetation present. Patches that appeared unoccupied were sampled with additional time and effort to be confident in their unoccupied classification. Additional information regarding field-survey protocol can be found in Beever et al. (2011; 2013) and Billman et al. (2021).

Our surveys spanned three scales to assess variability in species-climate relationships: inter-regional (among ecoregions), intra-regional (i.e., among mountain ranges within a single ecoregion), and between two sampling periods. The inter-regional comparison included data from the Southern Rockies (New Mexico), Central Rockies (Wyoming and Montana), and the Northern Rockies (Montana and Idaho) which represent three of six recognized lineages of O. princeps (Schmidt et al., in prep, Galbreath et al., 2010). Site-selection protocol differed slightly among ecoregions (e.g., elevational-transect design in Northern Rockies vs a near-census in Southern Rockies), so we provide an in-depth description of each region’s site-selection protocol in the Supplemental Information. The intra-regional dataset covers the Beaverhead Mountains, Italian Peaks, Lemhi Range, and the Tobacco Root Mountains nested within the Northern Rockies, which were surveyed in 2018–19. The temporal-analysis data were collected in the Central Rockies, with the “historical sampling period” spanning 2011–2016 and a “recent sampling period” from 2019–2021. These sites were all located within the Greater Yellowstone Ecosystem, including Yellowstone and Grand Teton National Park. All sites included in the temporal analysis were surveyed in both periods. However, several fewer sites were occupied in the latter period and excluded from the density analysis but retained (marked as unoccupied) in the occupancy analysis. Although these time periods seem relatively short when examining changes in occupancy, interannual patch turnover has been documented to be high for pikas in response to weather variability (e.g., Rodhouse et al., 2018). We excluded all talus patches without evidence of current or past pika occupancy, as these sites may not always have been structurally or otherwise suitable (following Billman et al., 2021). However, this may create minor differences between our study and others that have retained such no-evidence sites.

2.3 Climate data

Seasonal climate data were extracted from the standalone program, ClimateNA (Wang T et al., 2016). These data are locally downscaled from PRISM climate rasters using bilinear interpolation and account for site-specific elevation, latitude, and longitude. Summer (JJA) and winter conditions (DJF) are hypothesized to be the most stressful to pikas, so we extracted Summer: minimum, mean, and maximum temperature (respectively, Tmin, Tmean, Tmax), precipitation, relative humidity (RH), and climate moisture index (CMI), and Winter: Tmean and precipitation. We hypothesized each of these variables represents a mechanistic pathway of climatic influence on pikas. For example, winter precipitation (i.e., snowfall) has been shown to have a positive effect on site occupancy as higher snowpacks buffer pikas from sub-zero temperatures that occur above the snow (Beever et al., 2010; Johnston et al., 2019). Similarly, snowfall may be associated with higher population densities as heavy snowpacks provide a water source for high-quality forage in spring and summer in montane systems (Yandow et al., 2015). Full descriptions of the hypothesized mechanisms tied to each climate stressor can be found in Table 1 of the SI appendix.

For occupancy, we calculated the average climate conditions over the four years’ seasons prior to surveying. Winter metrics included the winter immediately preceding surveys and the previous three winters, whereas summer metrics comprised the previous four summers. The 4-year window minimized overlap in the temporal dataset in the Central Rockies (i.e., ~95% resurveys were 4+ years after initial surveys, mean = 5.2 years). We excluded current summer conditions for occupancy because a large proportion of surveys occurred in June, so current-season metrics would therefore represent future conditions (after surveys occurred) in many instances. For density predictors, we included only the winter conditions immediately preceding surveys, as well the previous year’s summer conditions (to avoid using data from after the date when a survey occurred in early summer).

2.4 Data analysis

For each region separately, climate metrics were scaled and centered using Z-scores to evaluate the relative importance of metrics that often differed by one to several orders of magnitude. Numerous studies have already sought to identify the top models and predictors of pika occupancy and density using multivariate modeling and AIC (e.g., Beever et al., 2011; Johnston et al., 2019; Billman et al., 2021), so we explicitly did not include that objective in this study. Rather, we examined the univariate effects of each climate variable on both responses. We used simple linear regression modeling for population density and logistic regression modeling for occupancy. As expected, density estimates showed a strong positive skew, so estimates were log-transformed to better meet the assumptions of parametric testing. We excluded unoccupied patches in the density analyses to differentiate the possible divergent drivers of occupancy and density. Given the species’ consistently high detectability (p ≥ 0.90, Beever et al., 2013; Ray et al., 2016), we used naive logistic regression models for occupancy, which assume perfect detectability (p = 1); naive models for this species have already been shown to be robust against parameter estimate bias in previous investigations (e.g., Rodhouse et al., 2010). We used the lm() function in base R for density analyses, and the glm() function in the lme4 package (Bates et al., 2015). All plots were created using ggplot2 (Wickham, 2011). Effect sizes and 95% confidence intervals were extracted from all model summaries. Coefficients for the logistic regression models are on the log-odds scale. Principal component analysis (PCA) plots were created using the ggfortify package (Tang et al., 2016), more specifically the prcomp() and autoplot() functions which scaled the eight climate variables prior to plotting (Figure 1). All analyses were conducted in R Statistical Software (v.4.2.2; R Core Team, 2021). Overall, we evaluated effect size strength (i.e., magnitude), directionality (+/−/non-significant if 95% CI included zero), and relative importance compared to other variables.

FIGURE 1

3 Results

There was marked heterogeneity in species–climate relationships across space and time. We surveyed and included a total of 1,867 talus patches for the inter-ecoregional dataset, 734 patches for the intra-regional dataset in the Northern Rockies, and 264 patches for the temporal dataset in the Central Rockies. Inter-regionally, approximately 30% of sites (N = 565) were historically occupied (i.e., possibly extirpated). Patches were distributed between 1,381 and 3,884 m above sea level, and breakdowns of patch occupancy by region and mountain range can be found in Table 2 of the SI Appendix. Patch sizes ranged from 1 pika home range (~315 m2, ~20-m-diameter circle) to patches in Grand Teton NP that spanned >1 km2. Lastly, patch-level population densities ranged from 0.23 to 4.0 pikas/50 m surveyed.

3.1 Inter-regional analysis

Across the three ecoregions, occupancy and density were all significantly negatively associated with the four temperatures metrics (mean β coefficient ± 1 SD; occupancy: −0.90 ± 0.21, density: −0.23 ± 0.08) and positively associated with the two moisture metrics (RH and CMI, mean β coefficient ± SD; occupancy: 0.62 ± 0.25; density: 0.25 ± 0.11). However, there was more nonstationarity in precipitation metrics across space than other metrics (Figure 2). Winter precipitation (i.e., snowfall) was positively associated with both responses in the Northern and Central Rockies, but unexpectedly not in the Southern Rockies. On the other hand, summer precipitation was only positively associated with both response metrics in the Northern Rockies, as well as occupancy in the Southern Rockies, but there were no relationships between summer precipitation and density in the Southern Rockies and both responses in the Central Rockies. The overall strength of species–climate relationships was stronger in the Northern and Central Rockies relative to the Southern Rockies. Notably, the precision of coefficient estimates did not consistently increase with greater sample sizes across regions.

FIGURE 2

3.2 Intra-regional analysis

All relationships coefficients were significant and consistent in directionality (within variables) across mountain ranges, with the exception of all Italian Peaks density variables, and Tmean in winter and relative humidity in summer for occupancy (Figure 2). Temperature was, again, consistently negative for both response metrics (when excluding the Italian Peaks, for density). Moisture and precipitation metrics were all positively associated with occupancy and density, as hypothesized (when excluding the Italian Peaks). For occupancy, the Beaverhead Mountains overwhelmingly displayed stronger species–climate relationships (i.e., larger coefficients) than did the other three mountain ranges, with only a few exceptions (Figure 2). Overall, coefficients showed a greater range across the various climate metrics in the intra-regional analysis relative to the inter-regional analysis (occupancy: inter-regional range = −1.18 to 1.47, intra-regional range = −1.76 to 3.30; density: inter-regional range = −0.39 to 0.36, intra-regional range = −0.41 to 0.64).

3.3 Temporal analysis

Similar to the inter- and intra-regional analyses, temperature consistently had a negative effect on both occupancy and density over time, with the exception of Tmax in summer for density (Figure 3). Variability in the strength of relationships to temperature was notably greater for density than for occupancy. Relative humidity was significant only in the historical sampling period for both response metrics. Climate moisture index was positively associated with occupancy in both periods and density in the historical period, but negatively associated with density in the current period. Winter precipitation (i.e., snowfall) consistently had a positive effect on both responses in both periods, whereas summer precipitation had no effect on occupancy in either period, and was only significant for density in the current period. However, this summer-precipitation effect on density was surprisingly negative (β = −0.38, 95% CI: −0.52 to −0.23). Overall, there was lower congruency in the strength and directionality than expected between the two time periods that were only several years apart (e.g., only ~69% of cases showed consistent directionality).

FIGURE 3

4 Discussion

Here, we examined how species–climate relationships vary across spatial and temporal scales using a readily detectable small mammal whose range spans approximately one-third of North America. Overarchingly, we found higher temperatures consistently had a negative effect on both pika occupancy and densities, while moisture metrics mostly had positive effects on both responses. However, the importance of precipitation in summer and winter showed the highest nonstationarity across space, with winter snowpack having no effect on either response in the Southern Rockies, and summer precipitation having no effect on either response in the Central Rockies. Moreover, the precision of effect sizes was generally the lowest for precipitation metrics, both at inter- and intra-regional scales. Our findings illustrate the existence of some nonstationarity of biotic responses to certain climate factors, cautioning extrapolation beyond single regions or time periods.

4.1 Among-ecoregion variability

At the inter-regional scale, which spans over twelve degrees of latitude and nearly ten degrees of longitude, there was general congruence in approximate coefficient size and directionality of temperature metrics for both responses. Overwhelmingly, warmer temperatures had a negative effect on both pika occupancy and density across the three vast ecoregions. In winter, warmer temperatures are typically correlated with lower snowpack, followed by lower snow retention into summer, which may explain why warmer winters have a negative effect. These results corroborate numerous other studies that highlight the sensitivity of pikas to warmer temperatures across seasons (Beever et al., 2003; Wilkening et al., 2011; Stewart et al., 2015). In contrast, several other investigations have found the effect of temperature on pikas to vary across study areas (Jeffress et al., 2013; Rodhouse et al., 2018). Nevertheless, high summertime temperatures have also been implicated in constraining dispersal abilities and corridors used by pikas (Henry et al., 2012; Castillo et al., 2014), suggesting that the previously occupied sites in hotter areas may be less likely to be recolonized. Given that all three ecoregions have warmed in recent decades, our results indicate that further temperature-related patch-level extirpations seem likely for the species, unless increases in precipitation possibly ameliorate thermal stress.

In addition to temperature, we found full congruence in the directionality of the two moisture metrics across ecoregions, suggesting that water balance constitutes a predictable selective pressure on pikas. However, the two seasonal-precipitation metrics were less consistently significant, showing several insignificant relationships depending on the ecoregion and response metric. Such differences appear to align with the amount of snowfall received and its duration, across our study regions. Downscaled data from ClimateNA are tightly correlated with weather-station measurements for temperature metrics (R2 ≥ 0.99) but not as much with precipitation metrics (R2 ≥ 0.64), leading to substantially lower predictive power for precipitation metrics (Wang T et al., 2016). Therefore, we cannot rule out possible climate-downscaling discrepancies as a contributor to insignificant relationships with precipitation. Still, we suspect the lack of consistent relationships with precipitation may be due to unaccounted-for interactions with temperature that may be biologically relevant. For example, another study found that high temperatures were indeed more influential on pika occupancy in drier national parks (Jeffress et al., 2013). Whereas Erb et al. (2011) found that annual precipitation best predicted pika occupancy across 69 sites in the Southern Rockies, we found that only precipitation in summer had a small positive effect on occupancy, whereas winter precipitation surprisingly had no effect in the Southern Rockies. West of the Rocky Mountains, in the Great Basin ecoregion, Millar et al. (2018) described the climate envelope of pika-extant sites as ranging from 7–83 mm for summer precipitation, whereas our pika-occupied sites had consistently more summer precipitation (Southern Rockies: 137–392 mm, Central: 68–247 mm, and Northern: 82–203 mm). Millar et al. (2018) suggested that Great Basin pikas can accommodate substantially more-arid summer conditions than what we have found in the Rocky Mountains, indicating that thresholds of pika sensitivity may vary across ecoregions. For example, 99.4% of our historically occupied (possibly extirpated) sites in the Rockies were still wetter in term of summer precipitation (Southern Rockies: 135–387 mm, Central: 57–239 mm, and Northern: 75–208 mm) than pika-occupied sites in the Great Basin (Millar et al., 2018). For summer Tmax, Millar et al. (2018) reported their pika-extant sites ranged from 11.4 to 29.9°C; our extant sites fell entirely within their ranges (14.4–26.8°C). However, our warmest (Tmax) historically occupied sites across our three regions only reached 27.4–29.1°C, suggesting a possible lower threshold for summer heat tolerance in the Rocky Mountains relative to the Great Basin.

Across temperature and precipitation metrics, our results continue to illustrate how climate-envelope thresholds may differ substantially across species’ ranges. Such differences may also, in part, reflect genetic, behavioral, habitat (e.g., microrefugia, forage quality), or other differences among populations that remain difficult to model. Our findings reiterate other studies’ calls for consideration of the context-dependencies of scale prior to formulating range-wide generalities (Schwalm et al., 2016). Nonstationarity in precipitation metrics especially warrants caution of their inclusion in broad-scale modeling, whereas temperature and moisture metrics appeared to play a more predictable role in shaping occupancy and population densities, herein. For precipitation metrics, the use of spatially varying coefficients (SVCs) in distribution modeling may be one effective solution to account for this large spatial variation in the relative and absolute importance of these metrics.

Investigations into spatial nonstationarity in species-environment relationships have been growing in recent years, with studies often documenting the greatest spatial nonstationarity in climate metrics specifically (Pease et al., 2022b). Identifying the key drivers of variability in these relationships is especially important for management decision-making and policy. Given that management actions are generally implemented at local- to regional-scales, results from one region applied elsewhere should be evaluated with particular care since species–climate relationships vary across spatial and temporal extents, over which environmental stressors can differ (Pease et al., 2022a). Yet, species-specific data limitations may constrain robust modeling efforts at local scales when resources and data are limited. This may be especially pronounced in rural areas, where species occurrence records are less common than near urban settings (Bowler et al., 2022). One solution to fill spatial gaps at local scales may involve heavier reliance on citizen-science data for additional occurrence and density records in some instances, but experts may still need to determine true absences for presence-absence models. In addition, there are widespread biases in occurrence records across species due to various reasons (e.g., body size, rarity, and group size, Callaghan et al., 2021), but intraspecific species–environment relationships often do not correlate to patterns with heterospecifics in the same community or in closely related species (e.g., Rossi et al., 2023). Accordingly, considerations may be necessary at the species level. Our results show that the level of species-specific nonstationarity depends on which climate metrics are evaluated in modeling (i.e., some metrics, like temperature, exhibited surprising levels of stationarity over large spatial domains), a pattern that is likely prevalent across taxa. We therefore caution that community-scale modeling that does not account for species-specific relationships with climate may be ineffective for predicting the future distributions of species, especially when extrapolating into novel conditions and the distant future.

4.2 Within-ecoregion variability

Relative to the inter-regional comparison, intra-regional coefficients were more variable, showing a larger range of strengths (occupancy: within-region: −1.76 to 3.30, among-regions: −1.18 to 1.47; density: within-region: −0.41 to 0.64, among-regions: −0.39 to 0.36). However, this may partially be due to smaller sample sizes within mountain ranges relative to across mountain ranges (Table 2, SI Appendix). Excluding the Italian Peaks, responses in the other three mountain ranges showed consistent directionality of each response to climate, with both summer and winter precipitation and moisture-related metrics having reliably positive effects. For these four variables, coefficients still varied widely in most cases (e.g., relative humidity for occupancy, coefficient range: 0.32–2.3 on the log-odds scale). Temperature had consistently negative effects on both occupancy and abundance, in alignment with the inter-regional analyses. The general consistency of response directionality among the three mountain ranges provides some promise that within ecoregions, results from one mountain range may be relevant (to some extent) to other nearby mountain ranges for both occupancy and density in a given time period.

Results from the Italian Peaks diverged from the other three ranges, whereby all climate relationships were insignificant for density, and winter Tmean and summer relative humidity were insignificant for occupancy. To note, the Italians Peaks were the only mountain range where pika population densities decreased slightly with elevation (see Figure S2 in Billman et al., 2021). The underlying geology of the Italian Peaks is predominantly limestone. Unlike the larger metamorphic and igneous talus found throughout the other mountain ranges, limestone talus in this mountain range was typically smaller in size (10- to 50-cm-diameter rocks, on average) and, consequently, shallower than larger-clast talus slopes. Although numerous studies have promoted optimism for the microclimate buffering capacity of talus as a long-term refugia for pikas (e.g., Millar et al., 2013), one recent study found that deep, subsurface temperatures in talus in Colorado (USA) had still warmed substantially between 1963–64 and 2009–21 (Monk and Ray, 2022). Accordingly, we suspect that the lack of pika–climate relationships for density in the Italian Peaks may largely be due to poor thermal buffering of the small, limestone talus; smaller rock diameters and shallower talus depth likely translate to warmer interstitial temperatures at low elevations and colder temperatures at high elevations. Thermally stressful conditions at the lowest and highest elevations in this mountain range would suggest a linear relationship for density would have no detectable effect; instead, the relationships may be quadratic, with the most favorable conditions for population densities at the mid elevations. On the contrary, the Beaverhead Mountains generally showed the strongest effect sizes, particularly for occupancy (Figure 2), which may have been due to topographic similarity of sampled watersheds in this mountain range. This range is especially linear and all watersheds on the west have comparable elevational spans of habitat and minimum elevations of occupancy, whereas the east-side watersheds retracted very little and also have comparable minimum elevations of current occupancy to one another (Billman et al., 2021). These elevational patterns of occupancy (vs extirpation) likely explain why the pika-climate relationships there were often the strongest compared to the other ranges. Lastly, the Beaverhead sites exhibited the broadest climate space in a PCA that included all eight climate variables (Figure S1). Overall, these context-dependencies highlight necessary caution when applying interpretations among mountain ranges with differing topography and geology, for example.

Although some empirical studies found that species–climate relationships become more stationary, and predictable, at broader scales (e.g., Miller and Hanham, 2011), others have found stationarity decreases at larger scales, such as when study regions are spaced far apart (e.g., McAlpine et al., 2008). Discrepancies here likely arise because far-apart study domains may or may not be similar in terms of climate, vegetation, and species composition, among other factors. Because the four mountain ranges in our Northern Rockies intra-regional analysis had similar climatic conditions and other similarities, like vegetation, we predicted that stationarity would be higher in terms of coefficient directionality and magnitude relative to the inter-regional analysis. However, we surprisingly found less congruence in effect-size magnitude and more consistency in climate-relationship directionality (when excluding the Italian Peaks). Interestingly, pika–climate relationships for precipitation and moisture metrics were consistently stronger (i.e., larger magnitudes) at local scales relative to the inter-regional analysis, suggesting that variability in range-wide modeling datasets may dampen the overall strength of species–climate relationships at local scales, which are more relevant to managers.

On the contrary, we found comparable coefficients for temperature metrics between intra- and inter-regional analyses for density (mean temperature β coefficient: −0.22 vs. −0.23, respectively), and occupancy (−1.06 vs. −0.90, respectively). Because eco-evolutionary processes are hierarchically structured in nature, patterns found at local scales often do not align with those at much broader scales (Osborne et al., 2007), making these findings particularly surprising. For species with large distributions that span immense climate gradients, species–climate relationships should generally be examined at the ecoregional scale when possible due to comparable topography, geology, vegetation, and climatic conditions (Smith et al., 2019). Within said ecoregions, we advise sampling efforts span the largest gradient of environmental conditions possible to evaluate variability in relationships over space. Given what we found intra-regionally, anomalous conditions in select mountain ranges (within an ecoregion) may warrant their exclusion from regional models, but they should also receive additional sampling efforts in the future would be beneficial to better understand and predict species–climate relationships therewithin.

In this case, pika–climate relationships for temperature metrics may be more practical than precipitation and moisture metrics for extrapolating to new regions or at different scales if necessary, but further exploration seem warranted. Our findings overall have important implications for species distribution modeling in the face of climate change, given that temperatures are rising across most of the planet, but changes in precipitation vary considerably and are less predictable (IPCC, 2022). Yet, precipitation and moisture are predicted to grow increasingly important for terrestrial species moving forward to maintain body cooling requirements with rising temperatures (Riddell et al., 2019). Provided that extirpations often occur at the hottest and driest sites within a species range, understanding spatial nonstationarity in precipitation and moisture metrics are critical for informing and modeling extinction risks.

4.3 Temporal variability

In general, we found high agreement of directionality within climate metrics over time, but less so for climate-density relationships than for climate-occupancy. For both responses, we identified only one case where the directionality of an effect changed sign through time; density was positively influenced by climate moisture index in the historical period, but CMI had a negative effect in the more-recent survey period. In contrast, the results from the occupancy analysis suggested moisture had a consistently positive effect in both periods, highlighting how the drivers of occupancy and density may not always align, a result documented in other studies as well (e.g., Dibner et al., 2017). Climate-density relationships are often less predictable than occupancy, since shorter-term processes generally drive population densities whereas changes in occupancy occur over longer periods and are often influenced by larger-scale processes like macroclimate (Beever et al., 2013; Schulz et al., 2019).

There were several cases where species–climate relationships were significant in one period but not the other in our analysis. This temporal nonstationarity appeared unpredictable and was found in temperature, moisture, and precipitation metrics (density: Tmax, relative humidity, and summer precipitation; occupancy: relative humidity). One possible explanation for temporal divergence is that correlations among climate metrics are often nonstationary themselves over time, whereby changes in temperature do not always correlate with changes in precipitation, for example (e.g., Bueno de Mesquita et al., 2021). In our case, summer relative humidity had a positive effect on occupancy in the historical period but not the more-recent period, whereas all other relationships with occupancy and climate were significant (except summer precipitation in both periods). This may indicate that temperature, winter precipitation, and CMI were more limiting in the more-recent survey period than relative humidity. Similarly, the relative importance of climate variables may have changed over time (Marcinkowski et al., 2015). Our mixed results regarding temporal nonstationarity somewhat stand in contrast to another pika study that found relatively high temporal transferability of SDMs in California (Stewart et al., 2015). Although the ecological explanations discussed above may indeed account for this change in variable importance over time, the range of climate values could differ between periods, where conditions may never have reached some stress threshold in the more-recent period. However, this possibility is unlikely in the case of relative humidity, given that the historical-period range was 44.5–70.5% and the recent-period range was 44.5–71.6% – nearly identical. In fact, the overall climate space of surveyed sites between periods had remarkable overlap when using a PCA that included all eight climate metrics (Figure S2). Overall, these patterns caution interpretations of results from single-time period surveys, including those in the two spatial analyses above. While biological interpretations of nonstationarity in species–climate relationships remain important, mathematical evaluation of whether the range of predictor values (prior to scaling and centering) might vary substantially over time over time would provide a valuable complement; this, in turn, could strongly alter the magnitude and directionality of predictors on the response across time periods.

A majority of studies on temporal stationarity in ecology thus far have focused on tree rings. These studies have frequently found that correlations between ring growth and climate fluctuate between significant and non-significant; such fluctuations have been attributed to changes in weather patterns through time. For example, one study found that the effect of winter precipitation on ring growth has decreased over time and become insignificant, whereas the effects of spring temperature and prior-year summer temperatures have become significant recently (Marcinkowski et al., 2015). This often unpredictable de-coupling of species–climate relationships on interannual and interdecadal timescales has wide-ranging implications for SDMs, which typically assume temporal stationarity (i.e., mean and variance change do not over time). However, assuming stationarity is unlikely to capture the complexity of underlying ecological processes and how variable they can be through time. Similar to previous tree-ring studies, our temporal-scale results highlight how assuming stationarity, even over a short time span, could also be misleading and generate incorrect interpretations of animal–climate relationships as well (Bueno de Mesquita et al., 2021). Accordingly, recent studies have recommended that temporal nonstationarity (rather than stationarity) should be the baseline assumption when modeling species–climate relationships through time (Tumajer et al., 2023).

Responses to weather and climate conditions are also highly variable across species (Rossi et al., 2023). Therefore, the degree of nonstationarity can be expected to differ markedly among species within the same community and site (Hereş et al., 2022). This adds an additional layer of complexity when modeling but must be considered prior to applying interpretations from one species to another. Collectively, a stronger emphasis should be placed on understanding temporal nonstationarity in animals to more effectively inform future SDMs, as well as management strategies that rely on timely, near-term forecasts. More broadly, in the face of deterministic climate change, directional selection across species may be an important evolutionary process for many populations over time (decades to centuries); however, increased frequency of extreme events may instead lead to destabilizing selection as well.

4.4 Important considerations and caveats

The primary goal of this study was not to identify top models for predicting occupancy and density, but rather to evaluate variability in species–climate relationships over space and time using univariate climatic factors. Accordingly, there are key interactions among climate factors that likely also explain important mechanisms of stress in this species. We only modeled linear relationships for density, which may or may not be the most descriptive model structure. We acknowledge these simple models constitute a great over-simplification of the full dynamics governing species occupancy and abundance. However, given 1) particular questions that we and our management partners sought to address in the various regions, and 2) our primary objective in this study to assess nonstationarity in species–climate relationships across space and time, we believe our current approach most directly addresses this objective. In addition, our study did not model relationships with non-climatic factors that are also important for occupancy and density, such as habitat suitability and connectivity, forage quality and abundance, interspecific interactions, microclimate temperatures, among others. As described elsewhere, metapopulation dynamics, Allee effects, human disturbances, and other transient processes also cloud the relationship between species and their environment (Holt, 2020). We acknowledge that, particularly relevant for our analyses of density, error and uncertainties are not propagated in other ways through our analysis; Bayesian approaches may better account for this noisiness, explicitly.

Given that American pikas have been considered a sentinel species for understanding metapopulation dynamics and the impact of climate change on alpine species, our findings of consistently strong species–climate relationships for most variables and scales are particularly exciting and highlight how strongly macro-climate appears to influence this species along with micro-climate conditions. Although our sampling design differed slightly among the Southern, Central, and Northern Rockies, we found many consistent relationships among regions (i.e., stationarity, see Tmax in summer) which suggests that methodological differences in study design may not always lead to divergent results and interpretations. Nevertheless, we do not encourage applying different survey designs over space when it’s possible to avoid. The presence of local adaptation across regions and the two pika lineage here was also beyond the scope of this study but may well be an important factor in species–climate relationships when comparing multiple ecoregions that have long had vastly different climatic conditions. These considerations are important when evaluating the applicability of SDMs for any given species, or otherwise projecting results from one study region or timescale to another. However, further studies examining species’ relationships with univariate climate factors would be beneficial to improving our understanding of species’ ecology, physiology, and conservation.

Nonstationarity, if unaccounted for, decreases model performance and reliability, but so does climate novelty when models aim to predict species’ responses beyond the range of the calibration data. Nonstationarity can reflect true biological or mechanistic changes in species–environment relationships over space or time; in contrast, model error can emerge simply due to inappropriate extrapolation into novel conditions. Whereas the focus of our investigation was to examine issues related to nonstationarity, ignoring climate novelty during SDM creation and application will also provide biologically unrealistic inferences and deteriorate model predictability. In addition, novel conditions and nonstationarity may act synergistically in influencing model performance and thus, overall utility. The development of increasingly flexible model structures should seek to account for both ecological nonstationarity and novel conditions when possible, as both are predicted to increase under continuing contemporary climate change (Radeloff et al., 2015).

5 Implications

The impending loss of global biodiversity underscores the need to evaluate species’ vulnerabilities to further shifts in climate, overexploitation, invasive species, and land-use patterns (Díaz et al., 2019). Species’ exposure and adaptive capacities are now recognized to vary widely across spatial scales (Beever et al., 2023) but given that ecological relationships are notoriously variable over multiple spatial and temporal dimensions, understanding patterns in new places, time periods, and under no-analog conditions can be difficult. However, there are growing opportunities to now capitalize on “big-data” systems that span many disciplines to better evaluate the importance of scale. Growing recognition of nonstationarity in ecological processes has highlighted remaining questions related to the appropriate scales for investigations, the limits of extrapolation, and the overall usefulness of spatially explicit models. Although our results were mostly comparable across three of the four ranges within the single ecoregion, range-wide modeling using climate envelopes from one region and time period would provide inaccurate and unrealistic results in most cases – particularly for precipitation, according to our results. Holistically, our results stress how a better understanding of the appropriate spatial and temporal scales of species–climate responses could create more effective, tailored climate-adaptation strategies under shifting conditions.

Analyses in recent years have begun exploring new mathematical methods in their models, like spatially varying coefficients (SVCs), to account for nonstationarity in species–climate relationships over space and time (e.g., Humphreys et al., 2022; Sultaire et al., 2022), justification for which is further supported by our analyses. Similarly, geographically and temporally weighted regression and machine learning in ecology should be considered in more instances, given their often-superior performance relative to other methods in forecasting processes in cases where nonstationarity is high (Feng et al., 2021). Nevertheless, model structure choice should ultimately be motivated by the study question at hand (Segurado et al., 2006). While we encourage additional investigations into nonstationarity in species–climate relationships across spatial domains, we anticipate future investigations into how these relationships change through time will yield the greatest insights in a rapidly changing world. Although the full complexity of ecological networks can never be fully modeled, we hope the results here provide further motivation to continue discussions surrounding the relationships of patterns and scale in order to improve biodiversity management and forecasts moving forward.

Data availability statement

The original contributions presented in the study are publicly available. This data can be found here: USGS Science Data Catalog (SDC), https://doi.org/10.5066/P9WB1EWC.

Ethics statement

Ethical review and approval was not required for the animal study because this study only included observational data.

Author contributions

PB, conceptualization, methodology, data collection and curation, formal analysis, investigation, writing–original draft, writing–review and editing, visualization, and supervision. EB, conceptualization, methodology, data collection and curation, investigation, writing–original draft, writing–review, and editing. MW and DR, data collection and curation, writing–original draft, and writing–review and editing. All authors contributed to the article and approved the submitted version.

Funding

This work was funded by numerous organizations and grants, including the American Society of Mammalogists GIAR (PB, MW), Northwest Scientific Association (PB), Sigma Xi (PB), American Alpine Club (PB), Western National Parks Association (EB), U.S. Geological Survey (EB), University of New Mexico Joseph Gaudin Scholarship (MW), and the New Mexico Research Grant (MW). Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fevo.2023.1202610/full#supplementary-material

References