Multiscale ecological niche modeling exhibits varying climate change impacts on habitat suitability of Madrean Pine-Oak trees

Anthropogenic climate change and increasing greenhouse gas emissions are expected to globally impact the biological function, community structure, and spatial distribution of biodiversity. Many existing studies explore the effect of climate change on biodiversity, generally at a single spatial scale. This study explores the potential effects of climate change on the habitat suitability of seven tree species at two distinct spatial scales: the Coronado National Forest (CNF), a local management area, and the Sierra Madre Occidental (SMO), an ecoregional extent. Habitat suitability was determined by extrapolating Ecological Niche Models (ENMs) based on citizen-science tree occurrence records into future climatic conditions using projected 30-year normals for two anthropogenic emissions scenarios through the end of the century. These ENMs, examined at a spatial resolution of 1 km2, are constructed using a mean average ensemble of three commonly used machine learning algorithms. The results show that habitat suitability is expected to decrease for all seven tree species at varying degrees. Results also show that climate-forcing scenario choice appears to be far less important for understanding changes in species habitat suitability than the spatial scale of modeling extent. Additionally, we observed non-linear changes in tree species habitat suitability within the SMO and CNF dependent on forest community type, latitude, and elevational gradient. The paper concludes with a discussion of the necessary steps to verify the estimated alters of these tree species under climate change. Most importantly, provides a framework for characterizing habitat suitability across spatial scales.


Introduction
Observed changes in Earth's physical and biological systems are directly related to the forcings of anthropogenic climate change (Rosenzweig et al., 2008;Hansen and Stone, 2016). Increasing greenhouse gas concentrations are expected to impact the world's climate through rising temperatures, prolonged drought, and increasing contrasts between wet and dry seasons (Yanahan and Moore, 2019). These human-driven alterations are negatively impacting biodiversity globally (Kannan and James, 2009). These observed impacts include the progressive decoupling of species interactions (i.e., plants and pollinators), extinctions of local populations along range boundaries, and the extension of species geographic limits toward the poles or higher elevations (Hughes, 2000). Due to global-scale impacts of anthropogenic climate change, it is likely that all species may encounter varying alterations to their livelihood, habitat, and geographic range limits.
For biodiversity occupying sensitive habitats (i.e., mountains), climate change poses a significant threat to species physiology, distribution, and community structure (Hao et al., 2007;Draper et al., 2019). Montane plants and plant communities, referred to as the sentinels of climate change (Malanson et al., 2019), often respond sooner to environmental changes than other species and communities. Additionally, temperature increases caused by climate change are more significant at higher elevations (Chen et al., 2011;Pepin et al., 2015). While some species perish, others survive by in situ adaptation or relocation to a higher or lower elevation (Colwell and Rangel, 2009). This observation varies among species as the ability of a species to survive environmental change is determined by their specific biological traits. Plant traits such as propagule production and dispersal, post-fire regeneration, and shade intolerance drive the relative ability to out-compete other species in any given location (Alexander et al., 2018). This trend may lead to decreased suitable habitat for high elevation plant species and communities. These communities may become increasingly fragmented as small, resilient patches of remnant vegetation and become stranded among new vegetative communities (Coe et al., 2012;Yanahan and Moore, 2019).
There are numerous studies, often model based, that examine how species habitat suitability is altering due to climatic change (Morin and Thuiller, 2009;Elith et al., 2010;Iverson and McKenzie, 2013). However, there is little reported research examining the importance of modeling the effects of climate change on plant species at various spatial scales. This study fills this gap in the literature by comparing the sensitivity of habitat suitability modeling to both changes in modeling domain and changes in the climate. The region of comparison is the Coronado National Forest (CNF), a topographically and ecologically complex management area.
While land managers of the CNF often focus primarily on the ecological occurrences within the forest boundaries, ecological function rarely follows geopolitical boundaries. Thus, in addition to examining species responses to anthropogenic climate change at the local scale, we also examine responses at the regional scale. The CNF represents the northernmost limits of the Sierra Madre Occidental (SMO), a mountain range that extends across the United States and Mexico, comprising the Madrean archipelago and mountains of Sonora, Mexico. By incorporating a regional scale for comparison, we develop a separate scenario and alternative evidence for alterations observed or not observed at the local scale. Regional models have been commonly linked to their larger counterparts, the continental and global scale (Gallien et al., 2012;Taucare-Ríos et al., 2016).
The primary objective of this study is to examine how anthropogenic climate change may alter the spatial distribution and quality of suitable habitat for seven common tree species (i.e., four pine, three oak) that provide the general structure of Madrean Pine-Oak Woodlands found throughout the SMO and CNF. To predict the geographic response of each dominant species to anthropogenic climate change, we employ Ecological Niche Modeling (ENM), a spatially explicit modeling approach to predict species habitat suitability across space using a mean average ensemble of three commonly used machine learning statistical models. We fit models with environmental variables that reflect current conditions along with future projections of climate under two anthropogenic emissions scenarios, formulated in 30-year climate normals for three time periods: 2040, 2070, and 2,100. These scenarios are based on social analogs, that range from moderate to severe in greenhouse gas emissions. Niche models are then developed for each tree species for current (2022) conditions, along with the three future climate conditions under the two emissions scenarios, at the two spatial scales of interest: Local and Regional. Additionally, we explore habitat suitability change for each tree species across elevation and latitude within the CNF using Generalized Additive Models (GAMs). We anticipate that habitat suitability will decrease in geographic range and in quality (characterized by the maximum values of habitat suitability) for all species examined in both climate change scenario. We anticipate our proposed framework will provide an important template for local-scale and regional-scale explorations of climateinduced habitat suitability change and should assist land managers who are anticipating the impacts of anthropogenic climate change.

