Projections of Future Suitable Bioclimatic Conditions of Parthenogenetic Whiptails

This paper highlights the results of bioclimatic-envelope modeling of whiptail lizards belonging to the Aspidoscelis tesselata species group and related species. We utilized five species distribution models (SDM) including Generalized Linear Model, Random Forest, Boosted Regression Tree, Maxent and Multivariate Adaptive Regression Splines to develop the present day distributions of the species based on climate-driven models alone. We then projected future distributions of whiptails using data from four climate models run according to two greenhouse gas concentration scenarios (RCP 4.5 and RCP 8.5). Results of A. tesselata species group suggested that climate change will negatively affect the bioclimatic habitat and distribution of some species, while projecting gains in suitability for others. Furthermore, when the species group was analyzed together, climate projections changed for some species compared to when they were analyzed alone, suggesting significant loss of syntopic areas where suitable climatic conditions for more than two species would persist. In other words, syntopy within members of the species group will be drastically reduced according to future bioclimatic suitability projections in this study.


Introduction
Climate change affects biodiversity by modifying species habitats [1,2] and has been documented as the primary cause of abundance and distribution losses in many species [3,4].The changes could be damaging [5], such as a decrease of prairie wetland habitat and a decline of future wetland areas, leading to reduction in species population as habitats disappear [6,7].Climate change can also lead to shifts and contractions in species distributions and composition [2,8,9] and decreased production of male offspring by temperature-dependent species [10].In addition, shifts toward warmer temperatures could influence disease dynamics [11] and trigger outbreaks [12].While some species could adapt to climatic change because of their ability to disperse [13,14], other species (e.g.reptiles) may not be mobile enough to adapt to local climate pressures [7].
Although most reptiles can tolerate warmer temperatures due to their scale-covered skin [15] and sufficient mobility to evade thermal stress [16], their primary habitats are still vulnerable.For example, Sinervo et al. [17] suggested that if global temperature continues to rise, global extinction could average 16% by 2050 and 30% by 2080, while equatorial extinctions could reach 23% by 2050 and 40% by 2080.A study in Texas for Canyon Lizards (Sceloporus merriami) indicated that a 2 • C rise in air temperature could diminish species movement, causing energy shortfalls and population size reductions [18].Furthermore, experimental studies have revealed that phenotypic variance observed in reptile hatchlings is highly induced by physical conditions during embryogenesis [19,20].For example, changes in temperature during incubation affect developmental rate, size of hatchling and, in some cases, determines gender of the offspring [19,21,22].Behavioral and other morphological changes have also been reported [19,23,24].
Many studies have considered the diverse reptile species documented in Arizona and New Mexico [25][26][27][28] and Texas [29][30][31].However, we are unaware of any studies that evaluated future climate scenarios on habitats and distribution of members of the Aspidoscelis tesselata complex and related species.Species distribution models (SDMs) and associated exploration of model parameters are ideal tools to better understand future steps for species management and policy [32].Our primary objective was to develop models of present-day and potential future distributions of suitable environmental conditions for members of the A. tesselata complex, progenitor and other related species.
In this study, Aspidoscelis (Teiidae) is used in replacement of Cnemidophorus referring to a clade describing all whiptail species native to North America [33].This genus includes the A. tesselata species group, which consists of the triploid Colorado checkered whiptail (A.neotesselata), and two diploid species, the New Mexico whiptail (A.neomexicana) and the common checkered whiptail (A.tesselata; [34]).Also included in this study is the gray-checkered whiptail (A.dixoni; [34]), a diploid parthenogen no longer considered as a species [35].Except for A. neomexicana, hybrid-derived obligatory parthenogens of the A. tesselata species group comprise the A. tesselata complex [36].The unisexual, common checkered whiptail resulted from a single hybridization event between a female of the marbled whiptail (A.marmorata) and a male of the spotted whiptail (A.gularis septemvittata) [37].No projections were conducted on the marbled whiptail because a U.S. Geological Survey (USGS)-defined species range was not available at the time of this report.In the past, the Plateau spotted whiptail (A.scalaris) has been described as a subspecies of the common spotted whiptail (A.gularis) by Maslin and Secoy [38].Therefore, this study included analysis of bioclimatic projections for the common spotted whiptail, the Plateau spotted whiptail, and the pervasive tiger whiptail (A.tigris).
We applied SDMs [32,39] to project the availability of suitable bioclimatic conditions using climate projections derived from four general circulation models (GCMs) and two representative concentration pathways (RCPs).We then contrasted the future projected climate envelope suitability results produced for two future time periods (years 2050 and 2070).Our objectives were: to develop models of present-day and potential future distributions of suitable bioclimatic conditions for whiptail species in the U.S. region only; and to compare how bioclimatic envelope suitability is projected to change from present to future.

Materials and Methods
We followed three main steps in this study: (1) Selection and processing of whiptail species, the species range and the bioclimatic variables; (2) selection of SDMs and evaluation of the current bioclimatic conditions; and (3) selection of GCMs, RCPs and the projection of the current conditions to future conditions.

