Climate Refugia for Pinus Spp. In Topographic and Bioclimatic Environments of the Madrean Sky Islands of México and the United States&nbsp;

doi:10.21203/rs.3.rs-391708/v1

Download PDF

Research Article

Climate Refugia for Pinus Spp. In Topographic and Bioclimatic Environments of the Madrean Sky Islands of México and the United States

https://doi.org/10.21203/rs.3.rs-391708/v1

This work is licensed under a CC BY 4.0 License

Journal Publication

published 01 Apr, 2022

Read the published version in Plant Ecology →

Version 1

posted

You are reading this latest preprint version

Context

Refugia are island-like habitats that are linked to environmental stability. Where topography acts as a deterministic control, microrefugia may continue to function as habitat under reduced rates of climate change. Continental island ecosystems provide propitious settings for identifying patterns of refugia at multiple scales and applying that knowledge to conservation.

Objectives

Our main objective was to identify microrefugia for pines where habitats are defined by topographic heterogeneity. Secondary aims were to describe climatic variation within microrefugia and examine how species’ response to seasonal climate alters spatial predictions of microrefugia.

Methods

We investigated how topography forms microrefugia in the Madrean sky islands, located in the borderlands of México and the United States. Our design incorporated pine species presence and absence field observations (P. strobiformis, P. engelmannii and P. chihuahuana, P. arizonica and P. discolor), modeled in relation to terrain, bioclimatic and remote sensing predictors.

Results

Terrain ruggedness, slope position and aspect defined microrefugia for pines within specific elevation ranges. Some species had narrow habitat preferences (e.g., P. chihuahuana); others exhibited a broader range of tolerance (e.g., P. arizonica). Hotspots of microrefugia were either limited to northern islands or occurred across central or southern latitudes. Response to seasonal climates shifted distributions of hotspots for species with open canopy structure and where regular fires occur.

Conclusions

Pine habitats with greater climate stability may provide holdouts and stepping-stones critical to species persistence with future change. Networks of refugia provide a promising focus for conservation, restoration, and fire management across a diverse, binational region.

Ecological Modeling

Madrean sky islands

climate change

conservation

microrefugia

Pinus spp.

remote sensing indices

The ability of refugia to facilitate species persistence under adverse conditions of the past signals their potential importance in protecting biodiversity with future change. Refugia are island-like habitats that are linked to long-term environmental stability (Keppel et al. 2018) and can be identified at broad spatial and temporal scales (Ashcroft 2010; Morelli et al. 2020). Within broadly favorable climates, microrefugia defined by topographic heterogeneity can facilitate species persistence outside their physiological limits by modification of regional climate (Boyko 1947; Dobrowski 2011; Hampe and Jump 2011; Keppel et al. 2012). Where topography acts as a deterministic control on the geographical distribution of microrefugia, these places may continue to function as habitat under redfuced rates of climate change (Mackey et al. 2012; Morelli et al. 2016). Furthermore, disturbances result in the formation of refugia where locations that are disturbed less severely or frequently than the surrounding landscape harbor biological legacies that contribute to ecosystem recovery (Krawchuk et al. 2016; Haire et al. 2017; Krawchuk et al. 2020). An important strategy for conservation is to identify refugia that exist within high priority landscapes (Lawler 2009).

Continental island ecosystems, or sky islands, provide propitious settings for identifying macro- and microrefugia and applying that knowledge to conservation. The isolation of mountain islands coupled with high climatic variation that characterize sky islands are key components of macrorefugia where the evolution and persistence of unique species and communities are supported at multiple scales (Hampe and Jump 2011; Ashcroft et al. 2012; Morelli et al. 2020). The configuration of mountain complexes, their origins and geology, and the presence of bridges and barriers to movement promote high species richness and endemism along topographic gradients (Warshall 1995). Topographically driven isolation increases species diversification across elevation, with greater effects at lower latitudes (Steinbauer et al. 2016). Sky island ecosystems occur across the globe: the Himalaya-Hengduan Mountains in southwest China (Luo et al. 2016); the Western Ghats in India (Krishnadas and Osuri 2020; Jins et al. 2021); and the sky islands of eastern Africa (Popp et al. 2008).

Important insights for biodiversity and conservation with forecasted climate change have emerged from studies of tepuis, the flat-topped 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 mountain species, adapted to a limited range of environmental conditions, concentrated in remaining habitat patches (Vegas-Vilarrúbia et al., 2012). Evidence from sky islands of the Guayana Highlands indicate that identifying potential refugia at regional scales can provide critical information for in situ conservation strategies in mountain environments (Scherrer and Körner 2011). However, complementary understanding of microrefugia that may increase the probability of species persistence at finer scales is needed (Vegas-Vilarrúbia et al. 2012).

Here, we investigate the potential for refugia to promote conservation in the Madrean sky islands, located in the borderlands of México and the United States. The Madrean sky islands are a stepping-stone archipelago comprised of montane habitats separated by lowland grassland, desert scrub and subtropical thornscrub and occur across steep and complex topography (Coblentz and Riiter 2004; Coblentz 2005). The region is a widely recognized biodiversity hotspot owing to its unique geography, specifically, the convergence of temperate and subtropical biomes between two cordilleras, the Rocky Mountains, and Sierra Madre Occidental (Warshall 1995; Spector 2002; Mittermeier et al. 2004). 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 (Pinus spp.) and oaks (Quercus spp.) drove vicariance and speciation, with movement fostered by low elevation passes between the two cordilleras and connectivity between the Chihuahuan and Sonoran deserts (Findley 1996).

Pines are an integral part of the diverse forests and woodlands that underpin ecological functions across the Madrean sky islands and are under extreme risk from habitat loss and increased isolation 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 including drought stress, disturbance, and limits to dispersal ability (Hampe and Petit 2005). Identification of pine refugia can provide a focus for conservation efforts for these keystone species as they exist within a wide range of communities across their distributions.

Within the Madrean sky islands, pines occur in higher elevation forests with spruce-fir (Picea engelmannii-Abies concolor) and mixed conifer species (e.g., Abies concolor, Populus tremuloides, Pseudotsuga menziesii) primarily on northern islands (Rehfeldt et al. 2012; Villarreal et al. 2020). At lower latitudes and elevations, pine-oak forests include a mixture of species including P. arizonica, P. engelmanii and Quercus spp. Pines form a lesser component of evergreen forests and woodlands where oaks (Quercus spp.) are predominant. Pinyon-juniper (P. edulis, P. discolor, Juniperus spp.) woodlands are common across middle to lower elevations and latitudes. Notably, Madrean pine-oak woodlands of the southwestern United States and México comprise one of the 36 global biodiversity hotspots endangered by climate change and human impacts (González-Elizondo et al 2012; Conservation International 2020).

Resilience of pine forests and woodlands to climatic change is also influenced by land-use practices and variation in environmental stewardship and management. The bi-national Madrean sky island region is culturally and economically diverse, and a wide range of land uses, jurisdictions and management goals create environmental pressures (Villarreal et al. 2019a). Urban and exurban growth on both sides of the United States-México border, as well as vegetation clearing for border wall construction are reducing habitat area, quality, and connectivity for a wide range of species (Flesch et al. 2010; Villarreal et al. 2013; Peters et al. 2018). Activities that prevent or encourage burning, combined with warmer drier climates, can alter fire regimes in ways that reduce fire-free intervals and limit forest recovery (Iniguez et al. 2016; Villarreal et al. 2020). Cooperative efforts between conservation groups and public-private partnerships can help harness scientific resources and address pressing conservation issues exacerbated by climate change in this transboundary region (Aburto-Oropeza et al. 2018; Villarreal et al. 2019a; Norman 2020). However, conservation planning in the Madrean sky islands has been hindered by the lack of high quality and consistent transboundary data sets (Pezzoli et al. 2014).