The local and regional scale
This study uniquely examines two separate spatial scales with the smaller scale being a subset of the larger. This purpose serves to examine how species are predicted to respond to climate change based on the extent and range of environmental variables you incorporate. The regional scale is represented by the SMO ( Figure 1A). The SMO is known for its high biological and cultural diversity and enormous environmental and economic importance (González-Elizondo et al., 2013). At approximately 290,000 km 2 , the SMO is home to a high diversity of flora, including 24 different species of Pinus, 54 species of Quercus, and many other trees species (Silva-Flores et al., 2014). Due to such a large geographic range of forested land, the SMO presents challenges in the management of species, and is often filtered down to distinct management zones.
The local-scale, CNF, extends to Southeastern Arizona and Southwestern New Mexico, USA and borders northern Mexico ( Figure 1B). With an area of approximately 7,200 km 2 , CNF comprises 12 geographically distinct management areas across 16 isolated mountain ranges, often called the "Sky Islands, " that form the northernmost extent of the Madrean Archipelago. Located in the Sonoran Desert, these ranges produce an orographic effect, creating cooler, wetter climatic conditions at higher elevations that reach upwards of 3,267 meters above the surrounding desert floor (Figure 1). These high elevation Sky Islands are biodiversity hotspots supporting over 7,000 floral and faunal species with many plant and animal communities usually only seen in more northern latitudes (Yanahan and Moore, 2019). The American Southwest climate, in which the CNF is located, is predicted to increase in temperature by [5][6][7][8], decrease in precipitation by 10%, and increase in the number and duration of periods of high temperature (Service U.F, 2013).The CNF, like the SMO, presents unique and challenging opportunities for land managers and climate change Frontiers in Ecology and Evolution 03 frontiersin.org research due to abrupt variations in topographic variables (i.e., elevation, slope, aspect) and the isolation between mountain ranges.

Tree species
The highly variable climate across the elevational gradient of the SMO and CNF (Figure 2A) allows for the development of distinct forest communities. These communities range from Broadleaf Evergreen Oak woodlands at lower elevations to Spruce-Fir forests at high elevations. Forest communities are distinct assemblages of trees at different successional stages or compositional maturity (Davis, 1981). While somewhat variable in occurrence across the landscape, the SMO and CNF houses major forest communities of interest (Service, 2013) Service, 2022). Due to limited data availability, the two highest elevation forest communities, Engelmann Spruce-Alpine Fir Forests and Bristlecone Pine-Limber Pine Forests located in the Pinaleño range, were excluded from this study.
Ecological community composition may range from few to many species, resulting in a complex system of species competing or relying on each other for resources. Forest communities, represented by dominant, co-dominant, and understory (sub-canopy) plants, rely on the dominant species, whose crowns commonly rise above the general canopy. Dominant trees provide various functions that guide overall forest community function and structure (Power et al., 1996). Observations suggest that dominant tree species direct the ecological processes that regulate microbial soil community assembly, which play a critical role in maintaining forest functional and structural stability (Ma et al., 2019). Dominant tree species are also driving factors in shaping spatial patterns of plant biodiversity, locally and regionally (Hao et al., 2007;Draper et al., 2019). Examining the spatial distribution of habitat suitability of dominant tree species under climate change is a relatively new approach to understanding ecosystem vulnerability, environmental change, and spatial conservation prioritization (Périé and de Blois, 2016;Wan et al., 2017;Dyderski et al., 2018). This study, while analogous to the overall topic, uniquely examines the seven dominant tree species of the SMO and CNF (Table 1).

