1. Introduction

Bioclimatic variables (biologically meaningful variables created from monthly temperature and precipitation) can have direct effects on birds leading to limits on their distributions ^[1] and changes in species habitat value ^[2] when these bioclimatic variables change (e.g., ^[3,4]). Several studies have already indicated how recent changes in climate have affected the geographical ranges of birds ^[5], their reproduction and migration ^[6]. This is one reason why numerous studies have based their predictions of the effects of climate change on species distributions on bioclimatic-envelope models (e.g., ^[7,8]).

The use of the bioclimatic variables alone to model the species' climate space has been met with objections by a handful of studies due to the absence of a broad range of climate change-related stresses in the model that could affect population ecology and physiology ^[9,10,11]. Studies like that of Jeschke and Strayer ^[12] and Dawson et al. ^[13] asserted that the consequences of climate change are a multifaceted problem that is not completely encompassed by assessing exposure of a focal species to climate change using bioclimatic-envelope models. A complete assessment of the impacts of climate change in a focal species includes other variables, such as land use ^[14], invasive species, and pollutants ^[12]. The interrelationships among these variables can become complex ^[15].

Despite the controversy over the limitations of bioclimatic-envelope models, many of their proponents have high praises for their predictive power ^{[16,17,18,19]}. These models have provided broad insights regarding the likely effects of climate change on species distribution ^[20] and biodiversity, especially the impacts on vulnerable species ^[11]. Beaumont et al. ^[2] and Elith et al. ^[21] have demonstrated that, even with just climate data alone, models were effective in establishing the current distributions of species, resulting in baseline models that could be used to predict the effect of climate on future distributions of conditions suitable for species. Further, in a modeling effort by Huntley et al. ^[22], bioclimatic variables were seen as the main range-limiting factor among environmental variables.

While other predictor variables such land cover is often used for modeling the present species distributions, it is not often utilized when projecting to future distributions ^[23]. In fact, other predictors, including land cover, were found to be not critical in species distribution modeling compared to climate variables ^[24]. Here, we evaluated bioclimatic-envelope models ^[16,25] in projecting availability of suitable bioclimatic conditions for nine bird species, identified as species of concern (SOC) in the South Central United States (U.S.), using various climate projections derived from general circulation models (GCMs) run according to Representative Concentration Pathways (RCPs) and post-processed via application of a simple statistical downscaling method. We compared future projected climate envelope suitability results produced by combinations of four GCMs and two greenhouse gas RCPs for two future time periods. Our objectives are (1) to develop models of present day and potential future distributions of suitable environmental conditions for multiple bird species of conservation concern in the South Central U.S. region, and (2) to compare how bioclimatic-envelope suitability is projected to change from present day to future conditions. Our hypothesis is that much of the area within the range of the focal species that currently possesses suitable bioclimatic conditions would be converted to unsuitable conditions in the future and that all birds considered here would be at high risk from climate change. Specifically, we expect to see bigger changes in some species due to their biology and vulnerability to climate variables.

---

2. Materials and Method

We employed a three-step procedure for this study. First, we gathered and processed datasets regarding the bird species occurrence, their species ranges, and the bioclimatic variables used in the bioclimatic-envelope modeling. Second, we developed models for current conditions, which included selecting SDMs and generating current bioclimatic conditions, and lastly, we modeled future conditions, including selecting species distribution models (GCMs) and RCPs, and projecting current to future conditions. The general flow of the study is shown in Figure 1.

---

2.1. Species data

After a rigorous selection process including feedback from valued stakeholders within the focal area for the South Central Climate Science Center, we identified nine birds as focal species. The nine bird species of interest-Northern/Masked Bobwhite Quail (Colinusvirginianus), Scaled Quail (Callipepla squamata), Pinyon Jay (Gymnorhinus cyanocephalus), Juniper Titmouse (Baeolophus ridgwayi), Mexican Spotted Owl (Strix occidentalis lucida), Cassin's Sparrow (Peucaea cassinii), Lesser Prairie-Chicken (Tympanuchus pallidicinctus), Montezuma Quail (Cyrtonyx montezumae), and White-tailed Ptarmigan (Lagopus leucurus)-are found within the South Central U.S. (Figure 2). We selected birds from this region because this region is known to be the driest region in the U.S. and changes of the climate could adversely affect diversity of animal species. For more information about the species selection process, see Appendix 1.

