Introduction

The ability of climate refugia to facilitate species persistence under adverse past conditions signals their potential importance in protecting biodiversity with future change (Keppel et al. 2012; Morelli et al. 2016). Refugia are characterized by environmental conditions conducive to the formation of relatively stable habitats (Keppel et al. 2018) and occur at multiple spatial and temporal scales (Ashcroft 2010; Morelli et al. 2020). Climate refugia of the past have been important in the formation of global biodiversity hotspots (Myers et al. 2000; Harrison and Noss 2017), and areas of relatively greater climate stability may provide holdouts and stepping-stones critical to species survival (Keppel and Wardell-Johnson 2015). Refugia are one of the key determinants of species’ persistence in the longer term (Keppel et al. 2012) and in complex shorter-term movements (Pecl et al. 2017). Environmental heterogeneity that defines refugia can facilitate range shift dynamics and buffering of habitats from a changing climate and associated disturbances, including drought, wildfire, and pests. The location and effectiveness of refugia will play an important role in both survival of individuals and populations, and ability to colonize more favorable habitats through time (Simons-Legaard et al. 2013; Andrews et al. 2018).

At local scales, microrefugia are defined by environmental heterogeneity often related to topography (Boyko 1947; Dobrowski 2011; Hampe and Jump 2011; Keppel et al. 2012). Topographically complex landscapes contain highly variable microclimates that can provide refugial habitats for specific taxa and serve as foundational populations when conditions become favorable in the broader landscape. Where topography acts as a deterministic control on the geographical distribution of habitat, these places may continue to function as microrefugia with reduced rates of climate change in comparison with the surrounding landscape (Mackey et al. 2012; Morelli et al. 2016). Understanding the topographic settings that form refugia can be instrumental in connecting species’ ecology with conservation and management of valued resources (Morelli et al. 2016).

Isolated mountain systems, or sky islands, are potentially important locations for climate refugia, due to outstanding species richness and endemism along topographic gradients (Warshall 1995; Coblentz and Riitters 2004; He and Jiang 2014; Mastretta-Yanes et al. 2018). Yet climate change poses a distinct threat to sky island refugia due to elevation constraints on forest species’ upward movements, leading to potential habitat loss and increased isolation between sky islands. For example, important insights for biodiversity and conservation with forecasted climate change have emerged from studies of tepuis, the flat-topped and isolated table mountains of the Guayana Highlands in South America (Vegas-Vilarrúbia et al. 2012; Nogué et al. 2013; Safont et al. 2016). Climate refugia were identified at broad scales where tepuis species, adapted to a limited range of environmental conditions, concentrated in remaining habitat patches (Vegas-Vilarrúbia et al. 2012). Identifying potential refugia at regional scales can provide critical information for in situ conservation strategies in mountain environments (Scherrer and Körner 2011), and complementary understanding of microrefugia that may increase the probability of species persistence at finer scales is also needed (Vegas-Vilarrúbia et al. 2012).

Here, we focus on describing climate refugia for pines (Pinus spp.) in the Madrean sky islands, located in the borderlands of México and the United States. The region is a widely recognized biodiversity hotspot owing to its unique geography, where the Rocky Mountains and Sierra Madre Occidental converge (Warshall 1995; Spector 2002; Mittermeier et al. 2004), bringing together temperate and subtropical biomes. The north–south configuration of these sky islands likely contributed to dispersal in the past, increasing biodiversity with fluctuating climates. Advancing and receding montane vegetation dominated by pines and oaks drove vicariance and speciation, with movement fostered by low-elevation passes connecting the two cordilleras as well as the Chihuahuan and Sonoran deserts (Findley 1996).

The region contains 17 native species of pines (Ferguson et al. 2013) which comprise an integral part of the diverse forests and woodlands that underpin ecological functions but are under risk from habitat loss and fragmentation due to climate change and wildfire (Barton and Poulos 2018; Shirk et al. 2018; Yanahan and Moore 2019). Many species of pines in the region occur at the margins of their distributions (Ferguson et al. 2013), increasing vulnerability to stressors such as drought, disturbance, and limits to dispersal ability (Hampe and Petit 2005). Notably, Madrean pine–oak woodlands of the southwestern United States and México comprise one of 36 global biodiversity hotspots endangered by climate change and anthropogenic impacts (González-Elizondo et al. 2012; Conservation International 2020).