In this study, we developed multi-scale species distribution models to identify and characterize refugia for pine species across the Madrean sky islands. The models serve to identify locations that function as microrefugia at fine scales where climatic change is expected to proceed more slowly due to heterogeneity in mountain topography (Hannah et al. 2015; Morelli et al. 2020). We then determined the range of variability in climate parameters associated with microrefugia, based on broader scale bioclimatic data. Last, we examined the variability in species response to seasonal climate to capture heterogeneity in species composition across forests and woodlands, and thus provide a perspective on shifts in microrefugia with current conditions.

Study region

The 74,788 km² 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 km² (Deyo et al. 2013). The international border divides the ecoregion near the center at 31°33 N.

The parallel valley and mountain configuration of the Madrean sky islands resulted from continental-scale deformation over the past several hundred million years in response to plate tectonics (Coblentz 2005). Valleys receded to expose older rocks derived from marine invasions, caldera eruptions and lava flows and metamorphic core complexes (Warshall 1995). Mountain ranges are surrounded by temperate desert grassland and both temperate and subtropical scrub communities formed at the junction of the Sonoran and Chihuahuan deserts.

Climate in the region is characterized by a bimodal precipitation pattern. Winter (November to March) frontal storms bring snow to higher elevations and rain to the lower elevations (Adams and Comrie 1997; Sheppard et al. 2002). In mid- to late-summer (July to August), moisture from the south creates a monsoonal weather pattern with thunderstorms. The arid foresummer period (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 at upper elevations.

Pine species descriptions

The five pine species vary in their association with particular environments (Table 1). Many studies reported relative affinities based on comparisons with other species, limiting interpretations of absolute limits or preferences. However, 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. 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. Affinities of the five pine species with the climatic environment, documented in recent literature. Geographic focus refers to data sources from Northern (mostly United States); Southern (primarily México) and Region-wide (sky islands and beyond). English and Spanish common names were reported in Felger et al. (2001, and verified by M.S. Gonzalez-Elizondo, pers. comm.). For some species, scientific names represent groupings of synonyms from Plants of the World Online (http://powo.science.kew.org) and International Plant Names Index (https://www.ipni.org): Pinus strobiformis; P. arizonica (syn. P. ponderosa var. arizonica); P. discolor (syn. P. cembroides var. bicolor); P. engelmannii; P. chihuahuana; (syn. P. leiophylla var. chihuahuana).

Species name	*Pinus strobiformis* Southwestern white pine Pino huiyoco		*Pinus arizonica* Arizona pine Pino blanco	*Pinus discolor* Border pinyon Piñón, pino piñonero		*Pinus engelmannii* Apache pine Pino real	*Pinus chihuahuana* Chihuahua pine Pino chino, pino chimonque
Geographic focus
Northern	Intermediate shade tolerance⁴ Broad environmental tolerance¹⁰ Associated with longer fire intervals^{3, 14, 15}		Moisture availability and fire intervals related to stand dynamics³ Thrives in disturbed sites with shade and exposed soil¹⁴ Tree establishment is associated with longer fire intervals^{3, 14, 15}	Water stress controls lower elevational limits; exhibits both drought resistance and tolerance of cold, more mesic conditions at higher elevations⁴		Water stress controls lower elevational limits; post-fire increase in light and temperature favorable to seedlings⁴ Sensitivity to water stress after wildfire can limit regeneration^{12, 13}	Water stress controls lower elevational limits; low light and deep litter control upper elevational limits; less shade tolerant than higher elevational pines⁴ Potential for adaptation to warmer temperatures, reduced soil productivity, and wildfires⁹ Sensitivity to water stress after wildfire can limit regeneration^{12, 13}
Southern	Linked to elevation and solar radiation, great temperature fluctuations across the year; large spatial distribution and niche breadth¹ Tolerant of semi-arid conditions⁶ Found in mesic sites in pine forests⁸ Inhabits northern aspect slopes and topographically sheltered sites¹⁰		Moderately sensitive to climate change⁶ Responds to changes in drought severity & cumulative water deficit⁷			Dry to moderately summer-warm open mountain slopes, mild climate, variable rainfall, elevation and maximum temperatures, open & mixed pine-oak² Sensitive to drought severity⁷ Found in semi-dry temperate areas in pine-oak forests⁸	Highly to moderately sensitive to climate change⁶ Sensitive to drought severity⁷ Found in semi-dry temperate areas in pine-oak forests⁸
Region-wide	Large range shift and contraction forecast; occurs in areas with high winter precipitation and low warmest month temperature¹¹		Inhabits mesic, north-facing slopes or canyon bottoms⁵	More drought tolerant than other pines⁵			More drought tolerant than other pines⁵


¹ Aguirre-Gutiérrez et al. 2015		⁷ González-Cásares et al. 2017			¹³ Poulos et al. 2020
²Ávila-Flores et al. 2016		⁸ González-Elizondo et al. 2013			¹⁴ Iniguez et al. 2008
³ Barton et al. 2001		⁹ Hess and Fule 2020			¹⁵ Iniguez et al. 2016
⁴ Barton 1993		¹⁰ Looney and Waring 2013
⁵ Ferguson et al. 2013		¹¹Shirk et al. 2018
⁶ Gómez-Mendoza and Arriaga 2007		¹² Barton and Poulos 2018

Field data and derivation of species’ presence/absence variables

We combined observations of pines from field studies conducted in the sky islands of México (Flesch et al. 2016; Flesch 2019) and the United States (Ganey et al. 1996; Sanderlin et al. 2013) to generate presence and absence data for each species (Fig. 1, Table 1). In México, vegetation data were gathered within 50 m of 1104 sampling points placed non-randomly in representative areas in major montane vegetation communities from 2009-2012. The majority of sky islands that support pines were sampled (Ferguson et al. 2013).

In the United States, data were collected within 318 plots between 2014 and 2015 (Sanderlin et al. 2014). Trees (> 10 cm diameter at 1.4 m above ground) distinguished by tree species were recorded within 0.1 ha circular plots (17.8-m radius). Plots were located along transects typically consisting of 12 points spaced at 300 m intervals (Sanderlin et al. 2013; Ganey et al. 2017), with transects being proportional to the approximate coverage of each montane vegetation type in the landscape.

For the response variable in species distribution models, species were considered present when that species was recorded as a tree within the sampled plot and absent when it was not recorded at the plot. Total sample sizes after combining observations from the México and United States field studies were n = 1422 for all species except P. discolor (n = 1204); a subset of plots did not include data for this species. Frequencies of occurrence 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).

Using field data provided distinct advantages over more commonly used data in species distribution modeling which are often comprised of gridded species’ occurrence records with coarse spatial resolution (Araújo and Guisan 2006). Field data 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).

Disadvantages of the sampling designs of the field studies include spatial bias (Araújo & Guisan, 2006) towards larger mountain ranges and potential for missing variability in mountain ranges where pines were assumed to be absent. Therefore, in some locations, predicted probabilities may reflect potential rather than actual distribution. Lower predictive ability may also result from factors not considered in the study design (i.e., island size, area of vegetation type) but should not obscure relationships to topography that define microrefugia.

Assembly of predictor variables

We assembled spatial layers of predictor variables across the geographic extent of each sky island complex (as defined by Deyo et al. 2013). We sampled the spatial layers at the field plot locations using the extract function in the R raster package (method=bilinear; Hijmans 2019). The dataset of predictor variables was appended to each species presence and absence database for use in modeling species distributions.

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 using in species distributions models, data sources, software used and definitions. Models were developed using 1) terrain metrics only; 2) top four predictors from the terrain model plus two remote sensing indices (NDVI or NDMI); 3) top four predictors from the terrain model plus two bioclimatic variables (TD and AHM).