Original species data included 27, 232 bird presence records, with 12, 260 for Northern/Masked Bobwhite Quail; 3629 for Scaled Quail; 3237 for Pinyon Jay; 2900 for Juniper Titmouse; 1763 for Mexican Spotted Owl; 1570 for Cassin's Sparrow; 993 for Lesser Prairie-Chicken; 485 for Montezuma Quail; and 395 for White-tailed Ptarmigan. We obtained species occurrence datasets from Natural Heritage programs in New Mexico, Arizona, and Texas. We only used presence data since there was not sufficient absence data available. Other online sources of presence data used in this study include: Biodiversity Information Serving Our Nation (BISON) ^[26], National Science Foundation's (NSF) bird specimen collection (ORNIS) ^[27], and NSF's biodiversity data portal (VertNet) ^[28]. For datasets that were provided in polygon format, we converted to point data, with a centroid (or point close to a centroid but contained within the polygon) generated for each polygon. We standardized dataset attributes across sources to match data from BISON and fixed errors in coordinates, as each data source has a different set of attributes for species record. Changing to the same projection was included in the process. In this study, we limited the datasets combined across sources to the same date range for which historical climate data are available (1950 to 2000). Also, we eliminated occurrences that fell outside the boundaries of species ranges obtained from the U.S. Geological Survey (USGS) National Gap Analysis Program (GAP) ^[29] to reduce the number of records for which the species had been misidentified or were otherwise geographic outliers. All spatial data analyses were conducted using ArcGIS for Desktop 10.3 (ESRI) ^[30].

---

2.2. Species distribution modeling

We obtained a set of 19 raster-based bioclimatic variables (Table 1) from among the WorldClim datasets ^[31] to describe present environmental conditions and explore the relationship between bioclimatic conditions and species distribution patterns. For each species, we removed one of each pair of highly correlated (r > 0.7) ^[32] environmental variables from the bioclimatic-envelope models to avoid collinearity among variables ^[33]. We chose between highly correlated variables based on the results of a species-specific literature search. In particular, we selected variables that were identified in one or more studies regarding each individual species of interest as having an effect on the species' range or population dynamics (Table 1). In cases where the results of the literature search could not differentiate between two highly correlated climatic variables, we used a qualitative assessment of the distribution of values of the variable at all presence points and of the relationship between the variable and species presence or pseudo-absence ^[34].

WorldClim provides climate projections statistically downscaled, using a "delta method" approach, to a spatial resolution of 30 arc-seconds, roughly 900 m at the equator. We only used the bioclimatic variables in the models because we assumed that change in climate could be a driving force of modifications to species habitats, especially for species of conservation concern and vulnerable species. If the climate continues to change, associated changes to species habitat could be damaging ^[11], and could lead to reduction in species population sizes as suitable habitats disappear ^[35].

