ORIGINAL RESEARCH article

Front. Ecol. Evol., 13 August 2026

Sec. Biogeography and Macroecology

Volume 14 - 2026 | https://doi.org/10.3389/fevo.2026.1886733

Climate change and conservation of caecilians (order Gymnophiona): prediction of distributional shifts and risk assessment in Chiapas, Mexico

  • 1. Laboratorio de Ecología Evolutiva, Instituto de Ciencias Biológicas, Universidad Autónoma de Ciencias y Artes de Chiapas., Tuxtla Gutiérrez, Mexico

  • 2. Biodiversidad, Conservación y Restauración A. C., San Cristóbal de Las Casas, Mexico

  • 3. Departamento de Sistemática y Ecología Acuática, El Colegio de la Frontera Sur Unidad Chetumal., Chetumal, Mexico

Abstract

Introduction:

Caecilians (Order Gymnophiona) remain poorly understood due to their fossorial lifestyle, which limits detection and ecological knowledge. This study is framed within the climate niche tracking theory, which posits that species may shift their geographic distributions to maintain suitable environmental conditions under climate change. We also evaluate the conservation mismatch hypothesis, which suggests that static protected natural areas become increasingly ineffective as species distributions shift beyond their boundaries In this study, we analyzed three species of caecilians (Dermophis mexicanus, Dermophis oaxacae and Gymnopis syntrema) in Chiapas, Mexico, to examine projected spatial and altitudinal changes and their representation within protected natural areas under future climate scenarios.

Methods:

We developed species distribution models under current and future conditions using two climate scenarios (optimistic and pessimistic) for 2060 and 2100, based on 235 occurrence records. From these models, we quantified changes in geographic range and elevation and estimated a protection index at both species-specific and overall levels.

Results:

We found marginally significant differences in predicted distribution area across scenarios (Friedman χ² = 11.571, p = 0.041). Altitudinal shifts were not significant (Friedman χ² = 3.41, p = 0.64), indicating no consistent pattern among species. Protected natural areas coverage declined from 58% under current conditions to 27% under severe scenarios by the end of the century.

Discussion:

Our results reveal an increasing mismatch between species distributions and protected natural areas under climate change, underscoring the need for more dynamic conservation strategies that account for spatial and elevational shifts.

1 Introduction

The order Gymnophiona, commonly known as caecilians, is poorly studied because their fossorial habits and cryptic nature hinder detection and field study. This has led to significant knowledge gaps regarding their biology, ecology, and distribution (; ). For example, a global analysis of emerging threats to amphibians assessed only 115 of the 206 described species, as 44% are data-deficient and 17% are threatened ().

On the other hand, studies on the potential distribution of caecilians are scarce or nonexistent and have primarily focused on range extensions based on new records (; ; ), rather than on predictive modeling. This limitation highlights the need to apply species distribution models under current and future Climate Change (hereafter CC) scenarios. Although the specific threats to caecilians remain unclear due to limited data, habitat loss is likely a major driver of population declines, as observed in other amphibians. However, other biotic and abiotic factors are also important, and although CC has not yet been demonstrated to be the direct cause of observed declines in caecilian populations, it remains an important factor to consider ().

Historically, the caecilian fauna of Mexico has been represented by very few species, primarily from the genus Dermophis, including Dermophis mexicanus and Dermophis oaxacae (). However, documented a range expansion of Gymnopis syntrema approximately 88 km westward from Guatemala, recording this species in Mexico for the first time from the Montes Azules Biosphere Reserve in the Selva Lacandona, Chiapas. Consequently, all three caecilian species known from Mexico occur in Chiapas, making it a critical region for the conservation of this amphibian group.

Given the restricted knowledge of their ecology and the potential effects of climate change on their distributions, understanding how these species may respond to future environmental changes is essential for developing effective conservation strategies. Accordingly, this study considered the climate niche tracking theory, which posits that species modify their distributional ranges to remain within favorable climatic conditions as environmental conditions change (; ). The conservation mismatch hypothesis was also examined, which predicts that static protected natural areas (hereafter PNAs) become progressively less effective as species shift beyond their boundaries under CC (; ). Therefore, our objective was to evaluate the three species of caecilians distributed in Chiapas and understand how their spatial and altitudinal distributions may change under two scenarios, as well as to determine the extent to which the network of PNAs will be able to protect a representative proportion of their distribution ranges. To achieve this, Species Distribution Models were used under current and future conditions, considering two scenarios (optimistic and pessimistic) for 2060 and 2100. Additionally, the rate of change in predicted area under different CC scenarios was calculated.

2 Methods