Environmental predictor variables
Ecological niche model requires the use of environmental predictor variables, and variable selection should reflect known The spatial extent of (A) the Sierra Madre Occidental (SMO) symbolized in dark green and, (B) the Coronado National Forest boundary in light blue overlaid.
Frontiers in Ecology and Evolution 04 frontiersin.org drivers of species distribution and habitat selection. The drivers of plant distribution are primarily climatic, edaphic (i.e., soils), and topographic. Thus, we gathered a collection of environmental variables that are commonly associated with plant distribution and habitat for examination. It is noted that many variables are highly correlated with others, which is known to obscure the contributions of individual variables in predictions (Zurell et al., 2020). In this study, the primary interest is in the final predictions, rather than any of the model coefficients and multicollinearity is known to not affect prediction quality, provided that predictions are made within the observed region of observations (Kutner et al., 2004). While multicollinearity does not affect predictions, there is still value in simplified model representations that achieve similar accuracy to models using many more variables. As such, we created a subset of variables by developing a matrix of all variable values extracted to each species occurrence. Then, pairs of variables in the matrix were examined, and if shared a correlation of more than 0.8 was present, we removed the second variable from consideration thereafter. Selecting variables that contribute to the understanding of species distribution and habitat suitability while reporting values beneath a threshold of 0.8 is a suggested rule of thumb in handling multicollinearity in statistical modeling (Grewal et al., 2004;Kim, 2019). This led to the collection of 22 variables used in subsequent ENM (Table 2).
In order to examine how species habitat suitability may change under future climate scenarios, we employed the use of bioclimatic variables commonly used in climate modeling research (Priti et al., 2016;Booth, 2018), produced from two distinct future anthropogenic emissions scenarios from the Coupled Model Intercomparison Project 6 (CMIP6) (Karger et al., 2018;Hamed et al., 2022), primary data modeling developed by the Geophysical Fluid Dynamics Laboratory (GFDL), downscaled to a 1-km 2 spatial resolution by the Climatologies at high resolution for the earth's land surface areas (CHELSA) research (Karger et al., 2017). The two scenarios, known as Shared Socioeconomic Pathway (SSP) 3-7.0 and SSP 5-8.5, are sectioned into three distinct 30-year normal: 2040 (2011-2040), 2070 (2041-2070), and 2,100 (2071-2,100). Each SSP, derived from anthropogenic emission scenarios, ranges from moderate to extreme climatic change through the end of the century (Alizadeh et al., 2022). SSP 3-7.0 infers that future global climate results in the medium to the high end of the range of future forcing pathways, leading to exceptionally high aerosol emissions and land use change (Meinshausen et al., 2020). , while similar in overall trajectory, represents the high end of the range of future pathways with a reference scenario in a high fossil fuel development world throughout the 21st century (Meinshausen et al., 2020).
While more than two scenarios exist, data availability and the feasibility of SSP 3-7.0 and SSP 5-8.5 becoming a reality drove selection over other available SSP used in this study. Thus, we developed seven ENMs that use the 19 bioclimatic variables associated with a specific time frame (Current, 2040(Current, , 2070 and SSP (3-7.0 and 5-8.5). Edaphic variables were selected from the

Ecological niche modeling
Ecological Niche Modeling (ENM) is the process of relating geographically referenced species occurrence records (i.e., presence/ absence records), environmental variables (i.e., elevation, soil moisture, annual mean temperature, etc.), and using statistical algorithms to predict habitat suitability for a given species across a region of interest (Pearson, 2007;Elith and Franklin, 2013). ENM is often considered synonymous with Species Distribution Modeling (SDM) and Habitat Suitability Modeling (HSM), though (Brown and Griscom, 2022) highlight some important differences between the terms in certain contexts. ENM is commonly applied to studies researching the distributions of plant species (Carrell et al., 2022), and consequences of future climatic change on their geographic distribution (Nakao et al., 2013;Moor et al., 2015). Our study uses ENM to investigate the geographic alterations in suitable habitat of tree species of the SMO and CNF under various anthropogenic climate change emissions scenarios (Figure 3). In ENM, no one model is constituted as "best. " Therefore, we implemented a mean ensemble approach, which simply averages probability outputs of predictive models, from three machine learning algorithms commonly applied in ENM. Ensemble models are examined to bolster model defensibility and account for possible alternate states of modeled system (Araújo and New, 2007;Forester et al., 2013). These models include: (1) Random Forest (Breiman, 2001); (2) Maximum Entropy (Elith et al., 2011); (3) Boosted Regression Trees (De'ath, 2007). The three models were chosen due to their observed robust performances in ENM studies (Fourcade et al., 2014;Mi et al., 2017;Yu et al., 2020). Random Forest, Maximum Entropy, and Boosted Regression Tree models were developed using default parameters found in the randomForest (RColorBrewer and Liaw, 2018), dismo (Hijmans et al., 2017), and gbm (Ridgeway and Ridgeway, 2004) packages in R, respectively. For model validation, Random Forest does not require a validation dataset by using an out-of-bag evaluation (Cutler et al., 2012), while Boosted Regression Trees and Maximum Entropy models were validated using a k-fold cross validation (Hastie et al., 2009).
Species presence data was gathered from commonly used citizenscience databases [e.g., Inaturalist (Matheson, 2014), Global Biodiversity Information Facility (Robertson et al., 2014), 'USGSs' Biodiversity Information Serving Our Nation (Guala, 2019), etc.] obtained through the Interface to Species Data Sources (spocc) R package (Chamberlain et al., 2021). The number of occurrence records used in both spatial scales are listed by species in Table 3. The use of species occurrence records collected through citizen-science methodology has been applied to a number of ecological modeling studies and is observed to improve ENM accuracy (Amano et al., 2016;Milanesi et al., 2020;Feldman et al., 2021;Matutini et al., 2021). Species pseudo-absences are generated at a 2-1 ratio for every observed presence in the dataset through random sampling of points across the modeling extent of the SMO and CNF, consistent with other pseudoabsence generated studies (Barbet-Massin et al., 2012;Liu et al., 2019).
Given that three different models were used and each have different conditions and modes of prediction, a maximum kappa classification threshold was implemented for consistency across statistical algorithm selection (Duan et al., 2014). We used the PresenceAbsence R package (Freeman and Moisen, 2008), which calculates the thresholds for any dataset using presence/absence response variables. By selecting "MaxKappa" as the threshold metric, model thresholds are optimized to the find the threshold that gives the maximum value of Kappa, which uses information from a confusion matrix to assess the improvement over chance prediction (Freeman and Moisen, 2008). Model accuracy was measured using the Area Under the Curve (AUC), Sensitivity, Specificity, and the True Skill Statistic (TSS). These metrics measure the ability of the models to correctly classify the presence and pseudoabsences in the training data. Specifically, AUC is an important index as it provides a measure of overall accuracy that is not dependent upon a particular threshold (DeLeo, 1993;Fielding and Bell, 1997). Sensitivity represents the proportion of correctly predicted presence records and quantifies omission errors, while specificity represents predicted absence records and quantifies commission errors (Shabani et al., 2018). TSS is a performance measurement and shows of SDMs in which predictions are expressed as presence-absence maps and shows significant correlation with those of the threshold-independent AUC statistic (Shabani et al., 2018). Each mean ensemble model predicts the [0.1] interval that represents the estimated probability of species occurrence (i.e., habitat suitability). Predictions are produced in a spatially explicit raster at a 1-km 2 spatial resolution.

Generalized additive models
As previously mentioned, habitat suitability was estimated using an ensemble of three machine learning algorithms at a 1 km 2 resolution at a local and regional scale. Post-ENM modeling, we examined both models in their current extent, and clipped both models to the geographic boundaries of the CNF ( Figure 1A) This process allows for the examination of local vs. regional scale models at 2 different scales and the examination of local vs. regional models when observed at the same scale. A shapefile of the plant communities of interest ( Figure 2B; USDA Forest Service, 2022), was used to partition the model grid cells into the distinct geographic boundaries of expertly defined plant communities. Grid membership in the partition was based on the geographic coordinates of the grid centroid. A GAM was used to relate the gridded estimates of habitat suitability to a gridded map of elevation and latitude obtained from a 1km 2 digital elevation model of North America (Portal, 2007) and a shapefile of latitude from ESRI (2011). The gridded elevation and latitude maps were resampled using bilinear interpolation to produce co-located elevation values for the grid of suitability predictions. In this study, the GAM includes only a single predictor: where y  represents the estimate of habitat suitability, x represents elevation, and α represents a global model intercept (Hastie et al., 2009). The function s x represents a smoothing spline, which is a series of polynomial curves connected at breakpoints (often called "knots") in the x dimension. The piecewise polynomials are required to be connected in such a way that the overall curve is both continuous and differentiable. The number of knots in the spline determines how much local variability is characterized: too many knots and the GAM will overfit the data, too few knots and the GAM will smooth over the true relationship between x and y . In the R implementation smoothing splines, the model smoothness is automatically selected by means of a penalized cross validation (Wood, 2004;Wood, 2017). It is the authors' experience using the default smoothing settings in many different contexts that the automatic smoothness selection rarely overfits the input data, especially when the sample size of inputs is large n ! 100 . The assumptions made by the GAM model are that the data vary smoothly (though no model form needs to be specified), and that the variance of the data about the true model form is constant. The satisfaction of these assumptions allows for the estimate of prediction intervals as in ordinary least squares (OLS) regression. The key difference between OLS and GAM modeling is that the GAM does not require pre-specifying the model form that describes the relationship between x and y . This allows great flexibility in determining smoothed representation of saturated scatterplots without having to guess a model form a priori. The popularity of GAMs for visualizing trends is evident in the default smoothing options of R's most popular data visualization package ggplot2 (Wickham et al., 2016).

Ecological niche models
Common modeling metrics (AUC, Sensitivity, Specificity, TSS) were generated producing the following values for all models developed at both the CNF and SMO scales and are reported as (Median, Mean) (Table 4). Statistically, models performed consistent with the standards of modeling accuracies and analogous research (Marmion et al., 2009;Williams et al., 2009). With AUC mean values above 0.90, it is inferred that the models have high predictive capability according to the predictor and response variables provided. Metrics for each specific model are found in the Supplementary Table S1.

End of century observations
Habitat suitability is observed to alter in both spatial distribution and quality (i.e., overall modeled probability) for all seven tree species of interest, regardless of the emissions scenario (Figure 4). For pine species, habitat suitability is observed to be negatively impacted by anthropogenic climate change in both quality and spatial distribution across the entire SMO extent. Pinus engelmannii (PIEN), is observed to be most negatively affected among the pine species. PIEN models produced a maximum probability prediction of 98% under current conditions with high values above 81% in both southern and northern extents of the SMO (Figure 4.2A). Under the predicted emission scenarios of SPP 3-7.0 and SSP 5-8.5 for the year 2,100, PIEN is predicted to decrease significantly, ranging in maximum probability of 61% (SPP 3-7.0) and 56% (SSP 5-8.5) (Figures 4.2B,C). The other 3 pine species exhibit characteristics that are complex to describe spatially. Pinus leiophylla (PILE) is observed to decrease in habitat suitability across the SMO over time but retain more pixels than others that could be interpreted as moderate suitability at 41-60% (Figure 4.3). Pinus arizonica (PIAR) (Figure 4.1) and Pinus strobiformis (PIST) (Figure 4.4) are observed to have more pixels that could be interpreted as low suitability (0-40%). However, PIAR and PIST contain pixel values that reach a maximum probability of 79 and 81%, respectively, ranging from 9 to 11% greater than the highest pixel predicted for PILE (70%). Oak species are predicted to fare much better than pine species, in terms of both habitat quality and spatial distribution within the extent of the SMO. The highest predicted probability values of oak species were those reported by Quercus arizonica (QUAR) models (Figure 4.5) for current and end-of-century conditions for both scenarios, which ranged from 96 to 93 to 94%, respectively. Quercus emoryi (QUEM) (Figure 4.6) and Quercus grisea (QUGR) (Figure 4.7), exhibit decreases in maximum probability prediction in both future scenarios to a greater, but not extreme extent than QUAR, with values ranging from 95-89% (QUEM) and 94-88% (QUGR). Increased suitability is observed for all oak species in more southernly extents of the SMO in future scenarios. Among these, QUEM may show the greatest increase in habitat suitability across the SMO in terms of habitat suitability, spatial exhibiting pixels to increase from a reported 0-40% suitability (Figure 4.6A) to greater than 80% (Figures 4.6B,C).
When compared visually, emissions scenarios show little variation on the predicted habitat suitability of species at the end of the century. Models reported contrasting, yet to be considered insignificant, maximum probability predictions for each species in SSP 3-7.0 and SSP 5-8.5. For example, PIAR shows a maximum probability of 79% at the end of the century for both future scenarios. The other 3 pines show a slight decrease of 61 to 56% (PIEN), 71 to 70% (PILE), and 82 to 81% (PIST) for SSP 3-7.0 and SSP 5-8.5, respectively. With exception of QUAR, which reports a higher maximum probability in SSP 5-8.5 (94%) than SSP 3-7.0 (93%), oaks follow suit with either constant or decreasing values among emissions scenarios at a fixed 89% (QUEM) or 89 to 88% (QUGR).

Spatial scale comparisons
A primary focus of this study was to not only observe species responses to anthropogenic climate change in terms of habitat suitability and spatial distribution, but to do so at both local and regional spatial scales. Model outputs were cropped to the boundaries of the CNF (i.e., local scale) and end-of-century predictions under SSP 5-8.5 (see Figure 5) were compared. The magnitude of differences in predicted habitat suitability for each species across the two model scales varied greatly depending on species. Some species, such as PIEN, showed similar predictions for suitability in the CNF using the two modeling Methodology flow-chart. scales. Other species, like QUEM, showed differences in model predictions that were larger than the projected changes in suitability due to climate change. The general trend observed among pine species is a decrease in habitat suitability when comparing SMO and CNF modeling extents. Variation from this trend is observed in PIAR and PIST models as some areas are shown to increase in habitat suitability in the SMO modeling extent. Another observation among pine species is that the SMO modeling extent may result in a spatially focused prediction, with less values predicted at 41-60% and more at 61-80%. Consistent with previous observations in section 3.2 (End of Century Observations), PIEN shows very little habitat suitability while PILE shows consistent values for both modeling extents at the end of the century for SSP 5-8.5.
Oak species show more consistency as compared to pine species among observed predictions in both modeling extents when compared in section 3.2 (End of Century Observations). Both modeling extents show higher estimates of habitat suitability when compared to most pine species, however, there are arguably more significant variations in predicted habitat suitability when comparing the two modeling extents. All three oak species visually show increased values predicted in habitat suitability in regional modeling extents than local modeling extents. QUEM, when cropped to the SMO modeling extent shows no pixels reporting values below 41%, with a large portion reporting greater than 61%. QUAR shows increases in the northern ranges of the CNF with values ranging from an observed 0-40% in CNF modeling extents to greater than 41% with most pixels reporting 81% and above in SMO modeling extents. Regardless of the specific pixel value being reported, visually, there is significant difference in the comparison between model outputs produced from the CNF modeling extent and SMO modeling extent.

Forest community suitability
All species showed variability in habitat suitability across community type and climate change scenario ( Figure 6). The most notable observation is the decreasing values in habitat suitability across time for all species, except for QUEM, which reports slight increase Habitat suitability for the seven tree species of interest across the Sierra Madre Occidental for current climatic conditions, end-of-century (2100) climatic conditions for SSP 3-7.0, and end-of-century (2100) climatic conditions for SSP 5-8.5. Habitat suitability probability is represented on a gradient of blue (low suitability) to red (high suitability). Predicted model probability range is in the upper right corner of each species map as [lower-upper%]. Figure column is assigned numerically while column is assigned alphabetically for rapid lookup (i.e., Pinus Strobiformis at end-of-century conditions for SSP 3-7.0 = 4.4.B).
Frontiers in Ecology and Evolution 09 frontiersin.org over time. Perhaps equally interesting is the lack of variation among modeling extents for the pine species in estimated suitability across time. When examining modeling extent, Oak species show quite significant variation in reported habitat suitability for all communities across time, with the SMO modeling extent estimating much higher suitability than the CNF. Additionally, pine species show consistent higher prediction values in habitat suitability for forest communities that are higher in elevation. Most interestingly, only one oak species (QUGR) follows the general trend observed among the pine species. This may suggest that higher elevations are becoming increasingly important for pine species survival, and QUAR and QUEM may be potential invaders to higher elevation forest communities.
It is important to note that these models only consider the suitability of individual species to their environment and do not account for competition among species. This means that even small changes in the bioclimatic suitability for one species could translate into a large and non-linear shift in the relative abundance of all species. Berger et al. (2008) highlights various approaches for modeling competition among plants. As a specific example, Seidl et al. (2014) demonstrates that live tree legacies play an important role in the ability of forests to recover from disturbances. Thus, the future of tree growth is partially a function of habitat suitability, and another part habitat legacy. However, much work needs to be done to explore the interaction of competition models in the face of a changing climate, which lies beyond the scope of this current effort.

Suitability and elevation
The two most natural characterizations of habitat suitability are over space and time. However, elevation drives climate variability over small geographical distances like the CNF (Randin et al., 2009). The GAMs produced visually compare the estimated changes in mean habitat suitability across elevation for both climate change scenarios and spatial modeling extents (Figure 7). GAM estimates across elevation allow us to visualize species responses to anthropogenic FIGURE 5 Habitat suitability models of all seven tree species of interest under the SPP 5-8.5 emissions scenario for 2,100. Each species model is cropped to the boundary of the Coronado National Forest under two distinct modeling scales: (A) Coronado National Forest and, (B) Sierra Madre Occidental.
Frontiers in Ecology and Evolution 10 frontiersin.org climate change that are difficult to detect accurately in mapping products like Figures 4, 5. All pine species are observed to show decreased habitat suitability for elevations below approximately 2,700 m, when compared to current conditions. Above approximately 2,700 m, all pines modeled to the extent of the SMO increase in habitat suitability, while the same species modeled at the CNF extent do not. Pine species also show a much higher level of agreement among the estimates for different emissions scenarios within modeling extents than the same emission scenarios across modeling extents. Generalized additive models estimates of oak species suitability across elevation varies drastically depending on the species. Most notably, models developed at the SMO modeling extent show higher suitability across all elevation than those developed at the CNF modeling extent. The commonly highest values in habitat suitability across elevation for pines were reported by current conditions while oaks report varying results, with each species exhibiting a different emission scenario as the highest. Oaks also report a dip in habitat suitability, at approximately 1,500 m, which is most likely the result of combinations in available land for habitat at that elevation. Figure 8 shows smoothed estimates of mean habitat suitability for all seven species across latitude within the SMO. Latitude is observed to be a driving factor in species response to a changing climate (Colwell et al., 2008) and represents and interesting variable to consider among multiple-scale modeling comparisons. Arguably, latitude would play little role in species response to climate change within the spatial extents of the CNF. The SMO, however, is a reasonably large ecoregion spanning more than 10 degrees in latitude. Habitat suitability of pine species was similarly reported to respond to climate change across latitude, with a small over decrease in suitability for PIEN. Oak species show similar responses across latitudes with general increases in habitat suitability under future emissions scenarios. Future emission scenarios for QUAR and QUGR eventually taper off at cross with current conditions at high latitudes, while QUEM exhibits higher values in habitat suitability across all latitudes. Interestingly all species respond to latitude greater than 32 degrees with an abrupt increase. This is likely due to the patchiness and sharp increases in elevation across space for mountain ranges occupying of the Madrean Archipelago and CNF.

Discussion
Our results indicate substantial alterations in habitat suitability among all tree species across all emissions scenarios and spatial scales. Pine species habitat suitability is generally observed to decrease in both spatial distribution and habitat quality. Changes in habitat Projected change in habitat suitability of the seven tree species of interest within four common forest community boundary designations in the Coronado National Forest. Lines represent the mean estimated suitability.
Frontiers in Ecology and Evolution 11 frontiersin.org suitability for oak species are more localized and species specific. These results align with anticipated shifts in species distribution from recent, analogous research (Du et al., 2021;Erfanian et al., 2021;Rödder et al., 2021). Our study also suggests that when modeling species distribution and habitat suitability under predicted anthropogenic emissions scenarios at local and regional scales, climate forcing scenario choice appears to be far less important for understanding drivers in species habitat suitability than the spatial scale of the models. This claim suggests that conservation efforts based on ENMs must consider and address their spatial scale of choice, and where possible, consider the sensitivity of model results to varying spatial extents. Ignoring options in the spatial extent at which models are used to train data overlooks the fact that all model results are conditional on a somewhat arbitrary choice of modeling extent. One shortcoming of models developed at smaller spatial scales is niche truncation (Chevalier et al., 2022). Local scale models incorporate a limited range of climate, omitting outside climatic influences that may ultimately affect the whole ecological niche of species (Peterson et al., 2018). Regional models overcome this shortcoming by providing a greater range of climate model inputs, thus avoiding inappropriate model extrapolation. That in mind, the SMO spans multiple countries with vastly different opinions of forest management. It thus becomes difficult to distinguish climate effects from geopolitical effects at the regional scale. This reinforces the need for researchers to carefully consider the consequences of their spatial modeling extent choice.
Ecological niche models techniques presented in this study ignore the reality of biotic relationships and that species are in competition with each other for limited resources (Fern et al., 2019). One goal addressed was to examine if alterations among species occupying a particular community (i.e., Broadleaf woodland, Coniferous Woodlands, etc.) would invade the boundaries of another community, Comparison of the smoothed relationship between mean habitat suitability and elevation for the seven tree species of interest in the Coronado National Forest for current conditions (black, solid) and 2,100 estimates (colored, dashed) under different emissions scenarios.
Frontiers in Ecology and Evolution 12 frontiersin.org which was examined in Figure 5. However, each species occupying a community generally neighbors the spatial boundaries of other plants. This means that even small changes in habitat suitability could lead to drastic and non-linear changes in the competitive balance between species, resulting in large alterations in community composition and distribution (Alexander et al., 2015(Alexander et al., , 2016. Future studies would benefit from more direct investigations of the change in the competition among species, rather than looking solely at the habitat suitability of species individually. Another consideration is the matter of spatial resolution. Smallscale variability is difficult to derive from course-grained, global climate data. Additionally, substantial changes in species habitat suitability projected in larger-scale studies may overlook small-scale areas like the CNF, which may present resilience to certain aspects of climate change. Local-specific features such as streams and springs could influence micro-climate forest communities to be more or less resilient (Wieser et al., 2003;Anderson et al., 2007). However, the smallest spatial resolution of climate data freely available at the time of this study was 1 km 2 , which would not capture fine-resolution microclimate. The heterogeneity that occurs with scale highlights the need for spatial models to match the scale of specific management boundaries (i.e., local, regional, global). This would allow models to be more applicable to the immediate climate adaptation science and land stewardship efforts.
Additionally, differences between modeled alterations and realized alterations may be the result of tipping points generated by extreme events (Saccone et al., 2009;Niu et al., 2014). Alters in typical climate conditions over a 30-year period may not be substantial enough to considerably change habitat suitability while increases in the yearly variance may supply a sufficient threat in forest resilience (Wang et al., 2016). For example, one exceptional drought could provide the catalyst for widespread forest fires or insect outbreaks Comparison of the smoothed relationship between mean habitat suitability and latitude for the seven tree species of interest in the Coronado National Forest for current conditions (black, solid) and 2,100 estimates (colored, dashed) under different emissions scenarios.
Frontiers in Ecology and Evolution 13 frontiersin.org that would noticeably alter the ecology of the CNF. Wildfires (expected to be more frequent and widespread in the future) could have large consequences due to the different fire resilience and postfire rejuvenation of each species (Barton and Poulos, 2018). Additionally, our results imply significant changes in species habitat suitability across elevational, community, and latitudinal gradients, but realized changes in species distribution and habitat suitability may be a result of temporal lags in individual response to climate change. Our study reinforces the point that modeling efforts may allow us to anticipate the future, but ultimately should be used in conjunction with observations in current and stochastic environmental occurrences.
Modeling the future is an exciting topic in the field of landscape ecology and management, however the true consequences of anthropogenic climate change will be realized in observation over time (Gonzalez et al., 2010;Tewari et al., 2017). The results of this study call for the long-term monitoring of forest species. Our models, reliant on geo-referenced species presence largely collected by citizen science, allow land managers to develop a more-focused resampling and ground-truthing program, to continually update and add to the reliability of habitat models. We also suggest a distinct monitoring of seedling and seedling habitat suitability, which present a significant role in forest change over time as tree regeneration is the driving factor in long-term survival (Ribbens et al., 1994;Warren et al., 2012). Additionally, we suggest that when considering the use of ENM for applied conservation efforts, to do so at multiple scales (Bradter et al., 2013;Hallman and Robinson, 2020). Doing so will allow for (1) the capture of a greater range of influencing environmental drivers on species habitat suitability and, (2) a greater confidence in model result and interpretation. Opportunities abound in the application of ENM to multiple spatial scales of interest as local, complex study areas like the CNF may misrepresent changes only observable at the regional scale, which in an era of global environmental change, may be costly to management and potentially fatal to biodiversity.

Data availability statement
The original contributions presented in the study are included in the article/Supplementary material, further inquiries can be directed to the corresponding author.