Araujo and New ^[36] emphasized that using a single SDM algorithm would not provide sufficient information regarding whether the selected technique gave the best predictive accuracy for the dataset. Here, we analyzed the bioclimatic condition by species distribution relationship using the following SDMs for each species: Generalized Linear Model (GLM), Random Forest (RF) ^[37,38], Boosted Regression Tree (BRT) ^[39], Maxent ^[40,41], and Multivariate Adaptive Regression Splines (MARS) ^[42]. We selected SDMs based on their performance with presence-only data ^[43]. We used the modeling tool Software for Assisted Habitat Modeling (SAHM) run within VisTrails ^[34,44] to create a workflow and develop bioclimatic-envelope models for present day conditions. When multiple species occurrences were present within a given pixel of the climatic data, a tool in SAHM consolidated them to a single occurrence per pixel. Since species lacked absence data, the tool randomly generated background points (i.e., pseudo-absences ^[41]) within a 95% minimum convex polygon defined by the presence data. Phillips et al. ^[45] suggested an improved method for background sampling that uses target group, however, the approach may not be applicable for our multisource datasets. Also, Phillips et al. ^[45] warned the use of target-group background for extrapolation studies that involves future climate conditions. We extended the species ranges obtained from the USGS National Gap Analysis Program ^[46] to generate a template boundary layer for each species. The length of the extension was based on the dispersal ability of the nine bird species ^[47,48]. In this study, we extended a maximum distance of 25 km beyond what USGS defined for the present species range. We used this template to restrict model development and future projection of the species. The choice of the extent of the species range was critical in the performance of the SDMs. Since SAHM assigns background points for each SDM, selecting a range that is too broad or too restricted could negatively affect the relationship between background and presence points ^[49]. To avoid the mismatch between species extent and presence/absence points, we used the unique geographic range of the species to define where background points should be assigned.

Further, we utilized the SAHM ensemble tool to produce ensemble maps for the present day distribution of suitable environmental conditions for each species. The ensemble map is a summation of binary maps generated from probability surfaces from each statistical modeling algorithm ^[50,51,52]. We used specificity = sensitivity as the threshold in discretizing the probability maps. This has previously been identified as the optimal threshold ^[53]. The final ensembles consisted of pixel values that showed the number of models in agreement that a particular pixel is suitable for the species. A pixel with a value of zero meant that none of the models identified bioclimatic suitability for the species at that location, while a value of five meant there was agreement across all five models.

We assessed confidence in individual model results in terms of concordance among the different distribution models. We had higher confidence that environmental conditions were suitable for a species when three or more (at least 60% of) algorithms were in agreement (e.g., ^[54]). We compiled information on various measures of model performance, including the Area Under the Receiver Operating Characteristic (ROC) Curve (AUC) for the test data, correct classification rate (%C) ^[55,56], the True Skill Statistic (TSS) ^[57], and the kappa statistic ^[58] for each algorithm by species combination. Swets ^[59] classified values of AUC as follows: those > 0.9 indicated high accuracy, from 0.7 to 0.9 indicated good accuracy, and those < 0.7 indicated low accuracy. Kappa values ranging from 0.2 to 0.5 were labeled poor, 0.6 to 0.8 useful, and those larger than 0.8 were excellent ^[60]. We checked other qualitative assessments of model performance, which included inspecting cross-calibration and deviance of residual plots. Calibration plots indicate whether models tend to over or under predict habitat suitability. Deviance of residual plots are used to identify individual data points that may require further inspection or whether there may be an important environmental layer missing from the model inputs ^[44]. The current version of SAHM ^[39] does not produce evaluation metrics for the ensembles, so we ran a code in R, with the same settings as the ones we used in SAHM, to generate ensemble statistics. We then compared the statistics of the five individual SDMs to their ensemble for all nine species. Finally, we used the current distributions estimated by the SDMs and projected each to climate conditions in the years 2050 and 2070.

---

2.3. Projecting models to future climatic conditions

Informed in part by previously published evaluations of model performance ^[61,62], we screened GCMs based upon their simulations of 20th century climate across the continental U.S. and regions that overlap the study area, as well as the general areas inhabited by the focal species (e.g., Central and Western North America). We used values showing bias of model output relative to observed historical data as one of several criterion to exclude GCMs. In particular, we excluded GCMs for which multiple variables had a relatively high bias (i.e., were more biased than two times the standard deviation of variation among biases of all models evaluated) or for which few evaluated variables were less biased (i.e., bias was less than half of the standard deviation of variation in bias among all models evaluated). Also, we excluded models with large values ( > 1 or < -1) for top of atmosphere energy imbalance (Wm-2) since these values may be an indication of long term drift in simulated climatic conditions (Forster et al. 2013). The final list of selected GCMs includes: Community Climate System Model version 4 (CCSM4) ^[63], Hadley Centre Global Environment Model version 2-Earth System (HadGEM2-ES) ^[64], Model for Interdisciplinary Research on Climate version 5 (MIROC5) ^[65], and Max Planck Institute Earth System Model, low resolution (MPI-ESM-LR) ^[66].