Study Area
Populations of the common checkered whiptail, tiger whiptail, Colorado checkered whiptail, Plateau spotted whiptail, common spotted whiptail, New Mexico whiptail and the gray-checkered whiptail occur within the southwest region of the United States (Figure 1).This region is known to have highly variable climate and high biodiversity.The Southwest is rich in cultural and natural resources along with expanding urbanization, which emphasizes the importance of predicting future climate and considering how changes may affect available suitable habitat [40].

Occurence Datasets
Species occurrence datasets were obtained from Natural Heritage programs in Arizona [41], Colorado [42], New Mexico [43] and Texas [44].There were a total of 15,663 occurrence points, distributed as follows: 1,262 for the common checkered whiptail; 11,002 for the tiger whiptail; 140 for

Occurence Datasets
Species occurrence datasets were obtained from Natural Heritage programs in Arizona [41], Colorado [42], New Mexico [43] and Texas [44].There were a total of 15,663 occurrence points, distributed as follows: 1,262 for the common checkered whiptail; 11,002 for the tiger whiptail; 140 for the Colorado checkered whiptail; 19 for the Plateau spotted whiptail; 2,469 for the common spotted whiptail; 664 for the New Mexico whiptail; and 107 for the gray-checkered whiptail.
Analysis of the A. tesselata species group consisted of a total of 2,173 occurrence points obtained by combining occurrence points from each species group member (A.tesselata, A. neotesselata, A. dixoni, and A. neomexicana).Separate range intersections between the tiger whiptail and common spotted whiptail, common checkered whiptail and Colorado checkered whiptail, tiger whiptail and New Mexico whiptail, common checkered whiptail and New Mexico whiptail, tiger whiptail, common spotted whiptail and common checkered whiptail, common checkered whiptail and New Mexico whiptail whiptails were determined based on phylogenetic relationship, documented geographic overlap.Only presence data was used, since no sufficient absence data was available.Other online sources of presence data used in this study included: Biodiversity Information Serving Our Nation (BISON) [45] and herpetological collections data [46].

Species Distribution Modeling
We utilized 19 raster-based bioclimatic variables that were obtained from among the WorldClim datasets [47].The set of variables was used to describe present environmental conditions and explore the relationship between bioclimatic conditions and species distribution patterns.WorldClim provides climate projections downscaled to 30 seconds, roughly 900 m at the equator.The entire list of climatic variables included in the analysis is shown in Table 1.The bioclimatic condition-species distribution relationship was analyzed using the following species distribution models/statistical algorithms for each species: Generalized Linear Model (GLM), Random Forest (RF) [48,49], Boosted Regression Tree (BRT) [50], Maxent [51,52] and Multivariate Adaptive Regression Splines (MARS) [53].These five SDMs were selected based on how they perform with presence-only data [53].The GLM is a linear regression adapted to binary count data.The method uses stepwise procedure to select covariates in the model.The MARS non-parametric algorithm build flexible models by fitting piecewise logistic regressions.Though it has similarities with GLM, MARS is better in accommodating non-linear responses to predictors and at the same time lessens the effects of outlying observations.The model RF uses decision trees through random grouping of the covariates.Random forest models both interactions of the variables and their nonlinear relationships and does not split the data into training and test as RF utilizes bootstrapping to fit individual trees.Like Random Forest, BRT also uses decision trees, but the method is robust to missing observations.BRT uses cross-validation by choosing models based on model comparisons of evaluation metrics [50].Maxent is best for presence-only modeling.While observed absence is valuable in modeling, data is oftentimes not available and using only presence data is unavoidable [54].The modeling tool Software for Assisted Habitat Modeling (SAHM) run within VisTrails [54,55] was used 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 [52]) within a 95% minimum convex polygon defined by the presence data.Species ranges provided by the USGS National Gap Analysis Program (GAP) [56] were used to generate a template layer for each species.This template restricted model development and projection to the present-day geographic range of the species based on 8-digit hydrologic unit codes.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 have negative effects on the relationship between background and presence points.To avoid the mismatch between species extent and presence/absence points, we used the USGS GAP geographic range of the species to define where background points should be assigned.
For each whiptail, one of each pair of highly correlated (r > 0.7) [57] variables was removed from the bioclimatic-envelope models to avoid collinearity among variables [58].Choices between variables were made based on the results of a species-specific literature search.In particular, variables that were identified in one or more studies regarding the species of interest as having an effect on the species' range or population dynamics were selected (see Appendix A).In cases where the results of the literature search could not differentiate between two highly correlated climatic variables, 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 was used [54].Through the SAHM ensemble tool, combination maps were produced for the current distribution for each species.A combined map is a summation of binary maps generated from probability maps output from each statistical modeling algorithm [59][60][61].The threshold was optimized by using specificity = sensitivity in discretizing the probability maps [62].Final map combinations consisted of pixel values that represented the number of models in agreement to indicate 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 5 meant there was agreement across all five models.
Confidence in individual model results was assessed in terms of concordance among the different distribution models.Confidence that bioclimatic conditions were suitable for a species was considered when three or more (at least 60% of) algorithms were in agreement [63].Information was compiled on various measures of model performance, including the Area Under the Receiver Operating Characteristic (ROC) Curve (AUC) for the test data, correct classification rate (Co%) [64,65] and the True Skill Statistic (TSS) [66] for each algorithm and species combination.The AUC value is the probability that the model would rank a randomly chosen presence observation higher than the randomly chosen absence observation.Swets [67] 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.The TSS indicates agreement between model predictions and true values of presence or absence for the species occurrence data.It is presented as an improved measure of model accuracy that, unlike the kappa statistics [68], is not dependent on species prevalence (i.e., proportion of occurrence points for which the species is present) [66].Other qualitative assessments of model performance such as calibration and deviance of residual plots were checked, which included the inspection of 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 layer missing from the model inputs [55].