However, conservation efforts in the Madrean sky islands are hindered by the lack of high-quality and consistent transboundary datasets (Pezzoli et al. 2014), including information on the location and attributes of climate refugia. Availability of data and models that support these efforts is especially important in areas where cultural and economic diversity, and a wide range of land uses, jurisdictions and management goals create extreme environmental pressures (Villarreal et al. 2019). The success of cooperative efforts between conservation groups and public–private partnerships to harness scientific resources and address pressing conservation issues exacerbated by climate change in the United States-México transboundary region depends on a foundation of well documented data and models (Aburto-Oropeza et al. 2018; Villarreal et al. 2019; Norman 2020).

Data needs for development of useful descriptions of climate refugia include physical, climatic, and ecological properties of the study landscape based on spatial data. Identifying physical characteristics of preferred habitats can provide insights into the microclimates that support plant survival and reproduction; for example, rugged terrain can afford protection from wind and desiccation, and slope position and aspect are linked to water accumulation and availability and evaporative demand (Dobrowski 2011; Hoylman et al. 2019). Bioclimatic variables integrate topographical and climatic variability; these variables are designed to capture climatic parameters that are meaningful to species’ biology and ecology, including mean annual temperature, annual precipitation, seasonality, and extreme or limiting environmental factors (Wang et al. 2016). In addition, remote sensing indices have distinct advantages in capturing refugial habitat characteristics, including canopy structure and greenness, and reflecting locally heterogeneous conditions and ecological context at high resolution (Cord et al. 2014; Villarreal and Petrakis 2014; Copenhaver-Parry et al. 2017). Seasonal reflectance characteristics can indicate changes in refugial habitat due to anthropogenic and natural disturbances like stand replacing wildfires, drought, insect outbreaks, or logging and are potentially useful for monitoring the condition of refugia through time.

In order to address existing data gaps regarding climate refugia for pines, we developed multi-scale species’ distribution models using field observations collected across the Madrean sky islands. Field data provide distinct advantages over data more commonly used in species’ distribution modeling which are often comprised of gridded species’ occurrence records with coarse spatial resolution (Araújo and Guisan 2006). The observations used in our analysis include reliable absence records that are especially important in estimating probabilities of occurrence by comparing environmental characteristics at sites where the species was observed with those at sites where it was not observed (Guisan and Thuiller 2005; Guillera-Arroita et al. 2015).

Our specific objectives were to (1) identify potential microrefugia for pine species at fine scales, defined by relationships to terrain, (2) determine the responses of pine species to the range of variability in climate parameters associated with topographically defined microrefugia, (3) examine variability in species’ occurrences to seasonal climate, and (4) identify hotspots and networks of pine microrefugia to inform regional and local conservation efforts.

Methods

Study region

The 74,788 km2 Madrean Archipelago Ecoregion (Omernik 1987), or Madrean sky islands, encompasses portions of the northern states of Sonora and Chihuahua in México, and southeastern Arizona and southwestern New Mexico in the United States (Fig. 1). Our study included 55 sky island mountain complexes within the ecoregion, with an area of approximately 45,000 km2 (Deyo et al. 2013). The international border divides the ecoregion near the center at 31° 33 N.

Fig. 1
figure 1

Elevation of the Madrean sky island study area with outline of island polygons (some island complexes are labeled for reference). Species’ distribution models were developed using field-based presence and absence data for pine species, recorded at sample locations (shown in red)

Full size image

Climate in the region is characterized by a bimodal precipitation pattern: frontal winter storms bring snow to higher elevations and rain to the lower elevations from November to March (Adams and Comrie 1997; Sheppard et al. 2002), and in mid- to late summer (July to August), moisture from the south creates a monsoonal weather pattern with thunderstorms. The arid foresummer (April to June), preceding the monsoonal moisture, is typically hot and dry and generally corresponds to the fire season. Annual precipitation varies greatly by elevation, with orographic uplift contributing to increased precipitation in upper elevations.

Pine ecosystems and study species