For future conditions, we used the downscaled data provided by WorldClim ^[31]. We downloaded raster data for two RCPs (2.6 and 8.5) available for all four selected GCMs and for two time periods (year 2050-average for 2041 to 2060 and year 2070-average for 2061 to 2080). We selected RCP 2.6 since it is the most aggressive among all RCPs in terms of greenhouse gas emissions reductions ^[67,68]. Also, near-term warming projected under RCP 2.6 is greater than under RCP 4.5, even though the greenhouse gas forcing is lower ^[69]. RCP 8.5 is the most extreme scenario in that it entails the highest projected increase in the concentration of multiple greenhouse gases in the atmosphere ^[70] and associated increases in global surface temperatures ^[71].

We projected the bioclimatic-envelope models developed using present day conditions to potential future climatic conditions as simulated by a total of 16 GCM by RCP by year combinations. To avoid generating hundreds of map results, again we used the SAHM ensemble tool to produce ensemble maps for the future bioclimatic-envelope suitability of each bird species. Once ensembles were generated, we further combined results across GCMs within each RCPs. This led us to having two sets of maps with information on projected suitable conditions for each species, one for each of the years (2050 and 2070) according to each RCP. Finally, after considering the agreement (overlap) of at least three species distribution models and two GCMs, we compared the current and future ensemble maps to determine areas of stability, gains, and losses in suitable bioclimatic conditions between present day and the two projected years.

---

3. Results

Six of the 19 variables-isothermality, minimum temperature of the coldest month, mean temperature of the wettest quarter, mean temperature of the driest quarter, and precipitation of the warmest quarter-were identified as relevant to the biology and distribution of all nine species, as reflected in the modeling results (Table 1). The bioclimatic variable, maximum temperature of the warmest month, was a driver of change for two quail species-the northern/masked bobwhite quail and the scaled quail.

We saw variability in of the performance of the five statistical models across species (Table 2). For instance, for the Pinyon Jay, the AUC scores were highest for the BRT model, followed by RF, Maxent, MARS, and GLM. For the White-tailed Ptarmigan, BRT also had the highest AUC, however, it was followed by Maxent and RF, and then MARS and GLM had equal values. AUC values generally indicated strong accuracy of the models, with most AUC values ≥ 0.70. Only GLM for the Pinyon Jay resulted to a slightly lower value (AUC = 0.68). Among all nine birds, only models for the White-tailed Ptarmigan all had AUC values ≥ 0.90. Comparing the five statistical models, BRT performed the best in terms of AUCs for all species (0.89 to 0.97). The values of percent of occurrences correctly classified for BRT were also higher than for the rest of the SDMs (80.3% to 90.5%). After BRT, Maxent and RF performed fairly well in terms of AUC, ranging from 0.78 to 0.94 for Maxent and 0.82 to 0.92 for RF. Low percentages of occurrence points correctly classified ( < 70%) were found for the GLM model. The ensemble models consistently had AUCs for the test data above 0.90 for all nine species, thus surpassing all individual SDMs. Table 3 shows the results of the True Skill Statistic (TSS) and Kappa statistic (K), with values varying across species and models. Similar to the results for AUC, the BRT and RF models performed fairly well in terms of the TSS. For the Cassin's Sparrow, for instance, BRT scored the highest (TSS = 0.92) followed by RF (TSS = 0.67), while GLM scored the lowest among the five models with TSS = 0.56. Overall, GLM had the lowest accuracy based on both the TSS and the K relative to the other four models.

