Literature Review

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.

Case Study

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).

Literature Review

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.

Fig. 1

Map of countries for which suitable climate has been projected under future climate scenarios for at least one amphibian species. Colors represent the number of studies in which a given country was included. A color version of this figure is available online.

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).

Table 1

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.

Fig. 2

Mean predicted change in suitable climate for amphibian species by citation and model parameters from the meta-analysis. Points represent means and bars represent 95% confidence intervals. Mean predicted change in suitable climate was calculated for each amphibian order within each study for each of the model settings, including Representative Concentration Pathway, Year, and Dispersal Limitation. Triangles indicate estimates from the case study, and circles indicate all other studies. A color version of this figure is available online.

Case Study

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.

Table 2

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).

Table 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).

Fig. 3

Percent change in climate suitability predicted for each species in the case study. For each species, predicted percent change in climate suitability is shown for the year 2050 under 13 general circulation models and two representative concentration pathways (2.6, 8.5) climate scenarios for both the default Maxent species distribution models (SDMs; regularization multiplier = 1; open red boxes) and the best-fit SDM (Akaike information criterion; ΔAICc = 0; shaded blue boxes) for three levels of dispersal: no dispersal, limited, or unlimited. Boxplots show median, interquartile range, and minimum and maximum values. A color version of this figure is available online.

Table 4

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.

Dispersal

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.

Model Selection

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.

Case Study

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.

Limitations

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.