Pine species occur in a variety of forest and woodland communities across the sky islands (Rehfeldt et al. 2012; Villarreal et al. 2020). Spruce-fir and mixed conifer forests are dominated by Picea engelmannii and Abies concolor but also include Pinus ponderosa, P. strobiformis, and Pseudotsuga menziesii. Pine–oak forests and Madrean evergreen forests and woodlands contain a diversity of pines (P. arizonica, P. engelmanii, and P. leiophylla) along with a variety of oaks (Quercus spp.). Pinyon-juniper woodlands are comprised of P. edulis, P. discolor, and Juniperus spp.

The five pine species included in this study vary in their association with particular environments (Table 1). Reports of broad environmental tolerance suggest that diverse microclimates associated with topographic variability could be utilized by species including P. strobiformis. Likewise, reports of drought tolerance suggest that cooler and wetter microclimates (e.g., mesic, north-facing slopes inhabited by P. arizonica) may provide refugia to some degree as warmer and drier regional climates affect local conditions (Table 1). For P. discolor, tolerance of cooler, more mesic conditions at higher elevations could accommodate upward migration, while drought resistance could promote survival in its current range. Sensitivity of P. engelmannii to drought severity could make microrefugia that are decoupled from regional climate especially important to this species. Although P. chihuahuana has shown potential to adapt to warmer temperature, poor soils, and wildfire, conditions conducive to regeneration may benefit from sheltered habitats found in microrefugia (Table 1).

Table 1 Affinities of the five pine species and climatic environments documented in recent literature
Full size table

Data collection and assembly

Field data

We generated presence/absence data for each species using observations of pines from two field studies, one in the sky islands of México (Flesch et al. 2016; Flesch 2019) and another in the United States (Ganey et al. 1996; Sanderlin et al. 2013) (Fig. 1, Table 1). Both studies sampled vegetation within circular plots that were spaced evenly at intervals (200 m in México and 300 m in United States) along transects that covered the major vegetation communities of the sky islands. Field sampling occurred between 2009 and 2012 across 14 sky islands in México (n = 1098 plots in 123 transects) and 2014–2015 across five sky islands in the United States (n = 318 plots in 28 transects; Sanderlin et al. 2014). Data were collected within 50 m radius of the plot center point in México transects, where the proportion of total vegetation volume composed of each dominant tree species was estimated visually in both the tall- and short-tree layers of each plot. In pine–oak woodlands, pines were classified as tall trees if they emerged from a canopy stratum of shorter understory trees that in these systems are typically evergreen oaks. In the United States study, trees (> 10 cm diameter at 1.4 m above ground) were recorded by species within 17.8-m radius circular plots (0.1 ha). For the response variable in species’ distribution models, presence/absence was recorded based on whether individual trees were recorded in a sampled plot. Total sample sizes after combining observations from the México and United States field studies were n = 1416 for all species except P. discolor (n = 1202); a subset of plots in México did not include data for this species. Occurrence frequencies for each species across datasets were P. strobiformis 14%; P. arizonica 35%; P. discolor 23%; P. engelmannii 30%; P. chihuahuana 17% (see maps of presence and absence data, Online Resource 2).

Spatial data

Terrain metrics were calculated from a 30-m Shuttle Radar Topography Mission (SRTM; Farr et al. 2007) digital elevation model (DEM) using SAGA-GIS (Conrad et al. 2015). Metrics included local morphometrics, characteristics of the catchment basin, relative position, and terrain ruggedness (Table 2). The slope position metrics describe locations relative to surrounding landscape, for example, Valley Depth is defined as the difference between elevation and an interpolated ridge level. Slope position metrics represent hillslope scale variations in moisture availability and microclimate (Hoylman et al. 2019).

Table 2 Predictor variables used in species’ distributions models, data sources and description, and interpretation for the purposes of this study
Full size table

We included two bioclimatic variables, Temperature Difference (TD) and Annual Heat Moisture Index (AHM), calculated for the normals period 1981–2010 at 1 km spatial resolution (AdaptWest Project 2015; Wang et al. 2016; Table 2). Both variables had high loadings in a principal components analysis of climate in the region (Villarreal et al. 2020). The AHM variable relates to temperature and water balance and its potential limits on plant reproduction and regeneration (Barton et al. 2001; Crimmins et al. 2011; Hess and Fulé 2020). TD represents variability in coastal influences on temperature extremes. In the context of the current warming trend, (i.e., low TD) less warming is generally expected to occur in locations that are nearer to coastlines (Ashcroft et al. 2009).

