Third-party data from ecoClimate was used for this study and are accessible at: <http://ecoclimate.org/downloads/>. Those interested would be able to access the data in the same manner as the authors. The authors had no special access privileges to this data.

Introduction {#sec001}
============

An accurate knowledge of species geographic distributions is critical to preserve biodiversity in a changing world. Despite the increasing availability of digitized information on biodiversity data and species occurrence, our knowledge about the geographic distribution of most species is regretfully incomplete \[[@pone.0205164.ref001]--[@pone.0205164.ref003]\]. This shortcoming, called the Wallacean shortfall, is more evident in highly diverse tropical ecosystems \[[@pone.0205164.ref004]--[@pone.0205164.ref006]\], which are paradoxically the most threatened ecosystems \[[@pone.0205164.ref007],[@pone.0205164.ref008]\]. Moreover, the lack of data on species distribution is frequently biased toward small, inconspicuous, and not easily detectable species \[[@pone.0205164.ref001],[@pone.0205164.ref009]\]. Even with the wide use of methodological tools to predict species ranges, large gaps still exist for these species, particularly because the primary data used in predictive models are scarce or incomplete, precluding a clear understanding of the threats from climate change and assessments of their conservation status.

Ongoing climate changes have already impacted organisms' distributions around the world in recent years \[[@pone.0205164.ref010]--[@pone.0205164.ref013]\]. Snakes are especially sensitive to climate changes \[[@pone.0205164.ref014]--[@pone.0205164.ref016]\] because most species maintain daily activities in a restricted range of temperatures \[[@pone.0205164.ref017]\] and because the reproductive strategies in tropical species depend on rainfall patterns and seasonality in precipitation \[[@pone.0205164.ref017]--[@pone.0205164.ref020]\]. The most important issue regarding the impact of climate change is that species that are unable to evolve rapid physiological adaptations or climatic tolerances (to new thermal conditions) may become extinct unless they have good dispersal abilities to track suitable habitats \[[@pone.0205164.ref021]--[@pone.0205164.ref024]\]. These impacts are particularly harmful for those species that are habitat specialists, which would probably fail to migrate following suitable climates \[[@pone.0205164.ref023],[@pone.0205164.ref025]\]. Despite this situation, with few exceptions, snakes have been neglected from studies assessing the impacts of climate change \[[@pone.0205164.ref024],[@pone.0205164.ref026],[@pone.0205164.ref027]\].

*Micrurus lemniscatus* Linnaeus 1758 comprises a taxonomic complex of six taxa: *M*. *l*. *lemniscatus*, *M*. *l*. *carvalhoi*, *M*. *diutius*, *M*. *l*. *helleri*, *M*. *potyguara* and *M*. *serranus*, which occur all across South America, except Chile and Uruguay. Despite the wide geographic distribution (encompassing such distinct biomes as the Amazon, Cerrado, Caatinga and Atlantic Forest) and medical importance due to their highly neurotoxic venoms, little is known about their ecology, taxonomic delimitation, or geographic distribution \[[@pone.0205164.ref028]--[@pone.0205164.ref032]\]. Some studies suggested that the *M*. *lemniscatus* complex is composed of two geographically distinct clades: one from Amazonia and another from the Caatinga and Atlantic forest biomes \[[@pone.0205164.ref029]\]. However, the limited number of individuals analyzed precluded a clear geographic delimitation of each form. Extensive compilations such as by Campbell & Lamar \[[@pone.0205164.ref033]\]reproduced these shortcomings and presented poor definitions of the species' geographic limits, which were based on previous studies without covering the entire scope of information from the specimens available in scientific collections. Thus, this group clearly exemplifies the necessity for new studies dealing with the Wallacean shortfall, as the knowledge about species distributions is very incomplete. Most importantly, similar to the observations for other coralsnakes \[[@pone.0205164.ref030],[@pone.0205164.ref034],[@pone.0205164.ref035]\], the high levels of habitat degradation in the biomes where they occur may represent a severe threat for their long-term survival. This degradation raises the point that we may be losing these species even before they are fully known.

In this paper, we address the Wallacean shortfalls for the *M*. *lemniscatus* complex by refining our knowledge about the geographic distribution of species based on ecological niche modeling tools and a set of thoroughly examined occurrence records from scientific collections. From the present ranges, we project the distributions into future scenarios of climate change and address the potential impacts for the long-term conservation of these species. Finally, we evaluate the efficiency of the current reserve network in maintaining suitable areas for the species in the present and future. Our findings add valuable new information to the geographic distribution of each taxon, pointing to the necessity of reconsidering the conservation efforts for these less conspicuous species at risk of extinction in the near future.

Materials and methods {#sec002}
=====================

Species and climate data {#sec003}
------------------------

Occurrence records were obtained from the scrutiny of all specimens belonging to the *M*. *lemniscatus* complex deposited in 33 museum collections in South America, the United States of America and Europe (see [S1 Table](#pone.0205164.s003){ref-type="supplementary-material"} and [S1 Appendix](#pone.0205164.s001){ref-type="supplementary-material"} for the list of specimens and the museum collections accessed). The specimens were carefully examined to confirm taxonomic identifications and avoid nomenclatural and georeferencing errors. A total of 768 occurrence records were mapped into a geographic grid encompassing South America, with cell boundaries following 0.5 × 0.5 degrees of latitude and longitude (the resolution of the climate data used for niche modeling, see below). Repeated records within a grid cell were excluded, thus reducing the spatial aggregation of occurrences and avoiding spatial autocorrelation effects in ecological niche modeling. The final number of records (unique records in each grid cell) was 75 for *M*. *l*. *lemniscatus*, 97 for *M*. *l*. *carvalhoi*, 33 for *M*. *diutius*, 49 for *M*. *l*. *helleri*, and three for *M*. *potyguara* ([S1 Fig](#pone.0205164.s009){ref-type="supplementary-material"} and [S2 Table](#pone.0205164.s004){ref-type="supplementary-material"}). Because there are problems in modeling data with few observations, we excluded *M*. *potyguara* from the analyses. We also did not use the data related to *M*. *serranus* owing to its unstable taxonomic position within the *M*. *lemniscatus* complex.

Because a clear spatial aggregation of records still remained for *M*. *l*. *carvalhoi*, probably representing a sampling bias in a highly populated region ([S1 Fig](#pone.0205164.s009){ref-type="supplementary-material"}), we applied the protocol proposed by Oliveira et al. \[[@pone.0205164.ref036]\] to reduce spatial aggregation and select geographically equidistant points for this species. We then compare the results obtained from using all records at 0.5°x0.5° resolution with the results after controlling for spatial aggregation. First, we calculated the Mahalanobis distances (D^2^) in environmental space, which was formed by the five bioclimatic variables used to build ecological niche models (ENMs, see below), between each one occurrence records and their centroid. Second, with the distances D^2^ and the geographical coordinates (i.e., latitude and longitude values) we fitted a simultaneous autoregressive (SAR) model to find the autoregressive coefficient *p* which measures the amount of spatial autocorrelation in data. Third, this coefficient *p* was used to calculate the effective degrees of freedom using the formula provided by Griffith \[[@pone.0205164.ref037]\]. Thus, the degrees of freedom were the effective number of independent records (presences and pseudo-absences). Finally, with the number of independent records (i.e., 60 points) we select in the environmental space the most equidistant ones using an algorithm that iteratively searches records which are most distant from each other. Since no considerable differences were observed between the results using all records (i.e., 97 points) and those after selecting equidistant records (i.e. 60 points) ([S2 Fig](#pone.0205164.s010){ref-type="supplementary-material"} and [S2 Appendix](#pone.0205164.s002){ref-type="supplementary-material"}), we presented here only the results using all records at 0.5°x0.5° resolution.

A set of bioclimatic variables at 0.5° resolution were obtained from the ecoClimate database ([http://ecoclimate.org](http://ecoclimate.org/), \[[@pone.0205164.ref038]\]). EcoClimate provides updated climatic simulations for several time periods (past, present, and future) derived from the Coupled Model Intercomparison Project--Phase5 (CMIP5) and Paleoclimate Modeling Intercomparison Project--Phase3 (PMIP3) \[[@pone.0205164.ref039],[@pone.0205164.ref040]\]. Given the variety of climatic simulations (as well as ecological niche models (ENMs hereafter), see below) currently available, we combined the output predictions from different climate models and ENMs following the ensemble approach of Araújo & New \[[@pone.0205164.ref041]\] (see details below). Thus, we used five coupled atmosphere-ocean general circulation models (AOGCMs: CCSM, CNRM, GISS, MIROC, and MRI) ([S3 Table](#pone.0205164.s005){ref-type="supplementary-material"}), with simulations for the present (represented by preindustrial data) and future (represented by mean values between 2080--2100), to derive 19 bioclimatic variables according to Hijmans et al. \[[@pone.0205164.ref042]\]. For the future, we used the emission scenario RCP4.5, which is an intermediate scenario between the lower (RCP2.6) and the higher (RCP8.0) emission scenarios \[[@pone.0205164.ref040]\].

To reduce collinearity between the bioclimatic variables in the ENMs, we applied a varimax-rotated factor analysis to the correlation matrix between pairs of variables and selected the variable with the highest loading in each one of the first five rotated factors ([S4 Table](#pone.0205164.s006){ref-type="supplementary-material"}). The selected variables were the mean annual temperature, annual temperature range, precipitation of the wettest month, precipitation of the driest month, and precipitation of the warmest quarter. Because the bioclimatic variables were highly correlated across the AOGCMs, we performed the factor analysis only with the AOGCM CCSM and then applied the selection for all other AOGCMs.

Ecological niche models and geographic distribution {#sec004}
---------------------------------------------------

The variety of alternative methods to model species niches have promoted a range of discussions about the individual performance of each method to estimate potential geographic distributions \[[@pone.0205164.ref043]--[@pone.0205164.ref045]\], although no agreement was reached regarding which is the best choice for particular aims. We opted to apply the ensemble forecasting approach \[[@pone.0205164.ref041],[@pone.0205164.ref046],[@pone.0205164.ref047]\] to generate consensus predictions about geographic distributions after combining the outputs from several different methods.

We used 12 ENMs, including presence-only, presence-background and presence-absence methods, ranging from simple bioclimatic envelope models (e.g., BIOCLIM) and distance-based methods (e.g., Euclidian distance) up to more complex regression models such as (GLM) and machine learning-based methods (artificial neural networks) \[[@pone.0205164.ref048],[@pone.0205164.ref049]\] ([S5 Table](#pone.0205164.s007){ref-type="supplementary-material"}). Models were built for the present using preindustrial climate data and projected into climate scenarios for the end of this century (mean simulations between 2080 and 2100).

The area used to adjust and project ENMs must correspond to regions that have been available for colonization to the species over relevant time periods \[[@pone.0205164.ref050]\]. The *Micrurus* genus is predominantly Neotropical, with more than 30 species occurring through the South America (the richest region for this clade) and only one species occurring in extreme south of North America \[[@pone.0205164.ref051]\]. Also, the taxa belonging to the *M*. *lemniscatus* complex are sibling species that were already considered as a single species distributed across most of the South American continent. Thus, we considered South America as the suitable area for model calibration and projection for this group.

For the modeling processes, we randomly divided species presence records (and pseudo-absences for those methods that use absences) into 75% for calibration and 25% for evaluation and repeated this process 50 times. Pseudo-absences were randomly selected in background regions (excluding cells with occurrence records) with the same proportion of species records (a prevalence of 0.5). The 50 repetitions in each method were converted in presence-absence maps according to thresholds established by the area under the ROC curve (AUC). All 600 models (50 repetitions × 12 modeling methods) were included in the consensus maps of each AOGCM, weighted by their model fit according to the true skill statistics--TSS \[[@pone.0205164.ref052]\] ([S6 Table](#pone.0205164.s008){ref-type="supplementary-material"}). The average across the ensemble outputs of each AOGCM resulted in the final consensus map of habitat suitability (varying from 0 to 1) for each species ([S1 Fig](#pone.0205164.s009){ref-type="supplementary-material"}). All models were generated using the computational platform BioEnsembles \[[@pone.0205164.ref046],[@pone.0205164.ref047]\].

The agreement among projected distributions from the different ENMs and AOGCMs was assessed through a hierarchical ANOVA \[[@pone.0205164.ref046],[@pone.0205164.ref047]\]. For this, the suitability of each cell from each combination of ENM and AOGCM was the response variable replicated (or nested) within the "time" component (present and future projections). We then used the sum of squares (SS) that can be attributed to each of these sources of variation (i.e., among ENMs and AOGCMs, both nested in time, and between the two time periods) to evaluate the uncertainties in species distribution. Low values of SS from ENMs and AOGCMs indicate high consistency in predictions from different methods and climate models. High values of SS from the time component indicate changes in suitability from the present to the future.

We used habitat suitability predictions to map the potential geographic distribution of each species in both time periods (present and future). Because ENMs do not consider biotic interactions and other local factors limiting species ranges (e.g., geographic barriers), their transferability (extrapolation in geographic space) may include areas that are inaccessible for the species. As we are dealing with poor dispersal species, we should be conservative to delineate more realistic and interpretable distributions. Thus, we applied the decision threshold based on the lowest presence threshold (LPT) \[[@pone.0205164.ref053]\], selecting the lower suitability value associated with the observed occurrence records. Following this, an area (or cell) was considered suitable if it had a certain value of suitability higher than the threshold. The final maps were compared with the expert opinions (N.J.S.Jr. and D.T.) to ensure that the predictions reflect the real potential distribution of each species.

Finally, we explored the relative importance of each bioclimatic variable used in niche models in determining the habitat suitability for each species by adjusting a multiple regression model between suitability and the bioclimatic variables. We then interpreted the standardized regression coefficients to set the contribution of each variable. Although the significance of these coefficients cannot be considered due to the implicit collinearity between suitability and the bioclimatic variables, this approach was useful for describing their relative contributions to the delimitation of species' geographic distributions and to infer how climate changes could impact species in the future.

Climate change and species representation in protected areas {#sec005}
------------------------------------------------------------

We evaluated the potential impacts of future climate change by comparing the number of currently occupied cells in relation to future distributions. We also performed a gap analysis \[[@pone.0205164.ref054]--[@pone.0205164.ref057]\] by measuring the level of representation of species distributions into the protected areas according to habitat suitability for the present and future. To delineate the geographic extension of this analysis for each species, we overlapped the predicted geographic distributions (after cutting by specific LPT thresholds) with the South American ecoregions obtained from The Nature Conservancy GIS database (<http://maps.tnc.org/about.html>). We then selected the ecoregions that intersect with the species distributions using QGIS v. 2.12 (<https://www.qgis.org/pt_BR/site/>) ([S3 Fig](#pone.0205164.s011){ref-type="supplementary-material"}) and calculated the proportion of protected areas in each cell grid overlapping the selected ecoregions. The spatial distribution of the protected areas was obtained from the World Database on Protected Areas of the IUCN (available at: <http://www.protectedplanet.net/search>).

From the proportion of each cell that is currently protected and the continuous values of habitat suitability, we applied the species representation index *SRI* proposed by Alagador et al. \[[@pone.0205164.ref057]\]to evaluate how the species suitability is distributed in relation to the protected areas in the present and how it would be in the future by considering climate change. We also assessed the degree to which reserve networks represent species better than expected by chance alone. For this assessment, we randomly selected through 1,000 permutations an equal number of cells as the number of protected grid cells and redistributed the suitability values within the grid cells. We then compared the *SRI* value obtained for each species using the real reserve network with those obtained using the corresponding random set of reserves. We applied two-tailed tests by counting the results in which the *SRI* value was lower or higher than the real *SRI*. In each case, we checked whether the real *SRI* was greater than the 95^th^ or lower than 5^th^ percentile of the *SRI* values obtained from the randomizations. The *SRI* index has the advantage of using continuous data regarding species distributions (habitat suitability) and reserve coverage (proportion of the reserve area by grid cell), thus reducing the uncertainty in choosing arbitrary thresholds to convert continuous suitability into species presence and absence and to convert a grid cell into being considered either protected or unprotected \[[@pone.0205164.ref057]\].

Results {#sec006}
=======

Niche models and species geographic distributions {#sec007}
-------------------------------------------------

The ecological niche models exhibited good predictive performance ([S6 Table](#pone.0205164.s008){ref-type="supplementary-material"}), with mean TSS values across the AOGCMs and ENMs higher than 0.5 for all species except *M*. *l*. *helleri* (TSS = 0.45). The current distribution of each species after applying the LPT threshold ([Table 1](#pone.0205164.t001){ref-type="table"}) was compatible with the known distribution from the observation records ([Fig 1](#pone.0205164.g001){ref-type="fig"}) and the expert opinions, except for *M*. *l*. *carvalhoi* (see below). Analyses of the uncertainties indicated that the climate models represented the lowest source of variation among the predictions ([Table 2](#pone.0205164.t002){ref-type="table"}). Niche models were the second highest source of variation for two subspecies (*M*. *l*. *lemniscatus* and *M*. *l*. *helleri*), and the differences through present and future predictions were higher than those for the other two taxa (*M*. *l*. *carvalhoi* and *M*. *diutius*) ([S4 Fig](#pone.0205164.s012){ref-type="supplementary-material"}).

![Habitat suitability and current geographic distribution.\
Consensus patterns of habitat suitability and potential present geographic distribution of the *Micrurus lemniscatus* species complex after applying specific decision thresholds (the lowest presence threshold). Black dots indicate occurrence records used in the models; A) *Micrurus l*. *lemniscatus*, B) *Micrurus l*. *carvalhoi*, C) *Micrurus diutius*, D) *Micrurus l*. *helleri*.](pone.0205164.g001){#pone.0205164.g001}

10.1371/journal.pone.0205164.t001

###### Changes in geographic distribution measured as the difference in the number of 0.5 × 0.5° cells in which species are predicted to occur in the present and in the future.

LPT: Lowest presence threshold used to cut the potential distributions.

![](pone.0205164.t001){#pone.0205164.t001g}

  Species                   [\*](#t001fn001){ref-type="table-fn"}Range size~P~   [\*](#t001fn001){ref-type="table-fn"}Range size~F~   \% of loss/gain   LTP
  ------------------------- ---------------------------------------------------- ---------------------------------------------------- ----------------- ------
  *M*. *l*. *lemniscatus*   1958                                                 2364                                                 20.7              0.36
  *M*. *l*. *carvalhoi*     2292                                                 1536                                                 -32.9             0.23
  *M*. *diutius*            616                                                  82                                                   -86.6             0.72
  *M*. *l*. *helleri*       1066                                                 600                                                  -43.7             0.63

\* Number of cells; P--present; F--future.

10.1371/journal.pone.0205164.t002

###### Summary of the hierarchical ANOVA output where species environmental suitability was used as response variable, and atmosphere-ocean general circulation models (AOGCM) and ecological niche models (ENM) were used as categorical variables nested into two time periods (the present and the end of the century).

The values represent the mean, maximum and minimum percentages of the sum of squares through all cells.

![](pone.0205164.t002){#pone.0205164.t002g}

                       *M*. *l*. *lemniscatus*   *M*. *l*. *carvalhoi*   *M*. *diutius*   *M*. *l*. *helleri*                                  
  -------------------- ------------------------- ----------------------- ---------------- --------------------- ------- -------------- ------- -------------
  **AOGCM \[Time\]**   5.75                      0.00--48.84             4.80             0.00--59.80           3.49    0.00--60.43    2.68    0.00--31.85
  **ENM \[Time\]**     25.29                     1.22--84.73             20.30            1.55--89.54           25.12   0.80--82.27    34.09   1.67--90.39
  **Time**             23.48                     2.26--79.39             29.19            0.85--79.29           27.69   1.69--80.76    21.59   1.54--67.44
  **Residual**         45.48                     9.56--79.43             45.70            3.26--79.54           43.70   12.48--77.19   41.64   6.39--72.43

The niche models indicated that *Micrurus l*. *lemniscatus* is widely distributed across the Amazon basin in northern Brazil, extreme northern Peru, southern Colombia, in small regions southern and western Venezuela, and across the Guianas ([Fig 1A](#pone.0205164.g001){ref-type="fig"}). In Brazil, the *M*. *l*. *lemniscatus* distribution also extends towards the northern part of the Cerrado and Caatinga biomes. *Micrurus l*. *carvalhoi* also showed a broad distribution, extending from the Atlantic rainforest in southeastern Brazil, through the Cerrado biome towards eastern Paraguay, and Misiones and Corrientes in Argentina, including part of Bolivia and Peru, reaching the southern limits of the Amazon rainforest ([Fig 1B](#pone.0205164.g001){ref-type="fig"}). The distribution of *M*. *l*. *carvalhoi* also includes most of the northeastern region of Brazil. *Micrurus diutius* exhibited a restricted distribution in the rainforest of the Guianas and northern Brazil ([Fig 1C](#pone.0205164.g001){ref-type="fig"}). The predicted distribution for *M*. *l*. *helleri* was congruent with the recognized occurrence in the western region of the Amazon rainforest, including northern Bolivia and Brazil, northern and eastern Peru, eastern Ecuador and central and southern Colombia ([Fig 1D](#pone.0205164.g001){ref-type="fig"}).

Extreme precipitation conditions and precipitation during the warmest season, respectively represented by the variables Bio13 and Bio18, had the highest standardized coefficients for all species except *M*. *diutius* ([Table 3](#pone.0205164.t003){ref-type="table"}). For this species, temperature variation was more important (Bio 7 --annual temperature range).

10.1371/journal.pone.0205164.t003

###### Standardized regression coefficients from the multiple regression models between habitat suitability and the bioclimatic variables.

Bold values indicate the two most important variables (i.e., with the higher coefficient values) for each species.

![](pone.0205164.t003){#pone.0205164.t003g}

  Species                   Bio1       Bio7        Bio13       Bio14   Bio18
  ------------------------- ---------- ----------- ----------- ------- ----------
  *M*. *l*. *lemniscatus*   **0.33**   -0.21       **0.45**    0.14    -0.24
  *M*. *l*. *carvalhoi*     0.27       0.06        **-0.36**   -0.29   **0.47**
  *M*. *diutius*            **0.27**   **-0.47**   0.21        0.12    -0.20
  *M*. *l*. *helleri*       **0.40**   -0.30       -0.16       0.08    **0.41**

Bio1 --Annual Mean Temperature, Bio7 --Temperature Annual Range, Bio13 --Precipitation of Wettest Month, Bio14 --Precipitation of Driest Month, Bio18 --Precipitation of Warmest Quarter.

Climate change and gap analysis {#sec008}
-------------------------------

Except for *M*. *l*. *lemniscatus*, all species were projected to contract their distributions in the future ([Fig 2](#pone.0205164.g002){ref-type="fig"}). *Micrurus diutius* exhibited the most severe contraction ([Fig 2C](#pone.0205164.g002){ref-type="fig"}), with a reduction of 86% in its suitable area ([Table 1](#pone.0205164.t001){ref-type="table"}). Substantial reductions were also projected for *M*. *l*. *carvalhoi* and *M*. *l*. *helleri*, with losses of 32 and 45%, respectively (Fig [2B and 2D](#pone.0205164.g002){ref-type="fig"}).

![Habitat suitability and future geographic distribution.\
Consensus patterns of habitat suitability and potential future geographic distributions of the *Micrurus lemniscatus* complex after applying specific decision thresholds (the lowest presence threshold); A) *Micrurus l*. *lemniscatus*, B) *Micrurus l*. *carvalhoi*, C) *Micrurus diutius*, D) *Micrurus l*. *helleri*.](pone.0205164.g002){#pone.0205164.g002}

The *SRI* estimates based on continuous values of habitat suitability and the proportion of reserve coverage in the grid cells indicated that for all species, less than 50% of the present suitable areas are protected in reserve networks ([Table 4](#pone.0205164.t004){ref-type="table"}). Even worse, all species will experience future suitability loss within the reserves, as less than 40% of the suitable areas are expected to be held in reserves at the end of the century. *M*. *l*. *carvalhoi* was the worst represented species in the protected areas (14% for the present and 13% for the future). The permutations ([Fig 3](#pone.0205164.g003){ref-type="fig"}) showed that for two subspecies (*M*. *l*. *lemniscatus* and *M*. *l*. *helleri*, [Fig 3A, 3B, 3G and 3H](#pone.0205164.g003){ref-type="fig"}), the current spatial configuration of the reserve network holds more suitable areas than those expected by their random distribution. However, for *M*. *l*. *carvalhoi* ([Fig 3C and 3D](#pone.0205164.g003){ref-type="fig"}), the network of reserves performed worse than expected to change, both for the present and future distributions. This poor performance means that the current distribution of the protected areas is insufficient to preserve *M*. *l*. *carvalhoi*, both at present and in the face future of climate changes. A similar result was observed for *M*. *diutius*, for which the loss of representation will result in a significant inadequacy of the reserves in the future ([Fig 3F](#pone.0205164.g003){ref-type="fig"}).

![Species representation index (*SRI*).\
Frequency distributions of 1000 permutations of species representation index by randomly selecting protected areas within the ecoregions encompassed by each species distribution at present and in the future. The arrow indicates the position of the real *SRI* value relative to the distribution of randomized values for each species; significance at a two-tailed test is also indicated as *p*-values; A and B) *Micrurus l*. *lemniscatus*, C and D) *Micrurus l*. *carvalhoi*, E and F) *Micrurus diutius*, G and H) *Micrurus l*. *helleri*.](pone.0205164.g003){#pone.0205164.g003}

10.1371/journal.pone.0205164.t004

###### Species representation index (*SRI*) and *p*-values from randomizations.

![](pone.0205164.t004){#pone.0205164.t004g}

  Species                       *SRI*~present~   *SRI*~future~   *p*~present~   *p*~future~
  ----------------------------- ---------------- --------------- -------------- -------------
  ***M*.** *l*. *lemniscatus*   0.39             0.40            \< 0.05        \< 0.05
  ***M***. *l*. *carvalhoi*     0.14             0.13            \< 0.05        \< 0.05
  ***M***. *diutius*            0.47             0.32            \< 0.05        \< 0.05
  ***M***. *l*. *helleri*       0.38             0.37            \< 0.05        \< 0.05

Discussion {#sec009}
==========

Reducing the Wallacean shortfalls {#sec010}
---------------------------------

Wallacean shortfall refers to our lack of knowledge of species distributions, both because we do not have adequate sampling efforts across multiple regions and lack straightforward methods to generalize from these data to obtain a clear picture of species distributions \[[@pone.0205164.ref001],[@pone.0205164.ref002]\]. Data compilations and syntheses from biological collections, coupled with ENMs, have been widely used in the last 20 years as a way to reduce Wallacean shortfalls. In this sense, if the knowledge of biodiversity improves in the data and theory/models, it is possible to better understand the potential impacts of climate changes on species\' geographic ranges \[[@pone.0205164.ref058],[@pone.0205164.ref059]\]. Here, we used this approach to generate robust predictions of the geographic ranges for the *M*. *lemniscatus* species complex, thus contributing to the guidance nonrandom field surveys and to the evaluation of the shifts in their distributions in response to climate change.

Regarding the current distribution estimated for each species, the potential distribution of *M*. *l*. *lemniscatus* is in agreement with the recent proposition of Silva Jr. et al. \[[@pone.0205164.ref031]\] and extends the distribution originally proposed by Campbell & Lamar \[[@pone.0205164.ref033]\] from a small area in northern Brazil to a wide area encompassing most of the Amazon biome. The map of habitat suitability indicates that potential areas for collection efforts are located in the western Amazonas and the Acre Brazilian states, south of Colombia and west of Venezuela. However, the distribution of *M*. *l*. *carvalhoi* largely differs from that suggested by Campbell & Lamar \[[@pone.0205164.ref033]\] and covers a larger area than that proposed by Silva Jr. et al. \[[@pone.0205164.ref031]\]. However, the wide area predicted by ENMs covering northern Argentina, through Bolivia, Peru and southwestern Amazonas must be considered with caution, given that no occurrence record was collected in these regions. These areas have low suitability values ([S1 Fig](#pone.0205164.s009){ref-type="supplementary-material"}), as indicated by the lowest presence threshold (0.23), probably as an effect of the high aggregation of occurrence records in the distribution core, thus resulting in lower suitability values in the distribution border where some occurrence records are placed. Despite this potential methodological artifact, our results clearly suggest that *M*. *l*. *carvalhoi* is allopatric with the Amazonian *lemniscatus*, probably indicating the evolution of distinct habitat preferences (humid forests in the case of *lemniscatus* and open and dry forests in the case of *carvalhoi*) \[[@pone.0205164.ref029]\].

The distribution of *Micrurus diutius* was slightly expanded compared to previous studies \[[@pone.0205164.ref031],[@pone.0205164.ref033]\], and a gap in the occurrence records in the northern border between Brazil, the Guianas and Suriname was suggested as an important region for further field surveys. *M*. *diutius* was recently reconsidered as a full species \[[@pone.0205164.ref060]\], but this status still warrants investigation. As *M*. *diutius* is sympatric with *M*. *l*. *lemniscatus* in most parts of its distribution, additional records from the undersampled but highly suitable areas could provide more accurate information to support splitting both species.

The distribution of *M*. *l*. *helleri* was more restricted than that proposed by Campbell & Lamar \[[@pone.0205164.ref033]\], and a larger area in northwestern Amazonia was indicated as having high suitability despite being undersampled. However, the status of this subspecies was recently reviewed by Silva Jr. et al. \[[@pone.0205164.ref031]\], who proposed its synonymy with *M*. *l*. *lemniscatus*. If so, the distribution of *M*. *l lemniscatus* will extend through most of Peru and Ecuador and probably reach the western and northern borders of the *M*. *l*. *carvalhoi* distribution.

Climate change and gap analysis {#sec011}
-------------------------------

Studies showing the impacts of climate change are common for vertebrate species worldwide \[[@pone.0205164.ref058],[@pone.0205164.ref061]--[@pone.0205164.ref065]\], but only a few of those studies have specifically addressed the potential effects on snakes \[[@pone.0205164.ref016],[@pone.0205164.ref024],[@pone.0205164.ref066],[@pone.0205164.ref067]\]. These few studies provide clear evidence that snakes are very sensitive to climate changes, with some species undergoing severe range shifts \[[@pone.0205164.ref066],[@pone.0205164.ref067]\], range contractions \[[@pone.0205164.ref068]\], or drastic population declines \[[@pone.0205164.ref016]\]. As the first study about the impacts of climate change on coralsnakes, we report here the alarming evidence that for three out the four species analyzed, the geographic distributions will probably shrink in the future.

Although more detailed (including experimental) studies are necessary to assess the susceptibility of coralsnakes to climate change, we can infer the potential impacts based on specific traits that make species more sensitive to climate change, as summarized by Foden et al. \[[@pone.0205164.ref069]\]. First, coralsnakes have specialized habitats and specific microhabitat requirements based on their semi-fossorial habit. For instance, drier soils caused by warmer temperatures and reduced precipitation may force these species to burrow deeper into the soil to find sufficient moisture or to hunt other fossorial prey \[[@pone.0205164.ref070]\]. If changes in soil humidity cause changes in the populations of other fossorial prey species such as amphisbaenians and caecilians, coralsnakes will probably also be affected because overall, snakes respond very strongly to prey densities \[[@pone.0205164.ref030],[@pone.0205164.ref070],[@pone.0205164.ref071]\].

Second, their physiology and ecology is strongly dependent on the specific ranges of climatic variables, such as temperature for thermoregulation, and precipitation and soil humidity for breeding cycles and foraging activities \[[@pone.0205164.ref030]\]. For instance, Marques et al. \[[@pone.0205164.ref017]\] showed that in triadal *Micrurus*, female ovulation occurs in spring, and males have a peak of spermatogenesis in the summer, synchronous with mating that occurs during this hot and humid season. Other similar studies corroborate the need for specific climatic conditions for the reproductive cycles of coralsnakes \[[@pone.0205164.ref030]\], and our exploratory analyses of the association between climate and habitat suitability confirms the importance of seasonal conditions (mainly precipitation) for these species. Therefore, rapid climate changes in the future probably will also require rapid evolutionary changes in phenology, which is rare (but see Moreno-Rueda et al. \[[@pone.0205164.ref068]\] and highly unlikely to occur in snakes \[[@pone.0205164.ref024]\].

Third, coralsnakes are also dependent on interspecific interactions, particularly regarding their diet specialization, consisting predominantly of amphisbaenians, caecilians and other snakes. Considering that *Micrurus* species are morphologically and ecologically conservative \[[@pone.0205164.ref029],[@pone.0205164.ref030],[@pone.0205164.ref033]\], and due to their high degree of diet specialization, they are unlikely to be able to switch to or substitute other food resources, increasing their susceptibility to disruption by climate changes.

And fourth, their poor ability to disperse to new suitable habitats because of their sedentary and subterraneous habitat prevents the species from tracking climate changes, which can be even worse due to the geographic barriers, given that some of these species are unable to overcome small rivers \[[@pone.0205164.ref030]\]. Thus, this negative combination of susceptible traits along with the climatically imposed reductions in suitable habitats is likely to highly threaten species in the future.

*Micrurus l*. *lemniscatus* was the only taxon for which future predictions indicate a gain in suitable areas (measured as the number of 0.5° cells) at the end of century, although part of this will occur in the western border of Peru, Ecuador and Colombia, i.e., outside the core of the distribution, where it is very unlikely for the species to colonize. In any case, a substantial increase in the suitable area was observed south of the current distribution (i.e., advancing towards the center of Brazil), and in central Colombia and southern Venezuela. As this species is dominant throughout the Amazon and is probably well adapted to the high temperatures and humidity characteristic of this biome \[[@pone.0205164.ref031],[@pone.0205164.ref033],[@pone.0205164.ref072]\], one could infer that it will be less affected by climate changes. However, we noted that as *M*. *l*. *lemniscatus* and *M*. *l*. *carvalhoi* are allopatric, the advance of the *M*. *l*. *lemniscatus* limits towards the *M*. *l*. *carvalhoi* distribution indicate the potential for novel interactions in the future, which can impose additional negative effects (e.g., competitive exclusion) on both species. Additional studies are necessary to evaluate the potential impacts of novel patterns of co-occurrence for these snakes, but the evidence provided here cannot be ignored.

Our study also indicates that for *M*. *l*. *carvalhoi*, *M*. *diutius* and *M*. *l*. *helleri*, the suitable habitats are poorly represented in the protected areas, both at present and for the future. This finding is particularly alarming for *M*. *l*. *carvalhoi*, which now has only 14% representation and 13% for the future, and whose distribution includes two highly threatened biomes---the Atlantic forest and the Cerrado. The Atlantic forest currently has only ca. 11% of its natural cover \[[@pone.0205164.ref073]\], and only ca. 7.2% of its remaining habitats are protected \[[@pone.0205164.ref074]\]. The Cerrado biome still retains ca. 45% of its natural cover \[[@pone.0205164.ref075]\], but its increasing transformation into pasture and cash-crop agriculture is noticeably rapidly reducing the natural areas. Even in the Amazon biome that still remains largely intact due to its great size, some species may be severely threatened. Thus, if the levels of habitat protection in the current reserve network are not sufficient, and species are being driven out of reserves due to climate change as we found here \[[@pone.0205164.ref076],[@pone.0205164.ref077]\], the remaining intact vegetation must be protected, and new areas in these biomes should be the focus of intense research and conservation actions. Most importantly, more attention should be paid to those areas considered as long-term climatic refugia \[[@pone.0205164.ref047],[@pone.0205164.ref078],[@pone.0205164.ref079]\], since they can function as buffers to the impacts of climate changes and provide a valuable alternative to preserve *in situ* species with poor dispersal abilities (i.e., without the need for migration or translocation).

*Micrurus diutius* was predicted to lose more than 80% of suitable area. In this case, it is important to note that the suitability value used as a threshold (defined from the LPT) was quite high (i.e., 0.72), resulting in a reduced distributional area. Use of the LPT, although appropriate in cases where a more conservative prediction is desired, can be influenced by the number of locality records, which increases as the sample size decreases \[[@pone.0205164.ref053]\]. *Micrurus diutius* had the lowest number of occurrence records (only 33 spatially unique records at the resolution of 0.5°), which can explain the high value of the LPT. Thus, the interpretation of the high decrease in its distributional area due to climate change, and consequently, its low representation in the current network of protected areas, should be made with caution. Nevertheless, the present distributional area obtained after applying the LPT provides valuable clues for identifying new sites that are at least as suitable as those where the species has been recorded.

In summary, through the application of niche model tools and carefully revised occurrence records, we were able to reduce, with a good deal of certainty, the Wallacean shortfall for the *M*. *lemniscatus* species complex. For the first time, we report the potential impacts of climate changes on the distribution of coralsnakes, with strong evidence that these organisms may be highly sensitive to the ongoing changes. Coralsnakes, similar to other snake species, are middle-order predators, and therefore, their extinction may have serious consequences for the functioning of ecosystems \[[@pone.0205164.ref080]\]. Thus, we undoubtedly need exhaustive studies focused on the unexplored aspects of ecology and evolution of coralsnakes, which might have far-reaching implications for understanding their responses to the changing climate and supporting more effective conservation management.

Supporting information {#sec012}
======================

###### List of museums.

List of museums from which specimens were examined.

(PDF)

###### 

Click here for additional data file.

###### Results of ENM for *M*. *l*. *carvalhoi*.

Results of ENM and analyses for *M*. *l*. *carvalhoi* after controlling for spatial aggregation in occurrence records.

(PDF)

###### 

Click here for additional data file.

###### List of occurrence records.

List of occurrence records with the collection IDs. Collection abbreviations correspond to the museums listed in the [S1 Appendix](#pone.0205164.s001){ref-type="supplementary-material"} above.

(PDF)

###### 

Click here for additional data file.

###### Number of occurrence records.

Number of occurrence records from museum collections and after mapping over the South American grid.

(PDF)

###### 

Click here for additional data file.

###### Climatic models.

Details of the climatic simulations (AOGCMs) used in the ecological niche modeling.

(PDF)

###### 

Click here for additional data file.

###### Output for factorial analysis.

Loadings of the bioclimatic variables in the first five axes of varimax rotated factor analysis, based on the CCSM AOGCM. Numbers in bold highlight the highest loading of the selected variable in each factor.

(PDF)

###### 

Click here for additional data file.

###### Ecological niche modeling methods.

Ecological niche modeling methods used to estimate species potential distributions.

(PDF)

###### 

Click here for additional data file.

###### Output of model evalution.

The mean true skill statistics (TSS) values across five AOGCMs (CCSM, CNRM, GISS, MIROC and MRI) for each ENM method.

(PDF)

###### 

Click here for additional data file.

###### Habitat suitability maps.

Consensus maps of habitat suitability derived from 12 niche modeling methods and 5 climate models. Hot colors indicate high habitat suitability; cool colors, low habitat suitability. Black dots indicate presence records used in the modeling processes; A) *Micrurus l*. *lemniscatus*, B) *Micrurus l*. *carvalhoi*, C) *Micrurus diutius*, D) *Micrurus l*. *helleri*.

(PDF)

###### 

Click here for additional data file.

###### Habitat suitability map for *M*. *l*. *cavalhoi*.

Consensus maps of habitat suitability from ENMs for *M*. *l*. *cavalhoi* for A) present, and B) future after controlling for spatial autocorrelation in occurrence records. Black dots are the most equidistant records select in the environmental space.

(PDF)

###### 

Click here for additional data file.

###### Terrestrial ecoregions.

Delimitation of the terrestrial ecoregions used to assess the level of species representation in protected areas. For each taxon, the selection of these areas was done by overlapping the map of habitat suitability (after applying specific decision thresholds) with the map of South American ecoregions. An ecoregion was selected if at least one cell overlapping it was predicted as present; A) *Micrurus l*. *lemniscatus*, B) *Micrurus l*. *carvalhoi*, C) *Micrurus diutius*, D) *Micrurus l*. *helleri*.

(PDF)

###### 

Click here for additional data file.

###### Spatial distribution of the sum of squares.

Proportion of the sum of squares (SS) from the ANOVA accounted for by the ecological niche models (ENMs), atmosphere-ocean general circulation models (AOGCMs) and time (the present and the end of the century). Dark red indicates high SS; light red, low SS. A) *Micrurus l*. *lemniscatus*, B) *Micrurus l*. *carvalhoi*, C) *Micrurus diutius*, D) *Micrurus l*. *helleri*.

(PDF)

###### 

Click here for additional data file.

The species data used in this study were derived from the research project "Phylogeography and toxinological diversity of *Micrurus* (Serpentes, Elapidae) in transitional areas of the Cerrado--Amazonia--Caatinga," integrated into the research network GENPAC (Geographical Genetics and Regional Planning for natural resources in the Brazilian Cerrado) supported by CNPq/CAPES/FAPEG. We thank PUC Goiás for supporting the "Workshop on Niche Models and New World Coralsnake Distributions" held on 27 June--01 July 2016 in Goiânia, Goiás. We thank Thiago F. Rangel for access to the computational platform BIOENSEMBLES used to generate the ecological niche models, which was also supported by the "Rede Clima" network of MCTI. We also acknowledge the World Climate Research Programme's Working Group on Coupled Modeling responsible for CMIP5 for producing and making available the climate models. Research by LCT, NJSJr and JAFD-F is continuously supported CNPq productivity grants. This paper was developed in the context of National Institutes for Science and Technology (INCT) in Ecology, Evolution and Biodiversity Conservation, supported by MCTIC/CNpq (proc. 465610/2014-5) and FAPEG.

[^1]: **Competing Interests:**The authors have declared that no competing interests exist.
