Climate change can lead to a geographic range shift of species in the future, which might challenge a species to maintain viable populations with lower dispersal abilities over time. Therefore, protecting stable habitats is important for the conservation of these species. Herein, we assess the effectiveness of the Brazilian protected areas to preserve the rare and threatened coral snake Micrurus brasiliensis and explore the occurrence of stable habitat areas through its geographic distribution at the end of 21st century. We used ecological niche modeling to generate the potential distribution of the species in the present, and then projected its distribution to past and future climatic scenarios. We assessed whether Brazilian reserves would encompass suitable habitats in the future and proposed areas in which conservation efforts could be directed based on habitat stability (refugia) over time. Our findings show that the potential distribution of M. brasiliensis have shifted over the time, and there is an expected decrease of more than 60% in the amount of suitable areas in the future. The protected areas will contain climatically less suitable areas in the future. We strongly suggest expanding the existing reserve network as well as the creation of corridors between protected areas, allowing the dispersal of M. brasiliensis, enhancing the opportunities for preserving viable populations in the long term.
Introduction
In the last century, anthropic activities raised the atmospheric mean temperature about 0.6℃ to 0.8℃, affecting the distribution and survival of many species worldwide (Colwell, Brehm, Cardelus, Gilman, & Longino, 2008; Parmesan, 2006; Parmesan & Yohe, 2003). By the end of the 21st century (2100), some more pessimistic scenarios of climate warming predict that global mean temperature will increase around 4.8℃, changing patterns of precipitation (IPCC, 2014). Hence, future climate changes are expected to affect individual species or entire ecosystems (Thomas et al., 2004; Walther et al., 2002).
Range shift is one of the most common responses of species in the face of climate change. It has been observed, for example, that some species are tracking suitable habitats and expanding their geographical ranges to higher latitudes and altitudes as a consequence of global warming since the industrial revolution (Araújo, Thuiller, & Pearson, 2006; Hughes, 2000; McCarty, 2001). Additionally, species may overcome stressful conditions through phenotypic plasticity or evolutionary adaptations to novel climatic conditions (Hoffmann & Sgrò, 2011; Williams, Shoo, Isaac, Hoffmann, & Langham, 2008). However, such adaptations are more likely to occur in species with short life cycles, large population sizes, and high genetic variance (Bradshaw, Holzapfel, & Crowder, 2006). Species with low dispersal abilities or lower reproductive and adaptive rates may not be able to overcome rapid climate change, and thus, they would be especially endangered under a scenario of future global warming (Walther et al., 2002).
Few studies have investigated the impacts of climate change on snakes, although some examples have provided evidence that, under a rapidly changing climate, most species would be committed to extinction (Lawing & Polly, 2011; Penman, Pike, Webb, & Shine, 2010). Coralsnakes (Family Elapidae) is a group particularly vulnerable to reduction in climatically suitable habitat due to its specialized habitat requirements and low dispersal abilities (Marques, Almeida-Santos, & Rodrigues, 2006), making it difficult for them to colonize suitable areas outside the existing range. Moreover, the high level of habitat fragmentation in tropical regions may impose additional restrictions for these species to colonize new areas, reducing their distribution to islands of natural areas embedded in a landscape of unsuitable habitats (Araujo, Cabeza, Thuiller, Hannah, & Williams, 2004; Ferro, Lemes, Melo, & Loyola, 2014). Within the coral snakes, the species Micrurus brasiliensis is probably highly susceptible to the combination of climate change and habitat fragmentation (Silva Jr., 2007). It is a small species (around 60 cm in length on average), characterized as having specialized fossorial or semi-fossorial habits and low dispersal ability (Almeida, Prudente, Curcio, & Rodrigues, 2016; Silva Jr., Pires, & Feitosa, 2016), which is distributed in the Brazilian savannas (also known as the Cerrado biome) as well as the transition zone between the Cerrado and Caatinga biomes. The few known individuals of M. brasiliensis were found in rugged terrain with open vegetation and sandy soils in areas known as dry forests; it is considered an endangered species due to the expansions of the soybean and sugar cane monocultures in the Central-West and Northeast geopolitical regions of Brazil (Silva Jr., 2007).
One opportunity for preserving species with low potential for dispersal is to preserve areas that are climatically or environmentally stable over time (Ashcroft, 2010; Oliveira et al., 2015; Terribile et al., 2012). Such environmental refugia allow viable populations to survive when conditions around them are unsuitable, enabling them to colonize adjacent areas when climatic conditions become favorable again (Ashcroft, 2010; Provan & Bennett, 2008; Tzedakis, Lawson, Frogley, Hewitt, & Preece, 2002). Moreover, refugia are extremely important because they harbor a high level of genetic diversity and variation between populations in different refugia (Carnaval, Hickerson, Haddad, Rodrigues, & Moritz, 2009; Provan & Bennett, 2008). In the case of species with lower dispersal abilities, refugia would allow the occurrence of populations in situ without needing to disperse long distances or be artificially translocated (Terribile et al., 2012). Thus, by considering the ongoing climate change and the expected high velocity of changes in the near future, identifying and evaluating the quality of refugia for such species is imperative so the impacts can be softened and extinction risks reduced (Dobrowski, 2011).
Given that not all areas potentially considered as refugia can be protected, since some of them have high potential for food productivity and are therefore useful for economic purposes, an evaluation for the potential of the existing protected areas to act as climatic refugia for vulnerable species is needed (Araujo et al., 2004). From this, additional refugia that complement the existing ones can be proposed, ensuring the long-term persistence of the species. Here, we used ecological niche modeling methods to find the geographic distribution and identify refugia for M. brasiliensis following the method of Terribile et al. (2012), and assessed whether the Brazilian reserves would be effective in protecting this coralsnake at present and at the end of the 21st century by considering climate change scenarios. Additionally, we identified additional areas based on the existence of populations within the continuous area of climatic refugia and the proximity to the current protected areas to guide future conservation efforts for this species.
Methods
Species Data and Climatic Variables
We obtained data on occurrence records for M. brasiliensis from Campbell and Lamar (2004), Silva Jr. (2007), and from field expeditions conducted by NJS Jr between 2007 and 2012. These data were carefully checked to validate species identity and locality data (Figure S1). Occurrence records were mapped on a grid with 0.5° × 0.5° resolution of both latitude and longitude that covered South America entirely.
We also obtained the climate data need for ecological niche modeling (see below) to model the species distribution at the present and project for the past and the future. Thus, we used climatic simulations of three time periods: from the preindustrial (representing current climate conditions), Last Glacial Maximum (LGM, 21,000 years ago—21ky BP), and future (mean data between 2080 and 2100 representing the climatic conditions for the end of 21st century), derived from five coupled Atmosphere-Ocean General Circulation Models (AOGCM)—Community Climate System Model (CCSM), Centre National de Recherches Météorologiques (CNRM), Goddard Institute for Space Studies (GISS), Model for Interdisciplinary Research on Climate (MIROC), and Meteorological Research Institute (MRI) (see Table S1, in Supplementary Material). For the future prediction, we used the RCP4.5 emission scenario, which is an intermediate scenario between the more optimist (RCP2.6) and the more pessimistic (RCP8.0) ones (Taylor, Stouffer, & Meehl, 2012). These data are available in the ecoClimate database at a 0.5° spatial resolution ( http://ecoclimate.org, Lima-Ribeiro et al., 2015, see also Table S1 in Supplementary Material for more details about the AOGCMs). The annual mean temperature, annual amplitude of temperature variation, precipitation of the wettest month, precipitation of the driest month, and precipitation of the hottest quarter variables were selected from 19 bioclimatic variables (according to Hijmans, Cameron, Parra, Jones, & Jarvis, 2005) through a factorial analysis with varimax rotation. The factorial analysis allows us to select orthogonal variables eliminating or decreasing the effects of the multicollinearity among the predictors in the modeling processes. We selected variables with the highest correlation to each of the five resultant factors (see details in Terribile et al., 2012).
Niche Modeling and Potential Distribution
We applied the ensemble forecasting approach (sensu Araújo & New, 2007; see also Diniz-Filho et al., 2009; Terribile et al., 2012) to generate consensus predictions about geographic distributions after combining the outputs from several different niche modeling methods and AOGCMs. Thus, the niche of the species was modeled using presence-only (BIOCLIM, Niche factor analysis, Euclidian distance, Gower distances, Mahalanobis distance), presence-background (maximum entropy—Maxent), and presence-absence methods (generalized linear models, generalized additive models, flexible discriminant analysis; multivariate adaptive regression splines; neural networks; and random forest; see Franklin, 2009 and Peterson et al., 2011 for a review of methods). As absence data are not available for this species, to satisfy those modeling methods based on presence-absence, we randomly selected pseudo-absences through the set of cells from which the species was not recorded, keeping prevalence equal to 0.5 (thus generating a dataset consisting of 50% presence and 50% pseudo-absence records). For each model, occurrence data were divided into two subsets: 75% of presence cells selected for calibration and 25% for testing the model’s predictive ability, repeating the sampling process 50 times. Thus, 3,000 models were generated through the combination of 12 modeling methods, 5 climatic models, and 50 repetitions (12 × 5 × 50) for the present, and projected into the past and future climatic scenarios. Each of the 50 models were converted into binary distribution (presence and absence, or 1 and 0, respectively) maps based on thresholds established by the area under the receiver operating characteristic curve, known as the area under the curve (Fielding & Bell, 1997). The frequency of presence of the species in each cell in these 50 models was used to generate the species habitat suitability maps, varying from zero (the species was recorded as absent in a cell in all the 50 results) to one (the species was recorded as present in all the 50 results; Figure S2). To generate the final maps of potential distribution through time, the habitat suitability maps from the combination of ecological niche models and AOGCMs were truncated by the lowest suitability in a presence record, which was 0.59 (lowest presence threshold; Pearson, Raxworthy, Nakamura, & Peterson, 2007).
Changes in Habitat Suitability Over Time and the Effectiveness of the Protected Areas
We assigned 0 (not protected) or 1 (protected) to each grid cell based on its overlap with the network of protected areas (hereafter PAs) from The World Database on Protected Areas (database available on http://www.protectedplanet.net). Then, we compared the mean habitat suitability by considering the following questions: (a) Does the current distribution of M. brasiliensis show higher suitability than the distribution projected for the future, regardless whether it is inside or outside of PAs? (b) Do the PAs exhibit higher suitability than the areas outside them for present and for future? To answer these questions, we first considered the variation in habitat suitability over time as the focus of the analysis and used a paired t-test with grid cells to compare suitability between present and future scenarios. We compared present versus future suitability inside and outside PAs, as well as across the entire South America (i.e., with no distinction between protected and nonprotected areas). Second, we considered the PAs as the focus of analysis and compared if protected areas attain higher suitability for M. brasiliensis than nonprotected areas. In this case, we compared suitability inside protected areas versus outside them separately for both present and future scenarios. All comparisons were performed using the package stats, with random resampling to set p-values in R (v. 3.2.3).
Areas of Habitat Stability
To identify areas of habitat stability (or refugia), we followed the above-described approach by converting the continuous frequency (suitability) maps from each time period into binary presence–absence predictions using the threshold equal to the lowest suitability value associated with an occurrence record (Aranda & Lobo, 2011; Pearson et al., 2007). Thus, cells with habitat suitability equal or higher than this threshold in the three time periods (past, present, and future) were identified as areas of habitat stability for the species.
Thinking of a worst case scenario, in which the species would lose suitable areas or would change its geographical range by tracking suitable conditions in the future, we proposed new areas where conservation efforts should be directed directed by expanding the currently protected areas as well as considering the creation of dispersal corridors within the areas of habitat stability. For this, we overlapped and compared the spatial distribution of the remaining Cerrado and Caatinga biome’s vegetation with the refugia to identify the areas that (a) have confirmed occurrence records for the species, (b) are near the protected areas, (c) are within the area of habitat stability, and (d) that has native remaining vegetation in adjacencies. By doing this, we identified the most important areas to protect current populations and promote dispersal through new areas within the refugia. The map of remaining vegetation for Cerrado and Caatinga biomes was obtained from the Centro de Sensoriamento Remoto of IBAMA (Instituto Brasileiro do Meio Ambiente e dos Recursos Naturais Renováveis, http://siscom.ibama.gov.br).
Results
The current species distribution predicted by models covers the central and northern region of the Cerrado, the western part of Caatinga, and a smaller region in the south of the Amazonia biomes (Figure 1). Predictions using the LGM revealed a more restricted distribution toward the centre and northwest region of Cerrado in the past in comparison with the present (Figure 1(a) and (b)). For the future, a displacement of suitable areas toward the Caatinga biome (Figure 1(c)) is expected. Moreover, a decrease in the suitable area of more than 60% (702 to 261 cells) was observed for the whole potential distribution between the present and future scenarios (Figure 1(b) and (c)).
Figure 1.
Potential geographic distribution of Micrurus brasiliensis in the past (a), present (b), and future (c), defined as the area with habitat suitability equal to or higher than 0.59 (i.e., the lowest suitability value associated with an occurrence record). Colors from yellow to red represent the habitat suitability varying from 0.59 to 1.00, respectively.
Besides range reduction, these areas will be, on average, less suitable in the future in comparison with present, meanpresent = 0.77, SD = 0.10; meanfuture = 0.68, SD = 0.05; t(875.29) = 19.19, p < .001, regardless whether inside, t(240.9) = 12.82, p < .001, or outside the PAs, t(634.03) =15.439, p < .001. We found no statistical difference in mean suitability when comparing the areas inside (mean = 0.77, SD = 0.11) and outside (mean = 0.78, SD = 0.10) the PAs in the present, t(446.37) = −0.62, p = 0.52. For the future, however, cells inside the PAs are expected to hold lower suitability (mean = 0.65, SD =0.04) than cells outside Pas, mean = 0.68, SD = 0.05; t(114.36) = −4.62, p < .001.
By considering the temporal and geographically continuous distribution of areas with high suitability, we found a potential refugia extending from the southern Caatinga to the northeast of the Cerrado biomes (Figure 2). Within this region, we identified three large areas in the Cerrado that include six populations near to small and fragmented PAs but that still have a high proportion of remaining vegetation in adjacencies. We suggest these areas as a priority for the expansion of the PAs (e.g., by creating dispersal corridors), thus reducing extinction risk by preventing isolation of populations.
Figure 2.
Area of habitat stability (light yellow) through the Cerrado (light gray) and Caatinga (dark gray) biomes; Areas indicated with dash-dotted lines indicate the areas where conservation efforts should be directed considering the presence of recorded populations (red dots); the fragmented and isolated protected areas are demarked with solid lines; and the remaining vegetation is shown in green.
Discussion
Our results clearly indicate a decrease in suitable areas for M. brasiliensis over time as a response to climate change. Our findings showed an expressive range reduction in the future, combined with a spatial displacement toward eastern Brazil and a decrease in the suitability in general. Similar results were also observed for other groups of organisms from the Cerrado biome (e.g., Collevatti et al., 2012, 2013; Terribile et al., 2012). This pattern probably reflects the effects of changes to the precipitation regimes (Schaller, Mahlstein, Cermak, & Knutti, 2011), which are expected to be more drastic mainly in the tropical regions of the world (in contrast to changes in temperature expected for temperate regions). Changes in the precipitation may directly affect M. brasiliensis, because the activity patterns of the genus Micrurus are strongly dependent of the seasonal rainfall regimes (Marques et al., 2006).
During the wet season (spring and summer), the availability of food for coralsnakes is higher, since these snakes feed on amphisbaenids and caecilians, which also live predominantly in the subsoil and emerge to the surface in wet seasons (Almeida et al., 2016). Moreover, chemical signals from prey are more evident in the soil during the rainy season, increasing the foraging activities of coralsnakes (Marques et al., 2006). Also, there is some evidence that oviparous snakes select wet locations for oviposition, since their offspring are more likely to survive in wet areas (Brown & Shine, 2004). Thus, these specialized ecological traits combined with the reduction in the suitable areas due to climate change reported in this study will probably have a negative effect on the activities of M. brasiliensis, and consequently, on their reproductive patterns and survival.
Furthermore, many studies have shown species range contractions in the future (e.g., Lemes & Loyola, 2013; Lemes, Melo, & Loyola, 2014; Loyola, Lemes, Faleiro, Trindade-Filho, & Machado, 2012), which indicate the necessity of identifying and protecting strategic areas to act as buffers under ongoing climate change. Unfortunately, the current network of protected areas is inefficient to safeguard most species from these changes (Araujo et al., 2004; Ferro et al., 2014; Urbina-Cardona & Loyola, 2008), and our study does not show a more optimistic scenario for M. brasiliensis. Furthermore, most of the protected areas are fragmented (e.g., Sobral-Souza, Francini, & Lima-Ribeiro, 2015), as we observed here. The decrease of habitat suitability for M. brasiliensis within PAs in the future enforces the need for identifying complementary areas, preferentially surrounding currently protected ones, which can also serve as dispersal corridors. Most importantly, such areas should be selected considering their stability over time, enhancing the chances of protecting the species while the climate is changing.
Implications for Conservation
Although our proposition for expanding the protected areas may sound quite general at first, it may be very useful as a general approximation on which areas, once connected, could enhance the viability of populations under climate change (Beier & Noss, 1998). Moreover, these remnants of native vegetation are in the core of the species’ stability area, which is a further reason to establish connections since there is very low uncertainty about where this species is probably going to be in the future. Indeed, isolated populations have a higher extinction risk (Diamond, 1975), mainly due to the effects of the demographic stochasticity, metapopulation dynamics, and inbreeding depression (Hanski, Pakkala, Kuussaari, & Lei, 1995; Keller & Waller, 2002; Opdam, 1991). Beyond that, fragmentation could increase the effects of the ongoing climate change (Opdam & Wascher, 2004), which is already causing extinction of snake populations worldwide (Reading et al., 2010). Thus, as fragmented habitats tend to hold fewer (sensitive) specialist species (Brown & Kodric-Brown, 1977; Diamond, 1975), corridors of suitable habitats might allow these species that are usually sedentary to disperse among patches, reducing the extinction risk due to rescue effect (Gilbert, Gonzalez, & Evans-Freke, 1998).
Finally, although our study was focused on a single species, its very particular habitat requirements and conservative ecological traits are shared for most species of Micrurus, which are probably also under severe threat due to climate change. Our approach exemplifies how conservation strategies can be oriented for species with low dispersal abilities, mainly based on preserving refugia in situ, and also reducing the uncertainties of where, all else being equal, conservation efforts should be focussed.
Acknowledgments
We thank “Rede Genética Geográfica e Planejamento Regional para Conservação de Recursos Naturais no Cerrado” (GENPAC), especially GENPAC 4 under the project “Modelos de nicho ecológico, distribuição potencial e o efeito de mudanças climáticas em espécies do Cerrado,” and GENPAC 14, under the project “Filogeografia e diversidade toxinológica de Micrurus (Serpentes, Elapidae) na área de contato Cerrado - Amazônia – Caatinga” (calling MCT/CNPq/FNDCT/FAPs/MEC/CAPES/PRO-CENTRO-OESTE No 031/2010). LCT and NJSJr’s research is supported by Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq) productivity grants (Processes 306418/2013-4 and 309443/2013-0, respectively). The research by CTC was supported by the Institutional Program of Scientific Initiation of UFG and CNPq.
Declaration of Conflicting Interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.
Funding
The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: The author received financial support from CNPq, CAPES and FAPEG for this research.