Using the Google Earth Engine cloud computing platform (Gorelick et al. 2017), we developed two median pixel top-of-atmosphere reflectance composites of cloud-free 30 m Landsat Thematic Mapper imagery representing two distinct seasons: the late-summer monsoonal rainy season (composite period: August 1, 2010–October 30, 2010, henceforth referred to as ‘summer’) and the spring growing season (composite period: March 1, 2011–May 30, 2011, henceforth referred to as ‘spring’) (Table 2; Fig. 2). For each seasonal composite, we calculated the Normalized Difference Vegetation Index (NDVI) as a proxy for canopy structure and chlorophyll content (Gamon et al. 1995) and Normalized Difference Moisture Index (NDMI) as a proxy for vegetation liquid water content (Gao 1996).

Fig. 2
figure 2

Summer (A) and Spring (B) Landsat Normalized Difference Vegetation Index (NDVI) and Summer (C) and Spring (D) Landsat Normalized Difference Moisture Index (NDMI) for the Madrean sky island study region. Minimum–maximum values are scaled to each composite image

Full size image

Data analysis

We developed species’ distribution models using boosted regression trees, a method which incorporates advantages of tree-based methods and improves performance by fitting many models and combining them for prediction (Elith et al. 2008). All models were parameterized using the Bernoulli family (for presence and absence data); the proportion of observations used in selecting variables was set to 0.5; complexity of individual trees was 5. To parameterize the learning rate, we started with η = 0.01 and adjusted to a slower rate (η = 0.005 or 0.001) to improve performance. We used tenfold cross-validation for model tuning and evaluation. Accuracy was evaluated using mean AUC (Area Under Curve) or ROC (Receiver Operating Characteristic) for the cross validations, which is most useful when the goal is predicting realized distributions and when data include known absences (Jiménez-Valverde 2012).

In order to identify locations that potentially function as microrefugia for pine species at fine scales (Objective 1), we started by constructing a terrain-only model (nine variables) for each species. We developed additional models using the top four terrain predictors to determine pine species response to the range of variability in climate parameters associated with topographically defined microrefugia (terrain-bioclimatic models including AHM and TD; Objective 2). We graphed partial dependence functions for the three variables with highest relative influence in the terrain-bioclimatic models and evaluated the trends in terms of local and broader-scale climate relationships. A final set of models was used to examine the variability in species response to seasonal climate and provide information on shifts in microrefugia with current conditions (Objective 3). For this purpose, we developed models that combined the top four terrain predictors with either NDMI (spring and summer) or NDVI (spring and summer).

Model results (probability of occurrence) for terrain-bioclimatic models and terrain-NDVI models were mapped to identify hotspots, or regions of high density separated by regions of lower density (Azzalini and Torelli 2007; Nelson and Boots 2008) for application in regional conservation efforts (Objective 4). Spatial predictions (1 km resolution) were smoothed using a Gaussian kernel in a 7 × 7 cell window (3 km radius). Smoothing the data can improve visualization of clusters of similar values (e.g., hotspots) while maintaining a continuous scale of predicted probability. We compared the distribution maps for each species to identify distinct geographic trends and information gained through inclusion of different variables in the models. Finally, we developed a map of species diversity hotspots by summing the smoothed probabilities, described above, from the five terrain-bioclimatic models and rescaling the values from 0 (low) to 1 (high). All analyses were conducted in R (R Core Team 2021).

Results

Model fit varied across the four models and for each of the five pine species (Table 3). The addition of bioclimatic and remote sensing variables to terrain predictors improved the fit of terrain-only models, especially for P. engelmannii and P. discolor based on cross-validation statistics (Table 3; see Online Resource 3 for ROC plots). Models for P. strobiformis had consistent values of AUC (≥ 0.80). The models for P. arizonica and P. engelmannii had slightly poorer overall fit (AUC mean < 0.80). Models that included terrain and bioclimatic variables had the best fit for all species except P. engelmannii which showed greater improvement with the addition of spring and summer NDVI to the best four terrain variables.

Table 3 Cross validation statistics for sky island pine species distribution models
Full size table

The influence of terrain on species’ distributions (Objective 1)

