Climatically suitable habitats under current and future scenarios for a potentially threatened snake Habitats climaticamente adequados para uma espécie

1 Programa de Pós-Graduação em Zoologia, Departamento de Zoologia, Instituto de Ciências Biológicas. Universidade Federal de Minas Gerais. Av. Antônio Carlos, 6627, Pampulha, 31270-901, Belo Horizonte, MG, Brasil. 2 Departamento de Ecologia, Centro de Biociências. Universidade Federal do Rio Grande do Norte. Campus Lagoa Nova, 59072-970, Natal, RN, Brasil. 3 Departamento de Ciências Biológicas. Instituto Federal Goiano, 75790-000, Urutaí, GO, Brasil. Henrique Caldeira Costa1 ccostah@gmail.com


Introduction
A major issue hampering the complete understanding of biodiversity is the lack of adequate distributional information for most taxa, the Wallacean shortfall (Cardoso et al., 2011;Whittaker et al., 2005).Since spatial distribution is among the criteria used to evaluate the conservation status of a species (e.g.IUCN, 2001), it is likely that lack of proper distributional knowledge on the species range has an unquestionable negative impact on the effectiveness of conservation practices (Bini et al., 2006).Wallacean shortfall is one of the main impediments towards making progress in species conservation (Whittaker et al., 2005), and contribute to the "biodiversity crisis" caused by an unsustainable use of natural resources (Baillie et al., 2010;Dubois, 2003).The current main drivers of such impact on biodiversity are land-use and climate change (Sala, 2000;Tylianakis et al., 2008).The effect of land conversion for agriculture is well recognized as a cause of decline and extinctions globally (Baillie et al., 2010;Mittermeier et al., 2004).Given climate change, species must either disperse to favorable environments or face altered climatic conditions by adjusting their behavior and physiology either through adaptation or natural selection (Parmesan, 2006;Skelly et al., 2007).Declines and extinction may happen to those not able to respond to these new environmental conditions (Sinervo et al., 2010).
Species responses to climate change will not be unimodal (all species will be affected), and multimodal responses (all species will respond differently to this environmental change).Since several species traits may be affected by climate changes, either positively or negatively (Hughes, 2000;Parmesan and Yohe, 2003;Parmesan, 2006), species responses to climate change will definitely vary.Nonetheless, specialist species, showing fine-tuned environmental requests, are expected to be more affected than generalist ones (e.g.Penman et al., 2010;Vasconcelos, 2014;Zank et al., 2014;Silva et al., 2015).
Tropical reptiles are among the taxa expected to be most affected by climate change since they are poor dispersers and generally exhibit lower physiological plasticity (Araújo and Pearson, 2005;Deutsch et al., 2008;Tewksbury et al., 2008;Aubret and Shine, 2010).Recent studies evaluating potential responses of reptiles to climate change have reinforced these predictions of shifts in range and population decline for both rare and common species (Sinervo et al., 2010), including snakes (Mesquita et al., 2013;Penman et al., 2010;Vasconcelos, 2014).
To evaluate what may occur to the distribution range of a taxon facing both climate and land-use changes, researchers often use species distribution models (SDM) (Guisan and Zimmermann, 2000;Guisan and Thuiller, 2005).By correlating the known records of a target species with their climatic data, SDM techniques allow projecting its Grinnellian niche (Peterson, 2001;Soberón, 2007) onto the geographic space, providing some support for practical conservation actions.These methods have been used to achieve a variety of objectives, such as: (i) to define the potential distribution of a given taxa, especially rare, threatened or invasive species (e.g., Martinelli and Moraes, 2013;Rödder et al., 2009;Silva et al., 2014a;Costa et al., 2015); (ii) to indicate suitable areas for future sampling (e.g.Marini et al., 2010); (iii) to test biogeographical and evolutionary hypotheses (Carnaval and Moritz, 2008;Stephens and Wiens, 2009;Porto et al., 2013;Silva et al., 2014b); (iv) to suggest the establishment of new conservation units (Diniz-Filho et al., 2004;Loyola et al., 2008;Nóbrega and De Marco Jr., 2011);and (v) to determine how species may respond to different climate change scenarios (Costa et al., 2012;Sinervo et al., 2010).
In this context, we used SDMs to assess how future scenarios of climate change may affect habitat suitability for Drymoluber brazili G , 1918 (Serpentes: Colubridae), a species occurring in undisturbed open areas mainly along the Brazilian Cerrado, with a few records in areas of Caatinga and Atlantic Forest (Costa et al., 2013).A previous study suggested that it may be the most threatened snake in central Brazil, due to a combination of a small number of offspring, diet and habitat specialization, terrestrial habit, rarity and sensitivity to habitat change (França and Araújo, 2006).However, it was not included in any species red list at national or regional level yet.Given the lack of a more refined large-scale data on the populations of D. brazili, rather than its geographic occurrences, the use of SDMs will be useful as a first assessment of the potential responses of the species to future climate changes.
In this paper, we used SDMs in order to answer the following questions: (i) Which areas are environmentally suitable for the occurrence of D. brazili, and where should new populations be sought?(ii) Considering current landuse changes, do areas with historical records of D. brazili still contain suitable conditions for its occurrence?(iii) How climate change in following decades may affect the distribution of this species, and which areas should receive priority for the conservation of D. brazili under future climate scenarios?
Considering the known records (Costa et al., 2013), we hypothesize based only on climatic variables that D. brazili may be able to occur in most Caatinga and Cerrado areas, especially in larger natural remnants.However, because of the supposed sensitivity of this species to habitat change (França and Araújo, 2006) and the small percentage of extant Cerrado areas (Mittermeier et al., 2004), localities with oldest records of D. brazili may no longer be suitable for its occurrence.Finally, based on studies with other Cerrado species (e.g., Marini et al., 2009;Siqueira and Peterson, 2003;Vasconcelos, 2014), we expected a future decrease on environmentally suitable areas for D. brazili.