Among all bird species, models for the White-tailed Ptarmigan resulted in a good level of accuracy for the TSS values ( ≥ 0.80). The ensemble models showed improved TSS (all values ≥ 0.85) and K (all values ≥ 0.80) for all species relative to the individual SDMs. Other qualitative assessments of model performance included an inspection of the calibration and deviance of residual plots. No models had calibration plots that were without some form of over or under prediction, with models under predicting lower probabilities of occurrence and over predicting higher probabilities.

---

3.1. Present day bioclimatic-envelope models

The present day maps of suitable climate envelopes differed among the models. As an example, all models of the Scaled Quail showed high suitable conditions in the western region of the species' range, in accordance with most of the reported species occurrence points in the region (Figure 3). Further, Figure 3 shows how the results of the five statistical models differ between each other even with using the same species occurrence points. Among the five models, only RF and BRT identified highly suitable bioclimatic conditions in areas where presence points were mostly concentrated and low suitability in areas with fewer points. These two suitability maps (Figures 3a and 3e) are associated with fairly high AUC values in Table 2 for RF (0.83) and BRT (0.92). MARS and GLM highlighted suitable conditions in the western region, and Maxent accentuated climate condition suitability more in the southern end of the range.

For each species, locations where all five models identified potential suitable bioclimatic conditions (score = 5) for the species usually matched where species occurrence points were found (Figures 4a,4b,4c,4f,4g,4h, and 4i). However, this was not entirely true for the Mexican Spotted Owl (Figure 4d) and the Lesser Prairie-Chicken (Figure 4e), where high agreements among models (score = 5) were also found in the northern region of the species ranges. Worth noting was the Mexican Spotted Owl, for which the models identified a large, suitable climate envelope area where no species occurrence points were recorded.

The suitable area (in km²) for each of the nine bird species' bioclimatic-envelopes, as modeled using the current climate conditions, as shown in Table 4. These suitable areas are contained within the specified USGS species ranges (Figure 2). The calculation of the suitable area was based on the agreement of at least three SDMs. The smallest area was for the Montezuma Quail (16, 089 km²), followed by the White-tailed Ptarmigan (43, 818 km²). These species also have the two smallest species ranges.

---

3.2. Projecting models to future climatic conditions

Figures 5,6,7, and 8 are maps of the projected bioclimatic-envelopes of the nine bird species generated through comparison of ensembles of suitable bioclimatic-envelopes between present and projected future climatic conditions. These maps show areas where ensembles agree between present day and future (stable), future ensemble projects newly suitable conditions (gain), present ensemble may be converted to unsuitable conditions in the future (loss), and areas where conditions are unsuitable now and in the future (unsuitable). For comparison purposes, we added the current day, suitable bioclimatic-envelopes in the first column, for which the area (km²) for each species is found in Table 4.

The total size of the area containing suitable bioclimatic conditions increased from present day to future for the years 2050 and 2070 for RCP 2.6 (Figures 5 and 9) for Juniper Titmouse (16.5% & 8.2%), Pinyon Jay (12.3% & 9.0%), and Cassin's Sparrow (41.0% & 78.6%). For the Juniper Titmouse and Pinyon Jay, the gains were mostly located in the northern part of the species ranges. For the Cassin's Sparrow, the gains were around the center of the range. For RCP 2.6, both years 2050 and 2070 showed a slight increase of suitable climatic conditions for the Mexican Spotted Owl (7.9% & 6.7%) and the Lesser Prairie-Chicken (23.3% & 15.0%). For both future time periods and RCP 2.6, major loss of suitable climatic habitat was projected for the Lesser Prairie-Chicken (48.8% average), while Cassin's Sparrow had the least loss (10.0% average; Figure 5).

The years 2050 and 2070 under RCP 8.5 (Figures 6 and 9) showed comparable trends to RCP 2.6 in terms of the amounts of gains and losses for the Lesser Prairie-Chicken. The only major difference observed was the lesser gain in 2070, with an associated loss of 63.9%. Also, about 56.1% of area of currently suitable climatic conditions for Pinyon Jay was projected to become unsuitable by 2070.