Projection to Future Conditions
GCMs were screened based on previously published evaluations of model performance [69,70] 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).Values for bias of model output relative to observed historical data as an important criterion to exclude GCMs were used.In particular, GCMs were excluded 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, large values (>1 or <−1) were used for top of atmosphere energy imbalance (W•m −2 ) (imbalance in energy inflow and outflow of the earth system at the top of the atmosphere) to exclude models since these values may be an indication of long-term drift in simulated climatic conditions.The final list of selected GCMs include: Community Climate System Model version 4 (CCSM4; [71]), Hadley Centre Global Environment Model version 2-Earth System (HadGEM2-ES; [72]), Model for Interdisciplinary Research on Climate version 5 (MIROC5; [73]) and Max Planck Institute Earth System Model, low resolution (MPI-ESM-LR; [74]).
For future conditions, the downscaled data provided by WorldClim [47] was used.Raster data was downloaded for two Representative Concentration Pathways (RCP 4.5 and RCP 8.5) available for all selected GCMs and for two time periods (year 2050-average for 2041 to 2060-and year 2070-average for 2061 to 2080).RCP 4.5 was selected over RCP 2.6 since it is more or less stable throughout the century among all RCPs in terms of greenhouse gas emissions reductions [75,76].RCP 2.6, in contrast to RCP 4.5, has a lower greenhouse gas forcing [77].For RCP 8.5, it is the most extreme scenario in that it entails the highest projected increase in the concentration of multiple greenhouse gases in the atmosphere [78] and associated increases in global surface temperatures [79].
Results of the current species distributions were projected to the future using the data from the four GCMs and according to the two RCPs.To avoid generating hundreds of map results, again the SAHM ensemble tool was used to produce combination maps for the future distribution of each species.Each RCP result from the four GCMs was combined.In the end, a set of projection maps for the year 2050 and another set for the year 2070 for each RCP and for each species was generated.Finally, the current and future combined maps were compared to determine stability, gains, and losses in suitable bioclimatic envelopes for the two projected years.

Individual Model Performance
Table 2 presents the performance of the five statistical models as applied to the seven whiptail species.Among the models, the BRT performed the best in terms of the AUC (with all AUC values ≥0.85) and the percentages of occurrence points correctly classified (with all %Co values ≥0.78), although it was only successful in modeling four out of seven species.Models like the RF and Maxent also performed well with most AUC values ≥0.75.The two other models, GLM and MARS, showed lower AUCs compared to the rest.AUCs for GLM ranged from 0.55 to 0.81, while for MARS ranged from 0.66 to 0.80.Although MARS was not the best performing model, it was most successful in modeling all seven whiptail species regardless of the sample size.Maxent, like MARS, also worked well with fewer occurrences.Table 3 shows the results of the True Skill Statistic (TSS), with values varying across species and models.Similar to the results of the AUC, the BRT model performed fairly well in terms of the TSS.Next to BRT, RF and Maxent models performed fairly well in terms of the TSS.