Occurrences of Drymoluber brazili
We gathered a total of 105 geospatial localities for Drymoluber brazili from specimens housed in herpetological collections (see Costa et al., 2013); literature records (Costa et al., 2013, additional records in Supplementary Material); online datasets from CRIA's (Centro de Referência em Informação Ambiental) Species Link (http://splink.cria.org.br;Supplementary Material).For new records (i.e., those not cited by Costa et al., 2013) missing exact geographical information, we used the Ornithological Gazetteer of Brazil (Paynter and Taylor, 1991), consulted researchers, and used Google Earth (Google Inc., 2015) to obtain proximate spatial references.For the occurrences lacking exact geographic coordinates, we used the coordinates of city center where samplings took place to obtain the georeference information.We discarded dubious occurrences or those for which we did not obtained reliable geographic information.From the initial 105 raw occurrences, considering the grid cell resolution we used in the modeling procedures (2.5 arc-min ≈ grid-cells with sizes of 0.041º in the tropics), 80 unique occurrences remained (Figure 1).

Environmental data, principal component variables, and modeling procedures
Considering the 19 bioclimatic variables available from WorldClim, (http://www.worldclim.org/current,Hijmans et al., 2005), we reduced variables collinearity by using a Principal Components Analysis (PCA) and deriving 19 new orthogonal variables (principal components; PCs hereon).We used the first seven PCs, which accounted for 98% of all original climatic variation (Table S1) as our final environmental layers.We obtained the same seven PCs for both pessimistic (A2a) and optimistic scenarios (B2a) modeled for the year 2080 from three different Atmosphere-Ocean General Circulation Models (AOGCMs; CCCMA-CGM2, CSIRO-MK2.0,and UKMO-HADCM3) from CIAT (http://ccafs-climate.org).We projected the linear combinations obtained for the 19 PCs obtained for the present conditions into both of the optimistic and pessimist future scenarios from each one of the AOGCMs we used.Then, we used the first seven PCs in our modeling procedures for the predictions of D. brazili future distributions.The use of PCs (or any other collinearity reduction method) is advised to avoid model overfitting that may result in biologically unreliable species potential distributions (Jiménez-Valverde et al., 2011;Serra et al., 2012;Silva et al., 2014a).
We divided the occurrences of D. brazili into 50 subsets of 70%-30% training-testing occurrences to, respectively, produce its potential distributions and evaluate them, considering both current and future climatic scenarios.Since D. brazili is a potential target species for future practical conservation actions (França and Araújo, 2006), we considered the ROC threshold, which balances for both omission and commission errors (Liu et al., 2005;Varela et al., 2014), to cut the modeled suitability matrices into presence/absence maps.Comparative studies have shown ROC threshold assures a higher prediction rates than others (Barbet-Massin et al., 2012;Liu et al., 2005).We used True Skilled Statistics (TSS hereafter; Allouche et al., 2006), a threshold-dependent statistics which varies from -1 to +1, to assess models' performance.According to this metric, negative or near to zero TSS values represent distributions no better than random, while values near +1 represent perfect agreement between the observed and the modeled spe-cies' distribution in current climatic scenarios (Allouche et al., 2006).Models with TSS values reaching at least 0.5 are considered as reliable and good representatives of the distribution range of the modeled species (Gallien et al., 2012).During the evaluation process, we allocated the pseudo-absences in the geographic space at random.A total of 1,400 presence/absence maps were produced for D. brazili, either considering its distribution in current (n=200) and future scenarios (n=1200).We made a strict consensus (sensu Marmion et al., 2009) of all 50 potential distributions obtained with each modeling algorithm and each scenario to represent the species distribution under that algorithm vs. climate scenario combination.For the species final distribution, we summed each strict consensus obtained for each combination of modeling algorithms and climatic scenario to obtain the final prediction for the species in all seven considered climate change scenarios.Later, considering the suitable areas in both present and future scenarios, we used the predictions from both A2a and B2a scenarios from all AOGCMs to highlight stable areas for the occurrence of D. brazili that could be important for its conservation.Finally, we overlaid the map of the Brazilian strictly protected areas obtained from the Brazilian Ministry of Environment (http://mapas.mma.gov.br/i3geo/datadownload. htm) onto those stable areas to discuss the protection status of D. brazili.We used Factorial ANOVAs to evaluate the range size differences for D. brazili among the different algorithms and each climatic scenario.In these analyses, both TSS values and the predicted areas for the species in each climatic scenarios were our response variables, while our modeling algorithms were the predictor variables.In case D. brazili was to be affected by the predicted future climate changes for South America, we expected its distribution to decrease in size or, eventually, that it would move southwards in order to pursuit the best environmental conditions for its development.A summary of all methods employed in this study is depicted in Figure 2.

Results
In general, all algorithms we used reached TSS values higher than 0.5 (MAX: 0.693+0.067;mean+standard deviation; SVM: 0.682+0.070;MAHAL: 0.667+0.071;EN-VSC: 0.682+0.067).Table 1 shows the mean TSS values obtained with each one of the 50 predicted distributions for D. brazili in each climate change scenario.The predictions for the present scenario were similar in all algorithms (Figure 3A).For the future scenarios, despite the fact the algorithms showed predictions with different sizes (Figure 3A), all of them predicted increases for the distribution of D. brazili when compared to the prediction for the present scenario (Figure 3A and B; Figure 4).Specifically, the predicted distribution for D. brazili under the CCCMA-B2a climatic scenario showed the highest increase in range size (Figure 3B).Predicted distribution ranges considering all modeling algorithms and current and future climatic scenarios are depicted in Figure 4.
Under the ROC threshold, the areas where D. brazili is predicted to occur in the present scenario are mainly found in southeastern and central Brazil (Figure 4), in ecotones between Cerrado and seasonal semideciduous forests where most of the confirmed records of this species came from.In the future scenarios, models show a general tendency of increase and, at the same time, fragmentation of suitable areas for D. brazili.Most algorithms show increased suitable areas in northeastern Brazil and towards the central-western South America, more specifically in Bolivia and central portion of the Andes (Figure 4).
Nonetheless, when considering the consensus maps of distribution models, we can notice that the strictly pro-tected areas correspond to just a small fraction of the areas predicted as suitable for D. brazili, in both current (Figure 5A) and future (Figure 5B-C) climatic scenarios.Additionally, when we overlap climatic suitable areas from all scenarios, only a few areas are predicted as stable (the darkest shades of red in Figure 5D).These stable areas basically consist of two main patches, one located in the northwest of the state of São Paulo and other in the east and northeast of the state of Minas Gerais.Other smaller patches are located between the bigger ones, along the states of Minas Gerais, Goiás, and Bahia.

Discussion
In this study, we modeled the potential distribution of D. brazili, a snake species known mainly from the Brazilian Cerrado biodiversity hotspot (Marris, 2005;  et al., 2004).Under current environmental conditions, D. brazili is expected to occur along the southeastern Cerrado and central portions of the Atlantic Forest.Known records from the Caatinga (probably linked to relictual savanna enclaves (Guedes et al., 2014)) and from the Paraguayan Cerrado were not predicted by the models as presence for the species, an indicative that these areas are sink populations for D. brazili (Pulliam, 2000).Under future climate scenarios, environmentally suitable areas will increase towards the northern Chaco, northwestern Cerrado and the Caatinga of northeastern Brazil.

Mittermeier
Other studies have used SDMs to predict the effects of future scenarios on different organisms inhabiting the Cerrado (e.g., Marini et al., 2009;Siqueira and Peterson, 2003;Vasconcelos, 2014).In general, these studies have found that environmentally suitable areas for sampled plants, birds, frogs and snakes are expected to decrease over time.Even widely distributed species showed potential future reduction (Mesquita et al., 2013).Contrary to those findings, our models predicted an increase in suitable areas for D. brazili under future climate scenarios.
Increase in suitable areas in the future have been found in models for most European reptiles (Araújo et al., 2006), but few studies have found this trend in South America (e.g., Fontanella et al., 2012).The future suitable areas predicted for D. brazili include most of the Cerrado, the Atlantic Forest of southeastern Brazil, the Caatinga of northeastern Brazil, and even western South American regions.This result is in accordance with predictive fu- ture expansions of the Cerrado, due to a 'savannization' of South American forests (Salazar et al., 2007) and a significant warming in tropical Andes (Urrutia and Vuille, 2009).
It is essential to investigate the dispersal abilities of D. brazili and the possible factors limiting its range (e.g., low adaptability to altered habitats) to better understand how this species could actually behave under the emergence of new suitable areas (Araújo and Pearson, 2005;Araújo et al., 2006).Additionally, it is important to note that the models present here are based solely on climatic variables, yet both biotic and historical processes determine the range of a species (Araújo et al., 2006;Soberón, 2007).There is little ecological data available for D. brazili, but it seems enough to suggest that this species is more sensitive to habitat change (França and Araújo, 2006) than to possible climate changes, given the predicted range increase for this species.Considering the rapid conversion of most of the Cerrado into agricultural land (Fearnside, 2002;Phalan et al., 2013;Zachos and Habel, 2011), this species could actually face a decrease in potential areas for its occurrence, despite an increase in climatic suitability.
Since many localities bearing historical records of D. brazili are now human-altered, the probability of the species occurrence in some regions has decreased, although climate still offers suitable conditions for its presence.One such example is the west of the state of São Paulo, stable for the occurrence of D. brazili under present and future climatic scenarios, but currently converted mainly into pastures and sugarcane plantations (Durigan et al., 2007).This high level of conversion of natural areas for livestock and agriculture may be the cause of absence of D. brazili records in São Paulo in the last 20 years.Additionally, the species is unknown from protected areas in that state, even from those that harbor preserved Cerrado remnants (Araújo and Almeida-Santos, 2011;Araújo et al., 2010;Sawaya et al., 2008).This situation led us to recommend D. brazili under some category of extinction risk in the state of São Paulo red list, contradicting Marques et al. (2009), who considered this species as 'least concern'.
Identifying areas bearing viable populations for D. brazili would be of great importance for conservation efforts.However, studies on population dynamics of snakes are nearly absent in South America (Guimarães et al., 2014), contrary to other regions such as the USA, Europe, and Australia (Shine and Bonnet, 2009).Since basic biology, ecology, and known occurrences are the primary data needed to these kind of analyses, we believe such knowledge shortfalls must be filled out, in order to improve conservation efforts and allow the results of SDMs to actively support practical conservation actions (Cardoso et al., 2011).
Available data (i.e.França and Araújo, 2006) suggest D. brazili is restricted to undisturbed open areas, especially Cerrado savannas.Costa et al. (2013) hypothesized an expansion in the distribution range of D. brazili in southeastern Brazil, especially in the lower Doce river region, in periods of Atlantic Forest shrinkage (21-6 thousand years ago [ka], see Carnaval and Moritz, 2008).Although the Cerrado has experienced an overall reduction of its area during the last interglacial (120 ka) and the last glacial maximum (21 ka), it probably has expanded along the Doce river basin between 21-6 ka (Werneck et al., 2012).This may give support to the hypothesis of Costa et al. (2013), but how the distribution of D. brazili has changed during past climatic events remains an open question for future analyses.

Conclusions
Our study provides a useful first approximation of the impact of ongoing climate change on the snake species Drymoluber brazili.Climatically suitable habitats for D. brazili will likely expand in the next decades.However, not only the climatic characteristics, but also the degree of habitat preservation must be taken into account to predict the occurrence of this species more accurately.We recommend that future sampling efforts should focus on putatively stable regions, especially those that match protected areas, like the Grande Sertão Veredas National Park, Mata Escura Biological Reserve, Grão Mogol and Rio Pardo State Parks, all in northern Minas Gerais.Our results suggest a good opportunity to reevaluate the conservation status of D. brazili, especially if new data on biological requirements are added to the interpretation of SDM predictions.

Figure 1 .
Figure 1.Spatial distribution of records of the colubrid snake Drymoluber brazili used in this study.The gray scale corresponds to the mean altitude in each of the grid cells represented in South America.Lighter areas have lower altitudes, while darker areas have higher altitudes.The grid cell size used in this study was equal to 0.041º degrees or 2.5 arc-min, what corresponds to ~4km in the tropics.See Costa et al. (2013) and the supplementary material for details.

Figure 2 .
Figure 2. Conceptual summary of all methods employed in this study to model climatically suitable areas for the occurrence of the colubrid snake Drymoluber brazili under current and future scenarios.

Figure 3 .
Figure 3. Results of model predictions for D. brazili under the different scenarios of climate change.(A) Range size of D. brazili considering environmentally suitable areas according to each different evaluated scenario and algorithm.(B) Range size of D. brazili considering only the effect of each climatic scenario.A2a refers to the pessimist scenario of each climate scenario, whereas B2a refers to the optimist scenarios in each of these scenarios as well.CC, CS, and HD refer to, respectively, CCCMA, CSIRO, and HAD climatic scenarios.Please refer to the methods section for a better understanding of the modeling procedures employed in the study.

Figure 4 .
Figure 4. Predicted environmentally suitable areas for the occurrence of Drymoluber brazili using four algorithms under current and future (year 2080) climatic scenarios.Future scenarios were based on a pessimistic (A2a) and an optimistic (B2a) emission projection for three different Atmosphere-Ocean General Circulation Models (CCCMA, CSIRO, and HAD).