Variable name	Source, description
Terrain metrics	30 m DEM, SAGA-GIS
Aspect	Local downslope direction with maximum rate of change (radians)
Slope	Local maximum rate of change in elevation (radians)
Elevation	Height above sea level (m)
Catchment Area	Size of the contributing area (m²)
Catchment Slope	Average slope of the contributing area (radians)
Mid-Slope Position	Intermediate slope position between high and low (0-1)
Slope Height	Vertical distance from the base of the slope to its crest (m)
Valley Depth	Vertical distance to a channel network base level (m)
Terrain Ruggedness Index (TRI)	Sum of the differences in elevation between a cell and its surrounding cells, calculated in 210 m radius (7 x 7 30-m cell) moving window
Remote Sensing Indices	Composites of 30 m Landsat imagery, Google Earth Engine
Normalized Difference Vegetation Index (NDVI)	Spring (2011-03-01 to 2011-05-30) and Monsoon (2010-08-01 to 2010-10-30)
Normalized Difference Moisture Index (NDMI)	Spring (2011-03-01 to 2011-05-30) and Monsoon (2010-08-01 to 2010-10-30)
Bioclimatic Variables	Normals period (1981-2010), AdaptWest
Temperature Difference (TD)	Difference between Mean Temperature of the Coldest Month (°C), and Mean Temperature of the Warmest Month (°C)
Annual Heat Moisture Index (AHM)	(Mean Annual Temperature + 10)/Mean Annual Precipitation/1000)

Using Google Earth Engine cloud computing platform (Gorelick et al. 2017), we developed two median pixel top-of-atmosphere (TOA) reflectance composites of cloud-free 30 m Landsat Thematic Mapper imagery representing two distinct seasons: the summer monsoon (composite period: August 1, 2010 – October 30, 2010) and the spring growing season (composite period: March 1, 2011- May 30, 2011) (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).

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). We selected AHM and TD based on their importance in a principal components analysis of climate in the region (Villarreal et al. 2020). Climate normals (1981-2010) for TD and AHM had low correlations (Pearson’s r < 0.50) with each other as well as terrain and remote sensing predictors. The 1981-2010 normals period overlapped or just preceded dates of field data collection.

The AHM variable is derived from Mean Annual Temperature (positive correlation) and Mean Annual Precipitation (negative correlation) (Table 2); temperature and water balance impose limits on plant reproduction and regeneration (Barton et al. 2001; Crimmins et al. 2011; Hess and Fule 2020). Differences in winter and summer temperatures determine TD values; TD represents variability in coastal influences on temperature extremes (Table 2). In the context of the current warming trend, less warming is generally expected to occur in locations that are nearer to coastlines (i.e., low TD; Ashcroft et al. 2009).

Model development and spatial predictions

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). We developed the models using the R packages gbm (Greenwell et al. 2019) and dismo (Hijmans et al. 2017). 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 10-fold cross validation for model tuning and evaluation. Accuracy was evaluated using mean area under the receiver operating characteristic (ROC) curve (AUC) 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).

We constructed a terrain-only model (nine variables) and three additional models using the top four terrain predictors: terrain-bioclimatic (AHM and TD); terrain-NDMI (spring and monsoon); and terrain-NDVI (spring and monsoon). We graphed partial dependence functions (R package pdp; Greenwell 2017) 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.

We mapped pine species distributions (probability of occurrence) at 1 km resolution based on the terrain-bioclimatic model and the terrain-NDVI model using the predict function in the R raster package (Hijmans 2019). Spatial predictions were smoothed using a Gaussian kernel in a 7 x 7 cell window (3-km radius) to improve visualization of clusters of similar values (e.g., hotspots) while maintaining a continuous scale of predicted probability. We compared the two distribution maps for each species to identify distinct geographic trends.

