Assessing the Impacts of Climate Change on the At-Risk Species Anaxyrus microscaphus (The Arizona Toad): A Local and Range-Wide Habitat Suitability Analysis

: Anaxyrus microscaphus (The Arizona Toad) is an at-risk species that is endemic to the southwestern United States. Despite conservation concerns, little is known about the ecological drivers of its distribution and habitat use. We investigated the potential distribution of A. microscaphus at the range-wide scale and local scales (i.e., Zion National Park), using MaxEnt to model habitat suitability under current and future climate scenarios. Our models incorporated 12 environmental variables, including climatic, geomorphological, and remotely sensed data. The results showed good model accuracy, with temperature and elevation being the top contributing variables. Currently, 42.6% of the park’s area provides a suitable habitat for A. microscaphus , but projections for 2050 and 2070 indicate a signiﬁcant reduction in suitable habitat across its range. Temperature was the most inﬂuential variable, with habitat suitability decreasing as the annual mean temperatures exceeded 10 ◦ C. Precipitation, vegetation, and topography variables also signiﬁcantly contributed to the models. The most suitable habitat within Zion National Park occurred along sloped rivers and streams and in valleys with sandy soils, emphasizing the importance of riparian habitat conservation for A. microscaphus survival and persistence. As climate change progresses, the species’ habitat is expected to become increasingly constrained across local and range-wide scales. Our models demonstrated a shift in the suitable habitat towards major river systems, indicating a potential reliance on larger permanent river systems as smaller, more ephemeral habitats decrease in size and abundance. Future management strategies should prioritize conserving and enhancing the resilience of these habitats. MaxEnt models can guide population survey efforts and facilitate the identiﬁcation of priority conservation areas, saving time and resources for species of concern such as A. microscaphus . Further research, including ﬁeld surveys and large-scale analyses, is necessary to further reﬁne our understanding of this species’ distribution and how it may be impacted by climate and habitat change.