Among the nine species of birds, the quails were the ones projected to be highly susceptible to climate change and appeared to be of most concern. The Northern/Masked Bobwhite Quail, Montezuma Quail, and White-tailed Ptarmigan all showed considerable losses of suitable climatic habitat for both modeled years and RCPs (Figures 7,8, and 10). The White-tailed Ptarmigan would lose 61.5% of its suitable habitat by 2050 and 66.7% by 2070, when both RCPs were averaged. Worse, gains of newly suitable regions for Montezuma Quail and White-tailed Ptarmigan were minimal. For the Scaled Quail, contraction would most likely occur in the northern part of the species range (Figures 7 and 8). A quantitative summary of the future projections shown in Figures 5 and 6 is provided in Figure 9, while for the future projections shown in Figures 7 and 8, the summary can be found in Figure 10. Bars in Figures 9 and 10 represent the percent increase (+) and decrease (-) in area of the present day suitable bioclimatic conditions for the focal species as projected according to RCPs 2.6 and 8.5 to the years 2050 and 2070. Decrease in area is related to loss, or area reduction, of the distribution of presently suitable conditions (Table 4). As an example, for Juniper Titmouse projected to 2050 (RCP 2.6), Figure 9 showed a percent increase of +16.5 and a percent decrease of-37.4, or a net of +20.9. This net positive value can be interpreted as a 20.9% increase of the present day suitable bioclimatic conditions for the Juniper Titmouse when projected to year 2050.

---

4. Discussion

---

4.1. Individual species modeling and ensemble performance.

The most important finding was the utility of the ensemble models, as produced by the SAHM ensemble function, for modeling the bioclimatic-envelopes. All ensemble efforts for all species resulted in AUCs > 0.90, showing more robust models than any individual SDMs. All five SDMs predicted different distributions of suitable bioclimatic conditions within species. Previous studies reported the same results ^{[52,72,73,74,75]}. Overall, values of most of the AUCs ( ≥ 0.70) indicated strong accuracy of the individual models.

Combining individual models to a single ensemble provided a more representative species distribution because it underscored the agreement among algorithms ^[36,52]. Although other studies (e.g., ^[75]) highlighted the limitations of the ensemble map in not outperforming individual models, in our case, the ensembles outperformed all five models. In other words, the fairly strong performance of the individual models led to a still stronger ensemble. Araújo and New ^[36] and Marmion et al. ^[76] agreed that the performance of ensemble models is dependent on the accuracy of the individual models.

Incorporating lower performing models (in our case the GLM and MARS) could make ensemble models less reliable than a single, high performing model ^[77]. Stohlgren et al. ^[52] emphasized that ensemble techniques are useful for newly arrived invasive species whose relationship to their environment has not yet been established. This study showed that ensembles could provide a valuable alternative to using a single SDM. In the future, we may test more SDMs, introduce ensemble thresholds to further exclude underperforming models, including those having low AUC values ^[77], and diversify the types of models selected to improve the performance of the ensembles ^[78].

---

4.2. Future projections: uncertainty and robustness

Although there were differences in results among the four climate models for both 2050 and 2070, all predictions entailed a significant reduction of the present day distribution of suitable bioclimatic conditions for most of the bird species, especially for quails. Suitability is projected to decline in the northern part of the Lesser Prairie-Chicken's species range. This projection goes against the northward movement often projected for species in response to climate change. This might be caused by the unusual pattern of precipitation in the area.