Terrain variables provided important information in pine species distribution models and additional predictors resulted in further improvement. Average AUC values were generally high for all five species (Table 3). Specifically, model fit statistics show that the addition of bioclimatic and remote sensing variables to terrain predictors improved the terrain-only models, especially for P. engelmannii and P. discolor (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 were lower overall, compared with models for other species (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 monsoon NDVI to the best four terrain variables.

Table 3. Cross validation statistics for sky island pine species distribution models. The AUC mean for the 10-fold cross validations from generalized boosted regression is given for each of the four models developed for the five pine species. Highest AUC values, indicating the best-fit model, are highlighted for each species.

Species	All terrain predictors	Top four terrain + AHM + TD	Top four terrain + two remote sensing variables
	AUC mean	AUC mean	AUC mean NDMI	AUC mean NDVI
P. strobiformis	0.815	0.844	0.828	0.831
P. arizonica	0.747	0.765	0.764	0.763
P. discolor	0.765	0.888	0.803	0.807
P. engelmannii	0.649	0.710	0.709	0.713
P. chihuahuana	0.745	0.803	0.760	0.767

The top four terrain predictors of pine species distributions (i.e., those included in Table 4 results) included terrain ruggedness (TRI) for P. strobiformis, P. engelmannii and P. chihuahuana and aspect for P. arizonica and P. discolor. Characteristics of the catchment basin 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 but were never top predictors in the terrain-bioclimatic models.

Table 4. Relative influence statistics as indicators of variable importance from generalized boosted regressions. Percent relative influence is given for the terrain-bioclimatic models (A) and the terrain-NDVI boosted regression models (B). TRI = Terrain Ruggedness Index; TD = Temperature Difference; AHM = Annual Heat Moisture Index. Best-fit models (from Table 3) are starred.

Pinus strobiformis		Pinus arizonica		Pinus discolor		Pinus engelmannii		Pinus chihuahuana
Terrain-bioclimatic model**		Terrain-bioclimatic model**		Terrain-bioclimatic model**		Terrain-bioclimatic model		Terrain-bioclimatic model**
AHM	39.64	Elevation	54.29	TD	33.51	Elevation	31.52	Elevation	27.59
Elevation	29.02	AHM	12.12	AHM	19.81	TD	23.35	AHM	23.61
TRI	9.39	Aspect	9.83	Elevation	16.43	AHM	18.83	TD	19.85
Slope Height	9.27	Slope Height	8.73	Aspect	10.57	Mid-slope Position	10.16	Slope Height	10.84
TD	6.44	TD	8.20	Catchment Slope	10.39	TRI	8.60	Valley Depth	10.76
Catchment Slope	6.20	Catchment Area	6.80	Slope Height	9.25	Slope Height	7.51	TRI	7.31

Pinus strobiformis		Pinus arizonica		Pinus discolor		Pinus engelmannii		Pinus chihuahuana
Terrain-NDVI model		Terrain-NDVI model		Terrain-NDVI model		Terrain-NDVI model**		Terrain-NDVI model
Elevation	52.26	Elevation	64.07	Elevation	27.39	NDVI monsoon	30.38	Elevation	43.15
NDVI spring	15.74	NDVI spring	12.66	NDVI spring	20.11	Elevation	28.87	Slope Height	17.87
Slope Height	9.58	NDVI monsoon	7.67	NDVI monsoon	16.14	NDVI spring	17.84	NDVI monsoon	11.46
TRI	8.09	Slope Height	5.59	Catchment Slope	14.25	TRI	8.97	NDVI spring	10.55
NDVI monsoon	7.58	Aspect	5.10	Slope Height	11.83	Mid-slope Position	8.47	Valley Depth	9.75
Catchment Slope	6.73	Catchment Area	4.88	Aspect	10.25	Slope Height	5.44	TRI	7.19

Elevation and one or both bioclimatic variables comprised top predictors of pine species distributions in the terrain-bioclimatic models, based on relative influence (Table 4-A). In the P. arizonica models, the predominance of elevation was especially notable (relative influence = 54%). 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.

Models that included NDVI variables showed different order of importance, e.g., NDVI monsoon was the top predictor for P. engelmannii distribution; all other species had elevation as the top predictor (Table 4-B). NDVI spring was the second most influential predictor in P. strobiformis, P. arizonica, and P. discolor models. Three species had both spring and monsoon NDVI as top three predictors: P. arizonica, P. discolor, and P. engelmannii. Terrain-NDMI results were not reported, due to similar AUC values compared with models that included seasonal NDVI.

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, i.e., P. engelmannii and P. chihuahuana differed in response pattern; P. engelmannii exhibited a more gradual decline as elevation increased and in contrast, P. chihuahuana response declined precipitously at approximately 2000 m.

Species common at higher elevations (P. strobiformis and P. arizonica) were associated with lower values of AHM (Fig. 3; left column). For P. strobiformis, predicted probability of occurrence increased with terrain ruggedness (TRI), and P. arizonica was more likely found on northern and northwesterly aspects (Fig. 3; middle column). In contrast, species common at mid- to lower elevation (P. discolor, P. engelmannii and P. chihuahuana) always included both bioclimatic variables in the top three set of predictors (Fig. 3; left and middle columns). Species presence was greatest in the middle to high range (P. discolor and P. chihuahuana) or lower range (P. engelmannii) of TD.

High probability of occurrence for the pine species had distinct geographic patterns with respect to latitude in both terrain-bioclimatic and terrain-NDVI models (Fig. 4 and Fig. 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 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.

General patterns in spatial predictions were similar, comparing the terrain-bioclimatic results (Fig. 4) with the terrain-NDVI results (Fig. 5) 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). Pine species associated with lower heat-moisture balance (P. strobiformis and P. arizonica; Fig. 4) also showed a strong relationship with Spring NDVI, conversely species that tolerated higher heat-moisture were associated with Monsoon NDVI (P. engelmannii and P. chihuahuana).

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 (Fig. 3 and 6). We identified a wide variety of settings for the formation of microrefugia across the region, described by complex topography within preferred elevational zones (Fig. 6). 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 clusters of high predicted probability, or hotspots, that circumscribe climate refugia essential to pine conservation: 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 (Fig. 4, 5 and 6).

Characteristics of pine microrefugia

Among the pines of the study region, resilience in forecasted warmer and drier climate may vary among co-occurring species (González-Cásares et al., 2017). The varied characteristics of microrefugia will provide important contributions to species resilience where topography slows local rates of climate change (Fig. 6). Each species preferred a specific range of elevation that corresponds to gradients in temperature, moisture and seasonality and regulate physiology and growth (Whittaker and Niering 1965; Barton 1993; Landsberg and Waring, 1997; Guisan and Thuiller 2005). The overlap in affinities across elevation reflects diversification of pine forests and woodlands where multiple species form communities (Fig. 3).

Across elevation, variability in terrain ruggedness, slope position and aspect were important factors in species distribution models (Fig. 3, Online Resource 4) that contribute to defining microrefugia. Specifically, local depressions associated with rugged terrain preferred by P. strobiformis (Fig. 3) are more sheltered from wind and desiccation (Dobrowski 2011). 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.

Distinct thresholds 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, P. chihuahuana probability of occurrence declined precipitously between ~1800 and 2100 m elevation, and notably in a narrowly defined range of temperature difference (Fig. 3). In contrast, P. arizonica exhibited a continuous, but variable, increase in response across elevations above ~1800 m and a steady decline with increase in heat moisture (AHM; Fig. 3). 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 and P. engelmanni was absent above ~2250 m (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, 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).

Microrefugia and seasonal climate response

Remote sensing indices have distinct advantages in capturing habitat characteristics, including differences in 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). In particular, seasonal reflectance characteristics can indicate changes in habitat due to anthropogenic and natural disturbances like stand replacing wildfires, drought, insect outbreaks or logging. Such recent and historical landscape changes influence the current distribution of species (vs potential) but were not otherwise accounted for in our models of microrefugia. Understanding seasonal changes in relation to microrefugia is especially important given the interactions between climate and forest disturbance (Villarreal et al. 2019b; Villarreal et al. 2020).

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 in the low to mid-elevation range exhibit spectral characteristics associated with a more open canopy structure, more specifically, greening of herbaceous understory or co-dominant plants in response to monsoon precipitation. The remote sensing variables were more important for species in some southern islands where fires occur regularly (Villarreal et al. 2020) and 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 monsoon 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 monsoon 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 located in more northerly habitats (Buermann et al. 2008).

Threats to pine refugia

The ability of sky island refugia to support pine species will be tested as climate continues to change. Climate warming is expected to increase drought stress reducing both forest productivity and growth while increasing tree mortality in semi-arid areas of the southwestern United States (Williams et al. 2013). Similarly, climate models predict increasing aridity and more frequent drought spells reducing long-term water availability in northwestern México (Cook and Seager 2013). Changes in fire frequency and severity attendant with climate change may produce conditions in which refugia could lose their functionality (Morelli et al. 2020).

The functioning of climate refugia for pine species in the sky islands is subject to the influence of dominant processes that vary across species distributions (Hampe and Petit 2005; Slaton 2015). Dispersal mechanisms, including long-distance dispersal and ability to tolerate changes in precipitation will influence persistence and movement upward into higher elevation zones for species such as P. chihuahuana and P. engelmannii (Table 1; Table 4). Local adaptations and ability to withstand drought stress will be important to success of northern populations of P. arizonica (Barton et al. 2001; Hess and Fule 2020; Fig. 4 and 5). Seed dispersal by animals offers potential advantages for long-distance dispersal in P. discolor (Copeland et al. 2018). In addition, the association of border pinyon with dry conditions (warmer slope aspects; Online Resource 4; Barton 1993) could infer resilience where drought is a dominant factor.

Regional studies predict shifts in species distributions with climate change across the sky islands. For example, comparisons of recent and historical data provide evidence that mature pine forest communities in the Madrean sky islands and Sierra Madre Occidental have been declining due to climate and land use change (Flesch 2019). Increased habitat loss and habitat patch isolation were forecast for montane species including P. strobiformis using comparisons of current and future climate change scenarios (Yanahan and Moore 2019). By 2080, a large range shift and contraction has been predicted for P. strobiformis across its distribution in the United States and México (Shirk et al. 2018). Increased frequency of high-severity fire has been noted in upper elevations of some northern sky islands (Villarreal et al. 2020), where pine species with limited distributions, e.g., P. strobiformis, favor mesic environments.

Conservation implications

Climate refugia of the past have been important to the formation of global biodiversity hotspots (Myers et al. 2000; Harrison and Noss 2017). In pine forests and woodlands, areas of relatively greater climate stability may provide holdouts and stepping-stones critical to species survival with future change (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 shorter-term movements that are more complex than simply uphill and poleward (Pecl et al. 2017). The large degree of environmental heterogeneity we observed across pine habitats shows promise for facilitating range shift dynamics and buffering of habitats from the negative effects of climate change.

Conservation and management strategies that incorporate refugia could mitigate the dire outcomes for pines forecast by previous studies. The location and effectiveness of refugia will play an important role in both survival of individuals and populations and the ability to redistribute into more favorable habitats through time (Simons-Legaard et al. 2013; Andrews et al. 2018). 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 infers local protection from fire (Iniguez et al. 2008) that allows cohorts of P. arizonica to establish within longer fire intervals (Iniguez et al. 2016).

Disturbance plays an important role in determining the ability of trees to disperse to favorable locations (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 survival of seedlings (Barton 1993). Species traits that promote seed dispersal, including by wind and animals, will be key to determining redistribution for particular species. 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).

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 for collaborative conservation efforts across ownership and political boundaries, based on areas with high probability of species occurrence identified in our study. Fire management has been identified as a unifying management theme in the Madrean sky islands (Villarreal et al. 2019a; Villarreal et al. 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. Moreover, activities to assist reforestation and migration are increasing in response to changing disturbance regimes that accompany climate change (Carbajal-Navarro et al. 2019). Spatial models of pine refugia will provide key data for prioritizing the allocation of resources for these efforts and thus increase the likelihood of success.

Funding was provided by the USGS Land Change Science Program.

Conflicts of interest/competing interests.

Not applicable.

Ethics approval.

Not applicable.

Consent to participate.

Not applicable.

Consent for publication.

Not applicable.

Availability of data and material (data transparency).

Data used in the analysis are available online (https://github.com/HaireLab/madrean-pines; DOI will be assigned via Zenodo upon publication of the manuscript). Spatial data layers of predicted probability and remote sensing indices will be published on ScienceBase and accessible through Data Basin (https://databasin.org/galleries/b27fbb9c42b849fe9029827095449d07/).

Code availability.

Computer (R) code is available online (https://github.com/HaireLab/madrean-pines).

Authors’ contributions.

All authors contributed to the study conception and design. Field data collection was done by Aaron D. Flesch, Jose M. Iniguez and Jamie S. Sanderlin. Data preparation was accomplished Sandra L. Haire and Miguel L. Villarreal; analysis was done by Sandra L. Haire. The first draft of the manuscript was written by Sandra L. Haire and Miguel L. Villarreal and all authors commented on previous versions of the manuscript. All authors read and approved the final manuscript.

Aburto-Oropeza, O, Johnson, A. F., Agha, M., Allen, E. B., Allen, M. F., González, J. A., … Taylor, J. E. (2018). Harnessing cross-border resources to confront climate change. Environmental Science and Policy, 87, 128–132. https://doi.org/10.1016/J.ENVSCI.2018.01.001
Adams, D. K., Comrie, A. C. (1997). The North American monsoon. Bulletin of the American Meteorological Society, 78(10), 2197-2214.
AdaptWest Project 2015. Gridded current and projected climate data for North America at 1km resolution, interpolated using the ClimateNA v5.10 software. https://adaptwest.databasin.org/pages/adaptwest-climatena
Aguirre-Gutiérrez, J., Serna-Chavez, H. M., Villalobos-Arambula, A. R., Pérez de la Rosa, J. A., Raes, N. (2015). Similar but not equivalent: Ecological niche comparison across closely-related Mexican white pines. Diversity and Distributions, 21(3), 245–257. https://doi.org/10.1111/ddi.12268
Andrews, C., Weiskittel, A., D’Amato, A. W., Simons-Legaard, E. (2018). Variation in the maximum stand density index and its linkage to climate in mixed species forests of the North American Acadian Region. Forest Ecology and Management, 417, 90-102.
Araújo, M. B., Guisan, A. (2006). Five (or so) challenges for species distribution modelling. Journal of Biogeography, 33(10), 1677–1688. https://doi.org/10.1111/j.1365-2699.2006.01584.x
Ashcroft, M. B. (2010). Identifying refugia from climate change. Journal of Biogeography, 37(8), 1407–1413. https://doi.org/10.1111/j.1365-2699.2010.02300.x
Ashcroft, M. B., Chisholm, L. A., French, K. O. (2009). Climate change at the landscape scale: predicting fine-grained spatial heterogeneity in warming and potential refugia for vegetation. Global Change Biology, 15(3), 656–667. https://doi.org/https://doi.org/10.1111/j.1365-2486.2008.01762.x
Ashcroft, M. B., Gollan, J. R., Warton, D. I., Ramp, D. (2012). A novel approach to quantify and locate potential microrefugia using topoclimate, climate stability, and isolation from the matrix. Global Change Biology, 18(6), 1866–1879. https://doi.org/10.1111/j.1365-2486.2012.02661.x
Ávila-Flores, I. J., Hernández-Díaz, J. C., González-, M. S., Prieto-Ruíz, J. Á., Wehenkel, C. (2016). Pinus engelmannii Carr. in Northwestern Mexico : A review. Pakistan Journal of Botany, 48(5), 2159–2166.
Barton, A. M. (1993). Factors controlling plant distributions: Drought, competition, and fire in montane pines in Arizona. Ecological Monographs, 63(4), 367–397. https://doi.org/10.2307/2937151
Barton, A. M., Poulos, H. M. (2018). Pine vs. oaks revisited: Conversion of Madrean pine-oak forest to oak shrubland after high-severity wildfire in the sky islands of Arizona. Forest Ecology and Management, 414, 28–40. https://doi.org/https://doi.org/10.1016/j.foreco.2018.02.011
Barton, A. M., Swetnam, T. W., Baisan, C. H. (2001). Arizona pine (Pinus arizonica) stand dynamics: local and regional factors in a fire-prone madrean gallery forest of Southeast Arizona, USA. Landscape Ecology, 16(4), 351–369. https://doi.org/10.1023/A:1011189408651
Boyko, H. (1947). On the role of plants as quantitative climate indicators and the geo-ecological law of distribution. Journal of Ecology 35: 138–157.
Bucholz, E. R., Waring, K. M., Kolb, T. E., Swenson, J. K., Whipple, A. V. (2020). Water relations and drought response of Pinus strobiformis. Canadian Journal of Forest Research, 50(9), 905–916. https://doi.org/10.1139/cjfr-2019-0423
Buermann, W., Saatchi, S., Smith, T. B., Zutta, B. R., Chaves, J. A., Milá, B., Graham, C. H. (2008). Predicting species distributions across the Amazonian and Andean regions using remote sensing data. Journal of Biogeography, 35(7), 1160–1176. https://doi.org/10.1111/j.1365-2699.2007.01858.x
Carbajal-Navarro, A., Navarro-Miranda, E., Blanco-García, A., Cruzado-Vargas, A. L., Gómez-Pineda, E., Zamora-Sánchez, C., … Sáenz-Romero, C. (2019). Ecological restoration of Abies religiosa forests using nurse plants and assisted migration in the Monarch Butterfly Biosphere Reserve, Mexico. Frontiers in Ecology and Evolution, 7, 421. https://doi.org/10.3389/fevo.2019.00421
Coblentz, D. (2005). The Tectonic Evolution of the Madrean Archipelago and Its Impact on the Geoecology of the Sky Islands. Connecting Mountain Islands and Desert Seas: Biodiversity and Management of the Madrean Archipelago II. USDA Forest Service Proceedings RMRS-P-36, 62–68.
Coblentz, D. D., & Riitters, K. H. (2004). Topographic controls on the regional-scale biodiversity of the south-western United States. Journal of Biogeography, 31(7), 1125–1138. https://doi.org/10.1111/j.1365-2699.2004.00981.x
Conrad, O., Bechtel, B., Bock, M., Dietrich, H., Fischer, E., Gerlitz, L., Wehberg, J., Wichmann, V., Böhner, J. (2015): System for Automated Geoscientific Analyses (SAGA) v. 2.1.4. Geosci. Model Dev., 8, 1991-2007, doi:10.5194/gmd-8-1991-2015. http://www.saga-gis.org/
Conservation International (2020). Biodiversity hotspots: Madrean pine-oak woodlands. Website: https://www.cepf.net/our-work/biodiversity-hotspots/madrean-pine-oak-woodlands
Cook, B. I., Seager, R. (2013). The response of the North American Monsoon to increased greenhouse gas forcing. Journal of Geophysical Research Atmospheres, 118(4), 1690–1699. https://doi.org/10.1002/jgrd.50111
Copenhaver-Parry, P. E., Shuman, B. N., Tinker, D. B. (2017). Toward an improved conceptual understanding of North American tree species distributions. Ecosphere, 8(6), e01853. https://doi.org/10.1002/ecs2.1853
Copeland, S. M., Bradford, J. B., Duniway, M. C., Butterfield, B. J. (2018). Life history characteristics may be as important as climate projections for defining range shifts: An example for common tree species in the intermountain western US. Diversity and Distributions, 24(12), 1844–1859. https://doi.org/10.1111/ddi.12813
Cord, A. F., Klein, D., Gernandt, D. S., de la Rosa, J. A. P., Dech, S. (2014). Remote sensing data can improve predictions of species richness by stacked species distribution models: a case study for Mexican pines. Journal of Biogeography, 41(4), 736–748. https://doi.org/10.1111/jbi.12225
Crimmins, S. M., Dobrowski, S. Z., Greenberg, J. A., Abatzoglou, J. T., Mynsberge, A. R. (2011). Changes in climatic water balance drive downhill shifts in plant species optimum elevations. Science, 331(6015), 324–327. https://doi.org/10.1126/science.1199040
Deyo, N. S., Van Devender, T. R., Smith, A., Gilbert, E. (2013). Documenting the biodiversity of the Madrean Archipelago: An analysis of a virtual flora and fauna. In: Gottfried, G.J., Ffolliott, P. F., Gebow, Brooke S., Eskew, L. G., Collins, L. C. Proceedings. RMRS-P-67. Fort Collins, CO: U.S. Department of Agriculture, Forest Service, Rocky Mountain Research Station, p. 292-299.
Dobrowski, S. Z. (2011). A climatic basis for microrefugia: the influence of terrain on climate. Global Change Biology, 17(2), 1022–1035. https://doi.org/10.1111/j.1365-2486.2010.02263.x
Elith, J., Leathwick, J. R., Hastie, T. (2008). A working guide to boosted regression trees. The Journal of Animal Ecology, 77(4), 802–813.
Farr, T. G., Rosen, P. A., Caro, E., Crippen, R., Duren, R., Hensley, S., … Alsdorf, D. (2007). The Shuttle Radar Topography Mission. Reviews of Geophysics, 45(2). https://doi.org/https://doi.org/10.1029/2005RG000183
Ferguson, G. M., Flesch, A. D., Van Devender, T. R. (2013). Biogeography and diversity of pines in the Madrean Archipelago. In: Gottfried, G.J., Ffolliott, P. F., Gebow, Brooke S., Eskew, L. G., Collins, L. C. Proceedings. RMRS-P-67. Fort Collins, CO: U.S. Department of Agriculture, Forest Service, Rocky Mountain Research Station, p. 197–203.
Findley, J.S. 1996. Mammalian biogeography in the American Southwest. Pp. 297-308 in Genoways, H.H., and R.J. Baker (editors), Festschrift for J. Knox Jones, Jr., Texas Tech University Press, Lubbock, Texas.
Flesch, A.D., C.W. Epps, J.W. Cain, M. Clark, P.R. Krausman, J.R. Morgart. 2010. Potential effects of the United States-Mexico border fence on wildlife. Conservation Biology 24:171-181.
Flesch AD (2019) Patterns and drivers of long-term changes in breeding bird communities in a global biodiversity hotspot in Mexico. Diversity and Distributions 00:1-15.
Flesch AD, Gonzalez Sanchez C, Valenzuela Amarillas J (2016) Abundance and habitat relationships of breeding birds in the sky islands and adjacent Sierra Madre Occidental of northwest Mexico. Journal of Field Ornithology 87:176–195. Doi: 10.1111/jofo.1215
Gamon, J., Field, C., Goulden, M., Griffin, K., Hartley, A., Joel, G., … Valentini, R. (1995). Relationships between NDVI, canopy structure, and photosynthesis in three Californian vegetation types. Ecological Applications, 5, 28–41. https://doi.org/10.2307/1942049
Ganey, J.L., Block, W.M., Boucher, P.F. (1996). Effects of fire on birds in Madrean forests and woodlands. In: Ffolliott, P.F.; DeBano, L.F.; Baker, M.B.; [and others], tech coords. Gen. Tech. Rep. RM-GTR-289. U.S. Department of Agriculture, Forest Service, Rocky Mountain Forest and Range Experiment Station. Fort Collins, CO, Pp. 146-154
Ganey, J. L, Iniguez, J. M., Sanderlin, J. S., Block, W. M. (2017). Developing a monitoring program for bird populations in the Chiricahua Mountains, Arizona, using citizen observers: Initial stages. Gen. Tech. Rep. RMRS-GTR-368. Fort Collins, CO: U.S, Department of Agriculture, Forest Service, Rocky Mountain Research Station. 30 p.
Gao B.C. (1996) NDWI—A normalized difference water index for remote sensing of vegetation liquid water from space. Remote Sensing of Environment. Dec 1;58(3):257-66.
Gómez-Mendoza, L., Arriaga, L. (2007). Modeling the effect of climate change on the distribution of oak and pine species of Mexico. Conservation Biology, 21(6), 1545–1555. https://doi.org/10.1111/j.1523-1739.2007.00814.x
González-Cásares, M., Pompa-García, M., Camarero, J. J. (2017). Differences in climate–growth relationship indicate diverse drought tolerances among five pine species coexisting in Northwestern Mexico. Trees – Structure and Function, 31(2), 531–544. https://doi.org/10.1007/s00468-016-1488-0
González-Elizondo, M. S., González-Elizondo, M., Tena-Flores, J. A., Ruacho-Gonzalez, L., Lopez-Enriquez, We. L. (2012). Vegetation of the Sierra Madre Occidental, Mexico: a synthesis. Acta Botanica Mexicana 100, 351-404.
Gorelick, N., Hancher, M., Dixon, M., Ilyushchenko, S., Thau, D., Moore, R. (2017). Google Earth Engine: planetary-scale geospatial analysis for everyone. Remote Sensing of Environment 202: 18–27 https://doi.org/10.1016/j.rse.2017.06.031.
Greenwell, B. (2017). Pdp: An R package for constructing partial dependence plots. The R Journal 9(1): 421-436. https://journal.r-project.org/archive/2017/RJ-2017-016/index.html
Greenwell, B., Boehmke, B., Cunningham J. and GBM Developers (2019). gbm: Generalized boosted regression models. R package version 2.1.5. https://CRAN.R-project.org/package=gbm
Guillera-Arroita, G., Lahoz-Monfort, J. J., Elith, J., Gordon, A., Kujala, H., Lentini, P. E., … Wintle, B. A. (2015). Is my species distribution model fit for purpose? Matching data and models to applications. Global Ecology and Biogeography, 24(3), 276–292. https://doi.org/10.1111/geb.12268
Guisan, A., Thuiller, W. (2005). Predicting species distribution: offering more than simple habitat models. Ecology Letters, 8(9), 993–1009. https://doi.org/10.1111/j.1461-0248.2005.00792.x
Haire, S. L., Coop, J. D., Miller, C. (2017). Characterizing spatial neighborhoods of refugia following large fires in northern New Mexico USA. Land, 6(1). https://doi.org/10.3390/land6010019
Hampe, A., Jump, A. S. (2011). Climate relicts: past, present, future. Annual Review of Ecology, Evolution and Systematics, 42, 313–333.
Hampe, A., Petit, R. J. (2005). Conserving biodiversity under climate change: the rear edge matters. Ecology Letters, 8(5), 461–467. https://doi.org/10.1111/j.1461-0248.2005.00739.x
Hannah, L., Flint, L., Syphard, A. D., Moritz, M. A., Buckley, L. B., McCullough, I. M. (2015). Place and process in conservation planning for climate change: a reply to Keppel & Wardell-Johnson. Trends in Ecology and Evolution, 30(5), 234–235. https://doi.org/10.1016/j.tree.2015.03.008
Harrison, S., Noss, R. (2017). Endemism hotspots are linked to stable climatic refugia. Annals of Botany, 119(2), 207–214. https://doi.org/10.1093/aob/mcw248
Hess, V. A., Fulé, P. Z. (2020). Is a Mexican pine species better adapted to the warming climate of the Southwestern USA? Frontiers in Forests and Global Change, 3, 60.
Hijmans, R.J. (2019). raster: Geographic data analysis and modeling. R package version 2.9-5. https://CRAN.R-project.org/package=raster
Hijmans, R.J., Phillips, S., Leathwick J., Elith J. (2017). dismo: Species distribution modeling. R package version 1.1-4. https://CRAN.R-project.org/package=dismo
Hoylman, Z. H., Jencso, K. G., Hu, J., Holden, Z. A., Allred, B., Dobrowski, S., … Seielstad, C. (2019). The topographic signature of ecosystem climate sensitivity in the western United States. Geophysical Research Letters, 46(24), 14508–14520. https://doi.org/https://doi.org/10.1029/2019GL085546
Iniguez, J.M., Swetnam, T.W., Yool, S. (2008). Topography affected landscape fire history patterns in southern Arizona, USA. Forest Ecology and Management 256:295-303.
Iniguez, J. M., Swetnam, T. W., Baisan, C. H. (2016). Fire history and moisture influences on historical forest age structure in the sky islands of southern Arizona, USA. Journal of Biogeography, 43(1), 85–95. https://doi.org/10.1111/jbi.12626
Jiménez-Valverde, A. (2012). Insights into the area under the receiver operating characteristic curve (AUC) as a discrimination measure in species distribution modelling. Global Ecology and Biogeography, 21(4), 498–507. https://doi.org/10.1111/j.1466-8238.2011.00683.x
Jins, V. J., Panigrahi, M., Jayapal, R., Bishop, T. R. (2021). Elevational gradients of reptile richness in the southern Western Ghats of India: Evaluating spatial and bioclimatic drivers. Biotropica. https://doi.org/10.1111/btp.12878
Keppel, G., Van Niel, K. P., Wardell-Johnson, G. W., Yates, C. J., Byrne, M., Mucina, L., … Franklin, S. E. (2012). Refugia: Identifying and understanding safe havens for biodiversity under climate change. Global Ecology and Biogeography, 21(4), 393–404. https://doi.org/10.1111/j.1466-8238.2011.00686.x
Keppel, G., Ottaviani, G., Harrison, S., Wardell-Johnson, G. W., Marcantonio, M., & Mucina, L. (2018). Towards an eco-evolutionary understanding of endemism hotspots and refugia. Annals of Botany, mcy173. https://doi.org/10.1093/aob/mcy173
Keppel, G., Wardell-Johnson, G. W. (2015). Refugial capacity defines holdouts, microrefugia and stepping-stones: a response to Hannah et al. Trends in Ecology & Evolution, 30(5), 233–234. https://doi.org/10.1016/j.tree.2015.01.008
Krawchuk, M. A., Haire, S. L., Coop, J., Parisien, M.-A., Whitman, E., Chong, G., Miller, C. (2016). Topographic and fire weather controls of fire refugia in forested ecosystems of northwestern North America. Ecosphere, 7(12). https://doi.org/10.1002/ecs2.1632
Krawchuk, M.A., Meigs, G.W., Cartwright, J.M., Coop, J.D., Davis, R., Holz, A, Kolden, C., Meddens, A.J.H. (2020). Disturbance refugia within mosaics of forest fire, drought, and insect outbreaks. Frontiers in Ecology and the Environment, 18(5):235-244.
Krishnadas, M., Osuri, A. M. (2020). Environment shapes the spatial organization of tree diversity in fragmented forests across a human-modified landscape. Ecological Applications. https://doi.org/10.1002/eap.2244
Landsberg J.J., Waring, R.H. (1997) A generalized model of forest productivity using simplified concepts of radiation-use efficiency, carbon balance and partitioning. Forest Ecology and Management, 95:209-228
Lawler, J. (2009). Climate change adaptation strategies for resource management and conservation planning. Annals of the New York Academy of Sciences, 1162, 79–98. https://doi.org/10.1111/j.1749-6632.2009.04147.x
Looney, C. E., Waring, K. M. (2013). Pinus strobiformis (southwestern white pine) stand dynamics, regeneration, and disturbance ecology: A review. Forest Ecology and Management, 287, 90–102. https://doi.org/10.1016/j.foreco.2012.09.020
Luo, D., Yue, J.-P., Sun, W.-G., Xu, B., Li, Z.-M., Comes, H. P., Sun, H. (2016). Evolutionary history of the subnival flora of the Himalaya-Hengduan Mountains: first insights from comparative phylogeography of four perennial herbs. Journal of Biogeography, 43(1), 31–43. https://doi.org/10.1111/jbi.12610
Mackey, B., Berry, S., Hugh, S., Ferrier, S., Harwood, T. D., Williams, K. J. (2012). Ecosystem greenspots: identifying potential drought, fire, and climate-change micro-refuges. Ecological Applications, 22(6), 1852–1864. Retrieved from http://www.ncbi.nlm.nih.gov/pubmed/23092021
Mittermeier, R. A., Robles Gil, P., Hoffman, M., Pilgrim, J., Brooks, T., Mittermeier, C. G., … da Fonseca, G. A. B. (2004). Hotspots revisited: Earth’s biologically richest and most endangered ecoregions. Mexico City, Mexico: CEMEX.
Morelli, T. L., Daly, C., Dobrowski, S. Z., Dulen, D. M., Ebersole, J. L., Jackson, S. T., … Beissinger, S. R. (2016). Managing climate change refugia for climate adaptation. PLoS ONE, 11(8), 1–17. https://doi.org/10.1371/journal.pone.0159909
Morelli, T. L., Barrows, C. W., Ramirez, A. R., Cartwright, J. M., Ackerly, D. D., Eaves, T. D., … Thorne, J. H. (2020). Climate-change refugia: biodiversity in the slow lane. Frontiers in Ecology and the Environment, 18(5), 228–234. https://doi.org/10.1002/fee.2189
Myers, N., Mittermeler, R. A., Mittermeler, C. G., Da Fonseca, G. A. B., Kent, J. (2000). Biodiversity hotspots for conservation priorities. Nature, 403(6772), 853–858. https://doi.org/10.1038/35002501
Nogue, S., Rull, V., Vegas-Vilarrubia, T. (2013). Elevational gradients in the neotropical table mountains: patterns of endemism and implications for conservation. Diversity and Distributions, 19(7), 676–687. https://doi.org/10.1111/ddi.12017
Norman, L. M. (2020). Ecosystem services of riparian restoration: A review of rock detention structures in the Madrean Archipelago Ecoregion. Air, Soil and Water Research, 13. https://doi.org/10.1177/1178622120946337
Omernik, James M. 1987. Ecoregions of the Conterminous United States. Annals of the Association of American Geographers 77 (1): 118–25. https://doi.org/10.1111/j.1467-8306.1987.tb00149.x.
Pecl, G. T., Araújo, M. B., Bell, J. D., Blanchard, J., Bonebrake, T. C., Chen, I.-C., … Williams, S. E. (2017). Biodiversity redistribution under climate change: Impacts on ecosystems and human well-being. Science, 355(6332), eaai9214. https://doi.org/10.1126/science.aai9214
Peters, R., M. Moskwik, W.J. Ripple, C. Wolf, G. Carreón-Arroyo, G. Ceballos, A. Córdova, R. Dirzo, P.R. Ehrlich, A.D. Flesch, R. List, T.E. Lovejoy, R.F. Noss, J. Pacheco, J.K. Sarukhán, M.E. Soulé, E.O. Wilson, J.R.B. Miller. (2018). Nature divided, scientists united: U.S.-Mexico border wall threatens biodiversity and binational conservation. BioScience. DOI:10.1093/biosci/biy063
Pezzoli, K., Kozo, J., Ferran, K., Wooten, W., Gomez, G. R., Al-Delaimy, W. K. (2014). One Bioregion/One Health: An integrative narrative for transboundary planning along the United States-Mexico Border. Global Society : Journal of Interdisciplinary International Relations, 28(4), 419–440. https://doi.org/10.1080/13600826.2014.951316
Popp, M., Gizaw, A., Nemomissa, S., Suda, J., Brochmann, C. (2008). Colonization and diversification in the African `sky islands’ by Eurasian Lychnis L. (Caryophyllaceae). Journal of Biogeography, 35(6), 1016–1029. https://doi.org/10.1111/j.1365-2699.2008.01902.x
Poulos, H. M., Barton, A. M., Berlyn, G. P., Schwilk, D. W., Faires, C. E., McCurdy, W. C. (2020). Differences in leaf physiology among juvenile pines and oaks following high-severity wildfire in an Arizona Sky Island Mountain range. Forest Ecology and Management, 457, 117704. https://doi.org/https://doi.org/10.1016/j.foreco.2019.117704
Rehfeldt, G. E., Crookston, N. L., Sáenz-Romero, C., Campbell, E. M. (2012). North American vegetation model for land-use planning in a changing climate: a solution to large classification problems. Ecological Applications, 22(1), 119–141. https://doi.org/10.1890/11-0495.1
Rodríguez-Trejo, D. A., Fulé, P. Z. (2003). Fire ecology of Mexican pines and a fire management proposal. International Journal of Wildland Fire, 12(1), 23–37. https://doi.org/10.1071/WF02040
Romero-Sanchez, M., Perez-Miranda, R., Gonzalez-Hernandez, A., Velasco-Garcia, M., Velasco-Bautista, E., Flores, A. (2018). Current and potential spatial distribution of six endangered pine species of Mexico: Towards a conservation strategy. Forests, 9(12), 767. https://doi.org/10.3390/f9120767
Romme, W. H., Everham, E. H., Frelich, L. E., Moritz, M. A., Sparks, R. E. (1998). Are Large, infrequent disturbances qualitatively different from small, frequent disturbances? Ecosystems, 1(6), 524–534. https://doi.org/10.1007/s100219900048
Safont, E., Vegas-Vilarrubia, T., Rull, V., Holst, B. K., Huber, O., Nozawa, S., … Silva, A. (2016). Plant communities and environmental factors in the Guayana Highlands: monitoring for conservation under future climate change. Systematics and Biodiversity, 14(4), 327–344. https://doi.org/10.1080/14772000.2015.1134700
Sanderlin JS, Block WM, Ganey JL, Iniguez JM, Cushman S (2014) Assessing large-scale effects of wildfire and climate change on avian communities and habitats in the sky islands, Arizona. Final Report for U.S. Fish & Wildlife’s Desert Landscape Conservation Cooperative.
Sanderlin, J. S., Block, W. M., Ganey, J. L., Iniguez, J. M. (2013). Preliminary Assessment of Species Richness and Avian Community Dynamics in the Madrean sky islands, Arizona. Merging Science and Management in a Rapidly Changing World: Biodiversity and Management of the Madrean Archipelago III. Proceedings. RMRS-P-67, 180–190.
Serra-Diaz, J. M., Scheller, R. M., Syphard, A. D., Franklin, J. (2015). Disturbance and climate microrefugia mediate tree range shifts during climate change. Landscape Ecology, 30(6), 1039–1053. https://doi.org/10.1007/s10980-015-0173-9
Scherrer, D., Körner, C., 2011. Topographically controlled thermal-habitat differentiation buffers alpine plant diversity against climate warming. Journal of Biogeography 38, 406–416.
Sheppard, P.R., Comrie, A. C., Packin, G. D., Angersbach, K., Hughes, M.K. (2002). The climate of the US Southwest. Climate Research, 21:219-238. doi:10.3354/cr021219.
Shirk, A. J., Cushman, S. A., Waring, K. M., Wehenkel, C. A., Leal-Sáenz, A., Toney, C., Lopez-Sanchez, C. A. (2018). Southwestern white pine (Pinus strobiformis) species distribution models project a large range shift and contraction due to regional climatic changes. Forest Ecology and Management, 411(October 2017), 176–186. https://doi.org/10.1016/j.foreco.2018.01.025
Simons-Legaard, E., D’Amato, A.W., Legaard, K., Sturtevant,B., & Weiskittel, A. (2013). Future distribution and productivity of spruce-fir forests under climate change: A comparison of the Northeast and the Lake States. Northeastern States Research Cooperative Final Report. https://nsrcforest.org/project/future-distribution-and-productivity-spruce-fir-under-climate-change
Slaton, M. (2015). The roles of disturbance, topography and climate in determining the leading and rear edges of population range limits. Journal of Biogeography, 42. https://doi.org/10.1111/jbi.12406
Spector, S., Spector, S. (2002). Biogeographic crossroads: priority areas for conservation. Conservation Biology, 16(6), 1480–1487.
Steinbauer, M. J., Field, R., Grytnes, J.-A., Trigas, P., Ah-Peng, C., Attorre, F., … Beierkuhnlein, C. (2016). Topography-driven isolation, speciation and a global increase of endemism with elevation. Global Ecology and Biogeography, 25(9), 1097–1107. https://doi.org/10.1111/geb.12469
Travers-Smith, H. Z., Lantz, T. C. (2020). Leading-edge disequilibrium in alder and spruce populations across the forest–tundra ecotone. Ecosphere, 11(7). https://doi.org/10.1002/ecs2.3118
Vegas-Vilarrubia, T., Nogue, S., Rull, V. (2012). Global warming, habitat shifts and potential refugia for biodiversity conservation in the neotropical Guayana Highlands. Biological Conservation, 152, 159–168. https://doi.org/10.1016/j.biocon.2012.03.036
Villarreal, M. L., Haire, S.L., Bravo, J. C., Norman, L.M. (2019a). A mosaic of land tenure and ownership creates challenges and opportunities for transboundary conservation in the United States-Mexico borderlands. Case Studies in the Environment, 3 (1): 1–10. Doi: https://doi.org/10.1525/cse.2019.002113
Villarreal, M. L., Haire, S. L., Iniguez, J. M., Cortés Montaño, C., Poitras, T. B. (2019b). Distant neighbors: recent wildfire patterns of the Madrean Sky Islands of southwestern United States and northwestern Mexico. Fire Ecology, 15(1), 2. https://doi.org/10.1186/s42408-018-0012-x
Villarreal, M. L., Iniguez, J. M., Flesch, A. D., Sanderlin, J. S., Montaño, C. C., Conrad, C. R., Haire, S. L. (2020). Contemporary fire regimes provide a critical perspective on restoration needs in the Mexico-United States Borderlands. Air, Soil and Water Research, 13(1).
Villarreal, M. L., Norman, L. M., Boykin, K. G., Wallace, C. S. (2013). Biodiversity losses and conservation trade-offs: assessing future urban growth scenarios for a North American trade corridor. International Journal of Biodiversity Science, Ecosystem Services & Management, 9(2), 90-103.
Villarreal, M. L., Petrakis, R. E. (2014). Conflation and aggregation of spatial data improve predictive models for species with limited habitats: A case of the threatened yellow-billed cuckoo in Arizona, USA. Applied Geography, 47, 57-69.
Wang, T., Hamann, A., Spittlehouse, D., Carroll, C. (2016). Locally downscaled and spatially customizable climate data for historical and future periods for North America. PLOS ONE, 11(6), 1–17. https://doi.org/10.1371/journal.pone.0156720.
Warshall, P. 1995. The Madrean Sky Island Archipelago: A Planetary Overview. In Biodiversity and Management of the Madrean Archipelago: The sky islands of Southwestern United States and Northwestern Mexico. General Technical Report RM-GTR-264. Tucson, Arizona: U.S. Forest Service.
Whittaker, R. H., Niering, W. A. (1965). Vegetation of the Santa Catalina Mountains, Arizona: a gradient analysis of the south slope. Ecology, 46(4), 429-452.
Williams, A. P., Allen, C. D., Macalady, A. K., Griffin, D., Woodhouse, C. A., Meko, D. M., … McDowell, N. G. (2013). Temperature as a potent driver of regional forest drought stress and tree mortality. Nature Climate Change, 3(3), 292–297. https://doi.org/10.1038/nclimate1693
Yanahan, A. D., Moore, W. (2019). Impacts of 21st‐century climate change on montane habitat in the Madrean Sky Island Archipelago. Diversity and Distributions 25: 1625-1638. https://doi.org/10.1111/ddi.12965.

OnlineResource1.docx

Download PDF

Journal Publication

published 01 Apr, 2022

Read the published version in Plant Ecology →

Version 1

posted

You are reading this latest preprint version

Climate Refugia for Pinus Spp. In Topographic and Bioclimatic Environments of the Madrean Sky Islands of México and the United States

Status:

Journal Publication

Version 1

Abstract

Figures

Introduction

Methods

Results

Discussion

Declarations

References

Supplementary Files

Status:

Journal Publication

Version 1