Introduction
Since the incorporation of new statistical methods and geographic information system (GIS) tools, the development of predictive species distribution models (SDMs) has greatly expanded the fields of ecology, biogeography, and conservation [1][2][3][4].SDMs describe how climatic and environmental factors relate to occurrence locations in geographic space, in order to delineate suitable habitats across local, regional, and global scales.Common applications for species modeling include forecasting changes in species distribution for current, past, and future climates, studying relationships between environmental parameters and species richness, mapping the range of invasive species, and conservation planning [5].
Remote sensing techniques such as species distribution modeling now allow scientists to determine the biotic and abiotic factors that drive landscape-level species distributions across the entirety of a species' range.As such, these techniques have greatly improved the ability of natural resource managers and conservationists to study and manage natural resources at the landscape level.For example, SDMs have been used in the context of conservation and management across many taxonomic groups and have proved useful in identifying and understanding species habitat and climatic requirements [6], hotspots and species co-occurrence patterns [7], patterns of endemism [8,9], threats to species persistence [10,11], predator-prey interactions [12], the spread of invasive species [13], and niche theory [14].
The complexity of the problems for which SDMs are employed to solve has increased as SDM capabilities have improved.Models that can handle data from multiple scales and dimensions are necessary to produce accurate and useful forecasts as anthropogenic pressures on ecosystems increase [15].Due to this increased demand, SDM approaches have been developed quite rapidly, combining new analytical frameworks and a variety of data sources.Recent developments have adapted SDMs to not only depict probable habitats, but also to infer the possible mechanisms underlying the presence or absence of specific species [16].Particularly for species that are sensitive to abrupt environmental changes, these contemporary models are better able to identify key areas of conservation concern or priority; this has substantially improved landscape-scale conservation action and precision.
Due to their narrow range of suitable environmental factors and relatively limited geographic extent, specialist species are ideal candidates for SDMs as their traits generally lead to a higher species model performance [17].Specialist species are also extremely susceptible to anthropogenic stressors such as climate change and land use change (e.g., agriculture, urbanization).These habitat alterations can fragment habitats, resulting in demographic isolation, population decline, or species extirpation [18,19].Habitat modeling provides a visual representation of the distribution of a species' fundamental niche and is thus often used as a key component in understanding environmental "hot spots"; it can help mitigate habitat fragmentation and allow resource managers to adequately plan for current and future climate scenarios.Therefore, SDMs for at-risk species and habitat specialists provide natural resource managers with a useful tool for conservation.
The Arizona toad (Anaxyrus microscaphus) occurs in the southwestern United States primarily near the Mogollon Plateau in western New Mexico, expanding through Arizona into far southwestern Utah and eastern Nevada along the Virgin and Colorado River basins and their tributaries [20,21].Historically, this species has occurred in the Agua Fria, Salt, Verde, Bill Williams, and Hassayampa Rivers in Arizona, and the Gila, Mimbres, and San Francisco Rivers in New Mexico [22,23].Anaxyrus microscaphus is a habitat specialist, utilizing sandy marginal zones or terraces and preferring a mixture of dense willow clumps and open flats or flood channels at elevations between 365 and 2700 m [24,25].Typically, adults are found within the vicinity of seasonal or permanent streams near slowly running shallow water, along permanent ponds, or within the proximity of larger streams in rocky canyons.Threats to the A. microscaphus habitat include urbanization, land alteration, agriculture, and impoundments [25].While this species has been categorized by the International Union for Conservation of Nature (IUCN) as a species of the least concern, A. microscaphus populations are largely characterized as declining [26], resulting in increased interest in the conservation of this species [23].Therefore, this highlights the suitability, as well as the need, for an SDM for this species and SDMs representing distributions under various future climate scenarios.The lack of these SDMs for A. microscaphus represents a major gap in the landscape-scale knowledge of this species' ecology and conservation.
One of the most unique and diverse landscapes in North America, Zion National Park (hereafter, Zion), provides refuge to many protected species within its boundaries, including A. microscaphus.It is thought that a suitable habitat for A. microscaphus currently remains within Zion, and sightings of this species are common along riparian zones and other seasonal water sources in the park.To protect sensitive habitats within the park, Zion complies with the National Environmental Policy Act (NEPA) in addition to other environmental regulations, particularly those involving the Endangered Species Act and the National Historic Preservation Act.However, park management often faces challenges regarding the decision-making processes of the NEPA with respect to the delineation of habitat ranges for protected and endangered species within the park.Analyzing species maps for potential alterations in the habitat of a threatened or endangered species is one of many steps in determining environmental impacts and the overall risk and vulnerability of the species.
Because A. microscaphus is characterized as a specialist species, it is imperative to provide a defensible procedure for developing species mapping within Zion and across its range to aid in decision-making processes.Therefore, our objectives were to (1) develop a species distribution model for the Arizona toad and delineate fine-scale species distribution maps within the park boundaries of Zion to facilitate conservation planning and decisionmaking in this critical area of protected habitat; to (2) develop an additional range-wide distribution model and forecast larger-scale changes in distribution under various predicted future climate scenarios; and to (3) test the importance of topographic and climatic variables in shaping A. microscaphus distributions at various spatial (i.e., within Zion National Park and range-wide) and temporal scales (i.e., current and future).

Study Area
Our study area included two scales of A. microscaphus distribution in North America.First, we studied A. microscaphus distribution on a localized landscape scale within Zion.Zion is located in the state of Utah in the western United States (Figure 1).Zion is located at the juncture of the Colorado Plateau, Mojave Desert, and Great Basin ecoregions.The elevation ranges from 2660 m at its highest point to 1117 m at its lowest point.Secondly, we studied A. microscaphus distribution on a range-wide scale across the full extent of its known and apparent range in North America.The ecoregions present at both scales included the Arizona-New Mexico Mountains, Great Basin, Mojave Desert, Colorado Plateau, and Utah High Plateaus ecoregions.The Arizona-New Mexico Mountains ecoregion is mostly semiarid with elevations ranging from 1371 to 3048 m and an average annual rainfall of 25-88 cm, which supports forests with perennial, intermittent, and ephemeral stream system microclimates.Vegetation below 1676 m consists of a surrounding desert habitat of juniper, yucca, stool, lechuguilla, prickly pear cactus, and mesquite.At elevations up to 2650 m, vegetation is dominated by semi-arid grasslands.The Mohave Desert ecoregion has a climate characterized by seasonal precipitation ranging from 8 to 25 cm and elevations consisting of majority desert plains and isolated mountains, ranging from −91 to 3344 m in elevation.The channel and stream characteristics consist mostly of seasonal runoff that filters through ephemeral drainage basins, with some of the runoff draining into the Colorado River [27].The Colorado Plateau ecoregion has an elevation ranging from 1500 to 2100 m and has average annual temperatures ranging from 4 to 13 • C. The climate is arid to semi-arid, with an average annual precipitation of roughly 14 to 40 cm.The ecoregion is characterized by geographic features including broad plateaus, ancient volcanoes, and deep canyons [28].In the Utah High Plateaus ecoregion, located to the far north of the study area, elevations range between 1500 to 4100 m.Precipitation is distributed throughout the year, ranging from 37.5 to 90 cm annually.An abundance of perennial streams drain into the Sevier, Virgin, and Colorado rivers [27].