2.1 Study area

The study was conducted in the state of Chiapas, located in southeastern Mexico. Chiapas covers an approximate area of 73,311 km², representing 3.7% of the national territory (). The state exhibits an altitudinal range from −6 to 4,060 m (). In addition, it presents high environmental heterogeneity, with 17 recognized vegetation types (). In terms of amphibian diversity, Chiapas ranks third nationwide, harboring 110 species, including three caecilian species (; ). Owing to its high biological richness, the state contains a network of 67 PNAs that protect 17.6% of its territory. These PNAs are distributed among different management categories: 18 Voluntarily Designated Conservation Areas, 29 state reserves, and 20 federal reserves. In turn, the federal PNAs are classified into six categories: Biosphere Reserves (7), National Parks (3), Flora and Fauna Protection Areas (5), Natural Monuments (2), Natural Resources Protection Areas (2), and Sanctuaries (1) (; ; ; Figure 1).

Figure 1

2.2 Conservation biology status classification and endemism

The three species of caecilians distributed in Chiapas (D. mexicanus, D. oaxacae, and G. syntrema) were considered. Finally, it was determined whether these species are included in any risk category according to the Mexican Official Standard NOM-059-SEMARNAT-2010 (hereafter NOM-059) and the International Union for Conservation of Nature Red List, and whether they are endemic.

2.3 Occurrence data

The occurrence database for each species was compiled from five sources: 1) the Global Biodiversity Information Facility () databases, 2) the National Commission for the Knowledge and Use of Biodiversity of Mexico (CONABIO), 3) herpetological museum collections, including the Smithsonian Institution National Museum of Natural History (USNM), the Royal Ontario Museum (ROM), the University of Michigan Museum of Zoology (UMMZ), and the Field Museum of Natural History (FMNH), 4) scientific publications, and 5) field records collected from 2008 to 2025. Subsequently, all obtained records underwent a quality control process. Records prior to 1960, as well as those located outside the study area or lacking geographic coordinates, were excluded. To avoid pseudoreplication and reduce spatial bias, the following criteria were applied: 1) removal of duplicate records, and 2) filtering clustered records located less than 1 km apart so that only one occurrence per pixel was retained. All data selection procedures were performed using the ellipsenm package and spThin package () in R ().

2.4 Selection of current and future environmental predictors

For each species, the accessible area (“M”) was estimated as the geographic region that the species could have occupied over its evolutionary history (). To achieve this, known occurrence records for each species were used to construct an Alpha Hull, to which a 10 km buffer area was applied. This procedure was carried out independently for each species. Within the generated “M”, 15 climatic variables from WorldClim version 2.1 () were obtained for the 1970–2000 period. Climatic variables eight, nine, eighteen, and nineteen combine precipitation and temperature, which present unusual spatial anomalies in the form of discontinuities between neighboring pixels (); therefore, they were excluded from the analysis.

Future environmental predictors were also calibrated within “M”. The selected Global Climate Model was the Model for Interdisciplinary Research on Climate, version 6 (MIROC6; ). Two scenarios were used: 1) optimistic and 2) pessimistic. The optimistic scenario corresponded to SSP1-2.6, a sustainable pathway. This scenario is characterized by a 2.6 W/m² increase in radiative forcing by the year 2100, representing a low-emissions trajectory aimed at limiting global warming to 2 °C. It assumes that significant measures to reduce greenhouse gas emissions will be implemented around 2050. The pessimistic scenario corresponded to SSP5-8.5, a fossil fuel-based development scenario. This scenario projects atmospheric CO2 concentrations reaching up to 1,135 ppm, with an additional radiative forcing of 8.5 W/m² by the year 2100 ().

Both current and future environmental variables were downloaded in raster format at an approximate resolution of ~1 km². The dimensionality of the climatic variables was reduced through Principal Component Analysis (PCA) using the raster package (), kuenm package (), and sf package () in R (). Components explaining at least 85% of the total variance were selected.

2.5 Current and future potential distribution modeling procedure

Only species with a minimum of 10 reliable records were used, in order to guarantee the robustness of the species distribution models. This threshold considers the effect of sample size on modeling performance and the uncertainty in predictions (; ). To select the best model for each species and evaluate its performance, three criteria recommended by were considered: 1) partial ROC; 2) omission rates ≤ 5%; and 3) models with delta AICc values ≤ 2. The best-performing models were converted into raster format in ArcMap 10.8 to generate probabilistic maps for each species. To obtain probabilistic maps for the state of Chiapas, the maps were clipped to include only the state’s political boundaries using the geoprocessing tools in ArcMap 10.8. Finally, a 10% minimum threshold was applied to transform probabilistic models into binary models (presence-absence maps). Current and future potential distribution models were generated using Maxent 3.4.4. (). It should be noted that the ecological niche models were initially calibrated using only climatic variables to characterize the climatic component of habitat suitability. However, major anthropogenic land transformations were subsequently incorporated when estimating current potential distributions; The procedures used to account for these factors are described next section.

2.6 Rate of change

To determine changes in distribution area, differences between the predicted areas of current and future potential distributions were calculated for each species (in km²). To evaluate differences in distribution area across scenarios, a non-parametric approach was used due to the small sample size and the dependence of the data. A Friedman test was applied, allowing comparisons among scenarios without assuming normality of the data. Subsequently, to identify which specific groups differed from one another, a post-hoc analysis was performed using Wilcoxon tests with Bonferroni adjustment. Likewise, altitudinal ranges were obtained for each species under both current and future potential distributions in order to detect shifts in altitudinal distribution. Only for the current potential distributions, agricultural and urban areas were removed using the Land Use and Vegetation Series VII (), and deforested areas were obtained from .

2.7 Species protection index

To calculate the Species Protection Index (SPI; ), the PNAs shapefile for the state of Chiapas obtained from was used, including federal, state, and Voluntarily Designated Conservation Areas. Overlaps among polygons belonging to different protection categories were identified. To avoid duplication within the protected surface area, redundant polygons were removed, prioritizing those corresponding to federal PNAs.

To apply the SPI, a conservation target based on each species’ distribution area was defined. These conservation targets were established according to the following criteria: 100% for species with distribution areas smaller than 10,000 km², 15% for species with areas larger than 250,000 km², and a log-linear relationship was used for species with intermediate-sized distributions. This approach recognizes that species with small geographic ranges require a greater proportion of their range to be protected to ensure their persistence.

Once conservation targets for each species were defined and the PNAs layer was refined, the percentage of suitable habitat within PNAs was calculated. To do this, the species distribution maps were spatially integrated with the PNAs cover layer using the sf package () in R (). The SPI was obtained using the following formula:

The resulting SPI is a continuous index from 0 to 100, where higher values reflect a greater degree of compliance with the species representation objectives. Subsequently, the global SPI for the state of Chiapas was estimated by averaging the protection scores obtained for each species. This procedure was applied to both current and projected distributions, assuming the network of PNAs will remain unchanged over time.

3 Results

A total of 235 occurrence records were obtained for the three caecilian species known to occur in Chiapas, which are vulnerable to CC according to the assessment by Larry D. (Table 1). Regarding their conservation status, two species are subject to special protection under NOM-059, while one is considered Near Threatened and two are classified as Least Concern according to the International Union for Conservation of Nature (IUCN) risk categories (Table 1).

Table 1

EspecieNo. registroNOM-059EndemismoUICNVulnerabilidad
Dermophis mexicanus191PrLCM
Dermophis oaxacae34PrMéxicoLCM
Gymnopis syntrema10NTH

Caecilian species distributed in Chiapas, conservation status, records and vulnerability to climate change according to .

NOM-59, Official Mexican Standard; Pr, Subject to special protection; IUCN, International Union for Conservation of Nature; LC, Least Concern; NT, Near Threatened. Vulnerability: M, Medium; H, High.

3.1 Change rate

3.1.1 Area

The species with the largest distribution area is D. mexicanus (72,469 km²), followed by D. oaxacae (8,522 km²), and finally G. syntrema (3,383 km²). When deforestation is accounted for, the distributions are reduced by nearly half (except for G. syntrema). No species showed an expansion of its distribution area under the different scenarios (SSP1-2.6 for 2060 and 2100, and SSP5-8.5 for 2060 and 2100); only reductions were observed. Finally, only G. syntrema is expected to experience local extinction in Chiapas under the SSP5-8.5 scenario by the end of the present century (Figure 2).

Figure 2

Significant differences were detected among scenarios in the percentage of distribution area (Friedman χ² = 11.571, p = 0.041). However, pairwise comparisons using Wilcoxon tests with Bonferroni adjustment did not detect significant differences between specific scenarios (p > 0.05 in all cases; Figure 3).

Figure 3

3.1.2 Altitude

Regarding the altitudinal profile obtained from the occurrence records for each species, D. mexicanus shows a broad distribution range, from low elevations to nearly 3,000 m. In contrast, D. oaxacae is restricted to intermediate elevations (≤ 2,000 m), whereas G. syntrema is distributed mainly in lowland areas (< 1,000 m; Figure 4).

Figure 4

On the other hand, based on the altitudinal ranges derived from current potential distributions and future scenarios, D. mexicanus exhibited the broadest altitudinal range, extending from lowland areas to nearly 4,000 m in its current distribution. However, under the SSP1-2.6 (2060 and 2100) and SSP5-8.5 (2060) scenarios, its upper altitudinal limit is reduced to approximately 3,000 m. In contrast, under the SSP5-8.5 (2100) scenario, a recovery of this range is observed, once again reaching altitudes close to current conditions. In contrast, D. oaxacae showed a more restricted distribution (approximately up to 2,600 m), with an expansion of its altitudinal range under the optimistic scenario (SSP1-2.6). However, a slight reduction of its upper altitudinal limit is observed under more severe climate scenarios (SSP5-8.5).

Meanwhile, G. syntrema was restricted to low elevations (below ~1,400 m), showing contraction in its altitudinal range across all scenarios, and it even disappears under the SSP5-8.5 scenario by the end of the century (Figure 5). Although changes in altitudinal limits were observed across scenarios, these were not statistically significant (Friedman χ² = 3.41, p = 0.64), suggesting no consistent altitudinal pattern among species across the evaluated scenarios (Figure 6).

Figure 5

Figure 6

3.2 Species protection index

The SPI indicates that D. mexicanus has less than 30% of its potential distribution within PNAs under all scenarios. In the case of D. oaxacae, the index shows a decreasing trend; currently, 53% of its distribution occurs within PNAs, decreasing to as low as 33% under some scenarios. However, under the SSP5-8.5 scenario, this value partially recovers to 48% by 2100. Despite this recovery, its total distribution is reduced by approximately half. Finally, G. syntrema maintains more than 70% of its distribution within PNAs in most scenarios, reaching nearly 100% under the SSP5-8.5 scenario in 2060, with its distribution becoming almost entirely restricted to the Montes Azules Biosphere Reserve (Figure 7). However, its total distribution decreases by more than 50% relative to its current range, and it disappears entirely under the SSP5-8.5 scenario by the end of the century (Figure 8).

Figure 7

Figure 8

Regarding the global SPI, the current potential distributions of caecilian species show a relatively high value (58%), which is maintained in three of the future scenarios, with a slight decrease reaching 51%. However, under the most catastrophic scenario, protection decreases considerably, falling to 27% by the end of the century (Figure 9).

Figure 9

4 Discussion

The occurrence of the three caecilian species known from Mexico in Chiapas may be associated with the state’s environmental conditions, particularly its high soil moisture and humid microhabitats, which are essential for caecilians because of their physiological dependence on moist environments (). Changes in temperature and precipitation regimes under future climate scenarios could alter these conditions and consequently affect the distribution and persistence of these species.

Considering the climate niche tracking theory, we observe that the potential distribution of caecilians is consistent with that reported for other amphibians, in which climate is a determining factor in geographic distribution. Several studies have shown that variations in temperature and precipitation alter habitat suitability, leading to spatial and altitudinal shifts in species (; ; ).

In this context, some species will be more affected than others. For example, species with a wide distribution, such as D. mexicanus, tend to exhibit greater stability in the face of CC scenarios due to their broader ecological tolerance and generalist habits (; ). Conversely, species with restricted distribution (such as D. oaxacae and G. syntrema) show greater sensitivity to environmental changes, increasing their vulnerability and the risk of contraction or total loss of their climatic niche ().

In the specific case of G. syntrema, the total loss of suitable conditions under the SSP5-8.5 scenario by 2100 suggests a possible collapse of its climatic niche. This is consistent with global predictions indicating that CC can lead to the local disappearance or even extinction of species unable to adapt or disperse quickly enough to new areas (). This situation in caecilians could be more critical since their fossorial habits and dependence on humid soil conditions limit their dispersal capacity (), making them more susceptible to rapid CCs.

Taken together, these results suggest that although there is a general trend toward change in spatial distribution (p= 0.041), the responses are not sufficiently differentiated between scenarios to yield statistically clear pairwise contrasts (p > 0.05). Moreover, altitudinal ranges may remain relatively stable over the short and medium term. This reinforces the idea that the effects of CC on these species may be progressive, cumulative, and context-dependent, rather than immediate or drastic.

Regarding the conservation mismatch hypothesis, the results obtained show a shift between the potential distribution of species and the coverage of PNAs. In the case of the genus Dermophis, although part of its current distribution is within PNAs, future models suggest that a significant proportion of its suitable climatic niche will fall outside these areas, and that this mismatch increases as CC scenarios progress, especially in the catastrophic scenario with high emissions at the end of the century (Figure 8). This pattern was also described by , who determined that, while protected areas maintain some capacity to support biodiversity, their overall effectiveness declines over time. They attribute this decline primarily to reduced habitat suitability and climate connectivity, especially under the high-emissions scenario. In this sense, the shift in the climatic niche results in many species being excluded from protected areas, reducing their effectiveness as a long-term conservation strategy.

In the case of G. syntrema, the pattern is even more critical. Although it generally maintains a high proportion of its current potential distribution (79.8%) across most scenarios (ranging from 85.9% to 99.1%), this does not necessarily indicate a favorable condition. Rather, it can be interpreted as an extreme contraction of its climatic niche, in which the species is restricted to a very small area (the Montes Azules Biosphere Reserve). Under the most severe scenario (SSP5-8.5), the disappearance of suitable conditions towards the end of the century suggests a possible process of local extinction, even within a PNAs. Furthermore, the global SPI shows this pattern more clearly, declining from approximately 58% protection in the current potential distribution to 27% by the end of the century under the SSP5-8.5 scenario. Therefore, PNAs do not guarantee the persistence of species under CC, especially for organisms with low dispersal capacity and high ecological specialization, such as amphibians (; ; ; ).

Taken together, these results suggest that the conservation of these species cannot depend exclusively on current PNAs, but requires dynamic strategies, such as the identification of climate refugia, landscape connectivity, and the integration of CC into conservation planning. Otherwise, even areas with high protection coverage could become insufficient to prevent biodiversity loss.

On the other hand, projections generated by species distribution models should be interpreted with caution, as they represent potential estimates of future climate suitability, rather than definitive predictions of species distribution. Since the models are based on correlations between presence records and environmental conditions, their results are inherently associated with uncertainties related to model assumptions, input data, and future climate scenarios (; ; ; ).

Finally, this work represents an initial effort to assess the potential distribution of caecilians in Chiapas and their response to CC. The results suggest that D. mexicanus could maintain relative stability, whereas species with restricted distributions in Chiapas, such as D. oaxacae and G. syntrema, would be more vulnerable, potentially losing their suitable climatic niche by the end of the century under the SSP5-8.5 scenario. Furthermore, the mismatch between the projected distribution and protected areas suggests that current strategies may be insufficient, highlighting the need for more dynamic conservation approaches.

Statements

Data availability statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author contributions

RC-H: Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Writing – original draft, Writing – review & editing. MG-T: Formal analysis, Writing – original draft, Writing – review & editing. JRC-V: Writing – original draft, Writing – review & editing. SL: Conceptualization, Formal analysis, Investigation, Methodology, Supervision, Validation, Writing – original draft, Writing – review & editing.

Funding

The author(s) declared that financial support was received for this work and/or its publication. This research was supported by the Secretariat of Science, Humanities, Technology and Innovation through the Estancias Posdoctorales por México Program, 2025 Call, awarded to RCH (CVU 550333).

Acknowledgments

The author RCH gratefully acknowledges support from the Secretaría de Ciencia, Humanidades, Tecnología e Innovación (SECIHTI) through the Postdoctoral Fellowships Program for Mexico, 2025. He also thanks the Universidad Autónoma de Ciencias y Artes de Chiapas and the Laboratory of Evolutionary Ecology for providing access and logistical support during the development of this research.

Conflict of interest

The author(s) declared that this work was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The author JC-V declared that they were an editorial board member of Frontiers, at the time of submission. This had no impact on the peer review process and the final decision.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Summary

Keywords

caecilians, climatic niche tracking, conservation planning, protected areas effectiveness, species distribution modeling

Citation

Cabrera-Hernández R, Gómez-Tolosa M, Cedeño-Vázquez JR and López S (2026) Climate change and conservation of caecilians (order Gymnophiona): prediction of distributional shifts and risk assessment in Chiapas, Mexico. Front. Ecol. Evol. 14:1886733. doi: 10.3389/fevo.2026.1886733

Received

20 May 2026

Revised

05 July 2026

Accepted

20 July 2026

Published

13 August 2026

Volume

14 - 2026

Edited by

Pablo Martinez, Federal University of Sergipe, Brazil

Reviewed by

Luis Amador, Auburn University, United States

Feiyun Tu, Hainan Normal University, China

Updates

Copyright

*Correspondence: Rodolfo Cabrera-Hernández, ; María Gómez-Tolosa, ; Sergio López,

Disclaimer

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.

Outline

Figures

Cite article

Copy to clipboard


Export citation file


Share article

Article metrics