Amphibian populations are threatened globally, and one of the hypotheses for these declines is climate change. Species distribution models are frequently used to predict changes in suitable habitat as a result of changing climates; however, these projections can be heavily influenced by choice of modeling approach. To evaluate global predictions for amphibians, we conducted a literature review of studies that utilize correlative species distribution models to project changes in climate suitability under various climate change scenarios. We paid particular attention to the use of model selection in choosing among candidate species distribution models (SDMs) so as to control for overparameterization of SDMs. In addition, we conducted a case study with three species of Slender Salamanders (Batrachoseps) to further investigate the impact of differences in modeling decisions on projected amphibian climate suitability. We found 83 studies (including the present case study) in which projections of future climate suitability were made for amphibian species. Of those studies, 36 included estimates of percent change in climate suitability and thus were included in our meta-analysis. These studies included projections for over 1000 species or species complexes, with the majority being Anurans (86%), and encompassed five continents with the most representation in South America, Europe, and North America. Across these studies, average projected change in climate suitability ranged from －70% to 167%, and these projected changes varied with dispersal assumptions representative concentration pathway (RCP) used, the projection year, and taxonomic order. Only three of the 36 studies reported the use of Akaike information criterion (AIC)-based model selection to choose a best-fit SDM. However, our case study demonstrated that predicted change in climate suitability varied whether the best-fit or default SDM was used for projections. Further, this result varied among species (ΔAIC = 0), suggesting that the impact of overparameterization differs across species. Our results illustrate that there is a pressing need to project climate suitability across more species and more geographic regions. In addition, we may need to revisit projections for previously investigated species to evaluate additional climate scenarios and whether overparameterization may have influenced projections. Our ability to accurately model future changes in climate suitability will be essential for successful conservation and management plans for amphibians.
The amphibian decline crisis is a complex issue (Houlahan et al. 2000) caused by a wide range of factors (Collins and Storfer 2003). Among these complex factors, climate change has been proposed as a major contributor to current (Smalling et al. 2019) and future declines in species (Carey and Alexander 2003). Amphibians may be particularly sensitive to the impacts of climate change relative to other terrestrial vertebrates because of their biological dependency on temperature and moisture for reproduction and survival as well as their limited dispersal capabilities (Donnelly and Crump 1998). To estimate the potential impact of climate change on amphibian declines, species distribution modeling (also known as ecological niche modeling; Elith and Leathwick 2009) has been used to project amphibian habitat suitability under future climate scenarios (Teixeira and Arntzen 2002; Parra-Olea et al. 2005). Numerous studies have predicted extensive losses of suitable habitat for a wide range of amphibian species under various climate projections (Milanovich et al. 2010; Barrett et al. 2014; Zank et al. 2014). However, a comprehensive evaluation of these studies has yet to be done and would provide a more global understanding of the regions of the world and species that may be most at risk to losses in suitable habitat. Furthermore, since the first use of species distribution modeling for projecting impacts of climate change on habitat suitability, the modeling process has advanced, with recent studies applying new recommendations (e.g., Warren and Seifert 2011; Phillips et al. 2017). These refinements to modeling approaches likely have significant impacts on the amount of suitable habitat projected for some species. Thus, it is crucial that we investigate how such modeling approaches influence projections to refine our predictions of the extent to which suitable habitat will be lost due to climate change.
Species distribution models (SDMs), which utilize georeferenced environmental data to predict the distributions of species, were first used to investigate environmental predictors of species occurrences and model habitat suitability (Elith and Leathwick 2009), but this approach was quickly adapted for conservation purposes (Rodríguez et al. 2007). To this end, SDMs can inform our understanding of baseline species distributions (Blank and Blaustein 2012), identify sites for habitat restoration or conservation (Giovannini et al. 2014), predict the spread of invasive species (Ficetola et al. 2010) or emergent pathogens (Yap et al. 2015; Richgels et al. 2016; Katz and Zellmer 2018), and predict future changes in habitat suitability as a result of global climate change (Parra-Olea et al. 2005). The development and growth of species distribution modeling in ecology (Elith and Leathwick 2009) occurred concurrently with increased attention on amphibian declines (Wake 1991); thus, it is no surprise that some of the earliest applications of this approach were for studying the effects of climate change on threatened or endangered amphibians (Teixeira and Arntzen 2002). Similar to other taxa (Peterson et al. 2002), distributional shifts are anticipated in response to climate change, with distributions predicted to shrink and shift into higher elevations and poleward (Dobson et al. 1989; Johnston and Schmitz 1997). As expected, some of the first projections of future climate suitability for amphibians predicted extensive losses for Golden-striped Salamanders, Chioglossa lusitanica (Teixeira and Arntzen 2002), and two Plethodon species (Parra-Olea et al. 2005) under climate change scenarios.
Yet, while numerous studies have projected losses in suitable habitat for amphibians, the extent of those losses remains unclear. Many choices in the modeling process affect the extent of projected suitable habitat (Thuiller 2004; Pearson et al. 2006; Buisson et al. 2010), such as which climate simulation (e.g., climate model intercomparison project [CMIP]), climate baseline, general circulation model (GCM), and representative concentration pathway (RCP) are used and which assumptions are made about dispersal (Real et al. 2010; Roubicek et al. 2010; Baker et al. 2016; Wright et al. 2016; Jarnevich and Young 2019; Thuiller et al. 2019). Even more primary, the modeling choices made in generating current-day SDMs can also influence projected changes in habitat suitability. The method or algorithm used to create SDMs (Beaumont et al. 2016; Sales et al. 2017; Thuiller et al. 2019), the thresholds chosen for habitat suitability (Liu et al. 2016), and the settings within different modeling frameworks, such as model tuning parameters (Warren and Seifert 2011) and identification of background localities (Rodda et al. 2011), can all lead to variations in current estimated habitat suitability which can then translate to differences in projected changes in habitat suitability under various climate change scenarios. As a result of the plethora of modeling decisions that influence projected habitat suitability, there remains much uncertainty about projections (Thuiller et al. 2019), making it difficult to make comparisons across different studies and, more importantly, affecting our ability to assess conservation risk (Wright et al. 2015).
In addition, the ease of use of some SDM tools has led to their use as black box tools, with little attention paid toward how different settings might influence the models created. Often, only the default settings are used in constructing SDMs (Merow et al. 2013; Morales et al. 2017; Katz and Zellmer 2018). Even when multiple approaches or settings are tested, the method by which the best-fit model is chosen among different candidate models can also lead to quantifiable differences in the amount of predicted suitable habitat (Warren and Seifert 2011). Specifically, it is important to control for overparameterization when comparing alternative candidate models (Warren and Seifert 2011). While differences between alternative models might only be slight when mapping current distributions, these differences may be amplified when spatially or temporally transferring the model (Warren et al. 2014). For instance, when predicting future habitat suitability, overparameterized models may result in less-suitable habitat because the over-fit model may make it harder to find analog environmental conditions across space or time. In fact, simulations suggest that overparameterized models perform worse in predicting future habitat suitability than do models that are less parameterized (Warren and Seifert 2011). More-complex parameterizations also led to greater predicted losses of suitable habitat for 34 European tree species (Brun et al. 2020) and 90 California species of concern (Warren et al. 2014). How the use of model selection to prevent overparameterization affects projected changes in future habitat suitability remains undertested relative to other modeling choices and yet may be central for interpreting differences in projections across studies.
Beyond potential effects of various modeling choices on projected outcomes, it is also important to consider the different classes of models and what each can tell us about future habitat suitability in the first place. Species distribution models can be constructed with correlative or mechanistic approaches, or a combination of the two. While these approaches work well under stationary conditions (Elith and Leathwick 2009), there are many scenarios in which conditions may change, such as with nonanalog climates (Dormann 2007; Fitzpatrick and Hargrove 2009), unmodeled gradients influencing species distributions, or species adaptations to climate change. Some have argued that alternative approaches to predicting habitat suitability, such as using mechanistic models which explicitly incorporate distribution-limiting processes, may provide better predictive power when dealing with changing climatic conditions (Kearney et al. 2010). However, mechanistic models have their limitations as well, requiring more-detailed ecological information (Dormann et al. 2012). Regardless, direct comparisons between correlative and mechanistic approaches have found congruent predictions (Kearney et al. 2010). Further, assessment of predictions made by correlative models using historical data indicates that this approach can generate reliable habitat suitability predictions under climate change (Morán-Ordóñez et al. 2017). While, ideally, multiple approaches should be used to validate predictions (Kearney et al. 2010; Diniz-Filho et al. 2019), fewer studies have utilized the mechanistic approach (Dormann et al. 2012). We thus focus in this study on evaluating correlative approaches for projecting habitat suitability under climate change.
Given the potential impact that each of these modeling choices could have on projecting changes in future habitat suitability of amphibians, it is timely to review the predictions that have been made using correlative models for future amphibian habitat suitability. Thus, we conducted a global assessment of predicted changes in climate suitability for amphibians and evaluated whether there are identifiable differences among studies that utilize alternative approaches to modeling species distributions. In addition, we investigated how model settings influence projected changes in suitable habitat for amphibians using one of the most widely used SDM tools, Maxent, which uses the machine learning approach maximum entropy (Phillips et al. 2017). While there are multiple settings that can be controlled within Maxent, we pay particular attention to the regularization multiplier (RM) which is used to control the extent to which additional parameters are added to the model (Phillips et al. 2006). To this end, we first conducted a literature review of studies that include SDM projections for amphibian species under global climate change scenarios. We then conducted a case study with three amphibian species to investigate the potential impacts of overparameterization on predicting habitat suitability declines. We created species distribution models for three species of Slender Salamanders (Batrachoseps) using the default Maxent settings and the best-fit alternative settings for 26 different future climate scenarios, consisting of 13 GCMs and two extreme RCPs (2.6, 8.5). We then calculated the proportion of habitat suitability predicted for each combination of GCM and RCP for both the default and best-fit models for all three species to quantify the difference in amount of predicted habitat. Finally, we discuss the outlook for the use of SDMs in studying amphibian declines as well as future considerations and research needs.
Materials and Methods
We conducted a Web of Science search for all articles on amphibians that predicted future distributions of species under climate change scenarios (September 2019). We used the search terms: TS = ([“niche model*” OR “distribution model*” OR “Maxent” OR “envelope model*”] AND [amphib* OR caudat* OR anura* OR caecilia* OR frog* OR newt* OR salamand* OR gymnophion* OR toad*] AND [“climate change” OR “climate scenario*” OR “ensemble forecast*”]). From these articles, we first identified the studies that used species distribution modeling to project changes in amphibian climate suitability. Here and throughout the paper, we define habitat suitability only in terms of climatic variables (e.g., temperature, precipitation). While many other variables, both biotic and abiotic, restrict the distributions of species, one of the primary concerns with climate change is availability of climatically suitable habitat. Further, bioclimatic variables are most frequently used in future projections due to their wide availability relative to other spatially explicit abiotic or biotic variable data layers. We thus restrict our analyses to bioclimatic variables only, which represent the realized climate space of organisms (Rodda et al. 2011).
For each of these articles, we recorded the species, genus, order, geographic range of the study, current climate models used, future climate models used (GCM), emissions scenarios considered (including either RCPs or the scenarios from the Special Report on Emissions Scenarios (SRES)), SDM approach used, threshold used for binary classification of suitable habitat, whether dispersal limitation was used (no dispersal, limited, unlimited), projection year, model selection approach used to compare candidate current SDMs, the RM if Maxent was used, and finally the predicted percent change in suitable habitat. Because some studies reported results for individual GCMs and others reported results for averages across GCMs, we used the average predicted change in suitable habitat for all GCMs considered for a single species within each study. Similarly, because some studies used RCP scenarios and others used SRES, we classified all RCP and SRES values into three general categories—pessimistic, moderate, and optimistic—based on similar temperature anomaly predictions (see Table S1 (Supplemenal_Material_S1.xlsx) in the Supplemental Materials available online). For studies that reported only average values across multiple RCP or multiple SRES, we used the averaged RCP value. We grouped projection year into 20-yr categories: 2030, 2050, 2070, and 2090. For studies that used Maxent, we categorized studies into those that used the default RM = 1 (default), those that used an increased RM > 1 without model selection (nondefault), and those that used Akaike information criterion (AIC) to select the best-fit RM (best-fit). We quantified summary statistics for the dataset to evaluate the frequency of all variables across studies.
To investigate which variables contributed to predicted changes in suitable climate, we conducted a meta-analysis via linear models using the metafor R package (Viechtbauer 2010). We first calculated effect sizes for each study by calculating the raw mean percentage change of future versus current suitable climate and standard deviation, with each species as a replicate. Thus, we only included studies with three or more species. We included as fixed variables: dispersal limitation, SDM approach, projection year, RCP, and taxonomic order. We were not able to include an RM category within the model because there were too few studies that utilized a best-fit RM value. Despite the known effect of both threshold approach and GCM on projections (Buisson et al. 2010; Liu et al. 2016), these variables were not included in our literature review because there were many variations of each, with insufficient repetition across studies to fully evaluate their effects. Similarly, we were unable to include which climate baselines were used (Roubicek et al. 2010; Baker et al. 2016; Jarnevich and Young 2019) because the climate baseline is relevant to the time period of observations used, which is infrequently reported in manuscripts. We compared 15 candidate models, including models with all combinations of the predictor variables and a null model. Because our sample size was limited, we did not consider any interactions among predictor variables. The candidate models were compared using AIC. The model with the lowest AIC score was considered the best-fit model. Models with ΔAIC < 2 were considered to have equal support whereas models with ΔAIC ≥ 2 were not considered to have support.
For our case study, we looked at projected changes in habitat suitability using default or best-fit SDMs for three species of Slender Salamanders. These species provided an interesting opportunity for evaluating the differential impacts of climate change, as they are similar in many biological traits and their ranges all met and slightly overlapped within a narrow area in southern California. However, despite their biological similarities and adjacent ranges, they inhabited distinct environments, with the southern-most species, B. major, primarily occupying sage scrub and oak woodlands in lowland coastal interior regions; the northern-most species, B. nigriventris, primarily occupying oak woodlands in moist mountain and foothills canyons near the coast; and B. gabrieli, with a very narrow distribution occupying talus slopes at the highest elevations in-between the other two species. While overlapping distributions could be problematic for modeling if interbreeding occurs, there are no known cases of hybridization among these species (Jockusch and Wake 2002; Jockusch et al. 2015). In addition to their distinct environmental needs, these species represent a spectrum of conservation needs, with B. major and B. nigriventris being of least conservation concern and B. gabrieli being data deficient with the potential for becoming a species of special concern (Thomson et al. 2016; Hansen 2017). Their unique but overlapping ranges, biology, and conservation status provided replicate systems for studying how modeling approaches influence suitable habitat projections.
Species locality data consisted of museum records (see Table S2 (Supplemenal_Material_S1.xlsx) in the Supplemental Materials available online) downloaded from Global Biodiversity Information Facility (GBIF; available at http://www.GBIF.org) using the rgbif R package (Chamberlain et al. 2019). Because B. gabrieli observations are much more rare, we supplemented the dataset with observations from the RASCals project on iNaturalist (Reptiles & Amphibians of Southern California, available at http://www.inaturalist.org/projects/rascals), N.A. HERP (North American Herpetological Education and Research Project, available at http://www.naherp.com), and our own observations. Only records with <1 km coordinate uncertainty were retained for the analysis. Localities were visualized to confirm that all observations are within the known range of each species. The majority of observations were recorded between 1950–2000. While some of the B. gabrieli observations were more recent, they were within the area of the species known range (Wake 1996). For each species, we spatially thinned the localities, removing any points within 5 km of one another for B. major and B. nigriventris and 1 km for B. gabrieli (due to limited range and observations). After thinning the dataset, there were 82 B. major, 36 B. nigriventris, and 23 B. gabrieli localities.
Environmental data were downloaded from WorldClim v1.4 (available at http://www.worldclim.org) at a resolution of 30 arc seconds and included 19 bioclimatic variables averaged from 1960–1990 (Hijmans et al. 2005). This version was selected because it best overlaps the years of most observations of each species from the museum record used in this study. We tested for correlations among the 19 bioclimatic variables at the observed localities for each species and removed highly correlated variables (r > |0.7|; Dormann et al. 2013; Petitpierre et al. 2017). For each species, we constructed 20 different candidate SDMs varying the RM from 1–20 at an interval of 1. The SDMs were created using Maxent v3.4 (Phillips et al. 2017) through the dismo R package (Hijmans et al. 2017) and implemented with the ENMeval package (Muscarella et al. 2014) to compare candidate SDMs. Maxent has been shown to perform well even with small sample sizes (n > 30), but for models built with fewer than 30 localities caution should be taken in interpreting the models (Wisz et al. 2008). We used the block partitioning method to create training and test datasets to build and validate each model. This method splits the dataset into four spatial blocks delineated by the latitude and longitude lines that divide localities and background locations as evenly as possible (Muscarella et al. 2014). For each candidate model, the training dataset consisted of three spatial blocks with the fourth remaining block reserved for test data, and this process was replicated across each block. To prevent overparameterization, we selected the best-fit model using corrected AIC (AICc; Warren and Seifert 2011). For each species, we compared the best-fit model (ΔAICc = 0) to the default Maxent model (RM = 1). We also assessed the area under the receiver operating characteristic curve (AUC) for both the training and test data, a method frequently used to evaluate model fit that does not account for number of model parameters. Comparing the difference in the AUC values between the training and test data can indicate how well the model performs for an independent set of data that were not used to build the model. However, only AICc was used to select the best-fit model for each species.
For both the best-fit and default models, we projected the SDM onto 13 different future climate (see Table S3 (Supplemenal_Material_S1.xlsx) in the Supplemental Materials available online) scenarios for the year 2050 (average of 2041–2060) for two different RCPs (2.6, 8.5) for a total of 26 future projections for each species for both models. Future climate projections were downloaded from WorldClim ( http://www.worldclim.org) and are downscaled CMIP5. For all projections we utilized the cloglog transformation, which is the new default output from Maxent v3.4 (Phillips et al. 2017). To calculate binary habitat suitability, we used the maximum training sensitivity plus specificity cloglog threshold (maxSS), which performs best for selecting a threshold with presence-only SDMs (Liu et al. 2016). We considered three different scenarios for dispersal: unlimited, limited, and no dispersal. For unlimited dispersal, habitat suitability was calculated across the entire extent of the study area; thus any cell could potentially become suitable in the future regardless of how distant it was from currently occupied cells. For dispersal limitation, we assumed a maximum total dispersal of 10 km over the projection period, as dispersal is thought to be highly limited for Plethodontid salamanders (Smith and Green 2005). For models with dispersal limitation, habitat suitability was calculated only within a 10-km buffer from the current predicted distribution for each species. For no dispersal, habitat suitability was calculated only within the extent of the current predicted distribution for each species. The percent change in habitat suitability was calculated as 100*([future suitability － current suitability] / current suitability; see Tables S1 (Supplemenal_Material_S1.xlsx) and S4 (Supplemenal_Material_S1.xlsx) in the Supplemental Materials available online).
To assess whether the percent change in habitat suitability was explained by whether the models were the best-fit versus default Maxent models, we created three candidate linear mixed-models with the lme4 R package (Bates et al. 2015). Percent change in habitat suitability was the response variable, and for the predictor variables we included the SDM type (best-fit vs default), species, and the interaction between SDM type and species. We compared this to two other candidate models: a model with only SDM type and species and a model with only species. Dispersal limitation, RCP (2.6 versus 8.5), and GCM were included as random variables in all three models. Model fit was evaluated with AICc using the R package AICcmodavg (Mazerolle 2019).
The literature search returned 201 citations, 82 of which included studies relevant to projecting future distributions of amphibian populations within their native ranges. Of these citations, 35 provided estimates of percent change in suitable habitat for at least one amphibian species. Combined with the case study we present below, there were 36 studies in our meta-analysis. The studies varied widely both in species and regions of the world represented. These studies involved 1092 species (or species complexes) with an average of 53 per study, covering 177 genera and representing 85.6% Anura, 13.6% Caudata, and 0.8% Gymnophiona (Table S1). Some studies investigated other taxa as well (e.g., mammals, reptiles; Ureta et al. 2018), although we do not report those results here. Collectively, the studies produced predicted distributions in 85 countries across 5 continents including Europe, Asia, North America, South America, and Africa (Fig. 1). Brazil was represented in the most studies (n = 9) followed by Spain, the United States, and Uruguay (n = 6 each), Portugal (n = 5), and Paraguay and Bolivia (n = 4 each). The majority of countries were represented in three or fewer studies.
In addition to the wide range of species and geographic regions covered, the studies also varied in methodological approaches. Multiple SDM approaches were used, including Maxent (n = 17), genetic algorithm for rule-set prediction (GARP; n = 1), generalized linear models (GLMs; n = 17), generalized boosted models (GBMs; n = 1), Bioclim (n = 1), and ensembles of multiple approaches (n = 14), with some studies comparing results across multiple methods (Salas et al. 2017). Further, studies varied in the climate models that were used, including different versions of current climate projections (CMIP3, CMIP5) as well as a wide variety of GCMs and RCP or SRES scenarios, ranging from worst case to best case scenarios (Table S1). Of the 36 studies, 3 studies used AIC model selection in choosing a best-fit SDM among candidate models.
Projections of amphibian future habitat suitability also varied in the extent to which dispersal was considered. The studies were classified into three broad categories that represented a continuum of ways of incorporating dispersal. (1) No dispersal: Future projected habitat suitability only considered within the current suitable habitat. Some even go beyond this and clip the current suitable habitat by the known distribution of a species. As a result, these studies only report percent loss of suitable habitat. Studies that do not consider dispersal ask the question “How much of the original habitat will remain under climate change?” (2) Limited dispersal: Future projected habitat suitability only considered within a buffer around current suitable habitat to where species could disperse. These studies usually use a cutoff based on expert opinion on possible dispersal distances for species or consider connectivity of suitable habitat (Schivo et al. 2019). Alternatively, some studies used a simulation to model whether the species could actually move to newly available habitat (e.g., Subba et al. 2018). (3) Unlimited dispersal: Percent of suitable habitat calculated across the whole study system (entire extent of the raster). With this approach, percent change in habitat suitability will be highly dependent on the total extent of the raster, particularly if a species is predicted to have increased habitat suitability. The greatest extents considered within this review were studies that calculated suitable habitat across entire continents (Europe, Pearman et al. 2010; Markovic et al. 2014; South America, Vasconcelos and Do Nascimento 2016). Along with dispersal, some studies also evaluated habitat suitability within protected areas (D'Amen et al. 2011; Ortega-Andrade et al. 2013; Popescu et al. 2013; Kafash et al. 2018).
To quantify predicted change in suitable habitat across studies, we conducted a meta-analysis of all studies that projected suitable habitat for at least three species (n = 24 studies). Mean percent change in habitat suitability of the studies ranged from －70.5% to 166.7%. Of the 15 candidate models describing the different approaches used across studies, and their association with change in projected suitable habitat, the model that included all predictor variables had the lowest AIC score (Table 1; Fig. 2). None of the other candidate models were supported (ΔAIC value > 11.8). Gymnophiona were predicted to have the greatest loss of suitable habitat followed by Caudata and then Anura (Table 1; Fig. 2). On average, studies that considered no dispersal did predict losses in suitable habitat whereas studies that considered unlimited dispersal predicted increases in suitable habitat (Table 1; Fig. 2). Both RCP and projection year were negatively associated with predicted change in suitable habitat (Table 1; Fig. 2).
Meta-analysis linear model results for variation in predicted change in habitat suitability across 36 studies in the literature review. Only the best-fit model results are shown (See Supplemental Table S2 (Supplemenal_Material_S1.xlsx) for all other models). Citation was included as a random variable while Dispersal Limitation, Representative Concentration Pathway (RCP), Year, and Order were included as fixed variables. Results include the lower (CI.L) and upper (CI.U) 95% confidence intervals. SE = standard error; AIC = Akaike information criterion.
For all three species, the best-fit SDM had an RM greater than 1 and, as expected, included fewer parameters than did the default SDM (RM = 1; Table 2). The default SDMs were not supported (ΔAICc > 9) for any of the species. Both the best-fit and default models showed good to moderate fit to the training data (AUC > 0.81) and test data (AUC > 0.71) for all three species.
Maxent species distribution model evaluation metrics for the three case study species (Batrachoseps). For each species, the best-fit model as selected by corrected Akaike information criterion (AICc) and the default Maxent model (regularization multiplier [RM] = 1) are shown. Results include area under the curve (AUC) for the training data and test data, variation in AUC for the test data, AICc score, ΔAICc relative to the best-fit model for each species, the Akaike weight (w) and number of parameters (k).
The best-fit model of predicted change in climate suitability was the model that included SDM type (best-fit versus default) and the interaction between SDM and Species (Table 3). Under no dispersal, the default SDM predicted greater losses in suitable habitat compared to the best-fit SDM (Fig. 3). In comparison, with limited or unlimited dispersal, the impact of SDM on predicted change in suitable habitat varied among species (Table 4; Fig. 3). Averaged across all climate scenarios, predicted habitat suitability for B. major changed less when the best-fit Maxent SDM was used (－9%, 54%) compared to the default SDM (－13%, 93%), ranging from slight losses in suitable habitat under no dispersal to moderate increases under unlimited dispersal (Fig. 3). Predicted changes in habitat suitability for B. gabrieli were similarly less extreme with the best-fit SDM (－16%, 176%) relative to the default SDM (－27%, 175%) in some dispersal scenarios, but under the unlimited dispersal scenario showed little difference among best-fit and default models (Fig. 3). In contrast, predicted habitat for B. nigriventris showed more-extensive predicted decreases with the default SDM in all three dispersal scenarios (－38%, 6%) and less of a predicted decrease or even an increase under the best-fit SDM (－15%, 63%). For all three species under the no-dispersal scenario, the default SDM frequently predicted extensive losses while the best-fit SDM predicted only slight losses or no change (Fig. 3).
Comparison of candidate models of predictors of percent change in climate suitability for three species of Batrachoseps salamanders. Percent change in climate suitability was quantified for each species using either the best-fit or default species distribution model projected onto 26 future climate scenarios (13 general circulation models [GCMs] and 2 [RCPs]) under three dispersal scenarios. The GCM, RCP, and dispersal were included as covariates in all candidate models. Results include log-likelihood (LL), corrected Akaike information criterion (AICc) score, ΔAICc, the Akaike weight (w), and number of parameters (k).
Confidence intervals for fixed effects in the best-fit model of percent change in suitable climate for Batrachoseps salamanders. Lower (CI.L) and upper (CI.U) 95% confidence intervals were calculated for the model with the lowest AICc score.
Predicting amphibian declines requires the integration of multiple tools, and species distribution modeling is one way to predict changes in suitable habitat for amphibian species. However, decisions in the modeling process can have significant implications for projections (Thuiller et al. 2019), and thus it is crucial to investigate how different modeling techniques affect predictions for better interpretation and understanding of the models generated. Through our literature review, we found extensive variation in modeling approaches used and subsequently extensive variation in predicted outcomes for amphibian bioclimatic suitability under climate change (Fig. 2). Studies predicting future bioclimatic suitability varied in the climate models and scenarios, model algorithms, year of projection, model selection approach, and consideration of dispersal among many other factors. Individual species predictions ranged from losses of 100% to gains well over 100%. While some of these differences appear to be due to species differences (Table 1; Fig. 3), variation in predicted suitability also was explained by RCP/SRES, dispersal limitation, and the projection year.
One of the primary predictors of percent change in climate suitability was the way in which dispersal was considered. For studies that considered no dispersal outside the current range, projected losses in suitability were the greatest, as expected (Table 1; Fig. 2). On average across all studies, there is a predicted change of －43% ± 2% (mean + standard error of the mean [SEM]) of bioclimatically suitable habitat within the current known or predicted ranges of amphibians. In contrast, studies that assumed unlimited dispersal across the entire extent of their study area frequently projected gains in suitability (Table 1; Fig. 2). Under an assumption of dispersal, there is a predicted change of －27% ± 5% (mean + SEM) of suitable habitat with limited dispersal and 18% ± 11% (mean + SEM) with unlimited dispersal. These results suggest that dispersal capabilities will be crucial for amphibians to be resilient to the changes in suitable habitat as a result of climate change (Carvalho et al. 2010).
How dispersal was modeled varied extensively across studies, with many authors indicating a lack of data to fully support any particular assumptions about dispersal. Although many studies evaluated the full range of possibilities, allowing for an analysis of the sensitivity of their estimates to assumptions about dispersal, some studies investigated only one potential dispersal scenario, which frequently was the no-dispersal scenario. While historically many amphibians were thought to be dispersal limited (Smith and Green 2005), there is mounting evidence that many species are capable of long-distance dispersal (da Fonte et al. 2019; Cayuela et al. In press). Even for species with highly limited dispersal, such as Batrachoseps and their close relatives (Smith and Green 2005), phylogeographic analyses suggest that extensive dispersal may have occurred during periods of historical global warming (Vieites et al. 2007). Thus, while some species may be at high risk of habitat loss due to climate change, others may experience increased habitat availability if they are able to disperse to newly available habitats. However, even if some species have extensive dispersal capabilities, habitat loss, fragmentation, and degradation due to additional human-mediated landscape modifications may restrict connectivity to climatically suitable habitats. Preservation of dispersal corridors will therefore be essential to assuring that species have access to climatically suitable habitat (e.g., Gonçalves et al. 2016; Subba et al. 2018). Regardless of the assumptions about the extent of dispersal, approaches that incorporate dispersal simulations (e.g., Subba et al. 2018) seem most promising for developing more-realistic projections for future habitat availability.
Because only three (one of which was this study) out of the 36 studies in our literature review used AIC to compare candidate SDMs (Wright et al. 2016; Kafash et al. 2018), we were unable to test the impact of this particular modeling decision on predicted change in suitability with our literature review. While many studies addressed model complexity by removing correlated variables, using a reduced set of predictor variables, conducting model averaging, or by using an increased RM (Struecker and Milanovich 2017; Subba et al. 2018), none of these approaches explicitly compare models with different numbers of parameters. As a result, such approaches could even lead to underparameterization. Although model selection has firmly become standard practice in most areas of ecology (Johnson and Omland 2004), this approach has seen limited application in SDMs despite evidence that methods such as AIC result in selection of models that improve our ability to correctly infer habitat suitability and make predictions across time and space (Warren and Seifert 2011; Katz and Zellmer 2018). Of the two studies that reported the best-fit RM values, both had RM values greater than the default value for each species. All three studies reported on average moderate declines (no more than －28%) and sometimes extensive increases (up to 129%) in projected suitable habitat. The lack of AIC-based methods in the literature perhaps results from the vast number of methods and approaches available for creating SDMs, leading to confusion and disagreement about best practices as well as the simplicity and ease of running Maxent as a black box tool (Morales et al. 2017).
We did, however, find support for our hypothesis in our case study, the results of which demonstrated that projected changes in suitable habitat vary depending on whether model selection is used to choose among a set of candidate SDMs. When comparing alternative SDMs, it is apparent that the impact of the RM on the projected suitable habitat varies across different species (Table 4; Fig. 3). For B. nigriventris, the best-fit model predicted smaller losses in suitable habitat whereas for B. gabrieli, the best-fit model predicted smaller gains in suitable habitat relative to the default SDMs (Fig. 2). In contrast, for B. major there was little difference in the predicted change in suitable habitat between the best-fit and default SDMs. For all three species, the best-fit SDMs had far fewer parameters than did the default SDMs (Table 2) and predicted a broader range of current habitat suitability. These results bolster previous research which suggests that modeling choices can result in extensive uncertainty over future habitat suitability projections across multiple taxa (Thuiller et al. 2019). Thus, researchers need to use caution in building SDMs with best practices (Morales et al. 2017), particularly before projecting models across time or space. Interestingly, in almost all cases, there was less variation in projected change in suitable habitat between different climate scenarios when the best-fit SDM was used as opposed to the default SDM. These results suggest that utilizing best practices for modeling may help to reduce our uncertainty in model projections.
It is unclear the extent to which previous models of future habitat suitability have been impacted by over- or underparameterization. Because the impact of the RM and overparameterization may vary by species, it may be necessary to revisit a number of these study systems and re-evaluate the extent to which species habitats are threatened due to future climate change. Although we are able to make estimates for global losses in predicted climate suitability through our meta-analysis, we caution that these estimates may be affected by overparameterization of models for some of the species. Joint estimates of projected changes in habitat suitability using stacked-SDMs may be of particular concern. If individual species models are overparameterized, then this bias may be amplified when combining across species, leading to underpredictions in community habitat suitability.
Beyond the impacts on predicting future habitat suitability, there are many other situations where overparameterized SDMs may be potentially underestimating habitat availability, limiting our ability to estimate effects on amphibian declines. For example, studies using current SDMs to determine conservation status may underestimate the amount of habitat currently suitable for species at risk. Studies designed to identify areas for habitat conservation or restoration, either currently or under future climate scenarios, could be affected as well. The effect of overparameterized models may also be particularly strong when transferring across space, such as when predicting potential ranges for invasive species. This issue may also arise when estimating the potential spread of emerging infectious diseases, such as Batrachochytrium dendrobatidis or Batrachochytrium salamandrivorans, a central question of concern for amphibian decline research.
Conversely, we must also consider the potential that by using a model selection approach such as AIC, it is possible that models may be too general and therefore will lead to overestimates of species habitat suitability. Which model selection approach is used may depend on the question at hand. In some cases, it may be preferable to have a more general and less parameterized model when very conservative results are preferred or when understanding how variables contribute to species distributions is most important. When conducting spatial or temporal projections it is especially important that models are not overparameterized. Alternatively, other studies may simply need very specific models for prediction of species presence within the known range of a species and be less concerned about how the model works, in which case a more highly parameterized model may be beneficial (Warren and Seifert 2011). Future research should carefully consider these potential effects and conduct comparative studies to assess the relative contributions of over- and underparameterization and the scenarios under which more or less parameterized models would be preferred.
The extreme differences in projected change in habitat suitability among the three species of Slender Salamanders in our case study warrant further discussion. On the one hand, these differences could arise from shifting climates becoming more advantageous to some species and less advantageous to others. Alternatively, the differences could reflect variation in quality of the models for each of the species. In fact, B. gabrieli, for which we unexpectedly predicted extensive expansion of habitat in the future, had the smallest sample size for locality data used to build the models. These results point to issues that can arise in modeling when using small sample sizes (Wisz et al. 2008) for species with limited ranges and narrow environmental tolerances (Saupe et al. 2012). Further, it is important to remember that these projections only indicate suitable bioclimatic habitat. For species with distinct microhabitat requirements, such as the talus slope associations for B. gabrieli, projections based solely on bioclimatic variables may vastly overestimate the availability of habitat. While some environmental constraints may be easy to build into predictive models, there are far fewer projections of changes in nonclimate-based environmental variables, limiting our ability to model the potential effects of other environmental changes. This limitation may be particularly important for species like B. major, which have projected bioclimatic suitability that overlaps regions that have been heavily urbanized or are likely to be converted to urban areas. Future studies need to focus on increased sampling for data deficient species, such as B. gabrieli, and on incorporating predicted changes in other environmental variables into suitable habitat projections.
It is important to note that our sample size for the literature review was limited in terms of the number of studies that reported predicted changes in habitat suitability for individual species. The majority of studies (n = 47) did not report percent change in suitable habitat for individual species and thus could not be included in our meta-analysis (e.g., Iosif et al. 2014). Similarly, many studies evaluated community-wide projected habitat changes but did not report percent change in suitable habitat for individual species (Hof et al. 2011; Ochoa-Ochoa et al. 2012; Lemes and Loyola 2013; Loyola et al. 2013). Twelve of the studies in our literature review only investigated one or two species (e.g., Toranza and Maneyro 2013). Although the studies we included in our meta-analysis assessed multiple species (n = 3–451), there were only 24 studies with at least three species, which made comparisons across studies difficult. Thus, comparisons of modeling approaches within studies may be subject to other confounding differences across studies. Due to wide variation among studies, lack of reporting, or small sample size, other variables not considered here include climate baseline, sampling bias, GCM, thresholding method, and amount of species range evaluated, among others. While we were not able to evaluate the effects of these modeling decisions, we account for the additional variation by including citation as a random variable in our meta-analysis. For any species of interest, a careful assessment of each of these modeling decisions needs to be considered before interpreting predicted changes in suitable habitat. Pairing predictions with field validation methods is one strategy for evaluating predictions (Courtois et al. 2016). Finally, most studies focus on changes in bioclimatic variables, but other environmental habitat changes can exacerbate the effects of a shifting climate (Rosenstock et al. 2015; Soares de Oliveira et al. 2016; Cobos and Bosch 2018). Regardless of these limitations, the small sample of studies observed here echoes previous concerns that amphibians are understudied relative to other taxa as to the potential impacts of climate change (Nabout et al. 2012).
Further, there is unequal representation of studies investigating predicted changes in suitable habitat across the three main lineages of amphibians. Although Anurans were the most widely studied with the greatest number of species for which there are projected changes in suitable habitat (followed by Caudates and finally Caecilians), Anurans and Caecilians were proportionally understudied relative to the total number of species within each of these taxonomic groups globally (see http://www.amphibiaweb.org). Few studies consider not just changes in predicted habitat for individual species but also changes in evolutionary history and ecosystem function (Loyola et al. 2014; de Pous et al. 2016; Lourenço-de-Moraes et al. 2019). In addition, the global representation of these studies is also lacking and does not fully represent patterns of taxonomic diversity around the world. Of the 36 studies included in our quantitative analyses, only one study was conducted within Africa (Ayebare et al. 2018), three within Asia (Duan et al. 2016; Subba et al. 2018; Ashrafzadeh et al. 2019), nine within North America, 10 within Europe, and 13 within South America. Most species were evaluated in only one study, and only a few studies replicated analyses for similar suites of species across the same geographic regions. Habitat suitability projections are needed for more species and across more geographic locations to better identify species and regions most at risk to the effects of climate change.
We investigated the potential impact of climate change on projected bioclimatic habitat suitability of amphibians, taking into account variation in approaches to correlative species distribution modeling. Our literature review highlights the wide variety of approaches used and assumptions made in projections of amphibian distributions. There is extensive variation in estimates of future climate suitability for amphibians. Our case study illustrates that some of this variation may be due to overparameterization of models. Based on these results, we are still unable to make clear predictions about changes in suitable habitat due to climate change for amphibian populations across the globe. Additional taxonomic and geographic coverage along with use of best practices in modeling are necessary for better estimates. Because previous predictions may be affected by overparameterization and because there is a general lack of studies predicting changes in suitable habitat for amphibians, we recommend revisiting projections for species already studied as well as creating novel projections for understudied species. Further investigation of alternative methods for species distribution modeling using rigorous, simulation-based approaches (Warren et al. 2020) will be necessary for evaluating the role that climate change may play in amphibian declines and to successfully plan for conservation of diverse amphibian lineages.
Supplemental material associated with this article can be found online at https://doi.org/10.1655/Herpetologica-D-1900066.S1.
Special thanks to the iNaturalist and N.A. Herp communities for observation data and to A. Clause and N. Hess for additional observation localities. We thank D.M. Green and the Redpath Museum (Montreal, Quebec, Canada) for organizing the symposium on amphibian declines.