Future Climatic Envelope
Figures 2-8 show the potential climatic envelope areas for the seven whiptail lizards as simulated by the four selected GCMs: A. tesselata (Figure 2), A. tigris (Figure 3), A. neotesselata (Figure 4), A. scalaris (Figure 5), A. gularis (Figure 6), A. neomexicana (Figure 7) and A. dixoni (Figure 8).Except for the projection for the Gray-Checkered Whiptail (Figure 8), all other future projections for the rest of the species identified potentially suitable bioclimatic conditions.
For the common checkered whiptail, most of the future loses (averaged loss of 29% and 33% for 2050 and 2070, respectively) were concentrated in the northern part of the species range (Figure 2).About 50% of the present suitable bioclimatic conditions would remain stable according to both year projections.A lesser percentage of area lost was observed for the tiger whiptail compared to the increase of suitable bioclimatic conditions (11% loss vs. 71% gain).The models in Figure 3 suggested that the tiger whiptail could spread northward within its species range.Whiptails that showed major losses of suitable bioclimatic envelopes were the Colorado checkered whiptail (Figure 4), Plateau spotted whiptail (Figure 5) and the gray-checkered whiptail (Figure 8).In fact, for the Colorado checkered whiptail, the models did not project any gain of suitable areas; instead, about 58% of the present suitable bioclimatic conditions would be lost by 2070.The gray-checkered whiptail was the only species to have lost all its present suitable bioclimatic envelopes; both 2050 and 2070 projections showed the same results with a minimal increase of area of about 9%.While 71% of the present suitable bioclimatic conditions would remain stable for the Plateau spotted whiptail, our models showed minimal to no increase of suitable bioclimatic conditions even within its range.
The species group, which includes A. tesselata, A. neotesselata, A. dixoni, and A. neomexicana in Figure 9, showed a different story for the species range of the Colorado checkered whiptail.More than 75% of the current suitable areas would remain stable in both 2050 and 2070 when species were combined.In comparison, individual model results for A. tesselata and A. neomexicana revealed the western region of the range would remain stable for both species.
Among the seven whiptails, two showed promising future bioclimatic suitability conditions: the common spotted whiptail; and the New Mexico whiptail.With only 9% of the present suitable bioclimatic conditions would be lost, the New Mexico whiptail could see 115% increase in suitable bioclimatic envelopes within its range in 2070.Models for the common spotted whiptail showed 87% increase.
Climate 2017, 5, 34 9 of 21 Whiptails that showed major losses of suitable bioclimatic envelopes were the Colorado checkered whiptail (Figure 4), Plateau spotted whiptail (Figure 5) and the gray-checkered whiptail (Figure 8).In fact, for the Colorado checkered whiptail, the models did not project any gain of suitable areas; instead, about 58% of the present suitable bioclimatic conditions would be lost by 2070.The gray-checkered whiptail was the only species to have lost all its present suitable bioclimatic envelopes; both 2050 and 2070 projections showed the same results with a minimal increase of area of about 9%.While 71% of the present suitable bioclimatic conditions would remain stable for the Plateau spotted whiptail, our models showed minimal to no increase of suitable bioclimatic conditions even within its range.
The species group, which includes A. tesselata, A. neotesselata, A. dixoni, and A. neomexicana in Figure 9, showed a different story for the species range of the Colorado checkered whiptail.More than 75% of the current suitable areas would remain stable in both 2050 and 2070 when species were combined.In comparison, individual model results for A. tesselata and A. neomexicana revealed the western region of the range would remain stable for both species.
Among the seven whiptails, two showed promising future bioclimatic suitability conditions: the common spotted whiptail; and the New Mexico whiptail.With only 9% of the present suitable bioclimatic conditions would be lost, the New Mexico whiptail could see 115% increase in suitable bioclimatic envelopes within its range in 2070.Models for the common spotted whiptail showed 87% increase.show areas where present and future projections agree (stable), future estimate projects new suitable conditions (gain), present estimate may be converted to unsuitable in the future (loss), and areas where conditions are unsuitable in the future (non).For the exact location of the species range, see Figure 1.In Figure 10, the Colorado checkered whiptail clearly showed the absence of any gains in both projected years, and some minute gains were observed for the Plateau spotted whiptail.The gainloss bar graphs for the common spotted whiptail exhibited quite similar amounts in all four climatic scenario combinations.Only the model for the New Mexico whiptail has shown gains of more than 100%.In Figure 10, the Colorado checkered whiptail clearly showed the absence of any gains in both projected years, and some minute gains were observed for the Plateau spotted whiptail.The gain-loss bar graphs for the common spotted whiptail exhibited quite similar amounts in all four climatic scenario combinations.Only the model for the New Mexico whiptail has shown gains of more than 100%.In Figure 10, the Colorado checkered whiptail clearly showed the absence of any gains in both projected years, and some minute gains were observed for the Plateau spotted whiptail.The gainloss bar graphs for the common spotted whiptail exhibited quite similar amounts in all four climatic scenario combinations.Only the model for the New Mexico whiptail has shown gains of more than 100%.

Discussion
Evidence in this study suggests members of the A. tesselata complex and related species have exclusive bioclimatic suitability parameters that will be drastically affected by climate change.Warmer global temperatures caused by increasing concentrations of greenhouse gases [80,81] have resulted in shifts in rainfall and temperature patterns across a large geographic scale [82,83].Effects of warmer temperatures have been documented in many reptiles [21,22] and include gender determination, mobility and dispersal restrictions and other behavioral and morphological changes [10,18,19,84].In addition, climate change may affect the availability of resources, such as prey and vegetation structure [85,86].Potentially, female ectotherms can ameliorate the effects of climate change on their offspring through basking behaviors in viviparous species and through nest-site