Though many robust features can be identified in future projections, multiple limitations associated with the projection of species distributions into the future under different climate scenarios have been documented. Three broad categories of uncertainties affecting the climate variables used to drive the SDMs include (a) uncertainties in future greenhouse gas concentrations ^[79], (b) limitations in the accuracy of GCM-simulated, large-scale physical climate responses to changing greenhouse gas levels ^[71], and (c) shortcomings and assumptions inherent to statistical downscaling methods used to refine GCM results to a finer level of spatial detail ^[80,81]. By utilizing data products derived from four GCMs and two RCPs, this study partially explores two of these three sources of climate variable uncertainty. Stoklosa et al. ^[82] specifically discuss approaches to account for some uncertainties in the climate variables used to drive SDMs. Furthermore, several authors have shown variability among future projections of suitable climatic conditions when using different climate models applied to the same species occurrence dataset ^[8,83]. This variability of results makes assessment of projections of future conditions a complex effort. First, there is no way to know which single SDM could provide the most accurate information for a species, although one could argue that the model with the highest accuracy in capturing the present day distribution of suitable climatic conditions may produce more accurate future projections. However, Thuiller ^[8] reasoned that even when a model has the highest AUC and K statistic, that does not mean the model provides the best estimate of the future distribution of suitable conditions, especially as every model is based on different assumptions. It is most fitting to use an ensemble model of future projections, as this ensemble represents the areas of agreement among individual model projections. The reliability of future conditions produced by ensembles may still be questioned, but ensemble results do incorporate the positive aspects of multiple models and provide a more conservative assessment of these conditions.

Our future projections showed that, although the nine bird species would be impacted by climate change, a vast area would still remain stable 50 years from now for most species. Even for the White-tailed Ptarmigan and Montezuma Quail, about 17% and 57%, respectively, of the distribution of presently suitable bioclimatic conditions would remain stable in 2070.

---

4.3. Caveats

