Assessing future habitat availability for coastal lowland anurans in the Brazilian Atlantic rainforest

ABSTRACT Global warming is expected to cause several modifications to physical environments, and sea level rise is a certain outcome. However, assessment of the potential impacts caused by sea level rise on biodiversity is still emerging. Therefore, we assessed the combined impact of global warming and sea level rise on the potential distribution of 19 coastal lowland anurans in the biodiversity hotspot Atlantic Forest. We applied a correlative species distribution model (SDM) (BIOCLIM) and GIS-based spatial analyses. We evaluated the extent of changes of potential distributions under extreme and moderate global warming scenarios as well as two extreme sea level rise scenarios. Our results suggest wide areas of suitable habitat for most species in the future. However, for 15% of these species the SDMs predict massive losses of range extent as a result of a combination of global warming and sea level rise. Such observations highlight an immediate need to consider the potential effects of sea level rise in conservation action plans. Since the current potential distribution of these anuran species is likely underestimated, we also analyzed their environmental niche under current conditions in order to provide a baseline for further field surveys. Considering this current state of knowledge for such species, species distribution modeling to help gather further information on unknown species is desirable.


Introduction
Global warming (GW) is a major threat to biodiversity worldwide (Brooks et al. 2006; Barnosky et al. 2011). The increase of atmospheric CO 2 in recent decades is expected to affect several environmental parameters in the short term, mainly air temperature (IPCC 2007), and it is vital that action is taken to mitigate predicted effects on biodiversity. Among several potential consequences, sea level rise (SLR) is particularly worrying for coastal terrestrial biota. In less than 100 years, water levels are expected to rise due to increased atmospheric warming (IPCC 2007;Bamber et al. 2009;Shepherd et al. 2012), changing worldwide costal shorelines (Menon et al. 2010). Although the Fourth Assessment Report (AR4) of the Intergovernmental Panel on Climate Change (IPCC) does not predict SLR to be above 60 cm (IPCC 2007), other studies focusing on SLR suggest much more alarming scenarios, with rising sea levels varying from 2 to 6 m by the year 2100 (Overpeck et al. 2006;Grinsted et al. 2010;Nicholls & Cazenave 2010). However, even though SLR is one of the most assured consequences of GW in the near future (Nicholls & Cazenave 2010), its effects on coastal terrestrial organisms remain largely understudied, hindering preventive conservation actions. To date, few studies have addressed potential consequences of SLR on terrestrial wildlife (e.g. Menon et al. 2010;Wetzel et al. 2012;Bellard et al. 2014;Oliveira et al. 2015). Consequently, a comprehensive assessment of the possible effects of SLR on terrestrial biodiversity is still needed.
Ecological communities and ecosystems related to shorelines (such as mangroves, saltmarshes, sandbanks, and estuaries) are the most threatened by SLR (Menon et al. 2010;Traill et al. 2011). SLR may additionally affect terrestrial ecosystems through water intrusions (e.g. Rios-López 2008) or permanent flooding (Menon et al. 2010;Bellard et al. 2014), jeopardizing biodiversity and ecosystem services. Menon et al. (2010) conducted one of the first studies assessing potential future impacts of SLR on coastal fauna at the global scale, and identified regions most likely to be affected. Recently, Bellard et al. (2014) investigated how insular hotspots and endemic species may be affected in the near future under distinct SLR scenarios.
A failure of conservation strategies is commonly observed as a result of insufficient baseline information or ineffective implementation (James et al. 1999;Nóbrega & De Marco 2011). The resources available for conservation in terms of time and money are often limited; prioritization of resources should be supported by theory to minimize biodiversity loss (Brooks et al. 2006). Thus, the scientific community should provide comprehensive guidance to facilitate resource allocation (Morais et al. 2013).
Although SLR is predicted to moderately affect the Brazilian shoreline (Menon et al. 2010), this area includes part of the Atlantic Forest, a biome that harbors astonishing biodiversity, rarity, and endemism levels (Ribeiro et al. 2011;Haddad et al. 2013;Toledo et al. 2014). However, humans have intensively disturbed the Atlantic Forest over the past 500 years (Ribeiro et al. 2009). The impact of such disturbance, together with the fact that anurans are the most threatened vertebrates on Earth (Hoffman et al. 2010) and quite sensitive to future GW (Foden et al. 2013;Oliveira et al. 2015), raises concern about their conservation in the Atlantic Forest. Thus, amphibian conservation assessment in this region is necessary and further studies are needed (e.g. Loyola et al. 2013).
Species distribution models (SDMs) have become a useful tool in ecology, evolution and conservation biology (e.g. Peterson & Nyári 2007;Capinha & Anastácio 2011;Thuiller et al. 2011). These models correlate information on a species' distribution with bioclimatic data in an attempt to predict areas providing suitable environmental conditions (i.e. which are part of the species' Grinnellian niche; see Soberón 2007). Owing to the many threats to biodiversity in the future, definition of conservation priorities is urgent. In this context, mitigation of biodiversity loss needs to be performed with the data at hand, given that the time needed for gathering exhaustive information is impractical and impeditive (e.g. Pearson et al. 2007;Pie et al. 2013).
Therefore, we evaluate the future availability of suitable habitat for terrestrial coastal lowland species under different scenarios of GW and SLR. We applied a comprehensive modeling and GIS approach to 19 coastal lowland anuran species serving as biological models. Since distribution information for these species may be underestimated, we also explore the species' potential distribution under present conditions in order to guide and foster further field surveys.