Discussion
Evidence in this study suggests members of the A. tesselata complex and related species have exclusive bioclimatic suitability parameters that will be drastically affected by climate change.Warmer global temperatures caused by increasing concentrations of greenhouse gases [80,81] have resulted in shifts in rainfall and temperature patterns across a large geographic scale [82,83].Effects of warmer temperatures have been documented in many reptiles [21,22] and include gender determination, mobility and dispersal restrictions and other behavioral and morphological changes [10,18,19,84].In addition, climate change may affect the availability of resources, such as prey and vegetation structure [85,86].Potentially, female ectotherms can ameliorate the effects of climate change on their offspring through basking behaviors in viviparous species and through nest-site selection in oviparous species [19].In this study, projections of bioclimatic suitability for the bisexual common spotted whiptail and the all-female New Mexico whiptail showed significant gains across their current range.Even though sexual reproduction would affect allelic frequencies responsible for specific offspring characteristics, variation in offspring phenotypes of parthenogenetic species results from alterations in maternal behavior, such as nest-site selection [19]; results suggest no disadvantage regarding reproductive mode.The effect of warming temperatures on nest-site locations and hatchlings of vulnerable whiptail species could be the focus of future research.
In this study, projections were strongly influenced by life history differences for each species, and in some cases, for each pattern class.For example, the common checkered whiptail, which occurs from southeastern Colorado [87] to extreme southeastern Arizona, New Mexico and western Texas [88,89], occupies a habitat characterized by sparsely vegetated areas such as canyon slopes, gullies, flatlands and bluffs [88].This species would experience significant bioclimatic suitability losses (29% and 33%) for the two future scenarios.Furthermore, evidence in previous studies suggests that different clones or pattern classes of this species occupy different habitats; in Presidio County, Texas, a clone occupies gravelly alluvial benches, while another occurs in sandy flood plains [90].Recently, several localities have showed declines due to agricultural and urban sprawl [87,91].Projected suitability losses are concentrated within the northern range for the species, potentially affecting only the clones occupying that area and further reducing interaction with parthenogens found there (i.e., Colorado checkered whiptail).
Contrastingly, analysis of future bioclimatic suitability revealed the tiger whiptail would expand northward within its current range.This generalist and bisexual species occurs in deserts, woodlands and other sparsely vegetated areas throughout Idaho, Oregon, west California and Texas into southern Baja California, Sinaloa and Coahuila, Mexico [92].Currently, there are no major threats identified for this species [93], and a subspecies, the Sonoran tiger whiptail (A.tigris punctinealis), has been considered to impact populations of the gray-checkered whiptail through hybridization, as the resulting hybrids are not capable of parthenogenetic cloning [94,95].Projected areas that show a gain in environmental suitability in the future, or that would remain stable for the tiger whiptail overlap with areas where climate suitability would be completely lost for the gray-checkered whiptail.
The Colorado checkered whiptail occurs only in southeastern Colorado at elevations below 2135 m [87,91,92] and occupies plains, grasslands and juniper woodland in arroyos, canyons, valleys, hillsides, but also parks and areas with human disturbance [91].This species is currently listed as Near Threatened due to loss of habitat and continues to be affected by fragmentation [93].In some areas, this species has been extirpated, while populations in already disturbed areas such as parks, buildings and rural landfills persist [91].Results of this study indicate there would be no gains in projected bioclimatic suitability, suggesting further fragmentation within its range for 2050 and 2070.
The Plateau spotted whiptail ranges from southwest Texas into northern Mexico [96] and occurs in sparsely vegetated areas including canyons and mountains [97,98].While no major threats have been identified [93], results of this study have revealed major losses of projected bioclimatic envelopes, and estimate that approximately 71% of the current suitable bioclimatic envelope would remain stable.In 2050, RCP 4.5 shows a major portion of its current range remaining stable, but experiencing major losses in the western portion of its current range by 2070.An earlier onset of loss of bioclimatic suitability is demonstrable by 2050 in RCP 8.5 projections.Several reports have suggested possible hybridization between A. scalaris and A. gularis [99,100].
The common spotted whiptail [35,101] is found from southern Oklahoma, Texas and southeastern New Mexico, while it remains unclear if populations in Aguascalientes, Queretaro and Veracruz, Mexico should be assigned to this taxon [88].Habitats include shortgrass prairie, desert grassland, plateaus, weedy areas, shrubby river bottoms, washes and rocky slopes, but it is also found in other sparsely vegetated areas characterized by sandy, rocky or gravelly soils [88,92,98].This whiptail is described as common within its range in Texas [89,98] and has been listed as Least Concern due wide range and number of stable populations [93].In this study, projected effects on bioclimatic suitability estimated only a 9% loss across its current range.Significant gains in climate suitability are projected for both years in the eastern region of its current distribution range, while experiencing some minor losses in the western region.The New Mexico whiptail, which occurs in New Mexico, northwestern Texas and neighboring Chihuahua, Mexico [88,92], is found in some protected areas and has no major threats [93].In general, habitat is perpetually disturbed and includes grasslands with scattered shrubs, river basins, arroyos, vacant lots and mesquite-creosote associations, but also desert and grassland ecotones and shrubby edges of desert playas [88].It is considered rare in higher elevations within pinyon-juniper woodlands where sandy alluvial benches are exposed [88,90].Populations in Conchas Lake and Fort Sumner in New Mexico are considered natural, while the populations at Petrified Forest National Park in Arizona are more likely introductions [93].Findings of this study suggest an increase of about 115% over the present suitable bioclimatic envelope for this species and only a small loss in the southwestern region of its current range.
According to our study, the entire present suitable bioclimatic envelopes would be lost for both 2050 and 2070 projections for the gray-checkered whiptail, no longer recognized as a species [35] but once listed as Near Threatened under criterion B1 [93].In the past, it was suggested that populations of this whiptail in Texas were benefitting from the grazing practices employed in that area [102].However, recent reports identified declines due to drought, mesquite invasion and exotic forage grasses [93].Similarly, threats in New Mexico include habitat conversion, BLM chemical brush control, overgrazing, mining activities and unregulated collecting [103,104].This whiptail was described from two separate areas: one in the lower ranges of the Chinati Mountains in southwestern Presidio County, Texas [89]; and the other near Antelope Pass in the Peloncillo Mountains, Hidalgo County, New Mexico [88].Habitat in Texas included dry washes, canyon bottoms, rocky plains and desert scrub [98], while in New Mexico, habitat is sandy or gravelly creosote bush flats [88].
The Aspidoscelis tesselata species group [34] has been the focus of numerous studies of allozymes, chromosomes, mitochondrial DNA and morphology [91].Each member includes three or more color pattern classes, or informal descriptions, first documented by Zweifel [105].Interestingly, each pattern class is found within well-defined geographic areas and can be distinguished by variable dorsal patterns.Several of these species and their respective pattern classes occur in syntopy in Arizona, New Mexico, Colorado and Texas [91,[106][107][108][109].Results of this study suggest that climate change will negatively affect the habitat and distribution of pattern classes within A. tesselata, while projecting gains in suitability for others.It is worth mentioning that the negative effects of the climate for RCP 4.5 scenario at year 2070 is quite similar to those in RCP 8.5 scenario for the year 2050 for almost all whiptail species.Furthermore, when the species group was analyzed together, climate projections changed for some species compared to when they were analyzed alone, suggesting significant loss of syntopic areas where suitable environmental conditions for more than two species would persist.In other words, syntopy within members of the species group and their pattern classes will be drastically reduced according to future environmental suitability projections in this study, but also as a result of the current threats of habitat loss and fragmentation for these taxa [93].

Caveats
While our results showed that some models outperformed the others, we purposely did not run an ensemble model that would have improved the bioclimatic condition selection.By combining individual models through an ensemble, thresholds could be introduced to exclude underperforming SDMs, those having low AUC values, and differentiate the types of models selected to improve the overall performance of the ensemble.The models we ran were based on the climate data alone and the projections were based on the occurrence records that were provided.We thought that unless there are dramatic changes to these non-climatic variables in 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.Apart from the lack of datasets projected according to the RCPs, scale is also an issue as vegetation and land-use datasets, for instance, are available at finer resolutions than the climate projections.However, we acknowledge the role that finer scale habitat changes, such as water availability, could play a significant role in species distributions.One important non-climatic variable that has shown to adjust the future distribution of suitable areas is the dispersal ability of the species [110].However, dispersal behavior is hard to incorporate in model projections as it could also change over time [111].Finally, we did not extend the study area to the south beyond the US-Mexico border and modeled only the U.S. portion of the species range.While the range of the species for the Mexican region could be available through the IUCN website [93], the occurrence points are nonexistent.Other third-party sources of occurrence points, such as the citizen science iNaturalist [112], have very limited data.For instance, the wide species range of the common checkered whiptail in Mexico has only one occurrence data available.The inclusion of this single data into our model could create accuracy problems, especially for the species group.The Plateau spotted whiptail has the same problem as well, with only five occurrence points available.Including the entire range of the species, including Mexico, would have enhanced our analysis.Although this limitation restricts the spatial extent of our projections, the inclusion of the Mexican region is beyond the scope of this research and should not disregard the extensive patterns of bioclimatic-envelope changes observed in the U.S. region of the species range.

Figure 1 .
Figure 1.In (a), specific ranges of the seven whiptail lizard species analyzed in this study: common checkered whiptail (Aspidoscelis tesselata), tiger whiptail (A.tigris), Colorado checkered whiptail (A.neotesselata), Plateau spotted whiptail (A.scalaris), common spotted whiptail (A.gularis), New Mexico whiptail (A.neomexicana) and the gray-checkered whiptail (A.dixoni)-as provided by USGS.Intersections of ranges are shown (in green) for (b) tiger whiptail and common spotted whiptail; (c) common checkered whiptail and Colorado checkered whiptail; (d) tiger whiptail and New Mexico whiptail; (e) common checkered whiptail and New Mexico whiptail; (f) tiger whiptail, common spotted whiptail, and common checkered whiptail; (g) common checkered whiptail and New Mexico whiptail.The ranges for the species group in (h) include A. tesselata, A. neotesselata, A. dixoni and A. neomexicana.

Figure 1 .
Figure 1.In (a), specific ranges of the seven whiptail lizard species analyzed in this study: common checkered whiptail (Aspidoscelis tesselata), tiger whiptail (A.tigris), Colorado checkered whiptail (A.neotesselata), Plateau spotted whiptail (A.scalaris), common spotted whiptail (A.gularis), New Mexico whiptail (A.neomexicana) and the gray-checkered whiptail (A.dixoni)-as provided by USGS.Intersections of ranges are shown (in green) for (b) tiger whiptail and common spotted whiptail; (c) common checkered whiptail and Colorado checkered whiptail; (d) tiger whiptail and New Mexico whiptail; (e) common checkered whiptail and New Mexico whiptail; (f) tiger whiptail, common spotted whiptail, and common checkered whiptail; (g) common checkered whiptail and New Mexico whiptail.The ranges for the species group in (h) include A. tesselata, A. neotesselata, A. dixoni and A. neomexicana.

Figure 2 .
Figure 2. Projected climatic envelope for the common checkered whiptail (Aspidoscelis tesselata) using (a) RCP 4.5, year 2050 (b) RCP 4.5, year 2070 (c) RCP 8.5, year 2050 (d) RCP 8.5, year 2070.Maps show areas where present and future projections agree (stable), future estimate projects new suitable conditions (gain), present estimate may be converted to unsuitable in the future (loss), and areas where conditions are unsuitable in the future (non).

Figure 3 .
Figure 3. Projected climatic envelope for the Tiger Whiptail (Aspidoscelis tigris) using (a) RCP 4.5, year 2050 (b) RCP 4.5, year 2070 (c) RCP 8.5, year 2050 (d) RCP 8.5, year 2070.Maps show areas where present and future projections agree (stable), future estimate projects new suitable conditions (gain), present estimate may be converted to unsuitable in the future (loss), and areas where conditions are unsuitable in the future (non).

Figure 2 .of 21 Figure 2 .
Figure 2. Projected climatic envelope for the common checkered whiptail (Aspidoscelis tesselata) using (a) RCP 4.5, year 2050 (b) RCP 4.5, year 2070 (c) RCP 8.5, year 2050 (d) RCP 8.5, year 2070.Maps show areas where present and future projections agree (stable), future estimate projects new suitable conditions (gain), present estimate may be converted to unsuitable in the future (loss), and areas where conditions are unsuitable in the future (non).

Figure 3 .Figure 3 .
Figure 3. Projected climatic envelope for the Tiger Whiptail (Aspidoscelis tigris) using (a) RCP 4.5, year 2050 (b) RCP 4.5, year 2070 (c) RCP 8.5, year 2050 (d) RCP 8.5, year 2070.Maps show areas where present and future projections agree (stable), future estimate projects new suitable conditions (gain), present estimate may be converted to unsuitable in the future (loss), and areas where conditions are unsuitable in the future (non).

Figure 4 .
Figure 4. Projected climatic envelope for the Colorado Checkered Whiptail (Aspidoscelis neotesselata) using (a) RCP 4.5, year 2050 (b) RCP 4.5, year 2070 (c) RCP 8.5, year 2050 (d) RCP 8.5, year 2070.Mapsshow areas where present and future projections agree (stable), future estimate projects new suitable conditions (gain), present estimate may be converted to unsuitable in the future (loss), and areas where conditions are unsuitable in the future (non).

Figure 4 .
Figure 4. Projected climatic envelope for the Colorado Checkered Whiptail (Aspidoscelis neotesselata) using (a) RCP 4.5, year 2050 (b) RCP 4.5, year 2070 (c) RCP 8.5, year 2050 (d) RCP 8.5, year 2070.Mapsshow areas where present and future projections agree (stable), future estimate projects new suitable conditions (gain), present estimate may be converted to unsuitable in the future (loss), and areas where conditions are unsuitable in the future (non).

Figure 5 .
Figure 5. Projected climatic envelope for the Plateau Spotted Whiptail (Aspidoscelis scalaris) using (a) RCP 4.5, year 2050 (b) RCP 4.5, year 2070 (c) RCP 8.5, year 2050 (d) RCP 8.5, year 2070.Maps show areas where present and future projections agree (stable), future estimate projects new suitable conditions (gain), present estimate may be converted to unsuitable in the future (loss), and areas where conditions are unsuitable in the future (non).

Figure 6 .Figure 5 .
Figure 6.Projected climatic envelope for the Common Spotted Whiptail (Aspidoscelis gularis) using (a) RCP 4.5, year 2050 (b) RCP 4.5, year 2070 (c) RCP 8.5, year 2050 (d) RCP 8.5, year 2070.Maps show areas where present and future projections agree (stable), future estimate projects new suitable conditions (gain), present estimate may be converted to unsuitable in the future (loss), and areas where conditions are unsuitable in the future (non).

Figure 5 .
Figure 5. Projected climatic envelope for the Plateau Spotted Whiptail (Aspidoscelis scalaris) using (a) RCP 4.5, year 2050 (b) RCP 4.5, year 2070 (c) RCP 8.5, year 2050 (d) RCP 8.5, year 2070.Maps show areas where present and future projections agree (stable), future estimate projects new suitable conditions (gain), present estimate may be converted to unsuitable in the future (loss), and areas where conditions are unsuitable in the future (non).

Figure 6 .Figure 6 .
Figure 6.Projected climatic envelope for the Common Spotted Whiptail (Aspidoscelis gularis) using (a) RCP 4.5, year 2050 (b) RCP 4.5, year 2070 (c) RCP 8.5, year 2050 (d) RCP 8.5, year 2070.Maps show areas where present and future projections agree (stable), future estimate projects new suitable conditions (gain), present estimate may be converted to unsuitable in the future (loss), and areas where conditions are unsuitable in the future (non).

Figure 7 .
Figure 7. Projected climatic envelope for the New Mexico Whiptail (Aspidoscelis neomexicana) using (a) RCP 4.5, year 2050 (b) RCP 4.5, year 2070 (c) RCP 8.5, year 2050 (d) RCP 8.5, year 2070.Maps show areas where present and future projections agree (stable), future estimate projects new suitable conditions (gain), present estimate may be converted to unsuitable in the future (loss), and areas where conditions are unsuitable in the future (non).

Figure 8 .Figure 7 . 21 Figure 7 .
Figure 8. Projected climatic envelope for the Gray-Checkered Whiptail (Aspidoscelis dixoni) using (a) RCP 4.5, year 2050 (b) RCP 4.5, year 2070 (c) RCP 8.5, year 2050 (d) RCP 8.5, year 2070.Maps showareas where present and future projections agree (stable), future estimate projects new suitable conditions (gain), present estimate may be converted to unsuitable in the future (loss), and areas where conditions are unsuitable in the future (non).For the exact location of the species range, see Figure1.

Figure 8 .Figure 8 .
Figure 8. Projected climatic envelope for the Gray-Checkered Whiptail (Aspidoscelis dixoni) using (a) RCP 4.5, year 2050 (b) RCP 4.5, year 2070 (c) RCP 8.5, year 2050 (d) RCP 8.5, year 2070.Maps showareas where present and future projections agree (stable), future estimate projects new suitable conditions (gain), present estimate may be converted to unsuitable in the future (loss), and areas where conditions are unsuitable in the future (non).For the exact location of the species range, see Figure1.

Figure 9 .
Figure 9. Projected climatic envelope for the Aspidoscelis tesselata species group using (a) RCP 4.5, year 2050 (b) RCP 4.5, year 2070 (c) RCP 8.5, year 2050 (d) RCP 8.5, year 2070.Maps show areas where present and future projections agree (stable), future estimate projects new suitable conditions (gain), present estimate may be converted to unsuitable in the future (loss), and areas where conditions are unsuitable in the future (non).

Climate 2017, 5 , 34 12 of 21 Figure 9 .
Figure 9. Projected climatic envelope for the Aspidoscelis tesselata species group using (a) RCP 4.5, year 2050 (b) RCP 4.5, year 2070 (c) RCP 8.5, year 2050 (d) RCP 8.5, year 2070.Maps show areas where present and future projections agree (stable), future estimate projects new suitable conditions (gain), present estimate may be converted to unsuitable in the future (loss), and areas where conditions are unsuitable in the future (non).

Figure 10 .
Figure10.This is a quantitative summary of the resulting projections shown in Figures2-9for all whiptail species including the species group according to these scenario combinations: RCP 4.5 to the year 2050, RCP 4.5 to the year 2070, RCP 8.5 to the year 2050, and RCP 8.5 to the year 2070.Increase in area (positive value) is related to the gain or area expansion of the present bioclimatic suitable conditions.Decrease in area (negative value) is related to loss or area reduction of the present bioclimatic suitable conditions.

Figure
Figure This is a quantitative summary of the resulting projections shown in Figures2-9for all whiptail species including the species group according to these scenario combinations: RCP 4.5 to the year 2050, RCP 4.5 to the year 2070, RCP 8.5 to the year 2050, and RCP 8.5 to the year 2070.Increase in area (positive value) is related to the gain or area expansion of the present bioclimatic suitable conditions.Decrease in area (negative value) is related to loss or area reduction of the present bioclimatic suitable conditions.

Table 1 .
List of 19 bioclimatic variables used in bioclimatic-envelope model development.Names and descriptions are based on WorldClim [47].

Table 2 .
The Areas under the Curve (AUC) associated with the test data and the percentages of occurrence points correctly classified (%Co) for the five different models.Model abbreviations are as follows: GLM = Generalized Linear Model, MARS = Multivariate Adaptive Regression Splines, BRT = Boosted Regression Tree, and RF = Random Forest.Species abbreviations are as follows: Common Checkered Whiptail (Aspidoscelis tesselata), Tiger Whiptail (Aspidoscelis tigris), Colorado Checkered Whiptail (Aspidoscelis neotesselata), Plateau Spotted Whiptail (Aspidoscelis scalaris), Common Spotted Whiptail (Aspidoscelis gularis), New Mexico Whiptail (Aspidoscelis neomexicana), and Gray-Checkered Whiptail (Aspidoscelis dixoni).NA values are due to the model not executing successfully due to an error related to sample size.