Terrain variables provided important information in pine species’ distribution models. Elevation was the most influential predictor in the majority of models (Table 4 A and B). In P. arizonica models, the predominance of elevation was especially notable (relative influence = 54%). Pine species exhibited specific affinities for elevation zones (Fig. 3; right column); P. strobiformis and P. arizonica were most common at highest elevations, reaching plateaus in predicted probability in areas greater than 2500 m. In contrast, P. discolor had the greatest positive response to elevation in a narrow, low-elevation zone; a more variable change in response was observed at middle elevations. Species associated with lower elevations, such as P. engelmannii and P. chihuahuana, differed in response pattern; P. engelmannii exhibited a more gradual decline as elevation increased; in contrast, the response of P. chihuahuana declined precipitously at approximately 2000 m.

Table 4 A, B Relative influence statistics as indicators of variable importance from generalized boosted regressions
Full size table
Fig. 3
figure 3

Partial dependence plots for the top three predictors from the terrain-bioclimatic model for each pine species included in the study (Table 4). Plots are arranged to aid comparisons across species in relation to Annual Heat Moisture Index (left column); either Temperature Difference or a top terrain predictor (center column); and elevation (right column). The y-axis scale represents the change in predicted probability relative to the predictor variable. Partial dependence plots for the top four terrain variables and partial dependence plots for the top three terrain-NDVI variables are located in Online Resources 4 and 5, respectively

Full size image

The top four terrain predictors of pine species distributions (Table 4 A and B) included terrain ruggedness (TRI) for P. strobiformis, P. engelmannii, and P. chihuahuana. Aspect was one of the best predictors for P. arizonica and P. discolor (Table 4 A and B). Characteristics of the catchment basin (area or slope) were important in P. strobiformis, P. arizonica, and P. discolor models. At least one of the relative height metrics (slope height, mid-slope position and valley depth) appeared in all models. Species exhibited preferences for rougher (e.g., P. strobiformis) or smoother terrain (e.g., P. chihuahuana); in general, cooler aspects were important for P. arizonica while warmer aspects were preferred by P. discolor (Online Resource 4).

Species’ response to broad-scale climatic factors (Objective 2)

One or both bioclimatic variables, along with elevation, comprised top predictors of pine species distributions in the terrain-bioclimatic models, based on relative influence (Table 4 A). Two climate variables exhibited high relative influence in P. discolor, P. engelmannii, and P. chihuahuana models (relative influence > 16%). In P. strobiformis and P. arizonica models, AHM and terrain metrics (i.e., TRI and aspect) had greater weight than TD.

Dramatic shifts in response were observed with minor changes in AHM for all species except P. arizonica which exhibited a more gradual change in response across AHM (Fig. 3). Similar threshold patterns were evident for P. discolor and P. chihuahuana in models with TD as a top predictor (Fig. 3). In contrast, the response of P. engelmannii to TD showed a more gradual decline as TD increased.

Species’ response to seasonal climate variables (Objective 3)

Remote sensing variables were important explanatory factors of occurrences of the five pine species. For high-elevation species, spring greenness (NDVI) was the second most influential predictor in terrain-NDVI models (i.e., P. strobiformis and P. arizonica; Table 4 B). Higher elevation species were associated with increased spring greenness (P. strobiformis and P. arizonica); P. arizonica predicted response decreased with summer greenness (Online Resource 5). Two additional species, P. discolor and P. engelmannii, exhibited strong and distinct relationships with both spring and summer greenness (Table 4 B). Two species were associated with a specific range of values for both spring and summer greenness (P. engelmannii, P. chihuahuana Online Resource 5). Other species exhibited characteristic relationships to summer (P. chihuahuana) or both spring and summer (P. discolor, P. engelmannii) greenness (Online Resource 5). Terrain-NDMI results were not reported, due to similar AUC values compared with models that included seasonal NDVI.

Identification of hotspots and networks of pine microrefugia (Objective 4)

High probability of occurrence for the pine species had distinct geographic patterns with respect to latitude in both terrain-bioclimatic and terrain-NDVI models (Figs. 4 and 5). Hotspots (highest probabilities) for high-elevation species (P. strobiformis and P. arizonica) were primarily limited to northern islands. In contrast, P. discolor high occurrence probabilities had a broad geographic range that stretched across the northern and central portions of the region. The only species with higher probabilities limited to southern parts of the region was P. engelmannii. Maximum occurrence probabilities for P. chihuahuana were scattered across the islands, with clusters of higher values in central and southern latitudes.

