Climate change and the future restructuring of Neotropical anuran biodiversity

222 –––––––––––––––––––––––––––––––––––––––– © 2019 The Authors. Ecography published by John Wiley & Sons Ltd on behalf of Nordic Society Oikos This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. Subject Editor: Ken Kozak Editor-in-Chief: Miguel Araújo Accepted 13 October 2019 43: 222–235, 2020 doi: 10.1111/ecog.04510 doi: 10.1111/ecog.04510 43 222–235


Introduction
Global climate change has the potential to affect profoundly the future distribution of biodiversity and the composition of ecological communities (Bellard et al. 2012). Potential effects of climate change include alterations in species range sizes, extinction probabilities (Lovejoy andHannah 2006, Ladle andWhittaker 2011) and local richness of species assemblages (Peterson et al. 2002, Thuiller et al. 2005, Lawler et al. 2009, Reu et al. 2014. Shifts in species' ranges will also alter the spatial composition of species assemblages across landscapes (i.e. their β-diversity; Tuomisto 2010a, b, Anderson et al. 2011). There is increasing interest in exploring changes in β-diversity with climate change (Buisson and Grenouillet 2009, Ochoa-Ochoa et al. 2012, Tisseuil et al. 2012, Molinos et al. 2016, Socolar et al. 2016, Zwiener et al. 2018) because shifts may influence both ecosystem functioning and provisioning of ecosystem services (Lovejoy andHannah 2006, Cardinale et al. 2012).
In general, species assemblages should become less diverse, and thus more homogeneous, where climate change precipitates the loss of geographically restricted, niche specialists and their replacement by more widespread, niche generalists (McKinney andLockwood 1999, Clavel et al. 2011). For example, climate-driven latitudinal and altitudinal shifts of species ranges have been associated with biotic homogenization in several taxonomic groups (Menéndez et al. 2006, Davey et al. 2012, Savage and Vellend 2015, generally without driving concomitant changes in species richness (Dornelas et al. 2014, Magurran et al. 2015. The ecological consequences of biotic homogenization may be significant, and can lead to reductions in overall community and ecosystem function, stability and adaptability (Olden et al. 2004).
Examples of climate-induced increases in β-diversity have also been observed (Moritz et al. 2008, Rahel 2010, Tingley and Beissinger 2013. However, the drivers and impacts of increased biotic heterogeneity are still unclear, but any shift in the natural structure of species assemblages may have extensive ecosystem consequences (Folke et al. 2004).
We suggest that regional topography should be a strong determinant of changes in β-diversity. Mountainous areas generally are home to many geographically restricted species, especially in the tropics, where high elevation gradients and rugged terrain favours the presence of local endemic species. Lowlands and plains, however, are frequently inhabited by more widespread generalists as there are relatively few geographic barriers to dispersal and they can represent large extents of similar habitat types. We predict that highland regions, currently dominated by habitat specialists, will likely see increases in generalists and reductions in habitat specialists, as climates associated with present day lowland habitats shift to higher altitudes (Parmesan 2006, Moritz et al. 2008, Chen et al. 2011. These shifts will lead to decreasing β-diversity and increasing biotic homogeneity (Menéndez et al. 2006, Davey et al. 2012. In contrast, climate warming in lowland regions is likely to induce range contractions, with previously widespread lowland species becoming restricted to 'islands' of remnant habitat. Local extirpations in lowlands may thus lead to reductions in α-diversity but increases in β-diversity, resulting in a regional biota that is increasingly heterogenous. To date, only few empirical studies (Burgmer et al. 2007, Davey et al. 2013, Savage and Vellend 2015 have evaluated directly the effects of climate change on patterns of regional β-diversity.
Climate model projections can be used to test predictions of regional trends in β-diversity. Here, we used information on anuran amphibians' (frogs and toads) present-day distributions and 21st century climate change scenarios to assess temporal changes in β-diversity across the Neotropical region. The Neotropics, encompassing the Middle American Cordillera Central and the South American Andes Mountains which run nearly the entire north-south extent of the Neotropics over ca 75° of latitude, the greatest latitudinal extent of any mountainous region in the world, provide an ideal case study. The Neotropics are also home to the greatest diversity of amphibians in the world, with approximately 3000 known species (Frost 2017). Amphibians, specially anurans, are important components of both terrestrial and aquatic assemblages. As omnivores when larvae, and small insectivores and carnivores when adults, they serve as critical links between the lowest and highest trophic levels (Duellman and Trueb 1994). Tropical species assemblages, which are characterised by both high species richness and high endemism, tend to flourish best in warm, humid climate, but may be especially sensitive to climate change (Stuart et al. 2004, Buckley andJetz 2008). Widespread declines have occurred among Neotropical amphibians in recent decades (Wake 1991, Pounds et al. 2006) and climate change has been implicated in these declines (Pounds et al. 1999, Kiesecker et al. 2001, Menéndez-Guerrero and Graham 2013. Here, we used β-diversity metrics to quantify the degree of anuran assemblage dissimilarity among sites (Olden and Rooney 2006, Buckley and Jetz 2008, Baselga 2010, Dobrovolski et al. 2012. We then tested whether climate induced changes in species richness, ecological generalism and phylogenetic diversity (PD) have potential to drive opposing changes in β-diversity within anuran assemblages in highlands versus lowlands.

Methods
We used species distribution modelling (SDM) to explore changes in spatial patterns of anuran β-diversity across the Neotropics. We focused our analysis on anurans, frogs and toads, which represent the vast majority of Neotropical amphibians (~90%), and which demonstrate remarkable endemism (Duellman 1999). Correlative models used in SDM generally ignore the effect of distribution-limiting factors, such as dispersal limitations and biotic interactions, and thus tend to overestimate the actual extent of species distributions and summed species richness Rahbek 2011, Mateo et al. 2017). However, our focus here is on broad-scale biogeographic patterns and determinants of regional shifts in β-diversity due to climate change, not localised presenceabsences. Our species distribution models (SDMs) delineate the environmental limits that define anuran distributions at relatively coarse scales. These models are not, therefore, intended to provide fine-scaled predictions of species distributions, where micro-climate, biotic interactions and local habitat variation will be increasingly important, reflecting the narrow environmental tolerances of anurans.

Biological and climate data
We extracted data on the current distribution of Neotropical anurans from IUCN digital range maps (available on line, IUCN 2015) for 2669 species of the ~3000 known amphibian species in the Neotropics (Supplementary material Appendix 1 Table A1.1). We followed the taxonomy of Löwenberg-Neto (2014), Morrone (2014) and Frost (2017) in the delimitation of the Neotropical realm, including the Andean region. Species with taxonomic discrepancies between the IUCN species names and the taxonomy of Frost (2017, n = 6), and species classified as introduced were excluded from the analyses.
Environmental data for baseline conditions were estimated using interpolated climatic variables from the WorlClim database (<www.worldclim.org>, Hijmans et al. 2005), with a temporal span of approximately 50 years, for the period 1950-2000, and a spatial resolution of 2.5 arcmin (i.e. ~4.5 km). We note this time period already encompasses recent climate change, but pre-climate change baseline data are lacking. Future climate scenarios, covering the period between 2041 and 2080 (hereafter referred to as 2070 scenarios), included three global climate models (GCMs) under the RCP 4.5 representative concentration pathways (Stocker et al. 2013), that assumes that global annual greenhouse emissions peak around 2040, with a slight improvement and stabilization of current emissions by 2100 (Van Vuuren et al. 2011a). We used the RCP 4.5 as a conservative greenhouse emission scenario, the lower-emission scenario -RCP 2.6 -is unrealistic (requires negative emission, Van Vuuren et al. 2011b). We evaluated outputs from the following three GCMs, which illustrate contrasting temperature and precipitation predicted changes: the HadGEM2-ES model that predicts high temperature increase and a reduction in precipitation (Jones et al. 2011, hereafter referred to as the hotter/drier scenario), the GISS-E2-R model that predicts moderate temperature increase and comparatively constant precipitation rates (Schmidt et al. 2014, moderate scenario), and the IPSL-CM5a-LR model that predicts warm temperatures and reduction in precipitation (Dufresne et al. 2013, warmer/drier scenario). These GCMs are frequently used in SDM studies, including in Neotropical areas (Tomillo et al. 2015, Ihlow et al. 2016.

Species distribution modelling
Species range maps were rasterized at the resolution matching the climate data, and a climate model was then fit to the observed current distribution of each species including 19 baseline (current) climate variables using Random Forest (Breiman 2001). Range maps have been previously used in amphibian biogeographical studies at large spatial scales (McKnight et al. 2007, Buckley and Jetz 2008, Lawler et al. 2009), and recent findings showed a strong agreement between SDMs trained with IUCN range maps and SDMs trained with true occurrences (Fourcade 2016). Random forests are machine learning algorithms based on classification of regression trees returning the class that is the mode of the classes output by individual trees (Breiman 2001); they have been widely used in SDMs, and proven to be effective for predictive modelling (Prasad et al. 2006, Mi et al. 2017. Random forest classifiers are robust to over-fitting but tend to better fit the observed data than other algorithms, and hence producing more accurate predictions of species' ranges (Lawler et al. 2006). Areas inside the extent of the species ranges were treated as presences, and 10 000 background points outside of the ranges were randomly selected as pseudo-absences. We modelled only species with at least three locations (following van Proosdij et al. 2016, Supplementary material Appendix 2). For each species, the calibration region (i.e. the mobility and accessible M area sensu Soberon and Peterson 2005) was restricted to occupied ecoregion(s) following the regionalisation of Morrone (2014), see Barve et al. (2011) and Cooper and Soberón (2018) (Supplementary material Appendix 3  Table A3.2).
All models were fit in R (<www.r-project.org>) using the package BIOMOD2 (Thuiller et al. 2014), with default settings (as suggested by Thuiller et al. 2009) and using 80% of the presences and pseudo-absences for calibration. The remaining 20% of the data were used to test model performance. Models were evaluated based on the true skill statistic (TSS, Allouche et al. 2006), that varies from −1 to +1, where 0 or less represent no discrimination and +1 represents perfect discrimination (Araújo and New 2007). Models with TSS scores of <0.70 were discarded from our analyses. A threshold set to the value at which TSS is maximized (TSSmax) was used to transform the model predictions (both current and future) into maps of presence/absence (Liu et al. 2011).
The fitted models were then projected under the three different GCMs (i.e. hotter/drier, moderate and warmer/ drier scenarios) to generate estimates of species distributional responses to climate change. For our purposes here, we assumed a full-dispersal scenario (i.e. perfect climate tracking), but see following section where we describe a relaxation of this assumption.

Uncertainties in forecasts of SDM
The projection of SDMs to predict shifts in the distribution of species under future climate change scenarios, where predictions cannot be cross validated, is challenging (Thuiller 2004, Elith et al. 2010, Peterson et al. 2018. We therefore evaluated some of the key factors that may affect the accuracy and transferability of the models to assess the robustness of our downstream analyses (Supplementary material Appendix 4): • Number and grain size of predictor variables. To evaluate the effect of these two common sources of uncertainty (Pearson and Dawson 2003, Braunisch et al. 2013, Manzoor et al. 2018), we reestimated baseline and future conditions using climatic variables with a spatial resolution of 10 arc-min (i.e. ~20 km). We then refit the SDMs on the eight least correlated variables (|rho| < 0.75) that have previously been identified as biological relevant (Menéndez-Guerrero and Graham 2013, Supplementary material Appendix 4). • Extent of calibration area (equivalent to M area). Several studies have highlighted significant effects of calibration area on SDM (Lobo et al. 2008, Barve et al. 2011. Restricting the calibration area to a hypothetically accessible area (M) relevant to the focal species is suggested to increase model performance (Saupe et al. 2012, Cooper andSoberón 2018). Our base models assume calibration areas based on the ecoregions of Morrone (2014), to evaluate model sensitivity, we also generated random pseudoabsences (or background data) at a maximal distance of two degrees from the species' range boundaries, following Barbet-Massin et al. (2012), Supplementary material Appendix 4. • Thresholding method: This is a procedure by which continuous SDM outputs are transformed into presence/absence maps, and represents another common source of uncertainty in SDM studies (Thuiller 2004, Liu et al. 2013).
Our base models threshold on the TSS, to further evaluate model sensitivity, we also selected the 10th percentile threshold (P10, Peterson et al. 2007). P10 is a conservative threshold that seems to be less sensitive to extreme environmental values (Radosavljevic and Anderson 2014). It includes 90% of the presence records with the highest probability of occurrence from the predicted presence range, and thus may be better suited for our input data (i.e. species' ranges, Supplementary material Appendix 4). • Dispersal scenarios: The importance of accounting for dispersal limitations when projecting SDMs is well recognised (Araújo and Guisan 2006, Midgley et al. 2007, Thuiller et al. 2008, especially for a group such as amphibians, whose dispersal capacity is lower than that observed in other vertebrates (Foden et al. 2013). Unfortunately, there is little data on dispersal for most Neotropical species (Stuart et al. 2008). Our full dispersal scenario can be considered as representing potential range shifts, to evaluate how this scenario might depart from realised range shifts, we also implemented a dispersal scenario of 16.9 km per decade, following the meta-analysis across multiple taxa -mainly plants and insects -by Chen et al. (2011) (Supplementary material Appendix 4). • Partitioning occurrence data for model evaluation.
Evaluating the accuracy of distribution models ideally requires testing them on independent data that accounts for spatial autocorrelation between testing and training localities, especially for studies involving transfer across time or space (Hijmans 2012, Wenger andOlden 2012).
To assess the effect of temporal transferability, we used two spatial-partitioning schemes (Supplementary material Appendix 4) for conducting spatially independent evaluations of our models (Muscarella et al. 2014).

Calculating β-diversity
To derive species composition of assemblages under current and future climatic conditions, we stacked (i.e. summed) the presences from the individual species models. We use a grid cell resolution of ~50-km (Rahbek 2005) as range map data can be unreliable at finer spatial scales (Hurlbert andJetz 2007, Hawkins et al. 2008), and because our focus was on the variation in the structure of biogeographical assemblages at broad spatial scales. We use the term species assemblage to refer to the set of species within a defined geographic area, following accepted usage in the fields of biogeography, macroecology and paleoecology (Elton 1927, Vellend 2016, and it does not necessarily imply coexistence at local scales, such as within the context of community ecology. A grid cell was considered occupied if it overlapped any part of a species projected distribution (McKnight et al. 2007).
We quantified spatial patterns of β-diversity using the Sørensen pairwise dissimilarity index (β sor , Sørensen 1948). This index can be partitioned into two additive components: turnover (β sim ), accounting for species replacement, and nestedness (β nes ), reflecting the loss (or gain) of species (see Baselga 2010 for a detailed explanation). Here, we focused our analyses on the turnover component, β sim , as it has been found to contribute more than nestedness to overall β-diversity among sites , Xu et al. 2015, and nestedness is captured in our separate analysis of species richness differences (Carvalho et al. 2012). We quantify turnover as: where a is the number of species in both cells, b is the number of species exclusive to the focal cell and c is the number of species exclusive to the adjacent cell. We calculated β sim for each cell as the average of the β-diversity values between a focal cell and each of the eight neighboring cells (Melo et al. 2009). Increasing the focal cell neighbourhood (e.g. 24 or 48 cells) has shown to not substantially alter β sim values (Melo et al. 2009). To evaluate changes in β sim with climate change (reflecting how differences in species spatial composition changes over time), we projected species distributions using the future climate scenarios and recalculated β sim , and then subtracted future projected β sim from current β sim . Negative values indicate a trend towards biotic homogenization, while positive values indicate a trend towards increasing biotic heterogeneity. All β-diversity metrics were calculated using functions in the betapart package implemented in R (Baselga and Orme 2012, <www.r-project.org>).

Assessing spatial patterns of diversity change under future climates
We assessed whether changes in β sim were associated with changes in species richness (species loss versus species gains), phylogenetic diversity (PD) and ecological generalism.
Changes in species richness (ΔSR). We estimated projected ΔSR by subtracting future projected richness from current richness. Negative values represented projected species losses, while positive values represent projected species gains.
Changes in PD (ΔPD). We used the phylogenetic tree from (Pyron 2014), the most complete amphibian phylogeny at the time of analysis, to calculate Faith's PD (the sum of the phylogenetic branch lengths connecting all taxa in a set, Faith 1992). Species missing from this tree (1315) were added as polytomies at the root node for the genus using the congeneric.merge function from the pez R library (Pearse et al. 2015). We used the standardized effect size of Faith's PD (SES.pd) to quantify PD corrected for species richness under both current and future climate scenarios. SES.PD was calculated by comparing observed PD values to null expectations (999 runs) generated through the independent-swap null model (Gotelli and Entsminger 2003). Negative SES values indicate species in assemblages are more related than expected by chance (clustered in the phylogenetic tree), while positive SES values indicate species in the assemblages are less related than expected by chance (overdispersed in the phylogenetic tree). PD analyses were conducted using functions in the R package picante (Kembel et al. 2010). Negative ΔPD indicate that assemblages are projected to become increasingly clustered, whereas positive ΔPD values indicated that assemblages are projected to become increasingly overdispersed.
Changes in ecological generalism (ΔEG). Species with wide environmental tolerances tend to have larger geographic ranges than specialists, with the latter being characterized by highly constrained niches and smaller geographic ranges (MacArthur 1972, Slatyer et al. 2013). Geographic range area may therefore provide a reasonable surrogate for the ecological generalism-specialism of species. We derived the mean of the current geographic range areas for the species present in each grid cell under both current and future climate scenarios, and calculated their difference. Here, positive values represent assemblages that are predicted to shift composition towards a greater proportion of generalists, whereas negative values represent assemblages that are predicted to shift composition towards a greater proportion of specialists. Range size provides just one proxy metric of ecological generalism, we therefore also estimated niche breadth assuming an alternative surrogate of ecological generalism (Brown 1984), using the inverse concentration metrics of Levins (1968). These metrics were calculated on the suitability score maps (which are functions of all climatic variables) generated by Random forest, as described above, using the ENMTools R package (Warren 2016, Supplementary material Appendix 5).

Statistical analyses
To describe changes in species assembly structure, we fitted generalized additive models (GAMs) with change in β sim as the response variable, and four predictor variables: ΔSR, ΔEG, ΔPD and elevation (derived from the shuttle radar topography mission [SRTM], Farr et al. 2007). GAMs were fit in the R library mgcv (Wood 2017) using default settings. All predictor variables were log transformed and standardised to a mean of 0 and a standard deviation to 1 prior to analysis. There was no significant multicollinearity among predictors (all VIFs <2, Quinn and Keough 2002). We generated the set of all possible additive models using the MuMIn R package (Barton 2018), with the intention to average the best models according to the Akaike information criteria (AIC, Burnham and Anderson 2003). Models with an AIC difference below 2 (ΔAIC <2) were considered to have significant support (Burnham and Anderson 2003). To evaluate statistical significance of individual predictors, we used likelihood ratio tests to contrast the best model with reduced models, removing one term at a time. Because there is some controversy regarding the fitting of GAMs (Hastie 2017) and interpratting model coefficients can be complicated, we also performed ordinary least squares (OLS) regression, with ΔSR, ΔPD, ΔEG, elevation and the spatial term (i.e. latitude × longitude) included as predictors.

Species distribution modelling
The average fit of the SDMs on present day climate was high ( X TSS = 0.78). Although choice of climate change scenario impacted the magnitude of projected gains or losses of species, overall trends were similar across scenarios. For brevity, we present results for the warmer/drier scenario IPSL-CM5a-LR in the main text because we believe amphibians are likely more sensitive to these climate trends, results for the alternative scenarios are presented in Supplementary material Appendix 6 Fig. A6.7, A6.8).
Our models predict large changes in anuran diversity across the Neotropics. Under 2070 projections, species retain, on average, just 58% of their current suitable area (42% range loss, Supplementary material Appendix 7 Fig. A7.9), but with large variation among species (1st quartile = 83% range loss, 3rd quartile = 24% range loss). Our results match closely to previous estimates of shifts in Neotropical amphibian alpha diversity under climate change projections (Lawler et al. 2009, Hof et al. 2011, including a slight tendency for SDMs to overestimate species richness, especially in the northern Andes (Lawler et al. 2009). The implications of these biodiversity changes have been discussed elsewhere, and we do not expand on them here.

Predicted changes of spatial patterns of β-diversity
We show that areas with high Neotropical amphibian β sim are currently concentrated in the highlands, primarily in the South American Andes (Fig. 1a). In contrast, β sim is relatively low over most of the lowland Neotropics, especially in the Amazonian, Guiana Shield and Cerrado region (Fig. 1a). These patterns are largely concordant with those of McKnight et al. (2007), who used the same metric to calculate anuran β-diversity for a smaller subset of species across the entire Western Hemisphere (see also Baselga et al. 2012).
On average, β sim is projected to increase by ~36% per cell by 2070, but there is large variation across space. In the Neotropical lowlands, where β sim is currently low, β sim is predicted to generally increase (i.e. there is a trend towards biotic differentiation, X β-diversity change = 39% ± 40.29 SD, Fig. 1a, 2a-b). This pattern of β-diversity change is projected to be greatest in the Amazon Basin, and some areas of the Cerrado and Chaco regions (Fig. 2b). Exceptions include several highland assemblages in Mesoamerica, such as Mexico's Sierra Madre Occidental mountain ranges, some parts of Sierra Madre Oriental and del Sur, and some highlands of Much of the area currently characterized as high β-diversity will remain proportionally high in β-diversity (i.e. yellow areas mainly in Central and South American Andes, except for northern Andes that are predicted to show a relative reduction in β-diversity [red areas]). In contrast, β-diversity is predicted to increase in some lowland areas (i.e. parts of the Amazon Basin, Cerrado and Chaco regions) where β-diversity is currently low (blue areas); whereas β-diversity will remain largely unchanged in other lowlands (e.g. Guiana Shield and parts of Amazonia; light green areas). We used a function generated by José Hidasi-Neto to generate the map (< http://rfunctions.blogspot. ca/2015/03/bivariate-maps-bivariatemap-function.html >). Guatemala (Fig. 2b). In contrast, β sim is predicted to decrease in approximately 75% of the Neotropical highlands, where current β sim is high (i.e. there is a trend towards biotic homogenization, X β-diversity change = −8.5% ± 6.66 SD, Fig. 1a,   2a-b), including most of the South American Andes and Guiana Shield highlands, some northern and central areas of Sierra Madre Occidental, and areas of Sierra Madre Oriental of the Mexican mountain system (Fig. 2b). Lowlands America showing predicted distribution of areas of biotic homogenization with respect to anuran species diversity (i.e. decrease in β-diversity), in red, versus biotic heterogenization (i.e. increase in β-diversity), in green. The majority of highland assemblages are predicted to become more homogeneous, whereas the majority of lowland assemblages are predicted to become more heterogenous (see violin plot). Analyses were performed under the assumption that a species can reach any area with suitable environmental conditions (universal dispersion). predicted to decrease in β sim include the Pacific lowlands of South America, especially the northern South America, some Caribbean areas of Colombia, and the lowlands from southern Argentina and Chile (Fig. 2b).

Predictors of change in β-diversity
The full GAM including ΔSR, ΔEG, ΔPD and elevation was the only model favoured by AIC (Table 1, ∆AIC <2). This model explained 65% of the total variance in β sim change ( Table 1). The results of the likelihood ratio tests showed that the inclusion of the each of individual predictors had a significant impact on the full model's predictive power (all p < 0.001, Table 1, Supplementary material Appendix 8  Table A8.3). Both ΔEG and ΔSR were negatively correlated with β sim change, while ΔPD was positively correlated with β sim change. The model including ΔEG as estimated by niche breath instead of range size, showed qualitatively similar associations (Supplementary material Appendix 5 Fig. A5.6).
The increasing heterogeneity of lowland anuran assemblages (Fig. 2) is primarily driven by range contraction and local extinction of generalist species (Fig. 3a-b, see also red areas in Fig. 1b-c). Range contraction and local extinction also tend to drive decreases in species richness, but result in higher phylogenetic dispersion -thus future anuran assemblages in these regions will tend to be species-poor and comprised of taxa that are more distantly related than current assemblages (Fig. 3c, see also blue areas in Fig. 1d). The increasing heterogeneity of some Mesoamerican highland assemblages is largely explained by local species extinctions, which tend to impact more generalist species (Supplementary material Appendix 9 Fig. A9.10).
The biotic homogenization predicted in the majority of highlands (Fig. 2), is driven by increases in generalist species (i.e. range expansion) and, to a lesser extent, loss of specialists ( Fig. 3a-b). This pattern is especially clear in the central and southern Andes (Fig. 4c-d). Range expansion will tend to increase local species richness, but coincides with the loss of specialist species. In this region, future anuran communities are predicted to become more phylogenetically clustered (formed of taxa that are closer together on the phylogenetic tree) than current assemblages (Fig. 3c). Anuran assemblages in the northern Andes, will also tend to become more homogeneous with projected climate change. However, here the shift in community structure appears to be driven more by local extinctions of specialists, resulting in assemblages with lower species richness and dominated by generalists (Fig. 4a-b). The predicted homogenization of northern assemblages is estimated to be stronger than that in the southern Andes ( X β-diversity change = −18.84% ± 7.31 SD versus X β-diversity change = −6.81% ± 4.41 SD, respectively). The local extinction of specialists is similarly predicted to drive greater homogenization of anuran assemblages in the northern Pacific lowlands of South America, whereas local species gains (especially of generalists) will tend to homogenize the southern Pacific lowlands and some areas of the southern Andes (Supplementary material Appendix 10 Fig. A10.11).
OLS models produced results analogous to those obtained in the GAMs, in which ΔSR and ΔEG were negatively associated with β-diversity change, and ΔPD was associated positively with β-diversity change (Supplementary material Appendix 12 Table A12.4).

Model sensitivity and robustness
We explored robustness of our findings to various sources of uncertainty. To encompass variability in future emission scenarios, we used three different and contrasting GCMs, and show largely congruent patterns of β sim changes (Supplementary material Appendix 6 Fig. A6.7, A6.8). Our predicted changes of β sim are clearly dependent on the assumption of unlimited dispersal ability of species under climate change. However, assuming a more realistic dispersal-scenario did not importantly affect the general spatial patterns of β-diversity shifts described in this study (Supplementary material Appendix 4). Given sufficient time, it is also possible that some species might be able to adapt to novel climatic conditions in their current locations (Bellard et al. 2012), in such cases our estimates of range contraction may be overestimates.
Although the IUCN digital range maps used here have been used widely in macroecological studies, the geographic distributions of some species are much better known than others. For example, the distribution models for anuran species in the 'data deficient' (DD) category under the IUCN Table 1. Results of the best (i.e. full) generalized additive model relating change in anuran β-diversity to change in species richness (ΔSR), elevation, change in the average level of ecological generalism (ΔEG), and change in phylogenetic diversity (ΔPD). Degrees of freedom (df), Akaike information criterion (AIC) and log likelihood (logLik) are presented. As measures of overall fit, the adjusted R 2 and the deviance explained are also presented. F-values (F) and their associated probabilities (p) for each variable are shown as well. Variables were standardised to a mean of 0 and standard deviation of 1 prior to analysis. See other models' AICs and logLiks in Supplementary material Appendix 8 230 Red List (<www.redlist.org/>) are likely less reliable than those for species in other IUCN categories. Nori and Loyola (2015) recently emphasized the need to consider DD amphibian species in broad-geographic scale studies, especially for the effectiveness of future conservation assessments. To evaluate sensitivity of our results, we recalculated our β-diversity metrics excluding DD species, and found again that results were qualitatively similar (Supplementary material Appendix 13 Fig. A13.13).

Discussion
Climate change is expected to affect patterns of species diversity and the structure of species assemblages, particularly in tropical regions (Buckley and Jetz 2008). Here, we show that, within Neotropical anuran amphibians, a trend towards biotic homogenization may be observed in mountainous areas, whereas a trend towards biotic heterogenization may be more likely in lowland and plains landscapes characterized by only modest variation in altitude. Our models indicate that biotic homogenization is mainly driven by an increase in generalist species and, to a lesser extent, loss of specialists, resulting in overall greater local species richness, but lower phylogenetic diversity (PD). In contrast, increasingly biotic heterogeneity is explained by range contraction and local extinctions of generalists, resulting in reduced species richness, but higher phylogenetic dispersion of remaining species (Supplementary material Appendix 14 Fig. A14.14 for changes in individual correlates of β-diversity). Such regional variation in shifts in β-diversity, which may have important consequences for ecosystem functioning (Lovejoy and Hannah 2006), indicate the importance of considering landscape-level variation in biodiversity responses to global change drivers.

Homogenization of Neotropical highlands
Our model predictions suggest biotic homogenization, especially along the Andean highlands, associated with a shift in assemblage composition towards higher species richness and a greater dominance by generalists, which also tend to be closely related.

231
and Vellend 2015, Zwiener et al. 2018, birds: Davey et al. 2013and fishes: Buisson and Grenouillet 2009, Tisseuil et al. 2012, Molinos et al. 2016). This has led some to suggest that we may be entering a new era -the 'Homogocene'. However, the relative responses of specialist versus generalist species have rarely been considered (but see Zwiener et al. 2018). Biotic homogenization may result not only as a consequence of the increasing spread of generalists, but also from the extirpation of more specialist, range-restricted, species (McKinney and Lockwood 1999, Rahel 2002, Olden and Poff 2003, Clavel et al. 2011. The Andes, especially the northern and central part, is one of the regions with the highest richness of range-restricted anurans in South America (Villalobos et al. 2013). Small range size obviously drives patterns of β-diversity, and range size has been shown to be a key predictor of anuran species vulnerability to extinction (Botts et al. 2013). In the north and southern Andes (and some of South American western lowlands), we observed an average increase in the mean geographic range size of anuran species within assemblages under climate projections, consistent with both an increase in large-ranged generalist species, and a local extinction of small ranged specialists. In the north, where there is a very high concentration of range-restricted species, extinctions are the primary driver of homogenization, as indicated by the decrease in local species richness projected under future climate change scenarios; whereas in the south the main driver is the range expansion of species, mainly generalists, as indicated by an increase in projected species richness. The more pronounced homogenization predicted in the northern versus southern Andes, suggests that climate-driven extinctions of specialists could have greater impact on patterns of β-diversity than the range expansions of generalists.

Heterogenization of Neotropical lowlands
Increasing biotic heterogeneity can also result from species extinctions (Rahel 2002, Olden and Poff 2003, Moritz et al. 2008, Molinos et al. 2016 and range shifts (Ochoa-Ochoa et al. 2012). We predict an increase in beta diversity Median is shown as red asterisk, and outliers are shown in red. The light gray polygon (i.e. the bag) contains 50% of the data points. (d) Biplot showing the predicted relationship between change in β-diversity and ΔEG in the southern Andes. Univariate generalized additive models (GAMs) were used to predict these relationships (light gray line), showing 95% confidence interval of the prediction shaded in dark gray and residuals (light gray points). Rugs on the x-axes (i.e. vertical lines above the x-axes) show the predictors values, and how they are distributed. Labels on the y-axes indicate the smoothed function (s) for the term of interest (ΔEG), and the estimated degrees of freedom (following the term). As measures of overall fit we present adjusted R 2 and the significance of predictors (p). of assemblages in Neotropical lowlands, notably within the Amazon Basin, and some parts of the Cerrado and Chaco. This increase in heterogeneity is predominantly driven by local extinctions, often among closely related species, but also by range reductions of generalists, which results in speciespoor assemblages dominated by a greater proportion of specialists. Other studies projecting future climate-induced shifts in amphibian species ranges have also identified the Amazon and the Cerrado regions as areas prone to high species losses, and with only few species gains (Lawler et al. 2009, Hof et al. 2011). However, our study is the first to translate such shifts into changes in β-diversity at a continental scale.

Predicting an uncertain future
In amphibians, it has been suggested that tropical communities may be particularly sensitive to climate change based on the tight linking between environmental conditions and β-diversity (Buckley andJetz 2008, da Silva et al. 2014). Here we have shown how such changes may reassemble amphibian communities across the entire Neotropical region.
Concerns about using SDM to infer impacts of climate change have been discussed widely elsewhere (Thuiller 2004, Hijmans and Graham 2006, Wiens et al. 2009, Elith et al. 2010, Heikkinen et al. 2012, Peterson et al. 2018. For instance, the uncertainty around global climate models (GCMs) could be larger than other sources of uncertainty in SDM (Diniz-Filho et al. 2009, Buisson et al. 2010, and it is difficult to identify a GCM that performs better than another (Martinez-Meyer 2005). We evaluated the sensitivity of our estimated shifts in β-diversity to several types of uncertainty: GCM's, grain size, number of predictors, extent of calibration area, thresholding, dispersion, partitioning of occurrence data for model evaluation, and inclusion of IUCN data deficient species. Predictions from these alternative models were generally highly congruent (Supplementary material Appendix 4 Fig. A6, A13), increasing our confidence in the main findings we present here. Nonetheless, we acknowledge that our results showing potential changes in the spatial configuration of anuran assemblages represent just a few of many possible future scenarios (Buisson and Grenouillet 2009, Supplementary material Appendix 4). However, perhaps the greatest uncertainty associated with our analyses surrounds the societal choices that will determine future emission scenarios.

Conservation implications
The distribution of β-diversity is a fundamental component of biodiversity relevant to conservation (Socolar et al. 2016).
Here we show that, in the Neotropical highlands, future anuran communities will become increasingly homogeneous. While many species might track climate changes by moving upwards in elevation, potentially increasing local species richness in some regions, greater biotic homogenization might reduce ecological functioning and stability of anuran assemblages (Olden et al. 2004). The projected increase in β-diversity of lowland assemblages might appear to be a positive outcome; however, it is largely driven by range contraction and localized extinction of more generalist species. The local extinction of such species might have lesser impacts on ecosystem functioning due to the high functional redundancy in species-rich amphibian assemblages (Naeem 1998, Strauß et al. 2010), but as their losses compound we might experience non-linear and potential catastrophic loss of function in the future. The changing β-diversity landscape, alongside the sharp decline of the amphibian richness around the world, represents another challenge for conservation biologists.

Data availability statement
IUCN digital range maps of Neotropical anurans are available on line at <www.iucnredlist.org>. Environmental data for baseline and future conditions are available at: <www.worldclim.org>. R code for SDM is available in Supplementary material Appendix 15.