Geographical extent
The Brazilian Atlantic Forest was one of the largest rainforests in the Americas, covering approximately 150 million ha in the past (Ribeiro et al. 2009). Now it has less than 12% of its original cover but it is still one of the world's biodiversity hotspots, owning a high number of threatened and endemic species (Ribeiro et al. 2011;Haddad et al. 2013). Longitudinal, latitudinal, and elevational variation is largely responsible for high habitat heterogeneity in Atlantic forests: there is a decrease of rainfall from the coast to the west causing distinct forest compositions; and the tropical forest in the north becomes subtropical toward the south (Ribeiro et al. 2009). Such complexity may explain the high amphibian diversity (Toledo & Batista 2012) and level of endemism, with about 90% of species in the Atlantic Forest are endemic (Haddad et al. 2013). In addition, the topographic complexity and elevation range (from sea level to about 3000 m) increases anuran rarity (Toledo et al. 2014). Thus, the Atlantic Forest is a priority for wildlife conservation (Zachos & Habel 2011).

Species data
We used anurans as biological models because they are highly sensitive to environmental shifts (Wake 2012;Loyola et al. 2013). Also, although this group is severely threatened (Hoffman et al. 2010), it is comparatively well known (Foden et al. 2013), which favors data access and the comprehension of potential effects.
Anuran species were selected according to the following criteria: (1) they are endemic to the Atlantic Forest; (2) information is available about altitudinal ranges; and (3) the distribution range is within coastal lowland areas. Although interesting for our aim, we exclude insular species because some islands have such a small extent that environmental information is not available. We selected exclusively Atlantic Forest lowland coastal anurans based on those amphibians listed by Haddad et al. (2013). We considered as "coastal lowland" those species with currently known altitudinal range distribution up to 60 m above sea level (asl) and no farther than 100 km from the coastline. In total we considered 19 species (Table 1). We accessed altitudinal range information for each species at the IUCN website (IUCN 2013) and in Frost (2016). Although incomplete and often biased for several reasons , occurrence points are valuable in SDM approaches (e.g. Capinha & Anastácio 2011;Nóbrega & De Marco 2011) and the use of range distribution maps to generate artificial points as surrogates may imply additional bias, since their accuracy needs to be reassessed in some regions, especially in South America (Ficetola et al. 2014). Additionally, there are no available range maps for some species present in our dataset; therefore, the use of point locations was the only option to assemble the models.
For each selected species we compiled occurrence records from the Specieslink online database (http:// www.splink.org.br/), which allows data to be accessed from the following collections: Museu de Zoologia "prof. Adão José Cardoso"; Coleção "Célio F. B. Haddad", Coleção Científica do Departamento de Zoologia e Botânica da Universidade Estadual Paulista (São José do Rio Preto), Coleção de Tecidos Animais do Departamento de Ciências Biológicas da Universidade Federal do Espírito Santo; and Coleção de Anfíbios do Museu de Biologia Mello Leitão. Since some selected species were only recently described, few or even no records were available for some of them. In such cases, and in order to increase the number of occurrence points as well, additional records were searched in GBIF (http://www.gbif.org/), in peer reviewed literature and in Frost's (2016) website. All occurrence records were verified for possible geographical or altitudinal errors in ArcGIS 10 (ESRI 2011; ESRI, Redlands, CA, USA) and Google Earth Pro TM (Version 7.1.2.2041; Google Inc., Mountain View, CA, USA).
Following our selection criteria our species dataset comprised a total of 19 species (Figure 1). These species are allocated into seven families and have a total of 91 occurrence records ( Table 1). Most of these species are not considered to be threatened, being classified as least concern (LC, n = 4), data deficient (DD, n = 6), or near threatened (NT, n = 1). Three of them are currently threatened, being classified as endangered (EN, n = 1) or vulnerable (VU, n = 2) (IUCN 2013). Additionally, five species are reported to have decreasing populations, and nine species have a population trend assigned as unknown. Finally, five species have no additional information in regard to their distribution range, occurrence points available in online database, assessments for threats or population trends.

Environmental data
We selected two extreme SLR scenarios of 3 and 6 m to simulate future potential loss of land of available habitat. We obtained SLR scenarios as raster files from the Center of Remote Sensing of Ice Sheets (https://www. cresis.ku.edu). Although according to IPCC's AR4 (IPCC 2007) rising sea levels above 60 cm are not expected, several uncertainties are involved in this topic and it is possible that SLR may reach several meters (Overpeck et al. 2006;Bamber et al. 2009;Grinsted et al. 2010) in different parts of the world (Gomez et al. 2010). The environmental variables used to generate the SDMs include 19 bioclimatic variables at a spatial resolution of 30 arc-second (approximately 1 km). For current conditions we obtained these climatic variables from the WorldClim homepage (http:// www.worldclim.org); and future climatic scenarios were accessed from the CIAT (http://www.ccafs-cli mate.org/). In order to compare the resulting outputs, we generate models for two future distinctive scenarios for the mid-point of a 30-year period (2071 to 2099), i.e. 2080; namely A2a, considered pessimistic, and B2a, optimistic (IPCC 2007). We used the three following climatic models projections for the year 2080: CCCMA-CGCM2, CSIRO-MK2, and HCCPR-HadCM3 developed by IPCC AR4 (IPCC 2007). We calculated the arithmetic mean of each climatic variable among the three climate models to reduce the uncertainty of forecasted climates that arises from variability among distinct climate models.
We avoided subjectively selecting variables since the specific bioclimatic variables differ in type, magnitude and significance for each analyzed species (Lawling & Polly 2011). However, since the variables describe means and extremes of temperatures and precipitation values, it is required to take into account the multicollinearity among the bioclimatic variables. To solve this issue we performed a principal component analysis (PCA), removing the collinearity and reducing the dimensionality of the climatic data. We used the four principal components (PCs) produced by PCA that encompassed the greatest cumulative variance (92%) through the data as input in each SDM.

Species distribution model
We assembled SDMs for the year 2080 to assess potential effects provided by future climate conditions and future SLR. We investigated expected potential effects of SLR on future habitat availability under GW by combining climatic scenarios and water levels in a logical manner. We explored a reasonable upper bound of a 3 m SLR (Hansen 2007) and an extreme scenario of SLR up to 6 m (Overpeck et al. 2006;Menon et al. 2010;Wetzel et al. 2012). We used this last extreme scenario in order to detect potential losses caused by sea level variations (e.g. lunar effects, marine inundation or storms), since the horizontal intrusions are difficult to predict (Wetzel et al. 2012;Bellard et al. 2014). We coupled the 3 m SLR with the B2a climatic scenario, and the 6 m with the A2a. Furthermore, we built SDMs for all species under current climatic conditions to provide additional information and also a starting point for future field survey in view to increase available information for such species.
Although there are a plethora of methods to assess potential distributions (Franklin 2009), we used the BIOCLIM algorithm implemented in the free software DIVA-GIS 7.5 (available at http://www.diva-gis.org/) to assemble the SDMs. BIOCLIM is a useful correlative presence-only modeling tool that summarizes a species' climatic envelope to predict its potential distribution (Beaumont et al. 2005;Araújo & Peterson 2012). It calculates the similarity between conditions by associating the values of climatic variables at any location of the study area to a percentile distribution of the values at previous known occurrences, ascribing suitability according to the distance to the 50th percentile (the median), i.e. the closer to the median (Hijmans & Elith 2013) (available at http://cran.r-project.org/web/ packages/dismo/vignettes/sdm.pdf).
We established distinct thresholds for each species in the SDMs, based on the ROC (receiver operating characteristic) curve plot-based approach to derive binary vector models (0/1). In this case, the threshold corresponds to the "most northwestern" point in the plot (Liu et al. 2005;Jiménez-Valverde et al. 2011). We opted for this method to minimize omission and commission errors, producing potential range maps with lower uncertainty while also reducing the risk of underestimating the extent of suitable sites. The same threshold criterion was used to generate binary range maps under current climatic conditions. We analyzed 95 models, three SDMs outputs (current, B2a and A2a) and two SLR scenarios implemented for 19 species. To evaluate the fit of SDMs we used the true skill statistics (TSS) method, which takes into account omission and commission errors to generate values that vary from −1 to +1 (Allouche et al. 2006). TSS values close to zero or −1 indicate model performance no better than chance, while the opposite situation, i.e. TSS values close to +1, suggest model with good discrimination capacity (Allouche et al. 2006).

Quantifying available suitability and the magnitude of potential impacts
We first calculated the total extent of predicted climate suitability (SDMs outputs) in number of grid cells and then converted them in square kilometers to assess the potential suitability and evaluate the magnitude of potential consequences of GW and SLR. After that, we extracted the potential loss of land by flood caused by SLR by overlapping and calculating the aggregate intersection in km 2 (Figure 2). Consequently, we could assess the combined impact of GW and SLR over potential suitable habitat. Then, we tested the significance of the dissimilarities observed in the range shifts among the two GW scenarios, as well as the potential effects provided by SLR on predicted suitable habitat availability, by the application of a simple paired Wilcoxon test, after the Shapiro-Wilk normality test.

Future suitable habitat availability
In general we obtained good model fit across SDMs, with TSS values varying from moderate (0.39) to high (1) for most species (global average ± SD = 0.86 ± 0.21) (Table 1). Therefore, all SDMs were used in the successive analyses.
Considering the final suitable area, i.e. the area estimated by modeling after taking into account the effects of GW and SLR, no cases of complete loss of future climatically suitable areas owing to either GW or SLR were detected. Instead, we observed divergent patterns showing final increased suitable areas for eight and nine species in B2a and A2a, respectively, whereas decreased final habitat size was detected for 11 species in B2a, and for 10 in A2a (Table 2). The B2a scenario combined with SLR of 3 m showed relatively low impacts with five species presenting approximately 10% and two more than 30% potentially suitable area loss (see an example in Figure 3). The A2a scenario coupled with SLR of 6 m revealed greater potential losses for the taxa in general, with 10 species presenting losses of approximately 10%, and four other species with over 20% suitable area loss (Table 2). For this scenario one species (Crossodactylus lutzorum) even loses more than 40% of the suitable area ( Figure 3D, E, F). Nevertheless, on average the A2a scenario suggested larger potential distributions (24,781.12 ± 30,226.37) than the B2a scenario (12,228.6 ± 14,751.8). We also found differences concerning the magnitude of changes in suitable area extent between GW scenarios (p = 0.018), with the A2a scenario producing greater changes than the B2a scenario ( Figure 4). Yet, outputs for the A2a scenario suggested larger potential distributions (26,662 km 2 ± 32,194) than those for B2a scenario (12,904 km 2 ± 15,410). In addition, SLR decreased future availability for suitable habitats in all cases. We observed differences between SLR scenarios (p = 0.0001) (Figure 4) as well, with the 6 m scenario affecting the potential distributions to a greater extent (1881 ± 2246) than the 3 m scenario (675 ± 770).

Discussion
Amphibians are facing high extinction risk (Hoffman et al. 2010;Wake 2012) that is likely to severely decrease biodiversity in the group at alarming rates (Wake 2012). Here we reveal an additional threat to be considered together with the other claimed threats. Understanding and quantifying future potential GW and SLR impacts on terrestrial fauna synergistically is still challenging (see Oliveira et al. 2015). For instance, Bernardo-Silva et al. (2012) identified priority areas for the conservation of two endangered red-bellied toads (Melanophryniscus dorsalis and M. montevidensis) in Brazil and Uruguay based on their climatic requirements. However, at least M. dorsalis, which is restricted to environments such as coastal dunes, is likely to suffer from SLR. Moreover, Oliveira et al. (2015) identified species from the Australasian biogeographic region likely to be highly sensitive to future loss of suitable habitats owing to GW and SLR coupled. Therefore, since available resources for implementing Conservation Units are limited, we argue that potential effects of SLR need to be considered in conservation strategies.
Although we identified 19 taxa likely to potentially suffer suitable area loss, we stress that the actual number of Atlantic Forest species under SLR threat is underestimated. For example, exclusively insular anurans were excluded (e.g. Cycloramphus faustoi, Scinax alcatraz, S. faivovichi, and S. peixotoi, all assigned as CR) owing to data limitation. Moreover, although dispersal overseas is a possibility for amphibians (Vences et al. 2003), we lack evidence showing that this is common. Yet, species isolated on islands possess environmental requirements and physiological constraints (Wells 2007), representing a special challenge for forecasting future distributions.
Our modeling approach may also improve categorization of endangerment of some species, as in the case of Haddadus plicifer and Leptodactylus hylodes. Morais et al. (2013) provided an overview on the risk of extinction and categorization in regard to data deficient (DD) species. These authors considered the extent of occurrence and the time passed since the description of the species to infer the threat category of DD species. According to their criteria, both species assigned as DD (IUCN 2013) are in fact CR.
Conversely, even after discounting losses caused by future SLR, L. hylodes is still predicted to have an increase in the extent of climatically suitable areas. However, H. plicifer is predicted to have losses caused by both GW and SLR, with only less than a half of its potential current suitable habitat remaining. These two species are among those that have the most restrictive habitat constraints, and based on the long time passed since their description, it is likely that they are rare (Morais et al. 2013).
The A2a scenario predicts a greater amount of atmospheric CO 2 and higher air temperatures in comparison to B2a (IPCC 2007). However, our results showed significant larger potential distributions under the A2a scenario, even with potential effects of SLR accounted. Patterns of expansion and contraction of suitable habitats are common when dealing with SDMs due two main reasons: the assumption that current climate is a good descriptor for species distribution ) and the non-consideration of the influence of biological interactions to determine species range (Soberón 2007). In fact, GW may provide increased habitat suitability for few taxa and decrease for others (Erasmus et al. 2002). For the species we analyzed, although several are poorly known, it seems that some of them may benefit from a warming climate, which is evidenced by a considerable increase in the extent of suitable habitat in the future, even with potential losses by SRL (e.g. Phyllodytes punctatus, Sphaenorhynchus palustris). Furthermore, there are species for which expansion of suitable areas (e.g. Aparasphenodon arapapa) or potential distribution towards inland (e.g. Chiasmocleis atlantica) may occur, diminishing potential effects of SLR (see supplementary material). Concerning the current habitat availability, our assessment provides a starting point for future field surveys. Because several species have only been recently described (e.g. Dendropsophus skuki and Melanophryniscus setiba), their real distributions are potentially broader than currently known. In addition, models revealed climatic suitability in mountainous regions for some species (e.g. Chiasmocleis atlantica, C. capixaba, and C. lacrimae), and such species may not be restricted to coastal areas. On the other hand, several species showed restricted potential habitat ranges (e.g. Crossodactylus lutzorum, Dendrophryniscus skuki, Haddadus plicifer, Leptodactylus hylodes, and Phyllodytes punctatus), including some described a long time ago (which suggest that they are truly rare) and, in such cases, development of further field surveys is necessary to improve the existing knowledge about their distribution range.
There are several factors that can bias SDMs, such as spatial autocorrelation and non-uniform data samples. Species with restrictive distributional range represent a special challenge to be overcome in SDM approaches, since few records spatially clustered tend to reduce statistical power (Wisz et al. 2008). Nonetheless, cases in which the number of occurrence records can be negatively correlated to the model performance (e.g. Girardello et al. 2009) highlight how complex it is to evaluate the models (Liu et al. 2005;Lobo et al. 2008). Although SDM predictions still require better evaluation approaches (Diniz-Filho et al. 2012), understanding factors that can jeopardize the conservation of species is also urgent, and cannot wait for the accumulation of complete knowledge (Pie et al. 2013). In this sense, our results stress potential events that may severely impact habitat availability for coastal lowland anuran species, and provide a starting point for further investigations in the Atlantic Forest. In agreement with the recommendations of the Brazilian Amphibian Conservation Action Plan (Verdade et al. 2012) and recent analysis of worldwide coastal anurans (Oliveira et al. 2015), we suggest long-term population monitoring, including of coastal lowland species, since this practice may enable the detection of population fluctuations, declines or local effects owing to GW and especially SLR.