We purposely excluded non-climatic variables, such as topography, vegetation cover, and land-use from our modeling effort. Our analysis did not disregard the fact that variables other than climate could add explanatory power in the model. However, we chose to determine the effects of the bioclimatic variables alone and did not want to restrict the distributions in present day and future by adding the non-climatic variables, especially for the more widely distributed species (e.g., Northern/Masked Bobwhite Quail, Pinyon Jay, Juniper Titmouse, and Cassin's Sparrow). We believe that bioclimatic modeling helps improve understanding of the vulnerability of the species of conservation concern considered here to climate change. Unless there are dramatic changes to non-climatic variables in the future (i.e., future changes in vegetation and land-use), the majority of any shifts in distribution of suitable conditions between present day and future would be driven by the climatic variables that we focused on. In other words, these non-climatic variables would not lead to dramatic changes in the future distributions of suitable conditions; they would simply limit the distribution of these suitable conditions further. Apart from the lack of datasets projected according to the RCPs, scale is also an issue as many vegetation and land-use datasets, for example, are available at finer resolutions than the climate projections. However, we acknowledge the role that finer-scale habitat changes, such as in water availability, could play in determining future distributions of species. Further studies on the impacts of climate change to these species of concern should address the interaction between climate and other non-climatic variables that may have profound effects on bird populations and their habitat. One important non-climatic variable that has been shown to adjust the future distribution of suitable areas is the dispersal ability of the species ^[84], which, unfortunately, we did not look cover in this research. Dispersal behavior is hard to incorporate in model projections as it could also change over time. By design, our models do not capture microclimate nor other local environmental characteristics such as soil moisture. However, they do provide a strong climate change context for future explorations of species population sustainability as affected by other drivers, including land use and water resources. There is a follow-up study after this that would look further beyond just bioclimatic effects on the bird species.

---

4.4. Biological relevance of projections

All five SDMs agreed that the distributions of the nine bird species within the areas where they currently occur contain the best climatic conditions for the species. However, a large area in the northern regions of the Mexican Spotted Owl and Lesser Prairie-Chicken ranges were also identified as currently containing suitable climatic conditions, even with few or no current species occurrences recorded. In the case of the Mexican Spotted Owl, the current climatic projection closely matches the species profile from the Environmental Conservation Online System (ECOS), mapping the home range within five States, including Colorado and Utah (both northern regions of the species' range). While it is unknown why Colorado and Utah support fewer owls ^[85], our projection indicated and confirmed the climatic suitability of these regions for the Mexican Spotted Owl.

As for the Lesser Prairie-Chicken, a large area in northwestern Kansas was identified by our models as currently containing suitable climatic conditions. The area is part of the Short Grass Prairie (SGPR), a type of habitat with which the Lesser Prairie-Chicken is associated ^[86]. Our results complement a recently developed habitat suitability model ^[87] that also showed a high level of suitability in the northwestern region of Kansas. The State of Kansas currently protects the most extensive remaining range of the lesser prairie-chicken among the five states-Kansas, Texas, New Mexico, Oklahoma, and Colorado-where the species occurs ^[88].

The climatic habitat loss predicted by our models for the White-tailed Ptarmigan was expected as the species are in alpine areas, which are likely to be affected by climate change under future conditions. The range of the ptarmigan is severely constrained by its dependence on the alpine areas that they occupy. As warmer temperatures continue, the species habitat is becoming unstable. Furthermore, warming also exacerbates the ecological instabilities caused by previous habitat degradation ^[89] that could further threaten the distributions of white-tailed ptarmigan.

Of all 19 bioclimatic variables used in the modeling, six were present in all of the 9 bird species. The minimum temperature of the coldest month, mean temperature of the wettest quarter, mean temperature of the driest quarter, and precipitation of the warmest quarter-were identified as relevant to the biology and distribution of all nine birds. Precipitation of the driest month was only relevant for the Mexican Spotted Owl and the Cassin's Sparrow, while precipitation seasonality seemed unnecessary in modeling the bioclimatic habitat for Scaled Quail, Cassin's Sparrow, and Lesser Prairie-Chicken.

---

5. Conclusion

Seven of the nine ensemble models for present day identified areas with potential suitable bioclimatic conditions that are currently occupied by the species. The other two models for the Mexican Spotted Owl and the Lesser Prairie-Chicken, after evaluation of existing studies, showed promising results as little gains were observed. Further, our results presented an opportunity for some of these other bioclimatically suitable areas to be evaluated as potential translocation sites for the Mexican Spotted Owl and the Lesser Prairie-Chicken. However, this has to be conducted with caution as there are many factors and environmental conditions at play that should be explored. In addition, locations that we identified as having suitable bioclimatic conditions could, following further evaluation and consideration of these other limiting factors, be given priority for conservation management.

We conclude that bioclimatic variables alone may not always correctly ascertain current or future suitable conditions. We recommend that suitable bioclimatic-envelopes be identified and that species-specific, non-climatic variables should be characterized to eliminate questionable modeled locations. Future research should incorporate temporal processes into the models to better capture the changing distributional patterns of the species distributions. For instance, bird species with a widespread geographical distribution may adapt to differences in local habitat availability or conditions across their range during movements in response to a change in climatic conditions. This process has to be accounted for when modeling species distributions, otherwise the SDM does not fully reveal the underlying processes that may drive expansions or shifts in species distributions.

---

Acknowledgments

The project described in this publication was supported by cooperative agreement G13AC00223 from the United States Geological Survey. Its contents are solely the responsibilities of the authors and do not necessarily represent the views of the South Central Climate Science Center or the USGS. This manuscript is submitted for publication with the understanding that the United States Government is authorized to reproduce and distribute reprints for governmental purposes. We thank the University of New Mexico, the Museum of Southwestern Biology, Esteban Muldavin, Rayo McCollough, and Mitch East from the Natural Heritage New Mexico for data and their review of this manuscript. We also acknowledge contributions from the Natural Heritage programs in Arizona and Texas; and all the other partner institutions of BISON and ORNIS (Appendix 2) for providing species occurrence datasets. We thank those that provided input on focal species from the Bureau of Land Management, George Miksch Sutton Avian Research Center, Natural Resources Conservation Service, NatureServe, New Mexico Department of Game and Fish, New Mexico State Land Office, Oklahoma Biological Survey, Oklahoma Department of Wildlife Conservation, Oklahoma State University, Playa Lakes Joint Venture, Texas Wildlife Department, The Nature Conservancy, U.S. Fish and Wildlife Service, and U.S. Forest Service. Science coordinators for the six LCCs that intersect the study area were also contacted.

---

Conflict of interest

The authors declare no conflict of interest.

---