Occurrence Data
Occurrence data for A. microscaphus were obtained using the Global Biodiversi formation Facility (GBIF; https://www.gbif.org/;accessed on 15 August 2020).To pr the spatial data for the model, the occurrence points were downloaded and cleaned the removal of duplicates and extreme outliers, and points located outside of the re able known range of the species).We filtered data to eliminate issues of spatial auto lation.Filtering was performed by removing localities that were within a 30 km rad each other [29].The spatial filtering step was performed using the SDMtoolbox wi tool 'Spatially Rarefy Occurrence Data'.A distance of 30 km was chosen based upo high spatial heterogeneity of the western mountains of the U.S. [30][31][32].

Data Acquisition 2.2.1. Occurrence Data
Occurrence data for A. microscaphus were obtained using the Global Biodiversity Information Facility (GBIF; https://www.gbif.org/;accessed on 15 August 2020).To prepare the spatial data for the model, the occurrence points were downloaded and cleaned (i.e., the removal of duplicates and extreme outliers, and points located outside of the reasonable known range of the species).We filtered data to eliminate issues of spatial autocorrelation.Filtering was performed by removing localities that were within a 30 km radius of each other [29].The spatial filtering step was performed using the SDMtoolbox with the tool 'Spatially Rarefy Occurrence Data'.A distance of 30 km was chosen based upon the high spatial heterogeneity of the western mountains of the U.S. [30][31][32].

Spatial Scale and Resolution
The study area extent for the focal species was determined by constructing a convex hull fitted to the spatially filtered occurrence points, followed by a 50 km buffer around the hull.The convex hull was generated in ArcMap 10.6.1 using the 'convex hull' tool.This technique is used to reduce the overfitting of a model [29,33].The hull was then used to clip the proceeding environmental variable rasters using the same extent, coordinate system, pixel count, and resolution.Each distribution (i.e., Zion-specific distribution and the complete distribution) was estimated based on data from their respective scales in order to improve the precision and accuracy of the models.
We used a 30 m grain size resolution for the environmental variable rasters.The grain size was determined primarily due to A. microscaphus being labeled as a habitat specialist and the need to delineate fine-scale habitat features (i.e., fragmentation, movement corridors) that may otherwise not be visible with a coarser grain resolution.The geomorphologic and remotely sensed covariates had a predetermined spatial resolution of a 30 m grain size and did not need further modifications to meet the grain size requirements.The climate data incorporated a two-step interpolation procedure based on the approach of Hutchinson and Xu [34] using weather station data (https://www.worldclim.org/data/worldclim21.html;accessed on 15 August 2020).A multiple linear regression model using Worldclim bioclimatic rasters and raster covariates followed by the thin plate spline (TPS) interpolation of the regression residuals was constructed in R [35] to create 30 m resolution climate rasters for the study extent of A. microscaphus [27,36,37].Following the methodology of Hutchinson and Xu [34], using the ANUSPLIN climate interpolation package, the R package, MACHISPLIN, was used as an alternate method for climate data interpolation.This resulted in all model covariates having a 30 m resolution and allowed the models to run with finer-scale resolutions.

Elevation and Remote Sensing Data
Digital elevation model (DEM) rasters were obtained from the NASA Earthdata website (https://earthdata.nasa.gov/;accessed on 15 August 2020) and mosaicked together using ESRI ArcMap 10.6.1 to form a master DEM of the study area.
Data acquisition for the remotely sensed environmental variables was performed using the USGS website (https://earthexplorer.usgs.gov/;accessed on 15 August 2020) and included atmospherically corrected level-two Landsat satellite imagery.The scenes collected used atmospherically corrected level-two Landsat satellite imagery and dated back to a time frame that correlated with the average rainfall and temperature, coinciding positively during the breeding seasons for A. microscaphus.By using imagery from an average rainfall and temperature year, the potential effects of substantial land cover change were minimized.The imagery for A. microscaphus was obtained using Landsat 7 imagery with <5% cloud cover with various dates from May 2002.The remote sensing covariates included the normalized difference vegetation index (NDVI), for characterizing various aspects of land cover, and the bare soil index (BSI), for characterizing areas of bare soil.

Environmental Variables
Other variables prepared for the modeling procedure were created using ArcMap 10.6.1.A slope raster was constructed using DEMs by calculating the maximum rate of change from the target cell and the eight surrounding neighbors in the raster.An aspect raster was also constructed using the DEM as the input data.We calculated a terrain ruggedness index (TRI) using the Vector Ruggedness Measure tool developed by Sappington et al. [38].The topographic position index (TPI) was calculated using the difference between a cell elevation value and the average elevation of the surrounding neighborhood of the cell.TPI variables were created using the Land Facet Corridor Analysis version 1.2.605 toolbox [39].The normalized difference vegetation index was calculated by obtaining relevant temporal Landsat imagery using the near-infrared and red color bands to generate an image displaying greenness.Lastly, the bare soil index (BSI) was created using the Raster Calculator tool [40].

Modeling Procedure
We modeled the potential distribution of A. microscaphus using a maximum entropy modeling approach (MaxEnt; [41]).We chose maximum entropy modeling over other available methods because it consistently results in better presence-only models [42,43].The use of variables in the constructed models using the a priori process was followed by the removal of highly correlated variables using Pearson's correlation coefficient.To adjust for multicollinearity, covariates displaying a high correlation above 0.9 were excluded in when preparing to build the model.We ran a total of 5 iterations for each model.

Digital Elevation Model
Elevation is an important variable in habitat selection in that A. microscaphus elevation rarely exceeds 2700 m [44].Elevation models were obtained from NASA's Earthdata website and clipped to the study area.

Slope
The slope and steepness of a region significantly influence runoff, especially in mountainous areas like Zion. A. microscaphus is commonly found in areas with little to no slope, although this species prefers breeding and egg deposition in lightly flowing water [44].

Aspect
Aspect was chosen as a variable display solar radiation for A. microscaphus because it prefers to live in valley bottoms or areas with high canyon walls surrounding streams and rivers [44,45].We hypothesize that this is important because this habitat could be influenced by solar exposure and a more southern-facing valley could constitute a preferred habitat.

Terrain Ruggedness Index
Because A. microscaphus prefers areas of relatively smoother topography, the terrain ruggedness index (TRI) should display a negative linear correlation within the MaxEnt models.

Topographic Position Index
The topographic position index was chosen for A. microscaphus to represent valley bottoms, which typically represent a high probability of habitat suitability due to the likelihood of alluvial accumulation and the presence of streams [44,45].

Normalized Difference Vegetation Index
Encroaching vegetation along the riparian habitat can potentially lead to the loss of a suitable habitat for A. microscaphus.We chose to utilize the NDVI because it quantifies vegetation along riparian areas that would negatively correlate with a suitable habitat choice for A. microscaphus.

Bare Soil Index
The bare soil index (BSI; [40]) was selected because A. microscaphus commonly occurs near areas displaying sandy soils or rocky slopes and sparse vegetation.The BSI relies on the short-wave infrared and red spectral bands to quantify soil mineral composition, while the blue and near-infrared bands display vegetative density.

Model Evaluation 2.5.1. Forecasting
Future model predictions were constructed using the bioclimatic variables from World-Clim (http://www.worldclim.org;accessed on 15 August 2020) for the years of 2050 and 2070.Yearly averages for 2050 were derived from the years (2041-2060), while 2070 projections were derived and averaged from (2061-2080).The Community Climate System Model (CCSM4) was used for forecasting.The spatial resolution for the future climate data was 30 s (~1 km 2 ).Best and worst-case future climate scenarios were used for 2050 and 2070, with representative concentration pathways (RCPs) of 2.6 for the best-case scenario and 8.5 for the worst-case scenario.MaxEnt was used to compare the differences between the current suitable habitat in comparison to the best and worst-case future scenarios.Forecasting within MaxEnt used the same variables the current model uses, excluding the NDVI and BSI.The NDVI and BSI were removed because of inaccuracies in the projected future measurements of vegetation in bare soil.

Model Assessment and Validation
The MaxEnt outputs for A. microscaphus were calculated statistically by using the area under the receiver operating characteristic curve (AUC) to quantify the strength of the model, as the AUC is one of the most common statistics used to characterize model performance [46].The AUC represents how well a model performs in predicting occurrences as opposed to a random selection of points.Model accuracy is considered excellent if the AUC is between 0.9 and 1, good if between 0.8 and 0.9, fair if between 0.7 and 0.8, poor if between 0.6 and 0.7, and failed if the AUC is between 0.5 and 0.6 [47].
The concept of the permutation importance of variables was used to assess the influence of predictor variables by measuring the drop in performance after randomly shuffling each variable's values.The use of this technique in SDMs helps to prioritize key environmental drivers, ensuring more informed conservation decisions and resource allocations.
Additionally, we used the bootstrap method and a random selection of 25% of the occurrence points as test data to cross-validate all models.Cross-validation is preferred with smaller sample sizes to avoid repeating the same selection of occurrence points for later iterations.

Results
Of the 327 occurrence points for A. microscaphus, 87 rarified occurrence localities were left for analysis and 17 (25%) were used for test data (Figure 2).The AUC for the A. microscaphus model was 0.853 for the training data and 0.810 for the test data.The variables selected for inclusion in the models included the annual mean temperature, mean diurnal range, isothermality, precipitation seasonality, precipitation of coldest quarter, aspect, bare soil index, elevation, normalized difference vegetation index, ruggedness, slope, and topographic position index.
Of the 327 occurrence points for A. microscaphus, 87 rarified occurrence localities were left for analysis and 17 (25%) were used for test data (Figure 2).The AUC for the A. microscaphus model was 0.853 for the training data and 0.810 for the test data.The variables selected for inclusion in the models included the annual mean temperature, mean diurnal range, isothermality, precipitation seasonality, precipitation of coldest quarter, aspect, bare soil index, elevation, normalized difference vegetation index, ruggedness, slope, and topographic position index.

Model Performance
The top five permutation importance variables for A. microscaphus were as follows: annual mean temperature (22.9%), elevation (22.7%), isothermality (11.5%),NDVI (11.1%), and mean diurnal range (11.0%) (Table 1).Response curves (Figure 3) display the range of suitable conditions for each environmental variable independent of all other variables.The AUC for A. microscaphus was recorded as 0.9 for the training data and 0.8 for the test data, indicating good model accuracy.

Model Performance
The top five permutation importance variables for A. microscaphus were as follows: annual mean temperature (22.9%), elevation (22.7%), isothermality (11.5%),NDVI (11.1%), and mean diurnal range (11.0%) (Table 1).Response curves (Figure 3) display the range of suitable conditions for each environmental variable independent of all other variables.The AUC for A. microscaphus was recorded as 0.9 for the training data and 0.8 for the test data, indicating good model accuracy.

Future Climate Trends
The current suitable habitat (>60% suitability) for A. microscaphus within Zion consists of 99,370 km 2 , which is 42.63% of the 233,078 km 2 (Figure 4).When interpolated into the future across the entire range of the species, the 2050 RCP of 2.6 W/m 2 habitat suitability decreased to 15,125 km 2 or 6.48% (Figure 5).Further, the RCP of 8.5 W/m 2 for the year 2050 increased to 7.59% suitability, with a surface area of 17,702 (km 2 ), an increase from the previous RCP of 17.04% (Figure 5).The 2070 RCP of 2.6 W/m 2 displayed an area of 17,942 km 2 , representing 7.7% suitability for the full range of A. microscaphus in North America (Figure 6).The suitable habitat area decreased considerably for the 2070 RCP of 8.5 W/m 2 down to an area of 11,976 km 2 , representing a suitable area of 5.14% and a 33.25% decrease from the 2070 RCP of 2.6 W/m 2 (Figure 6).
The current suitable habitat (>60% suitability) for A. microscaphus within Zion consists of 99,370 km 2 , which is 42.63% of the 233,078 km 2 (Figure 4).When interpolated into the future across the entire range of the species, the 2050 RCP of 2.6 W/m 2 habitat suitability decreased to 15,125 km 2 or 6.48% (Figure 5).Further, the RCP of 8.5 W/m 2 for the year 2050 increased to 7.59% suitability, with a surface area of 17,702 (km 2 ), an increase from the previous RCP of 17.04% (Figure 5).The 2070 RCP of 2.6 W/m 2 displayed an area of 17,942 km 2 , representing 7.7% suitability for the full range of A. microscaphus in North America (Figure 6).The suitable habitat area decreased considerably for the 2070 RCP of 8.5 W/m 2 down to an area of 11,976 km 2 , representing a suitable area of 5.14% and a 33.25% decrease from the 2070 RCP of 2.6 W/m 2 (Figure 6).

Discussion
Our models demonstrate the utility of remote sensing techniques in elucidating ecological relationships and improving our understanding of species ecology at the landscape scale.For example, the use of SDMs for the mapping of suitable habitats is a robust tool for spatially modeling species ecology in relation to geomorphological, remotely sensed, and climatic data to better understand landscape patterns that can then be used to frame management strategies or monitoring programs.Prior studies have undertaken habitat suitability analyses and the development of SDMs, demonstrating various methodologies and results [48,49].However, these tools have not been widely applied to A. microscaphus, especially in the context of climate change.Therefore, our results aid in understanding this species' habitat requirements at the landscape scale using remotely sensed data and a SDM framework, and how these may change under future global change scenarios.For example, our MaxEnt models generated a comprehensive map of the potential distribution of A. microscaphus, including Zion National Park, and reveal the unique contribution of the environmental variables to the prediction of habitat suitability both now and into the future.

Discussion
Our models demonstrate the utility of remote sensing techniques in elucidating ecological relationships and improving our understanding of species ecology at the landscape scale.For example, the use of SDMs for the mapping of suitable habitats is a robust tool for spatially modeling species ecology in relation to geomorphological, remotely sensed, and climatic data to better understand landscape patterns that can then be used to frame management strategies or monitoring programs.Prior studies have undertaken habitat suitability analyses and the development of SDMs, demonstrating various methodologies and results [48,49].However, these tools have not been widely applied to A. microscaphus, especially in the context of climate change.Therefore, our results aid in understanding this species' habitat requirements at the landscape scale using remotely sensed data and a SDM framework, and how these may change under future global change scenarios.For example, our MaxEnt models generated a comprehensive map of the potential distribution of A. microscaphus, including Zion National Park, and reveal the unique contribution of the environmental variables to the prediction of habitat suitability both now and into the future.
Our results predict a substantial range-wide reduction in habitat suitability across climate scenarios from 2020 to 2050 and from 2050 to 2070, which has significant conser-vation implications for A. microscaphus due to its current at-risk status.Our results are significant and help to fill a much-needed gap in the understanding of A. microscaphus population declines.Recent declines have been primarily attributed to factors other than habitat and climate, such as the introduction of invasive species and hybridization with closely related bufonid toads [44].Although these factors most definitely have had a significant impact on A. microscaphus populations and contemporary distributions, our results highlight the importance of habitat and climatic characteristics in shaping contemporary and future distributions of this species.Our results coincide with those of Suzart de Albuquerque et al. [45], which described the climatic factors (e.g., climate, topography, solar radiation) that significantly shape current and future (2081-2100) distributions of A. microscaphus.In their study, they also concluded that climate change may reduce habitat suitability and the distribution of A. microscaphus.Our model of the contemporary distribution of A. microscaphus includes not only the climatic and topographic variables proposed by Suzart de Albuquerque et al. [42], which were hypothesized to be the main drivers of the distribution of this species, but also includes the DEM (an absolute elevation value) and the NDVI (habitat vigor and health).The inclusion of the DEM and the NDVI improves the comprehensiveness of the representation of the absolute elevation and vegetation conditions, enhancing the predictive accuracy and ecological relevance of the model overall.Additionally, our forecasting models (i.e., predicted distributions in 2050, 2070) fill a temporal gap in predictions of distribution between now and 2100.The coupling of these models provides a time series of predictions of distributional changes in A. microscaphus over the next 80 years.
Zion is rich with microclimates that are affected by sharp landscape changes, including steep canyons, narrow and wide valleys, and steep slopes with differing aspects, which in addition to climatic variables, appear to influence the current and future distributions of A. microscaphus.Broadly, the identified geographic and climatic variables affect this microhabitat in ways that result in the increased or decreased identification of suitable habitats using fine-scale resolution models.Temperature contributed the most to model projections, with habitat suitability generally decreasing with an increase in the annual mean temperature above 10 • C. Because A. microscaphus is an ectotherm and relies heavily on the environment for thermoregulation, it is not surprising that temperature was influential in shaping this species' distribution.Additionally, precipitation (i.e., precipitation seasonality, precipitation of coldest quarter), vegetation (e.g., the normalized difference vegetation index), and topography-related variables also contributed greatly to the models.Variables related to precipitation are also not surprisingly important to A. microscaphus distribution because water and moisture are often the primary limiting factors for amphibians, especially in arid environments such as southwestern North America.Amphibians are commonly known to rely heavily on moisture and temperature for proper physiological function.Additionally, vegetation and topography are likely important because of their effects on the available heat and moisture, which this species requires for continued persistence.However, vegetation likely plays a prominent role with regard to providing proper thermal and escape cover.It is likely that these variables represent important ecological aspects that affect natural and life history characteristics, which in turn influence species distribution.For example, the smaller-scale model (i.e., only within Zion) demonstrates that the most suitable habitat for A. microscaphus occurs along sloped rivers and streams and in valleys that typically host sandy soils, which appear to be essential for this species' survival and persistence.This is further supported by the positive relationship between habitat suitability and the NDVI and its negative relationship with the Topographic Position Index (TPI).Specifically, a higher NDVI value indicates more green vegetation, which is indicative of higher water availability, more cover, and more preferable microclimates, which are characteristic of riparian habitats in this region.Additionally, a negative relationship with the TPI is indicative of the preference of A. microscaphus for valleys or lower-lying riparian areas.These are reasonable results considering that A. microscaphus is known to be strongly associated with riparian habitats throughout its range in the southwestern United States and that many of the individual variable relationships highlight important habitat characteristics that are common in this habitat type.Past studies have revealed similar ecological relationships.For example, Montgomery [44] found that broadscale habitat variables such as bioclimatic variables, the annual temperature range (lower probability with higher extreme temperatures), and precipitation (lower probability with decreasing rainfall) contributed to the probability of a habitat being occupied by A. microscaphus.Additionally, Montgomery [44] found that fine-scale variables such as shallow water, less woody debris, pebble substrate, and canopy cover contributed positively to A. microscaphus presence.This demonstrates the importance of riparian habitat conservation, which, based on our results, would likely benefit A. microscaphus populations greatly.As a result of these models and associated maps, land managers can visualize the probability of a suitable habitat for A. microscaphus and strategically allocate resources or construction projects in ways that benefit this species with more precise and accurate approaches.Moreover, the development of precise and accurate survey protocols, particularly those based on remotely sensed data, are especially important for A. microscaphus because of its status as an at-risk species and the limited knowledge of its habitat use at localized and landscape scales [50].
Based on the current climate projections, the A. microscaphus habitat will likely become increasingly more constrained across its range.For example, A. microscaphus habitat suitability decreased significantly as the models were projected further into the future with worst-case RCPs.A noticeable shift in habitat towards major river systems and away from smaller river systems can be observed in the 2070 forecasting maps, a concerning trend for this species considering its apparent association with riparian and riverine habitats.This may indicate a change in the amount of suitable riparian habitat available in the future as the climate warms and drought conditions increase in prevalence and duration.For instance, if smaller, more ephemeral streams in the landscape decrease as a result of rising temperatures and more severe and prolonged drought conditions, it is reasonable to expect populations of A. microscaphus to rely more heavily upon and persist near larger permanent river systems.This highlights a conservation and management target that can aid in bolstering or maintaining A. microscaphus populations in the North American landscape as conditions change.For example, we recommend that future management strategies should not only maintain larger, more persistent riparian habitats, but should also focus efforts on conserving smaller, more ephemeral habitats to increase their resilience during changing climatic conditions.This demonstrates the effectiveness of SDMs in deepening the understanding of species distributions across large spatial scales by elucidating the drivers that influence contemporary and future distributions, providing much-needed, and previously unattainable, information for the refinement of decision-making and management to combat the landscape-level effects of climate change.
Additionally, our models of the potential distribution of A. microscaphus can assist biologists and scientists in not only monitoring known populations, but also discovering populations by focusing surveys on areas of increased habitat suitability [2][3][4]51].For species of concern, such as A. microscaphus, priority conservation areas can be easily and cost-effectively identified by MaxEnt models [2][3][4][52][53][54].These models can also be used to guide population survey efforts for this species within Zion and across its range in the southwestern United States.For example, it has been found that as much as 70% of the time required to plan surveys for new species occurrences is saved when using SDMs to guide survey efforts [55].This is an important consideration when aiming to maximize the utility of conservation resources that may need to conduct surveys on and monitor multiple species [51].
Due to many complex ecological processes, areas predicted to be suitable may not be inhabited by A. microscaphus; therefore, our models are better described as representations of the species' fundamental, rather than realized, niche.Similar to many other distribution modeling studies, the ecological niche-based distribution models do have associated uncertainties and recognized weaknesses; therefore, our maps should only be used as estimations of the potential species distribution rather than definitive representations of the true distribution of this species [47].A major limiting factor to conservation and management activities is often a lack of information about the spatial distribution of target species and how this may change over time [2][3][4]50].As such, even theoretical habitat suitability and distribution models are highly useful for applied conservation, as they can outline distributions to a degree of certainty while also improving our understanding of the factors that shape species distributions, which allows for more thoughtful ecological insights [50].
Our results indicate that the current distribution of A. microscaphus in Zion National Park is a result of several different ecological and geological factors, but is also at risk of significant reduction across its range if current climate trends continue into the future.The future analysis and ground truthing of A. microscaphus distribution via field surveys is necessary to expand the model extent to the conterminous United States for a larger-scale analysis of the effects of climate and habitat change on suitable habitat.Also, expanding studies on local scales will help land managers and conservationists better understand how global and regional climate patterns influence local species populations for future needs.A noticeable lack of public literature regarding conservation efforts coupled with species models exemplifies the need to utilize SDMs based on their explicit conservation purposes.A standard protocol could objectively be assigned to private and government agencies to explore a middle ground to meet a standardized SDM approach for species on protected and private land alike.

Figure 1 .
Figure 1.Location of Zion National Park in southwestern Utah, U.S.

Figure 1 .
Figure 1.Location of Zion National Park in southwestern Utah, U.S.

Figure 2 .
Figure 2. The study extent and range of occurrence of localities for Anaxyrus microscaphus across the full extent of the known species range.

Figure 2 .
Figure 2. The study extent and range of occurrence of localities for Anaxyrus microscaphus across the full extent of the known species range.

Figure 3 .
Figure 3. Partial dependence plots displaying the marginal response of the 12 environmental variables selected for Anaxyrus microscaphus in the MaxEnt model.Each response curve demonstrates the range of suitability for each environmental variable if each variable were used to create a MaxEnt model independent of the outside influence of other variables.

Figure 4 .
Figure 4. MaxEnt model output for 30 m resolution habitat suitability maps for Anaxyrus microscaphus in Zion National Park, Utah, USA.Suitability classes describe the ranking of habitat preference for the species being modeled.Models were created using MaxEnt.Subfigure indicates the location (red square) of the suitability map relative to the state of Utah, USA.

Figure 4 .
Figure 4. MaxEnt model output for 30 m resolution habitat suitability maps for Anaxyrus microscaphus in Zion National Park, Utah, USA.Suitability classes describe the ranking of habitat preference for the species being modeled.Models were created using MaxEnt.Subfigure indicates the location (red square) of the suitability map relative to the state of Utah, USA.

Figure 5 .
Figure 5. Forecasting SDM of Anaxyrus microscaphus contrasting the suitable habitat for future best and worst-case climate scenarios with the representative concentration pathways that describe greenhouse gas concentrations of 2.6 W/m 2 and 8.5 W/m 2 for the year 2050.Subfigure indicates the location (red square) of the suitability map relative to the state of Utah, USA.

Figure 5 .
Figure 5. Forecasting SDM of Anaxyrus microscaphus contrasting the suitable habitat for future best and worst-case climate scenarios with the representative concentration pathways that describe greenhouse gas concentrations of 2.6 W/m 2 and 8.5 W/m 2 for the year 2050.Subfigure indicates the location (red square) of the suitability map relative to the state of Utah, USA.

Figure 6 .
Figure 6.Forecasting SDM of Anaxyrus microscaphus contrasting the suitable habitat for future best and worst-case climate scenarios with representative concentration pathways that describe greenhouse gas concentrations of 2.6 W/m 2 and 8.5 W/m 2 for the year 2070.Subfigure indicates the location (red square) of the suitability map relative to the state of Utah, USA.

Figure 6 .
Figure 6.Forecasting SDM of Anaxyrus microscaphus contrasting the suitable habitat for future best and worst-case climate scenarios with representative concentration pathways that describe greenhouse gas concentrations of 2.6 W/m 2 and 8.5 W/m 2 for the year 2070.Subfigure indicates the location (red square) of the suitability map relative to the state of Utah, USA.

Table 1 .
Permutation importance values for each bioclimatic variable within the MaxEnt model for the Arizona toad.The permutation value is determined by randomly permuting the values of each independent variables against the training points.Values are then normalized to provide percentages; higher values suggest a greater influence on the model.