Climate change is significantly altering species distributions worldwide, expanding the range of some and reducing suitable habitats for others. Concerning snakes, these changes have direct public health implications, as venomous snake species may invade new areas, including densely populated urban zones. This increases the need for appropriate antivenom distribution strategies, particularly in remote areas or those displaying limited healthcare infrastructure. While several studies have focused on climate change snake habitat effects, few have addressed how this affects snakebite treatment center logistics and organization. In this context, we aimed to assess if snake distribution changes driven by climate change require the restructuring of antivenom treatment center networks (ATCN) in Brazil. To this end, we applied occurrence records from herpetological collections alongside bioclimatic and land use variables to species distribution models employing five different algorithms considering projected 2061–2080 climate change scenarios under mitigated (SSP 2−4.5) and unchecked (SSP 5–8.5) conditions. The developed models indicate suitable snake area decreases, as well as species richness, particularly under the SSP 5–8.5 scenario. We identified a spatial mismatch, with suitable snake areas displaying up to 63% ATCN coverage deficit, particularly in extremely remote regions, especially for the Bothrops and Micrurus genera. In contrast, ATCNs carrying antivenom, particularly for the Lachesis and Crotalus genera, were noted in up to 92% of unsuitable snake areas. These findings highlight the need for constant ATCN location reassessments considering future climate change scenarios to ensure more effective snakebite incident responses, also emphasizing the importance of adopting an integrated health approach considering human, animal, and environmental health interactions to improve snakebite risk management in Brazil.
Many aspects of living organism are associated with their climatic niches, making climate change an important biodiversity pressure factor (Tozato et al., 2015). In this sense, changes in temperature and rainfall patterns and extreme weather events threaten not only ecosystems and species but also human survival (IPCC, 2022; Werneck et al., 2023). Furthermore, climate change also intensifies habitat degradation and losses, mainly driven by deforestation (Myers et al., 2000; Werneck et al., 2023).
About 28% of all animal species and plants are currently threatened, and over 1 million are at risk of extinction, with Latin America experiencing some of the steepest animal population declines since the 1970s (Myers et al., 2000; Werneck et al., 2023; IUCN, 2024). Due to their dependence on environmental temperature to regulate physiological functions, snakes and other ectothermic animals are particularly vulnerable to these changes, both directly, i.e., thermal stress or indirectly, due to disease outbreaks and decreased food availability (Winter et al., 2016; Lourenço-de-Moraes et al., 2019; Werneck et al., 2023).
Previous studies have demonstrated that climate change can drive substantial alterations in snake distributions worldwide (Martinez et al., 2024). In general, some venomous species are projected to shift or expand their climatically suitable ranges under future scenarios, potentially increasing interactions with human-modified landscapes (Werneck et al., 2023; Huang et al., 2024). In contrast, several forest-associated snake species in Brazil are projected to experience reductions in suitable habitat under climate change scenarios, including taxa occurring in the Atlantic Forest and Amazonian regions (Lourenço-de-Moraes et al. 2019), which may lead to population declines and increased extinction risks (Vasconcelos, 2014; Martinez et al., 2024).
Observed and predicted distribution shifts pose an additional public health challenge, as venomous snakes may expand their geographic ranges and occur in novel areas, including peri-urban environments and urban green spaces associated with remnant natural habitats. Such shifts may broaden the geographic scope of snakebite incidents, particularly in regions where human populations interface with fragmented natural landscapes (Gutiérrez et al., 2017; Matos and Ignotti, 2020). These changes have direct implications for the strategic management and equitable distribution of snake antivenom, particularly in remote or resource-limited regions where healthcare infrastructure is often inadequate (Yousefi et al., 2020; Guerra et al., 2023; Chowdhury et al., 2024).
The World Health Organization (WHO) estimates about 5.4 million snakebite incidents annually worldwide, causing 1.8–2.7 million envenomation cases due to venomous snakes and resulting in about 81,000–138,000 deaths (WHO, 2023). In Brazil, about 29,000 snakebite incidents take place each year, averaging 139 deaths, with about 86% of all cases not specifying the responsible snake genus (SINAN, 2024). These high incidence rates largely reflect human movement into natural habitats, increased human activities in previously sparsely populated areas, and greater reporting and access to healthcare, rather than shifts in the snake species ranges themselves.
Seventy-eight venomous snake species belonging to four genera occur in Brazil, namely Bothrops lato sensu (lanceheads/jararacas), responsible for 86.5% of all snakebite cases in the country, Crotalus (rattlesnakes/cascavéis), responsible for 9.7%, Lachesis (bushmasters/surucucus), accounting for 2.8%, and Micrurus (coral snakes/corais-verdadeiras), causing 1.1% of cases (Uetz et al., 2024; SINAN, 2024).
Envenomation symptoms and physiological effects can range from pain and swelling at the bite site to necrosis, hemorrhage, paralysis, and respiratory failure, depending on the snake genus and the amount of injected venom (Pinho et al., 2000; Pardal et al., 2007; Casais-e-Silva and Brazil, 2009; Luciano et al., 2009). Treatment effectiveness is directly associated to the quick administration of the correct antivenom (Lira-da-Silva et al., 2016).
Since 2014, Brazilian laboratories producing antivenoms have had to adjust their production facilities, limiting production capacity (ANVISA, RDC 17/2010; Schneider et al., 2021) and consequently reducing the number of served hospitals and municipalities, primarily affecting patient care response times (Salomão et al., 2018; Guerra et al., 2019).
Climate change is projected to affect the potential distributions of venomous snakes (Martinez et al., 2024), particularly in tropical zones, where models suggest higher risks of species losses by 2070 in South American and southern African countries. This is alarming, as venomous snakes play both essential ecological roles and their venoms hold enormous potential for pharmaceutical and biotechnological development (Martinez et al., 2024).
While most studies assess how climate change effects affect the expansion or reduction of suitable venomous snakes habitats (Guerra et al., 2023; Martinez et al., 2024) or how treatment center accessibility affects snakebite incident response times (Oliveira RAD et al., 2022), none have focused on investigating how current antivenom treatment center networks (ATCN) are spatially organized in relation to known present-day venomous snakes distributions, whether antivenom availability reflects current species ranges or historical distributions, or assessing how climate change may alter this spatial alignment.
In this sense, this study aimed to evaluate whether potential venomous snake distribution shifts due to climate change may require ATCN restructuring. We hypothesize that the current ATCN distribution in Brazil does not cover potential snakebite incident risk areas, both in the present and in the future, requiring spatial distribution ATCN reassessments.
Material and methodsOccurrence records and layersWe requested access to digital records from the following herpetological collections: National Museum of Rio de Janeiro (Museu Nacional do Rio de Janeiro, MNRJ; Rio de Janeiro), Institute of Biology, Federal University of Rio de Janeiro (Instituto de Biologia, Universidade Federal do Rio de Janeiro, ZUFRJ; Rio de Janeiro – transferred to MNRJ), National Institute of the Atlantic Forest (Instituto Nacional da Mata Atlântica, MBML-INMA; Espírito Santo), Laboratory of Vertebrate Zoology, Federal University of Ouro Preto (Laboratório de Zoologia de Vertebrados, Universidade Federal de Ouro Preto, LZV-UFOP; Minas Gerais), and the Miguel Lillo Foundation (Fundación Miguel Lillo, FML; Tucumán, Argentina). A remote requisition was required due to the 2021 and 2022 Covid-19 pandemic restrictions.
The database generated by Nogueira et al. (2019), who conducted a comprehensive mapping of all snake species, venomous and non-venomous, with at least one confirmed record in Brazil, was also included. All specimens examined by those authors were identified to the species level by specialists, ensuring identification accuracy. Some considered snake collections from widely recognized scientific community institutions, such as the Butantan Institute (Instituto Butantan, IBSP), Vital Brazil Institute (Instituto Vital Brazil, IVB), Museum of Zoology, University of São Paulo (Museu de Zoologia, Universidade de São Paulo, MZUSP), and the Emílio Goeldi Museum of Pará (Museu Paraense Emílio Goeldi, MPEG).
We excluded individuals with no known sampling location, as well as captive specimens, records presenting duplicate coordinates and/or with no location indications and those presenting identification errors from the analysis. The inability to directly evaluate Bothrops lato sensu and Micrurus specimens, which include several species presenting similar morphologies, required consideration of the identifier's reliability and the sampling location.
Regarding Lachesis, none of the compiled occurrence records included collection dates, preventing distinctions between recent and historical records and increasing the risk of temporal mismatch with contemporary environmental predictors. To address this limitation and minimize overprediction, we adopted a conservative data curation approach grounded in historical biogeography and previous distributional assessments, which consistently indicate that L. muta is restricted to continuous tropical forest formations, with southern and southwestern limits associated with the Amazon basin and the northern Atlantic Forest, and no evidence of stable populations in Cerrado-dominated regions or in the southern Atlantic Forest (e.g., Campbell and Lamar, 2004; Moura et al., 2016). Following the same procedure adopted by Citeli et al. (2020), we excluded occurrence records from southeastern Brazil, as well as from the states of Goiás and Mato Grosso outside the Amazon domain, from the analyses to reduce the inclusion of historical, mislocated, or ecologically implausible records that could inflate model predictions and generate false positives.
A complete list of all species included and excluded from the analyses, together with the number of records per species and relevant notes, is provided in the Supplementary Material (Table S1). To document the spatial and temporal structure of the occurrence dataset, a map of all raw records is provided in the Supplementary Material (Fig. S1), depicting occurrences following initial spreadsheet-based cleaning and prior to any spatial filtering or processing in the R environment. We categorized occurrence records according to the availability and period of collection dates, distinguishing records without collection year information, historical records (≤1949), and more recent records (≥1950), following the temporal criteria proposed by Guerra et al. (2023).
We obtained bioclimatic variables from the WorldClim database (Hijmans et al., 2005 – https://worldclim.org/data) at a 5 arcminute resolution, corresponding to about 10 × 10 km pixels, and land use variables from the Land Use Land Cover (LULC) database (www.usgs.gov/land-resources/eros/lulc), at a 25 km resolution. To ensure consistent dataset comparisons, we adjusted LULC variables to match the bioclimatic variable resolution employing the nearest neighbor method.
Given the 5 arc-minute (∼10 km) resolution of the environmental predictors, we applied a spatial thinning distance of 10 km to occurrence records to reduce spatial autocorrelation and avoid over-representation of densely sampled localities within single grid cells. This was not intended to fully correct for sampling bias related to uneven survey effort, but rather to minimize pseudo-replication and improve model stability and generalization at the spatial predictor resolution. While spatial thinning alone does not eliminate all sources of sampling bias, it is a widely adopted and conservative approach to reduce local clustering effects and prevent model overfitting, particularly when working with coarse-resolution environmental layers.
Regarding Bothrops and Micrurus genera modeling, we modeled each species individually, as more than one species occurs in Brazil, followed by a compilation of all species-level models into a single genus map. This was not necessary for Crotalus, as only one species occurs in the country.
Regarding Lachesis, a recent study recognized that populations with disjunct distributions that were previously considered a single species are actually two distinct species, one from the Amazon and the other endemic to the Atlantic Forest biome (Hamdan et al., 2024). However, the occurrence records employed herein were obtained prior to this discovery, meaning that the model was generated based on all occurrence Lachesis genus records from both Amazon and Atlantic Forest locations. A thorough taxonomic analysis of all individuals deposited in the accessed zoological collections would be necessary to separate the database and apply new models, requiring time and funding. Given that both species occur in forests, with no habitat overlap, and the fact that the model did not consider bioclimatic variables that limit the distribution ranges of each species, i.e., precipitation of the coldest quarter for the Amazon species and precipitation of the driest quarter for Atlantic Forest species (Hamdan et al., 2024), the model based on all Lachesis records was maintained. More importantly, this joint modeling approach has no implications for healthcare planning or antivenom access. Envenomations caused by both Lachesis species are clinically indistinguishable and are treated with the same genus-specific antivenom. Consequently, the relevant operational unit for this study is antivenom treatment demand rather than species-level taxonomic resolution. Thus, jointly modeling both species does not affect the ATCN coverage nor the public health risk interpretation.
We did not model Bothrops and Micrurus species with fewer than 80 initial occurrence records, as the application of spatial thinning substantially reduces the number of effective records available for calibration. A minimum initial sample size was therefore required to ensure sufficient post-thinning occurrences for reliable model fitting. We applied a 10-kilometer spatial thinning to all occurrence records to reduce spatial autocorrelation, minimize sampling bias, and improve model generalization to unsampled areas.
Our environmental variable selection followed a species-specific approach within a standardized modeling framework. For each species, we pre-selected an initial pool of bioclimatic and land-use variables based on prior ecological knowledge, including habitat associations, climatic tolerances, and known distributional constraints, following approaches commonly adopted in recent ecological niche modeling studies (e.g., Tardin et al., 2025). To reduce multicollinearity and potential model overfitting, we subsequently filtered predictor variables separately for each species employing Pearson correlation analyses and variance inflation factor (VIF) calculations implemented in the R Studio software v. 2023.12.1, package USDM (script available as a separate file in the Supplementary Material). We iteratively removed variables showing high pairwise correlation (|r| ≥ 0.7) or high multicollinearity (VIF ≥ 3). We selected these thresholds based on widely accepted recommendations in ecological modeling literature to ensure model stability and interpretability (Zuur et al., 2010; Dormann et al., 2013).
The inclusion of land-use and land-cover (LULC) variables in ecological niche modeling involves known methodological challenges, including the categorical nature of some predictors, potential increases in model complexity, and temporal or spatial mismatches with occurrence data. These limitations were explicitly considered during the model design. Nevertheless, LULC predictors were included because venomous snake species distributions are strongly mediated by habitat structure, forest cover, and anthropogenic landscape transformation, which are not fully captured by climatic variables alone. Given the applied objective of evaluating spatial mismatches between snake suitability and antivenom treatment center coverage, climate-only models would provide an incomplete representation of ecological and public health risk.
No predictors were exempt from the collinearity filtering process. Although we applied the same selection criteria across all species, the final set of retained environmental variables differed among species according to their ecological requirements and statistical relationships among predictors.
Modeling protocolWe employed five modeling algorithms implemented in the BIOMOD2 package v. 4.2-4 within the RStudio platform v. 2023.12.1, namely Generalized Linear Models (GLM), Generalized Additive Models (GAM), Random Forest (RF), Gradient Boosting Models (GBM), and Maximum Entropy (MaxEnt). Algorithms differed in their treatment of occurrence data, using presence–pseudoabsence approaches (GLM, GAM, RF, GBM) or presence–background data (MaxEnt) (Thuiller et al., 2020). Detailed information on model settings, tuning parameters, and calibration procedures for each algorithm is provided in the modeling script available as a separate file in the Supplementary Material.
We delimited the model calibration area (M) in a species-specific manner. For each species, we defined M based on its known geographic distribution polygon obtained from the International Union for Conservation of Nature (IUCN) Red List spatial database. To account for uncertainties in distribution boundaries and potential short-term range expansions, we applied a buffer of 1 ° around each species distribution polygon using ArcGIS v. 10.8 (ESRI). We then used these buffered polygons as the calibration extent for model fitting and for the random generation of background points and pseudo-absence data, thereby restricting model calibration to ecologically and biogeographically realistic areas while minimizing excessive extrapolation.
For MaxEnt, we randomly generated 10,000 background points within the species-specific calibration area (M). For GLM, GAM, RF, and GBM, we randomly generated 500 pseudo-absence points within the calibration area while enforcing a minimum distance of 30 km from known occurrence records, to reduce spatial autocorrelation and minimize the risk of selecting false absences (Barbet-Massin et al., 2012). This approach improves model discrimination by avoiding the inclusion of environmentally suitable locations near presences as absences, thereby contributing to a more robust model calibration. The minimum distance criterion also mitigates the effects of spatial sampling bias and reduces the risk of overfitting, particularly when occurrence records are spatially clustered. Although more complex bias-adjusted pseudo-absence strategies are available, random pseudo-absence generation within the accessible area has been shown to provide robust and comparable results across species, especially when combined with spatial constraints (Barbet-Massin et al., 2012). We consistently applied the same procedure to all modeled species.
We calibrated the models using a repeated 5-fold cross-validation procedure implemented in biomod2, in which the dataset was randomly partitioned into five subsets for iterative training and testing. This cross-validation procedure was repeated 10 times for each algorithm, resulting in multiple replicate model evaluations, with prevalence fixed at 0.5 (Guisan et al., 2017). We assessed model performance using both the area under the Receiver Operating Characteristic curve (AUC) and the True Skill Statistic (TSS), as implemented in the BIOMOD2 framework and calculated these complementary metrics for all individual models. We constructed ensemble models by retaining only those models that met minimum performance thresholds for both metrics (AUC ≥ 0.8 and TSS ≥ 0.7), ensuring robust and conservative predictions.
We generated both continuous suitability maps and binary ensemble predictions, where binarization was performed using the threshold that maximizes the sum of sensitivity and specificity (Liu et al., 2013), applied consistently across species to ensure comparability among models. We retained continuous suitability maps for complementary analyses to preserve information on ecological gradients and species-specific responses.
We then projected distributions of each snake species under two different climate change scenarios, namely SSP 2−4.5, which represents an intermediate mitigation scenario where current trends continue and moderate measures are applied to reduce greenhouse gas emissions, and SSP 5–8.5, characterized by high emissions, accelerated economic growth based on intensive fossil fuel use and limited mitigation efforts (Riahi et al., 2017). We applied the Model for Interdisciplinary Research on Climate (MIROC6) climate model covering from 2061 to 2080 to all projections.
Post-processingWe post-processed the binary maps using raster algebra operations to quantify habitat gains, losses, and genus-level species richness. Specifically, we subtracted binary suitability rasters generated from the combined effects of climatic and land-use and land-cover (LULC) variables (future minus current scenarios) to identify spatial patterns of habitat loss or gain. Subsequently, we summed suitable areas across species belonging to the same genus, for both current and future scenarios, to generate genus-level species richness maps. We implemented post-processing steps ed using the Raster Calculator tool in ArcGIS v. 10.8. However, all operations correspond to standard software-independent raster algebra procedures reproducible in any GsS or spatial analysis platform.
We assessed spatial mismatches between suitable snake genera areas and Brazilian ATCNs carrying antivenom, where the modelled spatial predictions, both current and future, were overlaid with a georeferenced and updated list available at the Brazilian Ministry of Health's virtual portal (Oliveira RAD et al., 2022). When necessary, we obtained geographic hospital coordinates at www.coordenadas-gps.net.This approach allowed for the quantification of (i) suitable areas for both the current and future climate scenarios (in km2) for venomous snakes where no ATCNs are available and (ii) the number of ATCNs located in unsuitable snake areas.
To estimate the ATCN coverage area within a maximum timeframe of 2 h, we applied routes provided by Google Maps between randomly selected municipalities based on different accessibility levels defined by the Brazilian Institute of Geography and Statistics (Instituto Brasileiro de Geografia e Estatística - https://www.ibge.gov.br/).
The accessibility levels defined by the Brazilian Institute of Geography and Statistics (IBGE) are based on the Geographic Accessibility Index, a national-scale territorial typology that classifies Brazilian municipalities according to their degree of geographic access to urban centers of higher hierarchical level. This index considers road and fluvial connectivity and distances between municipal seats and regional service hubs, reflecting the relative ease with which populations can reach essential services such as healthcare. Municipalities are categorized into four classes, namely very remote, remote, accessible, and highly accessible, representing increasing levels of geographic connectivity and infrastructure availability. Herein, we used these IBGE accessibility classes as a proxy for travel conditions to health facilities and they served as the basis for defining differentiated travel distances and buffers around antivenom treatment centers.
We calculated the distances and travel times used to estimate antivenom treatment center (ATCN) coverage areas using routes generated by Google Maps between randomly selected municipalities within each IBGE accessibility class. We based the travel time estimates d on the fastest available routes and relied primarily on the existing road network, as standardized and spatially explicit datasets integrating road and river transportation at the national scale are currently unavailable. Although fluvial transport is relevant in some remote regions of Brazil, particularly in the Amazon basin, its effects on accessibility were indirectly accounted for through the Geographic Accessibility Index defined by the Brazilian Institute of Geography and Statistics (IBGE), which incorporates territorial remoteness and connectivity constraints. We calculated average travel speeds for each accessibility class based on the relationship between estimated distances and travel times, which were subsequently used to define class-specific buffers around ATCNs corresponding to a maximum travel time of two hours.
From these calculations, we defined buffers around each ATCN reachable in at most two hours as follows: 40 km for very remote regions, 80 km for remote areas, 100 km for accessible locations, and 120 km for highly accessible regions. Once we estimated ATCN coverage areas for each accessibility level, we calculated suitable snake areas not covered by ATCNs, as well as unsuitable areas covered by ATCNs.
ResultsWe compiled a total of 17,727 occurrence records for the investigated venomous snake species, resulting in the modeling of 33 species across four genera: Bothrops (18 species), Micrurus (12 species), Lachesis (two species), and Crotalus (one species). All ensemble models displayed satisfactory predictive performance, indicating reliable discrimination between suitable and unsuitable areas (Table S2 – Supplementary Material).
Overall, the model projections revealed a consistent trend of reduction in suitable areas under future climate scenarios across all four genera, despite limited range expansion observed for a small number of species (Figures S2–S4 – Supplementary Material). These projected changes translated into pronounced spatial mismatches between climatically suitable snake areas and the current ATCN distribution, a pattern that is expected to intensify under worsening climate scenarios (Figures S5–S8 – Supplementary Material). Detailed information on the geographic distribution of antivenom types across Brazilian health units, including the availability of bothropic, crotaline, lachetic, and elapidic antivenoms, is provided in Table S3 (Supplementary Material). These patterns are synthesized in a set of multipanel figures (Figures S5–S8 – Supplementary Material), illustrating current and future mismatches between snake suitability and antivenom treatment centers across accessibility levels and genera.
At present, ATCNs are unevenly distributed across Brazil, with a strong concentration in densely populated regions (Figure S9 – Supplementary Material). This spatial configuration results in substantial gaps in coverage within suitable areas for medically important snakes, particularly in northern, northeastern, and central-western Brazil (Figure S5-S8 – Supplementary material). These regions consistently emerged as areas of elevated vulnerability across both current and future climate scenarios.
Currently, suitable Bothrops genus areas uncovered by ATCNs equipped with antivenom (SAB) range from 6% in very accessible areas to 63% in very remote areas, particularly in Northern Brazil areas located farther from access routes (Fig. 1), similar to the predictions for both future scenarios (SSP 2−4.5 and SSP 5–8.5). Likewise, the number of ATCNs carrying SAB located in unsuitable areas will gradually increase with decreasing accessibility, both in the current and future climate scenarios, especially in remote areas (Fig. 2).
Areas not covered by health units carrying bothropic antivenom (SAB) considering accessibility levels defined by IBGE (columns) and different climate scenarios (rows) for the genus Bothrops. The color gradient from green to red indicates the number of species occurring simultaneously in the same location. Darker red areas represent higher species richness. The maximum species limit varies according to scenario and Brazilian region. Gray areas represent unsuitable genus occurrence habitats. A reference map of Brazil depicting states and official macroregions is provided in the Supplementary Material (Figure S10).
Coverage by health units carrying bothropic antivenom (SAB) considering accessibility levels defined by IBGE (columns) and different climate scenarios (rows) for the genus Bothrops. The buffers around each health unit were defined as 40 km for very remote areas, 80 km for remote areas, 100 km for accessible areas, and 120 km for very accessible areas. The color gradient from green to red indicates the number of species occurring simultaneously in the same location. Darker red represent higher species richness. The maximum species limit varies according to scenario and Brazilian region. Gray areas represent unsuitable genus occurrence habitats.
Considering all accessibility levels, the highest average ATCN deficit is currently observed in suitable Micrurus areas, ranging from 11% in very accessible areas to over 55% in very remote regions. This remains consistent for the future climate scenarios, with a slight reduction projected for the SSP 5–8.5 scenario (Fig. 3). Although less frequent, accidents involving the Micrurus genus are noted throughout Brazil, with ATCNs carrying SAE located in unsuitable areas ranging from 2% to 44% in future climate scenarios. The projections indicate vast healthcare hub areas that will become unsuitable concerning Micrurus in the future, particularly in Northern, Northeastern, and Central-Western Brazil (Fig. 4).
Areas not covered by health units carrying elapid antivenom (SAE) considering accessibility levels defined by IBGE (columns) and different climate scenarios (rows) for the genus Micrurus. The color gradient from green to red indicates the number of species occurring simultaneously in the same location. Darker red areas represent higher species richness. The maximum species limit varies according to scenario and Brazilian region. Gray areas represent unsuitable genus occurrence habitats.
Coverage area of health units carrying elapid antivenom (SAE) considering accessibility levels defined by IBGE (columns) and different climate scenarios (rows) for the genus Micrurus. The buffers around each health unit were defined as 40 km for very remote areas, 80 km for remote areas, 100 km for accessible areas, and 120 km for very accessible areas. The color gradient from green to red indicates the number of species occurring simultaneously in the same location. Darker red ares represent higher species richness. The maximum species limit varies according to scenario and Brazilian region. Gray areas represent unsuitable genus occurrence habitats.
Suitable Lachesis genus areas uncovered by ATCN carrying SAL range from 4% in very accessible areas to 7% in very remote areas, following a similar general trend for the future SSP 2−4.5 and 5–8.5 climate scenarios (Fig. 5). Considering all accessibility levels, ATCNs carrying SAL located in unsuitable areas average 84%, with very remote areas being the least affected. This mismatch tends to intensify as climate scenarios become more severe (Fig. 6).
Areas not covered by health units carrying lachetic antivenom (SAL) considering accessibility levels defined by IBGE (columns) and different climate scenarios (rows) for the genus Lachesis. Suitable Lachesis occurrence areas are depicted in red. The maximum species limit varies according to scenario and Brazilian region. Gray areas represent unsuitable genus occurrence habitats.
Coverage area of health units carrying lachetic antivenom (SAL) considering accessibility levels defined by IBGE (columns) and different climate scenarios (rows) for the genus Lachesis. The buffers around each health unit were defined as 40 km for very remote areas, 80 km for remote areas, 100 km for accessible areas, and 120 km for very accessible areas. Suitable Lachesis occurrence areas are depicted in red. The maximum species limit varies according to scenario and Brazilian region. Gray areas represent unsuitable genus occurrence habitats.
On average, only 1% of suitable Crotalus areas are not covered by ATCNs carrying SAC considering all accessibility categories in both current and future climate scenarios (Fig. 7). In contrast, unsuitable areas covered by ATCNs represent the highest coverage, averaging 31% in the current scenario and reaching up to 37% in the most extreme climate scenario. This discrepancy is especially evident in remote areas, particularly in the Brazilian Amazon, Northeast, and Central-West (Fig. 8).
Areas not covered by health units carrying crotalic antivenom (SAC) considering accessibility levels defined by IBGE (columns) and different climate scenarios (rows) for the genus Crotalus. Suitable Crotalus occurrence areas are depicted in red. The maximum species limit varies according to scenario and Brazilian region. Gray areas represent unsuitable genus occurrence habitats.
Coverage area of health units carrying crotalic antivenom (SAC) considering accessibility levels defined by IBGE (columns) and different climate scenarios (rows) for the genus Crotalus. The buffers around each health unit were defined as 40 km for very remote areas, 80 km for remote areas, 100 km for accessible areas, and 120 km for very accessible areas. Suitable Crotalus occurrence areas are depicted in red. The maximum species limit varies according to scenario and Brazilian region. Gray areas represent unsuitable genus occurrence habitats.
Projections for future climate scenarios indicate that these mismatches will not only persist but intensify, with reductions in suitable areas occurring concurrently with an increasing proportion of ATCNs located in climatically unsuitable regions. This pattern suggests a growing misalignment between snakebite risk and healthcare infrastructure under climate change.
Variable importance analyses indicated that climatic factors played a stronger role than land-use variables in shaping projected suitability patterns. Isothermality and precipitation-related variables consistently ranked among the most influential predictors across species, whereas anthropized land-use variables showed moderate contributions, indicating a limited tolerance of several species to altered environments rather than widespread adaptation (Table S4 – Supplementary material).
DiscussionThis study reveals complex and multifactorial dynamics that extend beyond changes in venomous snake distributions driven by climate change alone. While projected losses of species richness and climatically suitable areas are consistent with previous ecological forecasts, the key contribution of this work lies in demonstrating how these shifts may reorganize snake distribution and exacerbate existing mismatches with the Brazilian ATCN network. By explicitly linking ecological suitability, geographic accessibility, and healthcare infrastructure, our findings expose a critical misalignment between antivenom supply and ecological risks, an imbalance that is projected to intensify under future climate scenarios.
Across both current and future conditions, insufficient ATCN coverage persists in climatically suitable snake areas, particularly in North, Northeast, and Central-West Brazil, where the Amazon and Cerrado biomes are predominant. At the same time, a substantial proportion of ATCNs remain in areas projected to become increasingly unsuitable for medically important snake genera, indicating public health resource allocation inefficiencies. Although Crotalus shows a relatively higher overlap between suitable areas and ATCN coverage, a considerable number of treatment centers remain positioned in unsuitable regions, highlighting that even apparently favorable scenarios mask structural distribution problems.
These spatial mismatches carry important public health consequences. The absence of ATCNs in suitable snake areas disproportionately affects vulnerable populations, especially in remote and socioeconomically disadvantaged regions. Longer travel times between snakebite incidents and medical care are likely to increase the risk of severe complications and mortality (Table S5 – Supplementary material). Populations relying on limited transportation networks and under-resourced healthcare systems face compounded barriers to timely treatment, while delayed care places additional strain on local health services through prolonged hospitalizations and intensive care needs (Oliveira RAD et al., 2022; Martinez et al., 2024).
From a public health perspective, an important consideration concerns how these modeled spatial mismatches relate to empirically recorded snakebite incidence. Although assessing spatial mismatches between ATCN and areas with higher snakebite envenoming incidences is highly relevant, such an analysis is constrained by data availability in Brazil. The most spatially refined officially available snakebite records report cases only at the municipality level, without providing geographic coordinates of actual accident locations. This limitation prevents fine-scale spatial analyses linking snakebite incidence directly to ATCN accessibility, particularly at sub-municipal resolutions required for travel-time–based assessments. Nevertheless, we compiled municipality-level snakebite incidence data and explored them as a complementary analysis, presented in the Supplementary Material (Figures S11-S14), allowing for broader spatial comparisons between snakebite occurrence patterns and ATCN distribution while acknowledging the inherent spatial uncertainty of these epidemiological records.
Together, these findings highlight the urgent need to reassess antivenom distribution strategies in Brazil, prioritizing regions of heightened ecological risk and socioeconomic vulnerability. Incorporating climate change projections into healthcare planning is, therefore, essential to prevent the widening of service gaps and to ensure that antivenom availability remains aligned with shifting snake distribution patterns and snakebite risks.
In this regard, the inadequate allocation of ATCN carrying antivenom in areas where medically important snakes do not occur constitutes not only a waste of public resources but also compromises antivenom availability where it is genuinely needed. According to the WHO, between 20% and 40% of health resources are wasted due to systemic inefficiencies that could otherwise be redirected to address critical healthcare delivery gaps (Mazon et al., 2015).
Given this scenario, it is paramount that antivenom allocation be guided by robust epidemiological evidence, incorporating updated information on the geographic distribution of venomous snake species, snakebite incidence patterns and healthcare facility response capacity. The implementation of data-driven public policies, combined with effective governance practices emphasizing transparency and accountability is of the utmost important in optimizing the use of public resources and enhancing the effectiveness of healthcare actions aimed at addressing snakebite envenomation incidents in Brazil (Santos and Rover, 2019).
Beyond immediate public health implications, the projected decline in snake populations driven by climate change and habitat loss may also disrupt ecological processes. Snakes play a key role in regulating small vertebrate populations, including rodents that act as disease vectors, and their decline may contribute to trophic imbalances and increased zoonotic risk (Martinez et al., 2024). Additionally, reductions in snake populations have implications for the pharmaceutical and biotechnology sectors, where venom is a fundamental resource for antivenom production and biomedical research (Oliveira AL et al., 2022; Martinez et al., 2024).
Importantly, snakebite risks do not map perfectly onto natural snake distributions. Underreporting in sparsely populated regions, as well as incidents occurring in artificial contexts, such as legal or illegal captive facilities, can further complicate the relationship between ecological suitability and observed incidence. These dynamics reinforce the need for integrative approaches that consider ecological, social, and infrastructural dimensions simultaneously.
An additional source of uncertainty relates to the temporal structure of the occurrence dataset used for species distribution modeling. As documented in the Supplementary Material (Table S1, Fig. S1), a substantial amount of occurrence records lacks precise collection dates. Excluding undated records would have drastically reduced sample sizes and prevented the modeling of most species, undermining the central objective of this study. Importantly, demonstrably historical records (≤1949) represent only a very small fraction of the dataset, whereas the majority of records either correspond to recent decades or lack temporal metadata altogether. In this context, and following a conservative SDM philosophy, we prioritized spatial representativeness over strict temporal filtering, while explicitly excluding areas with strong biogeographic evidence of species absence, as exemplified by the treatment of Lachesis occurrences. This approach balances data limitations with ecological realism and allows robust inference of broad-scale suitability patterns relevant to public health planning, while acknowledging that fine-scale local extirpations or recent range contractions may not be fully captured.
Some limitations of this study should be acknowledged. Not all Brazilian Bothrops and Micrurus species could be included due to insufficient occurrence records, and extreme climatic events were not explicitly modeled. These factors may further amplify habitat losses and risk patterns under future scenarios. Nevertheless, the broad spatial trends identified here are robust and reveal structural vulnerabilities in antivenom accessibility.
Geographic accessibility in a country as large and heterogeneous as Brazil poses challenges that extend well beyond the physical presence of ATCNs. In particular, Northern Brazil, although the least populated region in the country (Oliveira RAD et al., 2022), reports highest number of snakebite incidents (Table S5 – Supplementary material) and the greatest mobility challenges (Figure S15 – Supplementary material). These challenges are amplified by the particularities of each locality, such as the spatial distribution of the human population, the most common local snake bites, the local transport network and means of transport, as well as ATCN locations and travel times and costs. The prevailing cultural and socioeconomic conditions in each area should also be considered when planning effective and efficient ATCN access logistics (Oliveira RAD et al., 2022).
Ultimately, snakebite envenoming is a dynamic, multifactorial problem that cannot be addressed through static healthcare planning (Nori et al., 2014). Continuous monitoring, adaptive public policies, and cross-sector collaboration are essential. By integrating ecological niche modeling with accessibility metrics and public health infrastructure, this study provides a transferable framework that can be applied to other regions, taxa, and climate-sensitive health risks. Addressing snakebite envenoming through a One Health lens, recognizing the interconnectedness of human, animal, and environmental health, offers a pathway toward more equitable, resilient, and forward-looking healthcare systems under a climate change scenario (OHHLEP, 2022).
Author statementsGFCG, RT and RM conceived the ideas and designed the project; GFCG gathered the data; GFCG and RT performed the analyses; GFCG produced maps and figures; and all authors led the manuscript preparation, contributed critically to the drafts, and gave final approval for submission.
Compliance with ethical standardsThe authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
Declaration of Generative AI and AI-assisted technologies in the writing processDuring the preparation of this work the authors used Chat GPT in order to improve readability and language. After using this tool, the authors reviewed and edited the content as needed and takes full responsibility for the content of the publication.
GFCG would like to thank the postgraduate Program in Biodiversity and Health of the Oswaldo Cruz Institute (Fiocruz) and CAPES for the granted PhD fellowship (88887.614469/2021-00). RT would like to thank FAPERJ for the granted fellowship (E-26/200.238/2023). RM would like to thank CNPq (315184/2021–3) and FAPERJ (E-26/203.274/2017, E-26/210.254/2018, E-202.487/2018, E-26/200.967/2021, E-26/210.071/2023) for financial support. We also thank Ricardo Oliveira (FIOCRUZ) for providing the georeferenced list of treatment centers, and the curators of all herpetological collections who kindly shared specimen data. We are also grateful to Marcelo Weber (Santa Maria University), Rachel Davis (FIOCRUZ), Daniel Fernandes (UFRJ) and other researchers who directly or indirectly contributed with reviews and suggestions throughout the development of this work.
The following are Supplementary data to this article:
