Abstract

Our understanding of the impacts of ongoing global warming on terrestrial species has increased significantly during the last several years, but how climatic change has affected, and will affect, the distribution of earthworms remains largely unknown. We used climate niche modeling to model the current distribution of the giant earthworm Rhinodrilus alatus - an endemic species of the Cerrado Domain in Brazil, which is traditionally harvested and commercialized for fishing bait. R. alatus is sensitive to environmental changes because climate, in synergy with soil attributes, determine its annual reproductive cycle and distribution. The paleoclimatic reconstructions predict important geographical shifts from LGM (~21,000 yBp) to the present potential distribution of R. alatus: range expansion, fragmentation, and shrinkage of the current core area. Further, the 2070 scenarios predict substantial shrink and losses of stable areas. Our results indicate that climate change will not only affect the extent of the distribution, but will also promote significant fragmentation and a geographical shift to outside of the currently recognized geographical boundaries. In this context, populations of R. alatus would decline and traditional harvesting would collapse, requiring immediate implementation of management and conservation measures for the species and economically sustainable alternatives for the local community.

Key words
Climate change; ecological niche modeling; giant earthworm; harvesting; management; Rhinodrilus alatus

INTRODUCTION

Climatic change is causing important changes in the spatial occurrence and distribution of species worldwide and has been the focus of researchers seeking to anticipate how species, communities and ecosystems will respond (e.g., Parmesan 2006PARMESAN C. 2006. Ecological and evolutionary responses to recent climate change. Annu Rev Ecol Evol Syst 37: 637-669., Faleiro et al. 2013FALEIRO FV, MACHADO RB and LOYOLA RD. 2013. Defining spatial conservation priorities in the face of land-use and climate change. Biol Conserv 158: 248-257., Anadón et al. 2014ANADÓN JD, SALA OE and MAESTRE FT. 2014. Climate change will increase savannas at the expense of forests and treeless vegetation in tropical and subtropical Americas. J Ecol 102: 1363-1373., IPCC 2014IPCC. 2014. Climate Change: Impacts, Adaptation, and Vulnerability. Part A: Global and Sectoral Aspects. Contribution of Working Group II to the Fifth Assessment Report of the Intergovernamental Panel on Climate Chage. Cambridge.). Environmental factors (e.g., soil type, temperature) generally exhibit heterogeneous patterns, and if a species is affected by these factors, their geographical distribution should, theoretically, reflect environmental patterns associated with the constraints of their physiology and life history (Pulliam 2000PULLIAM HR. 2000. On the relationship between niche and distribution. Ecol Lett 3: 349-361., Sexton et al. 2009SEXTON JP, MCINTYRE PJ, ANGERT AL and RICE KJ. 2009. Evolution and ecology of species range limits. Annu Rev Ecol Evol Syst 40: 415-36.). In this context, strong ecological specialization and restricted distribution can make species susceptible to long-term climate change (Thackeray et al. 2010THACKERAY SJ et al. 2010. Trophic level asynchrony in rates of phenological change for marine, freshwater and terrestrial environments. Glob Chang Biol 16: 3304-3313., Barbraud et al. 2011BARBRAUD C et al. 2011. Contrasted demographic responses facing future climate change in Southern Ocean seabirds. J Anim Ecol 80: 89-100.). Furthermore, understanding how climate fluctuations of past glacial cycles could have influenced the current distribution of species can help anticipate their sensitivity to future climate scenarios (Sérsic et al. 2011SÉRSIC AN et al. 2011. Emerging phylogeographical patterns of plants and terrestrial vertebrates from Patagonia. Biol J Linn Soc 103: 475-494., Mellick et al. 2012MELLICK R et al. 2012. Palaeodistribution modelling and genetic evidence highlight differential post-glacial range shifts of a rain forest conifer distributed across a latitudinal gradient. J Biogeogr 39: 2292-2302., Forester et al. 2013FORESTER BR, DECHAINE EG and BUNN AG. 2013. Integrating ensemble species distribution modelling and statistical phylogeography to inform projections of climate change impacts on species distributions. Divers Distrib 19: 1480-1495.). In South America, particularly, the joint use of phylogeographic and paleoclimatic reconstructions has corroborated the sensitivity of distinct taxonomic groups to climate change based on their patterns of geographic expansion and contraction (see Turchetto-Zolet et al. 2013TURCHETTO-ZOLET AC, PINHEIRO F, SALGUEIRO F and PALMA-SILVA C. 2013. Phylogeographical patterns shed light on evolutionary process in South America. Mol Ecol 22: 1193-1213. for a review).

Inferring the geographical distribution of biological species has been done mainly using Ecological Niche Models (ENMs; also referred to as ‘Species Distribution Models’, Elith et al. 2011ELITH J et al. 2011. A statistical explanation of MaxEnt for ecologists. Divers Distrib 17: 43-57., Araújo and Peterson 2012ARAÚJO MB and PETERSON AT. 2012. Uses and misuses of bioclimatic envelope modeling. Ecology 93: 1527-1539.). These models have been used extensively to describe the current distributions of species and to project them into past and future climatic conditions (Araújo and Peterson 2012), based on the projection of occurrence records onto a spatially interpolated climate surface of a given area (i.e., environmental predictors, Guisan and Zimmermann 2000GUISAN A and ZIMMERMANN NE. 2000. Predictive habitat distribution models in ecology. Ecol Model 135: 147-186.). This modeling approach assumes that species are in quasi-equilibrium with their environmental requirements over the time period (under assumptions of niche conservatism, Warren et al. 2008WARREN DL, GLOR RE and TURELLI M. 2008. Environmental niche equivalency versus conservatism: quantitative approaches to niche evolution. Evolution 62: 2868-2883.) or that their distributions are determined mainly by the environment (Norberg et al. 2012). However, the inherent assumptions and simplifications of ENMs (i.e., not considering other important factors such as soil properties, species interactions, physiological and behavioral plasticity, and dispersal capability) may overestimate or underestimate the climatic tolerances of species (Norberg et al. 2012NORBERG J et al. 2012. Eco-evolutionary responses of biodiversity to climate change. Nat Clim Chang 2: 747-751., Fordham et al. 2013FORDHAM DA, AKÇAKAYA HR, ARAÚJO MB, KEITH DA and BROOK BW. 2013. Tools for integrating range change, extinction risk and climate change information into conservation management. Ecography 36: 001-009.). Therefore, this approach may be more effective with species that have climatically defined range limits, so that errors due to physiological adaptation or from estimates of paleoclimate have less affect on predicted spatial distributions (Hugall et al. 2002HUGALL A, MORITZ C, MOUSSALI A and STANISIC J. 2002. Reconciling paleodistribution models and comparative phylogeography in the Wet Tropics rainforest land snail Gnarosophia bellendenkerensis (Brazier 1875). PNAS 99: 6112-6117.). In spite of this, ENMs can be employed as a tool for forecasting species’ responses to climatic change taking on their limitations (Parmesan 2006PARMESAN C. 2006. Ecological and evolutionary responses to recent climate change. Annu Rev Ecol Evol Syst 37: 637-669., Mellick et al. 2012MELLICK R et al. 2012. Palaeodistribution modelling and genetic evidence highlight differential post-glacial range shifts of a rain forest conifer distributed across a latitudinal gradient. J Biogeogr 39: 2292-2302.).

Earthworms are ecosystem engineers and play a key role in the physical structure of soils, modulation of nutrient availability and modification of the rate at which soil microbes carry out important ecosystem processes such as N mineralization (Eijsackers 2011EIJSACKERS H. 2011. Earthworms as colonizers of natural and cultivated soil environments. J Hydrol 199: 360-377., Arnone et al. 2013ARNONE JA, ZALLER JG, HOFER G, SCHMID B and KÖRNER C. 2013. Loss of plant biodiversity eliminates stimulatory effect of elevated CO2 on earthworm activity in grasslands. Oecologia 171: 613-622., García-Palacios et al. 2014GARCÍA-PALACIO P, MAESTRE FT, BRADFORD MA and REYNOLDS JF. 2014. Earthworms modify plant biomass and nitrogen capture under conditions of soil nutrient heterogeneity and elevated atmospheric CO2 concentrations. Soil Biol Biochem 78: 182-188.). In this work, ENM was used as a tool to understand the environmental determinants that control the distribution of the Brazilian giant earthworm, Rhinodrilus alatus (Righi, 1971RIGHI G. 1971. Sobre a família Glossoscolecidae (Oligochaeta) no Brasil. Arquiv Zool 20: 1-96.), and its sensitivity to climatic change, by projecting its current distribution into future and paleoclimatic scenarios. Rhinodrilus alatus is an endemic of the savanna formation in southwestern Brazil (Cerrado Domain), which is restricted to the southern portion of the São Francisco River Basin. This species can be found in different savanna physiognomies, as well as in pastures and eucalyptus plantations (Drumond et al. 2013).

Different climate models with different greenhouse gas emission rates, great reductions in rainfall (35%) and increase in mean annual temperature (+ 4.0 °C), have been forecast for the São Francisco River Basin by 2080 (Marengo et al. 2012), and increases of 1 to 4 °C or 2 to 6 °C have been predicted for South America by the year 2100 (Nuñez et al. 2009NUÑEZ MN, SOLMAN SA and CABRÉ MF. 2009. Regional climate change experiments over southern South America. II: Climate change scenarios in the late twenty-first century. Clim Dyn 32: 1081-1095., IPCC 2014IPCC. 2014. Climate Change: Impacts, Adaptation, and Vulnerability. Part A: Global and Sectoral Aspects. Contribution of Working Group II to the Fifth Assessment Report of the Intergovernamental Panel on Climate Chage. Cambridge.). These scenarios may have a negative impact on the life cycle of R. alatus, affecting its aestivation, foraging and reproductive cycles with consequences for its survival and reproductive performance (Drumond et al. 2013) since, as in other earthworms, they are strongly influenced by the duration of dry and wet seasons, and intensity and regularity of rains (Jiménez et al. 2000JIMÉNEZ JJ, BROWN GG, DECÄENS T, FEIJOO A and LAVELLE P. 2000. Differences in the timing of diapause and patterns of aestivation in tropical earthworms. Pedobiologia 44: 677-694., Barbraud et al. 2011BARBRAUD C et al. 2011. Contrasted demographic responses facing future climate change in Southern Ocean seabirds. J Anim Ecol 80: 89-100., Arnone et al. 2013ARNONE JA, ZALLER JG, HOFER G, SCHMID B and KÖRNER C. 2013. Loss of plant biodiversity eliminates stimulatory effect of elevated CO2 on earthworm activity in grasslands. Oecologia 171: 613-622.). In addition, ectothermic animals are especially prone to be affected by climate change as a result of their particular ecological requirements (Sinervo et al. 2010SINERVO B et al. 2010. Erosion of lizard diversity by climate change and altered thermal niches. Science 328: 894-899., Sahlean et al. 2014SAHLEAN TC et al. 2014. Refining climate change projections for organisms with low dispersal abilities: a case study of the Caspian whip snake. PLoS One 9: e91994.), limited acclimation ability, low tolerance for warming, especially if their dispersal abilities are limited (Williams et al. 2008WILLIAMS SE, SHOO LP, ISAAC JL, HOFFMANN AA and LANGHAM G. 2008. Towards an integrated framework for assessing the vulnerability of species to climate change. PLoS Biol 6: 2621-2626., Deutsch et al. 2008DEUTSCH CA et al. 2008. Impacts of climate warming on terrestrial ectotherms across latitude. PNAS 105: 6668-6672., Travis et al. 2013TRAVIS JMJ et al. 2013. Dispersal and species’ responses to climate change. Oikos 122: 1532-1540.).

Rhinodrilus alatus is traditionally used as bait for amateur fishing, and in the last 80 years its harvest, concentrated in the dry season, is the main source of employment and income for some rural communities, particularly those of afrodescent (“Quilombola”, Drumond et al. 2015). However, this historical practice of exploitation can make this giant earthworm, and the communities that depend on their extraction (ca. 3,000 people), vulnerable to climate changes (Drumond et al. 2015) in the absence of long-term management that takes into account the possibility of population declines and changes in areas of occurrence due to climate change (Marengo et al. 2012MARENGO JA et al. 2012. Development of regional future climate change scenarios in South America using the Eta CPTEC/HadCM3 climate change projections: climatology and regional analyses for the Amazon, São Francisco and the Paraná River basins. Clim Dyn 38: 1829-1848., IPCC 2014IPCC. 2014. Climate Change: Impacts, Adaptation, and Vulnerability. Part A: Global and Sectoral Aspects. Contribution of Working Group II to the Fifth Assessment Report of the Intergovernamental Panel on Climate Chage. Cambridge.).

Our objectives in this study were: (i) to evaluate the environmental factors that determine the current distribution of the giant earthworm Rhinodrilus alatus, (ii) assess its potential geographic distribution until the Last Glacial Maximum (LGM), and (iii) assess the impact of future climate change on its geographic distribution. To achieve these objectives, we modeled the spatial distribution of R. alatus using climate data and employed the latest climate-change scenarios provided by the 5th Assessment Report of the IPCC (Taylor et al. 2012TAYLOR KE, STOUFFER RJ and MEEHL GA. 2012. An overview of CMIP5 and the experiment design. Bull Am Meteorol Soc 93: 485-498.) to reconstruction of past scenarios until the LGM and to forecast how the distribution of this species will change by the year 2070.

MATERIALS AND METHODS

MODELING THE POTENTIAL DISTRIBUTION OF RHINODRILUS ALATUS

Occurrence records of R. alatus (n = 670) were obtained by collecting individuals directly from the soil in all the know distribution area, during April 2006 to July 2011 (Figure 1 and 2, Table SIII). All samples were collected by a team of biologists and local extractors who hold the traditional technique of manual giant earthworms extraction. Permission to collect on private property was obtained verbally, while permission to work with wild fauna was granted by the Authorization and Information System in Biodiversity (SisBio), from Chico Mendes Institute for Biodiversity Conservation (ICMBio, authorization number: 30164-4). We minimized the main sources of error in ENM modeling: (i) taxonomic uncertainty – the specimen identification was based on morphologic criteria proposed by Righi (1971) and phylogenetic analysis (COI and ITS1-5.8S-ITS2 regions, see Siqueira et al. 2013SIQUEIRA FF et al. 2013. Genetic diversity and population genetic structure in giant earthworm Rhinodrilus alatus (Annelida: Clitellata: Glossoscolecidae). Pedobiologia 56: 15-21. for more details), (ii) inaccurate georeferenced data for the species – we employed only occurrence records obtained in situ, and (iii) biased sample – records that cover the total distribution and altitudinal range of the species were obtained. Another source of significant uncertainty in ENM projections is the choice of a single general circulation model (GCM) for climate forecasting (Buisson et al. 2010BUISSON L, THUILLER W, CASAJUS N, LEK S and GRENOUILLET G. 2010. Uncertainty in ensemble forecasting of species distribution. Glob Chang Biol 16: 1147-1157.). For this reason, different combinations of GCMs that represent alternate possible states of the system being modeled should be used (Araújo and New 2007ARAÚJO MB and NEW M. 2007. Ensemble forecasting of species distributions. Trends Ecol Evol 22: 42-47.). When these are combined, using consensus methods (see below), they can form a more accurate projection that outperforms single GCM models (Forester et al. 2013FORESTER BR, DECHAINE EG and BUNN AG. 2013. Integrating ensemble species distribution modelling and statistical phylogeography to inform projections of climate change impacts on species distributions. Divers Distrib 19: 1480-1495.).

Figure 1
- a) Earthworm with more than one meter in length; b) Detail of Rhinodrilus alatus, with laterally projected clitellum; c) Earthworms mating - arrows indicate the two individuals; d) a cocoon within the chamber; e) a recently born earthworm; f) an immature; g) an adult (arrow shows the projected clitellum); h) earthworm removed from aestivation chamber (larger hole) - the arrow indicates a gallery that leads to the chamber; i) gallery detail; j) earthworm with a mucus web excreted during the dry season.