Fig. 4
figure 4

Spatial predictions of pine species distributions based on the terrain-bioclimatic model for each pine species included in the study (1 km resolution with kernel smooth). The range of predicted probabilities for each species was stretched across a color ramp to emphasize neighborhoods that differed in low, medium, and high probability. Areas of high probability, or hotspots, for P. strobiformis and P. arizonica, were clustered in northern islands. Hotspots for P. discolor were evident on western islands and continued across toward southeast and central locations. The only species with a distinctly south-southwestern distribution of high probability was P. engelmannii. Only one hotspot was evident for P. chihuahuana, located near the United States-México border on Peloncillo-Pan Duro (see Fig. 1)

Full size image
Fig. 5
figure 5

Spatial predictions of pine species distributions based on the terrain-NDVI model for each pine species included in the study (1 km resolution with kernel smooth). The range of predicted probabilities for each species was stretched across a color ramp to emphasize neighborhoods that differed in low, medium, and high probability. Similar to predictions from terrain-bioclimatic models (Fig. 4), areas of high probability, or hotspots, for P. strobiformis and P. arizonica, were clustered in northern islands. Hotspots for P. discolor were more prominent on western and northern islands and scattered across central islands. Areas of high probability for P. engelmannii were found at lower latitudes. Distinct hotspots for P. chihuahuana were absent, but compared to terrain-bioclimatic models, its distribution had shifted northward (see Fig. 4)

Full size image

General patterns in spatial predictions of the terrain-bioclimatic results (Fig. 4) and the terrain-NDVI results (Fig. 5) were similar, but the magnitude of predicted probability of occurrence differed. Differences between spatial predictions for P. strobiformis and P. arizonica distributions were minor; however, greater differences in probability were observed for other species. The terrain-bioclimatic models resulted in higher probabilities overall for P. discolor than the terrain-NDVI models with the exception of near zero probabilities in the far south. The P. engelmannii probability values were higher on some western islands, based on terrain-bioclimatic models; terrain-NDVI models predicted higher maximum probabilities in the southern-most islands. Terrain-bioclimatic models for P. chihuahuana exhibited higher probabilities across middle and lower latitudes within its elevational range, with one notable hotspot on Peloncillo-Pan Duro (Fig. 4). In contrast, terrain-NDVI models for P. chihuahuana included some areas with probabilities in the medium range in more northern islands (Fig. 5).

Hotspots of diversity for the five pine species occurred across a broad swath of central latitudes and in large portions of some northern islands (Fig. 6). The gradient in summed probabilities (derived from the terrain-bioclimatic models, Fig. 4) illustrated that climate refugia with the greatest diversity of pine species (colored orange to red) were associated with smaller, more isolated islands to the north, including the Pinaleños, Santa Catalinas, and Chiricahuas (see Fig. 1 for reference). Areas with climate refugia for a diverse array of pines were also located in El Gato and Peloncillo-Pan Duro which border the Sierra Madre Occidental to the south and east. Climate refugia with lowest species diversity were most evident on southwestern and northeastern islands (i.e., represented by colors ranging from blue to green).

Fig. 6
figure 6

Climate refugia diversity hotspots illustrate areas where greatest sum of predicted probabilities (scaled 0 to 1) for the five species occur across the region (shaded orange to red). Intermediate scaled sums are colored blue to green, and lowest sums are blue. Values were calculated from smoothed terrain-bioclimatic model predictions (Fig. 4)

Full size image

Discussion

The Madrean sky islands are a stepping-stone archipelago that has served to facilitate movement of species between the Sierra Madre Occidental and the Rocky Mountains during past climatic changes (Warshall 1995). The presence of potential microrefugia on contemporary sky island landscapes may provide key locations to maintain local populations, or to support shifts in species distributions under current and future conditions. Climate refugia for pines were evident within the highly variable environments found on isolated mountains of the Madrean sky islands where topography provides deterministic control on habitat and potentially reduces the rate of climate change in microrefugia (Figs. 3 and 7). We identified a wide variety of settings for the formation of microrefugia across the region, defined by complex topography within preferred elevational zones. The availability of presence and absence data from field studies and a fine-scale digital elevation model enabled us to identify microrefugia as a primary focus. At a regional scale, we identified occurrence hotspots that define climate refugia (Figs. 4 and 5), which are essential foci for conservation of pines and a diverse group of taxa that depend on them for habitat. This included a northern range for P. strobiformis and P. arizonica, north-central locations for P. chihuahuana and P. discolor, and across more southern islands for P. engelmannii. Climate refugia harboring the greatest species diversity occurred both on smaller, more northern islands and across middle latitudes farther south (Fig. 6).