Figure 2
- Geographic distribution of the Rhinodrilus alatus (black dots, n = 670) and potential distribution model under Current scenario using the MaxEnt algorithm (Phillips and Dudík 2008). Warmer colours (red/yellow) of the logistic output correspond to regions with a higher probability of occurrence. The cubic convolution interpolation technique was employed. Drainage of the southern portion of the São Francisco River (purple) and tributaries (pink), and Cerrado, Atlantic Forest and Caatinga Domains in Minas Gerais State, Brazil are shown.

The environmental data consisted of 19 predictor for further modeling from WorldClim dataset (Hijmans et al. 2005HIJMANS R, CAMERON S, PARRA J, JONES P and JARVIS A. 2005. Very high resolution interpolated climate surfaces for global land areas. Int J Clim 25: 1965-1978.), trimmed to the state of Minas Gerais. We selected five variables with the least autocorrelation (pairwise r<0.60, correlation analysis employed ENMTools v.1.4.3 program (Warren et al. 2008). This procedure was done to avoid overparametrization of our modelling with redundant variables. In addition, the resulting selection was ecologically meaningful for R. alatus (Drumond et al. 2013). For example, duration of dry and wet seasons, and regularity of rains co-determine the life cycle of R. alatus, as in other earthworms, with a strong effect for its survival and reproductive performance (Jiménez et al. 2000JIMÉNEZ JJ, BROWN GG, DECÄENS T, FEIJOO A and LAVELLE P. 2000. Differences in the timing of diapause and patterns of aestivation in tropical earthworms. Pedobiologia 44: 677-694., Barbraud et al. 2011BARBRAUD C et al. 2011. Contrasted demographic responses facing future climate change in Southern Ocean seabirds. J Anim Ecol 80: 89-100., Arnone et al. 2013ARNONE JA, ZALLER JG, HOFER G, SCHMID B and KÖRNER C. 2013. Loss of plant biodiversity eliminates stimulatory effect of elevated CO2 on earthworm activity in grasslands. Oecologia 171: 613-622., Drumond et al. 2013DRUMOND MA et al. 2013. Life history, distribution and abundance of the giant earthworm Rhinodrilus alatus (Righi 1971): conservation and management implications. Braz J Biol 73: 699-708.). The final selected variables were: mean annual temperature (bio 1), mean diurnal temperature range (bio 2), temperature seasonality (bio 4), annual precipitation (bio 12) and precipitation seasonality (CV, bio 15). All layers were downloaded with a spatial resolution of 2.5 arc-min, or roughly 4x4 km².

We used the program Maxent version 3.3.3a (Phillips and Dudik 2008PHILLIPS SJ and DUDÍK M. 2008. Modeling of species distributions with Maxent: new extensions and a comprehensive evaluation. Ecography. 31: 161-175.) for modeling the geographic distributions of R. alatus. The selection of Maxent algorithm was done because the predictive performance is consistently competitive with the highest-performing methods presence-absence data (see Elith et al. 2011). In addition, Maxent’s can be modelled with presence-only data (Phillips and Dudik 2008). One presence record per pixel was used for a total of 670 records. We used Maxent’s auto-features: (i) default regularization multiplier parameter (1.0), (ii) convergence threshold (10^-5) and (iii) maximum number of iterations (500). We also used 70% of the localities for model training with replicate sample sets chosen by sampling with replacement (bootstrap, 100 replicates), and contrasted them with 10,000 background points sampled from the remainder of the study area to determine the Maxent distribution and test model performance. Maxent’s logic is based on the assumption that the potential geographic range of a species can be determined by a set of spatialized covariates (such as mean annual temperature and rainfall, and duration of the dry season) that define its fundamental niche, under assumptions of niche conservatism (Norberg et al. 2012). Species occurrences appear as high probabilistic clusters against a climatic/geographic background that has the maximum entropy permitted by the information not available and associated with randomly sampled localities (see Phillips and Dudík 2008PHILLIPS SJ and DUDÍK M. 2008. Modeling of species distributions with Maxent: new extensions and a comprehensive evaluation. Ecography. 31: 161-175., Elith et al. 2011ELITH J et al. 2011. A statistical explanation of MaxEnt for ecologists. Divers Distrib 17: 43-57.). We used area under the curve (AUC) to infer the accuracy of the models, which has been employed extensively with ENMs (Elith et al. 2011ELITH J et al. 2011. A statistical explanation of MaxEnt for ecologists. Divers Distrib 17: 43-57., Jiménez-Valverde 2012JIMÉNEZ-VALVERDE A. 2012. Insights into the area under the receiver operating characteristic curve (AUC) as a discrimination measure in species distribution modeling. Glob Ecol Biogeogr 21: 498-507.).

Models for R. alatus were projected onto current conditions, paleoclimatic conditions and the conditions of future climate-change scenarios (^{See Supplementary Material - Table SI} SUPPLEMENTARY MATERIAL SUPPLEMENTARY MATERIAL Table SI - list of the 17 global climate models downscaled and calibrated Coupled Model Intercomparison Project Phase 5 (CMIP5) for 2070 under the RCP4.5 and RCP8.5 scenarios. Table SII - AUC values for all models (n=100 replicates) building with Training and Test data. In addition, AUC values without one variable and with only one variable are present. Table SIII - Occurrence registers of Rhinodrilus allatus employed to interpolation with scenopoetic predictors (spatial resolution of 2.5 x 2.5 arcmin, or roughly 4x4 km²) by MaxEnt algorithm. Longitude (Long) and Latitude (Lat) coordinates in decimals degrees. Figure S1 - Logistic output of the probability of presence for all occurrence points of R. alatus. Figure S2 - Stable surface from current to future for the RCP4.5 scenario. Light and dark green represent stable areas above and below of 900 meters a.s.l., respectively. Satellite images of some stable areas are displayed. Scale bar equals 15 km. ). Current distributions were based on mean climatology from the 1960 to 1990 baseline period (Hijmans et al. 2005). Both the past and the future models were calibrated (bias corrected) using WorldClim 1.4 (release 3) as baseline of ‘current’ climate. For the reconstruction of past scenarios, we used three models of the Last Glacial Maximum (LGM, ~21,000 years before present - yBp) and nine of the Mid-Holocene (MH, ~ 6,000 yBp) from downscaled global climate models (Coupled Model Intercomparison Project Phase 5 - CMIP5). These scenarios represent cool-dry (LGM) and cool-wet conditions (MH) of the Quaternary climate (Clapperton 1993CLAPPERTON CM. 1993. Nature of environmental changes in South America at the Last Glacial Maximum. Palaeogeogr Palaeoclimatol Palaeoecol 101: 189-208., Behling 2002BEHLING H. 2002. South and southeast Brazilian grasslands during Late Quaternary times: a synthesis. Palaeogeogr Palaeoclimatol Palaeoecol 177: 19-27., Mayle et al. 2004MAYLE FE. 2004. Assessment of the Neotropical dry forest refugia hypothesis in the light of palaeoecological data and vegetation model simulations. J Quat Sci 19: 713-720.). Our rationale behind the selection of the distinct paleoclimatic models was to assess the sensitivity of R. alatus to climate variation and its correspondence with observed haplotype variation (see Siqueira et al. 2013 for details). For future scenarios, we used 17 models downscaled from CMIP5 (See Table SI). For our projections we selected the Representative Concentration Pathway 4.5 (RCP4.5) and 8.5 (RCP8.5) for the year 2070. Within the Fifth IPCC Assessment Report (Taylor et al. 2012), RCP4.5 represents the scenario with the medium-low concentration of greenhouse gases, with a predictive radiative forcing of + 4.5 W m^-2 and 1.8 °C (1.1 to 2.6 °C) global mean temperature (closest SRES-B1 equivalent), while RCP8.5 represents the scenario with the highest concentration of greenhouse gases, with a predictive radiative forcing of + 8.5 W m-², and 3.7 °C (2.6 to 4.8 °C) global mean temperature (closest SRES-A2 equivalent) (see van Vuuren et al. 2011VAN VUUREN D et al. 2011. The representative concentration pathways: an overview. Clim Chang 109: 5-31.). We not analyzed data for RCP2.6 scenario because recent studies show that GCMs for RCP2.6 (so called “negative emissions”) is unlikely (i.e., with Global CO₂ emissions of ~4 PgCyr^-1 at 2050, Peters et al. 2013PETERS GP et al. 2013. The challenge to keep global warming below 2 °C. Nature Clim Chang 3: 4-6.) when compared with the most current estimates (~10 PgCyr^-1 at 2017, Friedlingstein et al. 2014FRIEDLINGSTEIN et al. 2014. Persistent growth of CO2 emissions and implications for reaching climate targets. Nat Geosci 7: 709-715., Le Quéré et al. 2018LE QUÉRÉ C et al. 2018. Global Carbon Budget 2017. Earth Syst Sci Data. 10: 405-448.). The selection of the best (unlikely) and worst (more plausible) scenario is to examine whether areas presently occupied by R. alatus will be climatically suitable in the future, as well as to capture the directions of changes (Anadón et al. 2014). Average and standard deviation of annual rainfall and annual mean temperature from downscaled CMIP5 were congruent with regional estimates (see Marengo et al. 2012MARENGO JA et al. 2012. Development of regional future climate change scenarios in South America using the Eta CPTEC/HadCM3 climate change projections: climatology and regional analyses for the Amazon, São Francisco and the Paraná River basins. Clim Dyn 38: 1829-1848., Alvares et al. 2013ALVARES C, STAPE JL, SENTELHAS PC, GONÇALVES JLM and SPAROVEK G. 2013. Köppen’s climate classification map for Brazil. Meteorol Z 22: 711-728.).

MODELING THE EXTENT AND DISTRIBUTION UNDER CLIMATE-CHANGE SCENARIOS

Models for R. alatus that indicate the presence of suitable habitat were projected across current and past climates and probabilistic output was converted to presence/absence using minimum training presence logistic threshold (MinTP=0.191) as a cutting method for the models (Liu et al. 2005LIU C, BERRY PM, DAWSON TP and PEARSON RG. 2005. Selecting thresholds of occurrence in the prediction of species distributions. Ecography 28: 385-393.). We sum the 13 presence/absence maps to assess areas of historical stability and possible climate refugia during the climatic oscillations of the Quaternary (i.e., pixels that are classified as suitable across Current, MH (9 models, ~ 6,000 yBp) and LGM (3 models, ~21,000 yBp) scenarios). These intersections of all models or with higher congruence among models resulting in a consensual refugial map (see Carnaval and Moritz 2008CARNAVAL AC and MORITZ C. 2008. Historical climate modelling predicts patterns of current biodiversity in the Brazilian Atlantic forest. J Biogeogr 35: 1187-1201., and Werneck et al. 2012 for examples). A similar approach was employed in the building of the future consensus probability maps to RCP4.5 and RCP8.5 scenarios, summing all 17 presence/absence scenarios. Each scenario was combined with slices of 20% or more of the models (see Forester et al. 2013FORESTER BR, DECHAINE EG and BUNN AG. 2013. Integrating ensemble species distribution modelling and statistical phylogeography to inform projections of climate change impacts on species distributions. Divers Distrib 19: 1480-1495.).

Quantification of current and future stable, expansion and retraction areas of R. alatus was analyzed. Stable areas are those that currently correspond to R. alatus distribution and will continue to do so in the future. Expansion areas are those that do not correspond to R. alatus distribution in the present but could potentially be colonized in the future due to climate change. Retraction areas are those that currently support R. alatus, but will shrink in the future. We performed all geographic information system (GIS) analyses in ArcGIS v. 10 (ESRI 2011).

RESULTS

ENVIRONMENTAL PREDICTORS AND MODEL EVALUATION

ENMs estimated under current climate conditions had mean AUC values >0.98 (± 0.003 SD), indicating a very good ability to discriminate between locations of presence and absence of R. alatus. This metric was >0.83 for all models with any of the five predictors (Table SII). The resulting predictive distribution map is shown in Figure 2. Our model showed high accuracy (correctly classified 97% of presences). Based on Maxent jackknife scores, annual precipitation (42.8% of variance), temperature seasonality (19%) and mean annual temperature (19%) were the most important predictors of R. alatus distribution. R. alatus occurs in environments with narrow ranges of temperature seasonality and annual rainfall (Figure S1).

FROM LGM TO CURRENT CONDITIONS

Our models predict important geographical shifts in the potential distribution of R. alatus, from dryer and cooler conditions during LGM and cool-wet conditions in the MH to the current moist and warm conditions (Figure 3). A fragmented narrow-range potential distribution during LGM (Figure 3a) gave place to a more widespread, large range, but also a fragmented distribution in the moderate cool-wet Mid-Holocene (Figure 3b), extending to the northwest of the State. From the Mid-Holocene to the present, the potential distribution of R. alatus contracted towards a central region of the state of Minas Gerais, with medium size fragmented ranges in the mid-north; the northwest range disappeared. Some areas in the core distribution of R alatus showed a high probability of climatic stability from LGM to current scenarios (Figure 3c). These historically stable areas most southwards, scattered in Minas Gerais state, showed higher haplotypic diversity (Figure 3c).

Figure 3
3 - Predicted occupancy across Minas Gerais State, Brazil of Rhinodrilus alatus during: (a) Last Glacial Maximum (LGM) 21,000 yBp (cool-dry) and (b) Mid-Holocene (MH) 6,000 yBp (cool-wet). Consensus refugial projection indicating stability surface for Rhinodrilus alatus across the Quaternary (c) produced by overlapping predicted logistic outputs under LGM, MH and current scenarios.

FROM CURRENT CONDITIONS TO 2070

All 17 GCMs for RCP4.5 (r = 0.400 to 0.914; n = 91 comparisons) and all RCP8.5 (r = 0.356 to 0.831) were positively and significantly correlated (P<0.001). On the other hand, projections for the RCP4.5 (Figure 4a) and RCP8.5 (Figure 4b) scenarios did not differ substantially. In this case, the concordance among multiple GCMs indicates their convergence to a common spatiotemporal pattern, with both showing dramatic losses of suitable habitats (>90%) from the present range to 2070 (RCP4.5, Figure 4c; RCP8.5, Figure 4d). These projections show contractions and fragmentation of suitable habitat, and a southward latitudinal shift and altitudinal shift of suitable locations to outside of its currently recognized geographical boundaries (Figure 4 and Figure S2).

Figure 4
- Distribution models of (a) RCP4.5 and (b) RCP8.5 future scenarios for R. alatus. Red color represents areas with higher probabilities of occurrence (current) and stability surface based on a cut-off of 20% of the models (future). Change in suitable habitat from current to future for the (c) RCP4.5 and (d) RCP8.5 scenarios.

DISCUSSION

The high sampling effort allowed a confident characterization of the geographical and ecological range of R. alatus, enabling a congruent ENM, which confirms R. alatus as a savanna species associated with the Cerrado Domain (after Ab’Saber 2003AB’SÁBER A. 2003. Os domínios de natureza no Brasil: potencialidades paisagísticas. Ateliê Editorial, São Paulo.). In spite of its geographical range, well-defined temperature boundaries and seasonality in annual rainfall determine the geographical distribution of R. alatus and limit it from spreading (see Wiens and Graham 2005WIENS JJ and GRAHAM CH. 2005. Niche conservatism: integrating evolution, ecology, and conservation biology. Annu Rev Ecol Evol Syst 36: 519-539.). Narratives of the occurrence of this species by earthworm extractors (Drumond et al. 2008DRUMOND MA, CAMPOS SHC, GUIMARÃES AQ and NUNES JT. 2008. Ecologia e uso do minhocuçu Rhinodrilus alatus. MG-Biota 1: 5-23.), correspond to some potential areas outside its core distribution, corroborating the predictive power of the model and reinforcing the importance of local knowledge in improving biological data, as also stated by other researchers (Berkes and Folke 1998BERKES F and FOLKE C. 1998. Linking social and ecological systems: management practices and social mechanisms for building resilience. Cambridge: Cambridge University Press., Bender et al. 2014BENDER MG et al. 2014. Local ecological knowledge and scientific data reveal overexploitation by multiyear artisanal fisheries in the South-western Atlantic. PLoS One. 9(10): e110332.).

PAST SCENARIOS: RANGE EXPANSION AFTER THE LGM

Distribution patterns result from an amalgam of historical and current processes (Peterson et al. 2011PETERSON AT et al. 2011. Ecological niches and geographic distributions. Princeton University Press.). Long term fluctuations in precipitation and temperature during glacial periods have led to cycles of contraction, fragmentation, expansion and recolonization of the geographical ranges of vegetation (see Behling 2002BEHLING H. 2002. South and southeast Brazilian grasslands during Late Quaternary times: a synthesis. Palaeogeogr Palaeoclimatol Palaeoecol 177: 19-27., Ledru et al. 2006LEDRU MP et al. 2006. Millenial-scale climatic and vegetation changes in a northern Cerrado (Northeast, Brazil) since the Last Glacial Maximum. Quat Sci Rev 25: 1110-1126., Carnaval and Moritz 2008CARNAVAL AC and MORITZ C. 2008. Historical climate modelling predicts patterns of current biodiversity in the Brazilian Atlantic forest. J Biogeogr 35: 1187-1201., Werneck et al. 2012WERNECK FP et al. 2012. Climatic stability in the Brazilian Cerrado: Implications for biogeographical connections of South American savannas, species richness, and conservation in a biodiversity hotspot. J Biogeogr 39: 1695-1706., Bueno et al. 2017BUENO ML et al. 2017. Effects of Quaternary climatic fluctuations on the distribution of Neotropical savanna tree species. Ecography 40: 403-414., Costa et al. 2017COSTA GC et al. 2017. Biome stability in South America over the last 30 kyr: Inferences from long-term vegetation dynamics and habitat modeling. Glob Ecol Biogeogr 27: 285-297. for changes in Brazilian vegetation), followed by animal species that dwell in them in the northern and southern hemispheres. While the geographical range of forest species was reduced, the response of species of open domains was more variable, but with a marked tendency to expand during dry periods (Turchetto-Zolet et al. 2013TURCHETTO-ZOLET AC, PINHEIRO F, SALGUEIRO F and PALMA-SILVA C. 2013. Phylogeographical patterns shed light on evolutionary process in South America. Mol Ecol 22: 1193-1213.).

In spite of being a savanna species, our projections indicate that only a very small proportion of the extant range of R. alatus, near the periphery of the two actual core areas, would have been climatically suitable in the cool-dry conditions of the Last Glacial Maximum. However, with the advent of moist and warmer conditions of the Mid-Holocene (approx. 6 000 yBp), R. alatus could expand its distribution towards discrete habitats, mainly to the north and northwest of the state of Minas Gerais. In addition, R. alatus could select and/or displace its distribution towards the boundaries of riverine forests that constitute mesic enclaves within the Cerrado, even during the driest phases of the last glacial period (Redford and Fonseca 1986REDFORD KH and FONSECA GAB. 1986. The role of Gallery forests in the zoogeography of the Cerrado’s non-volant mammalian. Biotropica 18: 126-135.). As a result, its distribution would become increasingly patchy even within potential refugia. Fragmented distribution, variable environments and low dispersive capability compose the scenario for adaptive divergence and reproductive isolation (Schluter 2009SCHLUTER D. 2009. Evidence for ecological speciation and its alternative. Science. 323: 737-741.), giving rise to several genetic lineages of R. alatus. The match between refugia identified by the ENMs and the phylogeographic analysis of Siqueira et al. (2013) supports divergence in allopatry and niche conservatism in R. alatus across the modeled time span. In this sense, our paleo-ENM has provided the explicit spatial argument needed to bring together the complex phylogenetic population structure and the paleogeographic history of this giant earthworm. R. alatus was not detected in the northwest of its potential range (Drumond et al. 2013), suggesting that the São Francisco River may be a physical barrier to dispersal.

FUTURE SCENARIOS: GAIN AND LOSSES OF AREAS

The 2070 scenarios predict substantial losses of climatically suitable habitat across the range of R. alatus. However, consensus projections anticipated a gain in range to the south toward the Atlantic Forest Domain, a more densely populated and economically developed region. Although R. alatus occurs in disturbed areas, such as pastures and eucalyptus plantations, it is not probable that this savanna species could expand its distribution into agricultural land due to its sensibility to pesticides (Drumond et al. 2015DRUMOND MA et al. 2013. Life history, distribution and abundance of the giant earthworm Rhinodrilus alatus (Righi 1971): conservation and management implications. Braz J Biol 73: 699-708.). In addition, the annual life cicle of the R. alatus is strong associated with seasonal precipitation condition, and its survival is negatively affected by soils under high moist level and with less than three consecutive dry months per years (i.e., typical of Atlantic Tropical Forest; Drumond et al. 2013). Consensus projections also show a gain of suitable habitats along the eastern sides of the Espinhaço Mountain Range (where rupestrian grasslands over quartzite bedrock predominate) and the Iron Quadrangle (where rupestrian grasslands over ironstone predominate). However, displacements over long distances would be unlikely, due to the low dispersal capacity of terrestrial oligochaetes (Eijsackers 2011EIJSACKERS H. 2011. Earthworms as colonizers of natural and cultivated soil environments. J Hydrol 199: 360-377.), especially considering the short time span to 2070. In addition, taking into account current records of R. alatus in latosols and no records above 900 m, it is unlikely that this species can adapt to the harsh soil conditions of these two mountainous areas, where shallow, hard, coarse-textured, dystrophic soils derived from the nutrient-poor quartzitic and ironstone bedrock (the latter with high levels of heavy metals) predominate (Benites et al. 2007BENITES VM, SHAEFER CEGR, SIMAS FNB and SANTOS HG. 2007. Soils associated with rock outcrops in the Brazilian mountain ranges Mantiqueira and Espinhaço. Brazil J Bot 4: 569-577., Jacobi et al. 2007JACOBI CM, CARMO FF, VINCENT RC and STEHMANN JR. 2007. Plant communities on ironstone outcrops: a diverse and endangered Brazilian ecosystem. Biodiver Conserv 16: 2185-2200., Schaefer et al 2016SCHAEFER CEGR et al. 2016. The Physical Environment of Rupestrian Grasslands (Campos Rupestres) in Brazil: Geological, Geomorphological and Pedological Characteristics, and Interplays. In Fernandes GW (Ed), Ecology and Conservation of Mountain-Top Grasslands in Brazil. Springer, International Publishing, Switzerland, p. 15-70.). These soils, also characterized as being low in clay and organic matter content, probably would challenge the building of aestivation chambers and reduce their capability for water retention, hence reducing R. alatus survival, mainly in the dry season (Drumond et al. 2013DRUMOND MA et al. 2013. Life history, distribution and abundance of the giant earthworm Rhinodrilus alatus (Righi 1971): conservation and management implications. Braz J Biol 73: 699-708.).

FUTURE SCENARIOS: POPULATION DECLINE AND REFUGIA

The concensus of ENMs suggest that R. alatus will be a very threatened species in the coming decades, declining towards rarity or even becoming extinct, even in the area of its current core distribution. However, this worm may partially endure the enhancement of rainfall decline and uncertainty through behavioral adaptations that could change these forecasts. In fact, local knowledge suggests that the depth of the aestivation chambers increase with soil dryness and temperature (Drumond et al. 2015DRUMOND MA, GUIMARÃES AQ and SILVA RHP. 2015. The hole of local knowledge and traditional extraction practices in the management of giant earthworms in Brazil. PLoS One 10(4): e0123913.). By building deeper chambers, where soil humidity is higher and temperature lower, this worm may reduce chances of death by desiccation. In more vegetated woody areas within the Cerrado, soil humidity may be higher, due to hydraulic uplifting (Moreira et al. 2003MOREIRA MZ et al. 2003. Hydraulic lift in a neotropical savanna. Funct Ecol 17: 573-581., Oliveira et al. 2005OLIVEIRA RS et al. 2005. Deep root function in soil water dynamics in Cerrado savannas of central Brazil. Fuct Ecol 19: 574-581., Scholz et al. 2008SCHOLZ FG et al. 2008. Biophysical and life-history determinants of hydraulic lift in Neotropical savanna trees. Funct Ecol 22: 773-786.), which, we suppose, would benefit R. alatus in these patchy areas. However, large Cerrado areas have been deforested, lowering water tables and reducing or even drying water courses (Scott and Lesch 1997SCOTT DF and LESCH W. 1997. Streamflow responses to afforestation with Eucalyptus grandis and Pinus putula and to felling in the Mokobulaan experimental catchments, South Africa. Appl Soil Ecol 50: 1-13., Dye 2000DYE PJ. 2000. Water use efficiency in South African Eucalyptus plantations: a review. S Afr Forest J 189: 17-26.). On the other hand, some suitable habitats predicted for 2070 (RCP4.5), are congruent with relict climate localities that would have been maintained from the mid-Holocene to the present. These stable areas would fit the necessary conditions for the conservation of this species but they are not legally protected in the present.

Whatever the case, the traditional use of R. alatus would be affected at different levels. The predict changes in the distribution of R. alatus would cause an intense decline in the practice of extractivism. If this species adapts to changes by building deeper chambers, extraction would become harder, requiring increased effort and causing increased impact on the environment. In order to minimize these impacts, “ex situ” conservation measures should be implemented, although all present attempts for this, based on captive breeding, had not been successful (Drumond et al. 2008).

CONCLUSIONS

This study approach offers a climate framework on the regional scale for evaluating the sensitivity of R. allatus to paleoclimatic conditions, and the impacts of ongoing global warming. We confirmed that geographic-ecology distribution of R. allatus is well-defined and related to mean annual temperature and precipitation seasonality. In addition, some suitable habitats predicted for 2070, are congruent with relict climate localities during the Quaternary climatic fluctuations, and should be considered as priority areas for the establishment of protected areas. On the basis of projected current and future potential distribution, this approach illustrates to which extent these giant earthworm could potentially decline, and we alert to necessity of immediate implementation of management and conservation measures for the species and economically sustainable alternatives for the local community.

ACKNOWLEDGMENTS

The authors are deeply grateful to giant earthworms extractors for help with the fieldwork and for sharing their knowledge. Two anonymous reviewers provided numerous comments that improved a previous version of this manuscript. We would also like to thank the staff of Floresta Nacional de Paraopeba for their support during all stages of this project. This study was funded by Programa de Extensão Universitária (ProExt) of Ministério da Educação (Brazil) and Pró-Reitoria de Extensão (www2.ufmg.br/proex) of Universidade Federal de Minas Gerais (www.ufmg.br); Fundação de Amparo à Pesquisa do Estado de Minas Gerais (www.fapemig.br).

REFERENCES

Publication Dates

History