Fig. 7
figure 7

The Madrean sky islands form an archipelago of isolated mountain ranges between the Sierra Madre Occidental to the south and the Colorado Plateau to the north (map A). The sky islands are characterized by steep and rugged topography (photographs B, C, and D). We identified a wide variety of settings that create pine microrefugia across the region, facilitated by complex topography within preferred elevational zones (map A insets: zoomed in predicted probability surfaces and text boxes). Map insets of spatial predictions illustrate hotspots of high probability in dark green, and corresponding text boxes describe preferred topographic and climate conditions for each species: we observed a northerly distribution of high probabilities for P. strobiformis (map A top left) and P. arizonica (map A top right) north-central distribution of P. chihuahuana (map A bottom right) and P. discolor (map A middle left) and more southerly distribution of P. engelmannii (map A bottom left). Photo locations: Santa Catalina Mountains (B) by J.M. Iniguez, US Forest Service; Sierra Aconchi (C) and Sierra el Pinito (D) by A.D. Flesch, University of Arizona

Full size image

Characteristics of pine microrefugia

As the climate in the Madrean sky islands region becomes warmer and drier, resilience is likely to vary by species (González-Cásares et al. 2017). The varied characteristics of microrefugia represent important contributions to species persistence where topography slows local rates of climate change (Fig. 7). Each species preferred a specific range of elevation that corresponds to gradients in temperature, moisture, and seasonality, and regulates physiology and growth (Whittaker and Niering 1965; Barton 1993; Landsberg and Waring 1997; Guisan and Thuiller 2005). Overlap in affinities across elevation reflects diversification of pine forests and woodlands where multiple species aggregate as communities (Figs. 3, 6).

Across elevation, variability in terrain ruggedness, slope position, and aspect were important factors in species’ distribution models that contributed to defining microrefugia (Fig. 3, Online Resource 4). Slope and aspect are closely related to solar insolation and its influence on growing season length, moisture availability, and temperature stress. Species’ response to slope height indicates preference for particular microclimates associated with positions along hillslope gradients (Hoylman et al. 2019). Local depressions found in rugged terrain preferred by P. strobiformis (Fig. 3) are more sheltered from wind and desiccation (Dobrowski 2011).

Distinct thresholds in occurrences defined microrefugia for some species, where small changes in elevation and bioclimatic variables resulted in a dramatic shift in species response. Thresholds occur where a physical limit is reached or shift in functional response occurs (Romme et al. 1998). For example, probability of occurrence of P. chihuahuana declined precipitously between ~ 1800 and 2100 m elevation, and notably in a narrowly defined temperature range (Fig. 3). Species that occur in a narrow range of climate variability may experience an isolation effect that acts as a barrier to competition and negative influences such as pathogens (Ashcroft et al. 2012).

In contrast, P. arizonica exhibited a continuous, but variable, increase in response at elevations above ~ 1800 m and a steady decline with increases in heat moisture (AHM; Fig. 3). Species that exhibit tolerance to a wider range of conditions may be able to take advantage of variability along a climate continuum, dependent on species interactions and other factors. Species-specific trait responses, such as tight control of transpiration via stomatal closure exhibited by P. strobiformis, will play a role in acclimation to drier and warmer climate (Bucholz et al. 2020).

In some cases, pine species were entirely absent from specific elevation zones, evidenced by plateaus in response at upper or lower extremes in elevation. For example, P. strobiformis was absent below ~ 2000 m while P. engelmanni was absent above ~ 2250 m (Fig. 3). Microrefugia for P. engelmannii were associated with a low range of temperature difference (i.e., greater coastal influence), suggesting that in its optimal range of predicted probability, habitat for this species will be subject to slower rates of climate change than locations in the broader landscape (Ashcroft et al. 2009).

Microrefugia and seasonal climate response

Seasonal reflectance characteristics capture ecosystem responses to recent and historical landscape changes caused by anthropogenic and natural disturbances. Thus, models including remote sensing variables provide better predictions of the current (vs potential) distribution of species with information not otherwise accounted for in our microrefugia models. Our results indicate that seasonal remote sensing data are especially important for characterizing microrefugia for P. engelmannii and P. chihuahuana, two species that are widely distributed across a low to mid-elevation range (Fig. 3 and Online Resource 5). Pine habitats within this elevation range exhibit spectral characteristics associated with a more open canopy structure, which we associate with greening of herbaceous understory or co-dominant plants in response to summer precipitation.

The remote sensing variables were therefore more important for species in some southern islands where fires occur regularly (Villarreal et al. 2020), as well as where climate and disturbance conditions are favorable for seedling establishment and regeneration. For example, P. engelmannii exhibits fire-adapted traits such as thick bark, self-pruning of dead branches and "grass-stage" seedlings which enhance fire tolerance (Rodríguez-Trejo and Fulé 2003). Models for P. engelmannii showed the most improvement with the inclusion of summer greenness, and greater local spatial variability of hotspot predictions was apparent in the remote sensing model (Fig. 5).

Spatial hotspots for P. chihuahuana exhibited a distinct northward shift when seasonal greenness was considered (Fig. 5), compared with terrain-bioclimatic model predictions (Fig. 4). The affinity of P. chihuahuana with more open conditions that have a stronger ecosystem response to summer moisture contrasts with higher elevation pines that have more closed canopies and more subtle differences in canopy reflectance response during spring greening. These contrasting traits provide information that aides in distinguishing conifers at higher elevations with similar climatic conditions but are located in more northerly habitats (Buermann et al. 2008).

Management recommendations

Specific conservation measures will need to be developed to effectively conserve pine populations across the bi-national Madrean sky islands in light of existing priorities within places that vary in land ownership and tenure. Networks of potential refugia could become a useful focus to prioritize collaborative conservation efforts across ownership and political boundaries, based on areas with high probability of species’ occurrence identified in our study. It is critical that mapping and validation of refugia and monitoring their role in maintaining desired conditions be central components of regional management and conservation strategies (Morelli et al. 2016; Balantic et al. 2021). These assessments of the location and functions of refugia can take advantage of remote sensing data monitored in near-real time to evaluate the seasonal greenness and wetness conditions that are favorable to pines. Remote monitoring of temporal dynamics relative to species’ distribution thresholds may help to identify potential early warning indicators of climate change impacts to microrefugia. Such strategies are critical given the limited time and funds available for conservation efforts.

Fire management has been identified as a unifying management theme in the bi-national Madrean sky islands (Villarreal et al. 2019, 2020). In a restoration context, fine-scale information on pine refugia reported here could help identify areas for restoration and also protect refugia from risk of severe fire effects and shortened fire intervals (Barrows et al. 2020). Moreover, disturbance plays an important role in driving dispersal and recolonization into new habitats (Serra-Diaz et al. 2015), especially at the margins of species’ distributions (Slaton 2015; Travers-Smith and Lantz 2020). Fire will likely play a similar role where pines occur at their northern margins within the Madrean ecoregion as it affects soil and light conditions that are essential for growth and seedling survival (Barton 1993). Species that are reliant on animal-dispersed seed (i.e., P. strobiformis and P. discolor), have the ability to resprout (i.e., P. chihuahuana), or are fire-tolerant as grass-stage seedlings (i.e., P. engelmannii) may have an advantage given increased fire observed during recent timeframes on northern sky islands (Villarreal et al. 2020). Importantly, fire can result in the formation of refugia where locations that are disturbed less severely or frequently than the surrounding landscape harbor biological legacies and contribute to ecosystem recovery (Krawchuk et al. 2016, 2020; Haire et al. 2017). Recent studies indicate microrefugia played a role in survival of pines (i.e., P. leiophylla var chihuahuana and P. engelmannii) following a recent, large fire in the Chiricahua Mountains (Barton and Poulos 2018). In the Santa Catalina mountains, heterogeneous topography inferred local protection from fire (Iniguez et al. 2008) that allowed cohorts of P. arizonica to establish within longer fire intervals (Iniguez et al. 2016). Thus, managers and conservationists can benefit from understanding the relationships between microrefugia and species’ occurrence identified here with likely ties to fire patterns both in the past and in future.