Climate-driven changes in zoonotic risk of arenaviral hemorrhagic fevers in South America

Long-term changes in climate and extreme weather events are known to influence the spread of rodent-borne diseases^1,2,3. Attributed to direct and indirect effects of climate change, over 3000 mammalian species, including rodents, are projected to shift their distribution by 2070⁴, potentially amplifying the transmission risk of rodent-borne diseases. Shifts in the distribution of rodent reservoirs of these zoonotic diseases are driven by changes in food availability, migration of rodent populations, climate suitability, and human population dynamics. Prior research shows that bioclimatic variables such as temperature, precipitation, and changes in land use have been instrumental in modeling the changing risk of high-consequence rodent-borne diseases such as Lassa fever, Argentine Hemorrhagic fever^5,6,7,8 and Hantavirus⁹. Specifically, arenaviral spillover is strongly influenced by human-rodent contact and the distribution of rodent reservoirs^10,11. Sarute and Ross estimated an increase in human–rodent interactions influenced by anthropogenic drivers and climate change¹². It is estimated that the distribution of Calomys musculinus, the rodent reservoir of Junin virus (JUNV), which causes Argentine Hemorrhagic Fever (AHF), will undergo substantial changes in the future as a response to climate change^1,6.

South American New World Arenaviruses (NWAs) include Guanarito virus (GTOV) in Venezuela and Colombia, Machupo virus (MACV) in Bolivia and Paraguay, and JUNV in Argentina. These viruses have caused multiple human outbreaks, with case fatality ranging from 5 to 30%¹³. Despite their inclusion on priority pathogen lists by global health agencies^10,14,15, aspects such as disease dynamics, changes in reservoir distributions and factors associated with distribution changes have not been studied extensively¹⁴. In contrast, some of the Old World Arenaviruses (OWA), such as Lassa fever virus in Africa have been extensively studied in terms of disease dynamics and spillover risk^8,16,17. The distribution of Mastomys natalensis, the known reservoir of Lassa fever, has been shown to be sensitive to factors such as temperature range, precipitation seasonality, co-presence of humans and rodent populations, and land-use patterns^8,17. An extensive meta-analysis of literature on Lassa fever showed increased transmission risk in peri-urban settings that was associated with movement of rodents, humans, and climate change¹⁸. On the other hand, the impact of climate change on NWA spillover from rodents to humans remains poorly understood.

These climatic shifts, combined with anthropogenic land-use changes, could profoundly affect rodent populations and elevate the risk of arenaviral spillover¹⁹. In particular, increased wildfire frequency and intensity (driven by hot-dry conditions and land-use disruptions) create a post-fire environment that favors rapid repopulation by known rodent reservoirs of NWAs, a pattern linked to increased risk of human cases. Our previous study on risk of JUNV indicated that factors associated with changes in temperature, rainfall, and land-use could alter the habitat suitability of the rodent reservoir⁶. Rodent habitat suitability combined with human population density can shift the spatial and temporal patterns of human risk from NWAs²⁰. Given the climate change paradigm, changing patterns of land-use and other anthropogenic activities, the spillover risk of NWAs needs to be quantified based on changes in rodent-human interactions and force-of-infection.

We hypothesized that the spillover risk of rodent-borne NWAs is closely linked to the species distribution of their rodent reservoirs, and that projected climatic changes in South America may facilitate their spread to previously non-endemic regions. We further expect that climate-driven changes in the distribution of the rodent reservoirs will be predictive of changes in NWA spillover risk to humans. In this study, we estimated a proxy for the spillover risk by quantifying the FOI, defined as the probability of successful transmission to humans, based on human-rodent interactions under different climate change scenarios. FOI was mechanistically modeled using species distribution patterns and estimated population densities of six known rodent reservoir: Zygodontomys brevicauda (GTOV), Sigmodon alstoni (GTOV), Calomys callosus (MACV), Calomys musculinus (JUNV), Calomys laucha (JUNV) and Oligoryzomys flavescens (JUNV). This integrated approach combines species distribution models (SDMs) with mechanistic FOI modeling under three climate scenarios: (i) the current climate, and two future climate scenarios based on SSPs²¹ (ii) SSP2-4.5 (moderate climate change scenario) and (iii) SSP5-8.5 (extreme climate change scenario) for the period 2041–2060.

Hotspots for spillover risk of New World Arenaviruses will increase in the future

Projected spatial risk profiles under both future climate change scenarios (SSP2-4.5 and SSP5-8.5) indicated more widespread and elevated spillover risk, represented by force-of-infection (FOI) as a proxy metric, henceforth) for all three NWAs compared to projections for current conditions (Fig. 1). Elevated FOI was also predicted in some regions where these NWAs are already endemic (Fig. 1). Overall, the two climate scenarios produced similar FOI patterns, with SSP5‑8.5 generally yielding smaller deviations from the current climate than SSP2-4.5 (Figs. 1, 2).

In both climate scenarios, GTOV elevated FOI areas expanded beyond current endemic zones in Venezuela, with higher projected FOI in non‑endemic areas of Guyana, Suriname and Brazil (Fig. 1A). In contrast, endemic regions along Venezuela’s northern coast showed reduced FOI, while the interior habitats of GTOV reservoirs exhibited little change under either scenario (Fig. 1A). Some coastal endemic areas showed unchanged GTOV FOI under SSP2‑4.5 and slightly lower FOI under SSP5‑8.5 (Supplementary Fig. 1.4). For MACV, FOI was higher overall under both scenarios, with declines in endemic zones along the Bolivian Andes foothills and increased FOI in non‑endemic grasslands of Bolivia and Paraguay (Fig. 1B). JUNV showed reduced FOI in Argentina’s Pampas grasslands and higher FOI in the surrounding non‑endemic zones under both scenarios, with the magnitude of JUNV FOI change being greater under SSP5‑8.5 than under SSP2‑4.5 (Fig. 1C). Some endemic areas in Buenos Aires province showed little change, whereas FOI increased around the capital region, and non‑endemic regions along the borders of Argentina, Paraguay and Bolivia showed higher FOI in both scenarios, overlapping spatially with MACV high‑FOI zones (Fig. 1B, C).

For all three NWAs, average FOI increased under both SSP2‑4.5 and SSP5‑8.5 compared with the current climate (Fig. 2A, B). For GTOV and MACV, the increase in FOI was higher in eastern tropical portions of reservoir ranges than in the drier western regions separated by the Andes mountain range and other hilly areas (Fig. 2A). For MACV, this reflected a shift in FOI from the Andean foothills toward interior grasslands, yielding larger average increases in eastern Bolivia than in the higher‑altitude western regions (Fig. 2A). GTOV FOI increased by 5% and 8% (median values) under SSP2‑4.5 and SSP5‑8.5, respectively. For MACV, median FOI increased by 19% and 18% under SSP2‑4.5 and SSP5‑8.5, respectively, whereas JUNV FOI increased by 3% under both scenarios (Fig. 2B). Potential FOI hotspots (areas within 90th percentile of estimated FOI) for NWA spillover identified under current conditions persisted under both SSP scenarios (Supplementary Fig. 1.3).

Changes in the NWA force-of-infection (FOI) were found to be associated with climate and land-use changes in the reservoir habitats

For all three NWAs, the predicted increase in FOI under both future climate change scenarios corresponded to higher temperature-related features (annual range, seasonality), reduction in precipitation-related features (precipitation in warmest quarter or in wettest quarter), expansion of urban and crop land, and contraction in forest land (Fig. 3A; Supplementary Table 3.1).

For GTOV, changes in FOI from the current to the future scenarios were found to be associated with features related to temperature, such as annual range and seasonality (Fig. 3A). For the top three features in the SSP 2-4.5 scenario model, the estimated increase in FOI corresponded with an increase in temperature seasonality and the presence of crop land. In contrast, the estimated decrease in FOI corresponded with an increase in precipitation of the wettest quarter of the year (Fig. 3B). For the SSP 5-8.5 scenario, the predicted change in FOI showed a positive association with changes in annual range and seasonal patterns of temperature. In contrast, an increase in the presence of forest land corresponded with a decrease in FOI (Fig. 3B).

For MACV, the changes in FOI were associated with changes in precipitation and temperature features. Considering the top three features for the SSP 2-4.5 scenario, the increase in FOI corresponded to a decrease in isothermality and annual mean temperature, and a slight increase in the minimum temperature of the coldest month (Fig. 3A). Under the SSP5-8.5 scenario, changes in FOI were positively associated with precipitation of the wettest quarter. A decrease in FOI corresponded with a decrease in the mean temperature of the wettest quarter.

For JUNV, the change in FOI corresponded to changes in temperature, precipitation, and presence of urban or crop lands for both scenarios (Fig. 3A). Under the SSP2-4.5 scenario, FOI was negatively associated with temperature seasonality and precipitation in the warmest quarter and positively associated with urban land cover (Fig. 3B). Under the SSP5-8.5 scenario, FOI was slightly negatively associated with precipitation of the wettest quarter and negatively associated with mean temperature of the coldest quarter, while urban land showed a positive association (Fig. 3B).

The models for investigating the changes in FOI versus features performed moderately well. The models under SSP 2-4.5 scenario for GTOV, JUNV and MACV explained 77.4%, 77.8%, 90.2% of variance, respectively (Pseudo R-squared estimate). Similarly, models under SSP 5-8.5 explained 75.2%, 72.4% and 87.5% of the variance in the changes in FOI, respectively.

Predicted shifts in species distribution patterns of NWA rodent reservoirs

We predicted the species distribution patterns for the six rodent reservoirs of the three NWAs for (i) the current, and the future in years 2041–2060 represented by (ii) SSP 2-4.5 scenario and (iii) SSP 5-8.5 scenario, using a Species Distribution Modeling (SDM) framework. In general, the projected future species distributions showed different spatial patterns and magnitudes for presence probabilities for all rodent reservoir species compared to the current distributions, predicting a radical change in the habitats of these rodent reservoirs in the future. However, the differences in the probabilities of presence between the SSP 2–4.5 and SSP 5-8.5 scenarios, were more subtle (Fig. 4). We also predicted more widespread range of distribution for the NWA rodent reservoirs compared to their historically reported occurrences.

We projected the distribution of Z. brevicauda away from the Caribbean coast and towards forest and rural areas (Fig. 4A), with lower probabilities in existing metropolitan and highly populated areas. We also projected a westward shift of S. alstoni distribution for both SSP 2-4.5 and 5-8.5 scenarios. Our models did not predict similar shift patterns for the species distribution of C. callosus. We projected minor changes for C. callosus, with increased probability in northern parts of Paraguay and western Brazil, and lower probabilities in central Bolivia and northern Argentina (Fig. 4B). We predicted a decrease in distribution probabilities in the central region of Argentina for C. musculinus (Fig. 4C). Conversely, we predicted increased probabilities for C. laucha and O. flavescens in the same region, thereby showing changes in the distribution range of rodent species in this region (Fig. 4C). We projected the distribution of O. flavescens to shift northward, whereas that of C. musculinus to shift southward. The presence of C. laucha was projected to shift inland into the central region of Argentina. Detailed visual observations on maps of the geographical changes for each species distribution are shown in Supplementary Table 2.3.

The distribution of Z. brevicauda was mostly influenced by crop land and temperature seasonality. For S. alstoni the most important features (prediction drivers) were annual precipitation and temperature seasonality. For C. callosus, annual and diurnal ranges of temperature were the most important features. For C. musculinus, the most important features were presence of urban land, and annual precipitation. For C. laucha, crop land, maximum temperature in the warmest month and warmest quarter were the most important features. For O. flavescens, crop and urban land were equally important features followed by annual precipitation.

The fitted SDMs performed moderately well with a cross-validation (CV) accuracy between 77% and 87% and a test accuracy of 75% to 98% depending on the reservoir species (for details on precision, recall, F1 scores and confusion matrices, see Supplementary Information Section B 2.3).

Our models predicted that NWAs could expand into currently non-endemic areas as reservoir distributions shift under climate change, potentially increasing the risk of human spillover. Using an epidemiological perspective, we linked NWA rodent reservoir habitat patterns to potential human spillover risk. In our study, a higher FOI (proxy metric for the spillover risk) indicated a greater likelihood of reservoir presence and a greater likelihood of epidemiological contact between infectious reservoirs and humans. This approach differed from prior studies, which have used species distribution models for determining habitat suitability of disease vectors and reservoirs²².

In this study, the FOI was calculated under a density‑dependent framework for rodent–human transmission. We used a transmission rate parameter informed by Lassa fever (OWA) data, as a pragmatic proxy in the absence of NWA‑specific estimates. Although we lacked direct data on reservoir population densities, approximating rodent abundance from SDM‑derived occurrence probabilities offers a consistent, model‑based way to scale FOI across species and regions, while identifying population density as a key target for future field validation²³. Thus, disease risk might be overestimated. Salazar et al. found that the geographic extent of rodent reservoir habitat was substantially larger than the areas where NWA outbreaks had been recorded²⁴. However, it is highly probable that NWAs could spill over into non-endemic regions due to changes in human movement patterns, as well as changes and fragmentation of rodent habitats that may increase proximity to human settlements. This is in line with IPCC AR6 findings of increasing human-animal contact leading to the emergence of zoonotic spillovers, such as those of NWAs²⁵.

Our model framework aimed to reduce cumulative uncertainty in future projections that systematically occurred due to modeling choices and assumptions. By fitting ensembles of SDMs and FOI models across iterative data resampling and parameterization, we summarized the spread of resulting risk instead of relying on a singular model outcome of changing spillover risks. This approach is consistent with recommendations for spatial uncertainty propagation and ensemble SDM modeling²⁶. The uncertainty in FOI parameters and model structure was further reduced by the use of ensemble tree-based models with repeated cross-validation and recursive feature elimination, so that variation in fitted relationships between changes in the predictors and that of FOI was reflected in the distribution of out-of-sample predictions rather than a single deterministic surface (for details see “Methodological Discussion”).

Our outcomes conform with the findings of previous studies, which have shown that ecological traits of rodent reservoir species in the New World are sensitive to changes in temperature range and/or seasonality, as well as precipitation patterns and rainy season^27,28,29. Secondary effects of climate change, such as alterations in anthropogenic land use and subsequent effects on disease dynamics have also been reported previously³⁰. Considering that all six of the NWA reservoirs are highly adaptable species, habitat changes can be expected without significant loss in population size³¹. GTOV reservoirs (Z. brevicauda, S. alstoni) can be found in shrub and grassland areas with wet conditions. Hence, their sensitivity to seasonality in temperature and precipitation was expected³¹. In the case of MACV, we theorize that the apparent sensitivity to changes in temperature and precipitation might be ascribed to the habitat preferences of C. callosus, along with the expected increase in human migration towards non-urban areas, since C. callosus is known to inhabit areas with dry, semi-arid climate conditions with open vegetation^32,33. Moreover, we predicted overlapping areas with high FOI of BHF and AHF, caused by MACV and JUNV, respectively. The persistence of potential hotspots in endemic regions is consistent with findings from Portuguesa State, Venezuela, where Silva-Ramos et al. documented increased and persistent hotspot zones for VHF outbreaks³⁴.

Our findings are in line with prior studies showing that shifts in the habitat of rodent reservoirs of JUNV intersect with areas of human activity^6,35. Crucially, we predicted widespread habitat distribution across all six rodent species that may signal a possible fragmentation of rodent habitats in the future. Particularly for MACV and JUNV, habitat fragmentation might be expected in the future in all climate change scenarios for Andean foothill biomes as well as grasslands in northern Argentina, Bolivia, and regions in neighboring Paraguay and Chile³⁶. By predicting these fragmented habitats accurately, high-risk “bridge” areas between current endemic habitats could be identified to integrate into additional surveillance efforts.

Given that the current global climate trajectory aligns with the moderate climate change scenario of SSP 2-4.5^36,37, our predictions may be particularly salient for monitoring potential outbreak hotspots and optimizing the use of scarce public health resources. Furthermore, our methodological framework for these diseases can be adapted to other neglected diseases, as well as other climate scenarios. The outcomes can serve as a template for global health agencies to prioritize at-risk communities and countries facing disproportionate health burdens due to outbreak risk, enabling proactive rather than reactive responses.

Rodent reservoirs for the New World Arenaviruses (NWA)

We selected six rodent species that are known to act as reservoirs for NWA and contribute to zoonotic spillover events (Table 1). The criterion for selection was that the reservoir species had to be seropositive for the virus during the time of the outbreak, with at least two studies confirming their reservoir status during the time of the outbreak^{24,28,32,37,38,39,40,41}. The six rodent species, namely, Zygodontomys brevicauda and Sigmodon alstoni for GTOV, Calomys callosus for MACV, Calomys musculinus, Calomys laucha, and Oligoryzomys flavescens for JUNV, have been commonly associated with outbreaks of their respective NWA in human populations (Table 1)^{42,43,44,45,46}.

Species Distribution Model (SDM)

Modeling frameworks such as species distribution models (SDMs), which can incorporate bioclimatic and environmental factors to estimate the distribution of rodent reservoirs, can be effectively employed in predicting transmission risk for associated diseases⁴⁷. Since the modeled association between climate change and species distribution was not expected to be monotonic⁴⁸ but rather complex and subtle, modifying the SDM framework to employ ensemble tree-based techniques such as random forest (rf), extra trees (et), extreme gradient boost (xgb), and light gradient boosting machine (lgbm) is warranted⁴⁹. We followed a similar protocol to the one we previously used for modeling the species distribution of Calomys musculinus in Argentina⁶. The protocol was adapted to each of the rodent species selected for this study (Table 1). Here is a summary of the SDM protocol (Fig. 5) used to model the rodent reservoirs of NWA following the guidelines established for reporting Ecological Niche Models or SDMs to maximize reproducibility as closely as possible pertaining to our study⁵⁰.

Occurrence data on all six rodent species were sourced from Global Biodiversity Information Facility (GBIF) database (presence-only)⁵¹ using R-package dismo⁵² (Derived datasets⁵³). These data were cleaned to remove duplicates and restricted to years 1990 onward to match the period of climate/land use features. All occurrences in the South American continent were retained regardless of the country. Presence-only data were preprocessed to be converted to 1:1 presence-absence data using habitat suitability analysis with R-package USE::paSampling() utilizing a uniform approach⁵⁴. This was done to avoid sample location bias that commonly occurs with random sampling of a geographical extent without considering the environmental conditions, as discussed by Da Re et al. (see Supplementary Section C3.4 for details)⁵⁵. Environmental and ecological data were sourced from open-source datasets, such as 19 Bioclimatic variables from World Clim (Bioclim), Normalized Differential Vegetation Index (NDVI) and Digital Elevation Model (DEM) from Moderate-Resolution Imaging Spectroradiometer (MODIS) satellite database and five Land Use datasets from Land-Cover and Land-Use Change (LCLUC) program maintained jointly by NASA and University of Maryland^{56,57,58,59,60}. These datasets were resampled to the highest common resolution and clipped to the extent of species occurrence data mentioned above for each individual rodent reservoir. In total, three sets of 24–26 rasters and six occurrence datasets were used for the SDMs (Supplementary Table 2.1).

Similarly, for predicting the species distribution in the future, analogous raster data for Bioclimactic variables⁶¹, NDVI (data from 2020 for NDVI⁵⁷), DEM⁵⁸, and land use⁶² were downloaded for Shared Socioeconomic Pathways 6 (SSP) scenarios of Moderate (SSP 2-4.5) and Extreme (SSP 5–8.5) predictions for climate change in the future as established by Coupled Model Intercomparison Project v6 (CMIP6) for years 2041 to 2060²¹. These datasets were also preprocessed and resampled to the highest common resolution and averaged over the 20-year time period.

SDMs were developed by fitting four ensemble tree-based classifier algorithms, namely, random forest (rf), extra trees (et), extreme gradient boost (xgb) and light gradient boosting machine (lgbm). The occurrence datasets for the six rodent species were considerably smaller compared to the standard size of presence-absence datasets in similar studies, limiting our capacity to perform spatial block cross-validation which is a standard practice for ecological models⁶³. Instead, data was split randomly into train-test ratio of 4:1 (80% train; 20% test). The missing values of each predictor (raster features) from both the training and test sets were imputed using SimpleImputer() from scikit-learn python library⁶⁴ with the mean of the existing values of each predictor. This was necessary for fitting “rf” and “et” algorithms to the training data. Each algorithm was trained with fivefold cross-validation and the predictions for training and the test set were generated to extract cross-validation accuracy, precision, recall and the F1 score based on (i) fivefold cross-validation and (ii) test set confusion matrices. This process was iterated 100 times to improve the robustness of the models using 100 different presence-absence resamples from the occurrence data. To retain the appropriate number of features necessary to predict the distribution probabilities, we used a recursive feature elimination (RFE) algorithm with stratified cross-validation using the RFECV algorithm⁶⁴. Final model predictions were derived by averaging results across all iterations and the four classifier algorithms.

Since hyper-tuning did not yield significant improvements to any of the four classifier algorithms (details in Supplementary Section B: Hypertuning Results and Parameter Selection), the default hyperparameters were used. We expected a high level of collinearity between the rasters based on spatial correlation analysis (results not shown), given the nature of geospatial data. To counter this, RFE was performed at each iteration to improve the accuracy of the model without having potential issues with multicollinearity. In each iteration, 10 features were selected per classifier algorithm based on stratified cross validation of RFE performed before the iterative model training using RFECV() from scikit-learn⁶⁴. Final fitted models were taken as the average of all the iterations and of the four classifier algorithms.

Based on the fitted models, the current distribution probabilities of each rodent species were imputed to the geographical extent of their occurrence data using pyimpute library in python⁶⁵. Similarly, the future distributions of the rodents in response to SSP 2-4.5 and SSP 5-8.5 scenarios of CMIP6 climate change were imputed based on the fitted models. Changes in the distribution probabilities were mapped onto the raster for each of the scenarios by subtracting current probabilities from the future probabilities of species distribution. The process is illustrated in Fig. 5.

Zoonotic risk for spillover and force of infection

Based on the distribution probabilities of the SDM for each rodent reservoir, a risk profile for human outbreak, i.e., zoonotic spillover of NWA, was modeled using infection dynamics simulation. The current and CMIP6 scenario-based projections of human population density were overlaid on the SDM probabilities to generate a risk profile based on the force of infection for NWA spillover. Force-of-infection (FOI) was calculated based on the contact rate between humans and the rodent reservoirs and the possibility of the rodent testing positive as reservoir for the NWA.

The mechanistic model used for estimating FOI was a density dependent contact rate model with binomial sampling as detailed in Eq. 1.

Where, ({FOI}) is the force of infection defined as the contact rate between susceptible humans and infectious rodents resulting in the successful transmission of infection, (beta) is the transmission rate parameter derived from review of analogous viral transmission dynamics studies (see Appendix: Force of Infection), (H) is the population of susceptible humans set at 0.95 times that of total human population in the same geospatial coordinates of SDMs (author’s expertise and from study performed on Lassa Fever in Nigeria⁶⁶) and (R) is the infectious proportion of the rodent population based on the binomial sampling between 1 and 15 rodents per grid cell, adjusted with the probability of presence of rodents in the given geospatial coordinates based of the SDMs. Given the lack of data on rodent populations, we estimated the number of rodents present in each grid cell based on the estimation of rodent family size found in one location by Han et al.⁶⁷ and binomial sampling. The denominator (left(H+Rright)) represents the total density of the interacting populations of human hosts and rodent reservoirs. Equation (1) was adapted from a generalized formula used in similar transmission studies of vector and rodent-borne infections with and without the effect of climate change^68,69,70.

A similar methodology was followed using the projected SDM probabilities and the projected human population under the SSP 2-4.5 and SSP 5-8.5 scenarios of CMIP6⁷¹. For generating these maps, the FOI of each reservoir species of the same virus were combined ({FO}{I}_{{virus}}^{t}=mathop{sum }nolimits_{r=1}^{n}{FO}{I}_{r}^{t}) where r was the rodent reservoir and (tepsilon {{current},{SSP}2-4.5,{SSP}5-8.5}).

A similar process was done for ({FO}{I}^{{SSP}2-4.5}{andFO}{I}^{{SSP}5-8.5}).

The differences ((delta)) between the future climate change scenarios (SSP 2-4.5 and SSP 5-8.5) and the current scenario of FOI (({delta }_{{foi}}^{{ssp}2-4.5}{and}{delta }_{{foi}}^{{ssp}5.85})) were calculated and scaled to be from −100% to + 100%.

The map depicting FOI for each of the three viruses under study was further reclassified in 10% quantiles and the geospatial zones that had FOI values in the top 10th percentile (≥90% of the FOI estimated range across the geographical bounds of the raster) were depicted as potential hotspot zones for outbreak of zoonotic arenaviruses.

Association of changing risk with climate change

We modeled the associations between the changes in FOI and the changes in 24 bioclimatic and land-use features that were used to predict the species distribution patterns of the reservoir species. NDVI and DEM features were dropped from these analyses due to a high number of missing values. The difference (({boldsymbol{Delta }})) between future climate change scenarios (SSP 2-4.5 and SSP 5-8.5) and the current scenario of FOI estimates (Eq. 2a), along with the corresponding features (Eq. 2b), were used as inputs in a random forest regressor model for each of the three viruses. The ({{boldsymbol{Delta }}}_{{boldsymbol{f}}}^{{boldsymbol{ssp}}}) and ({{boldsymbol{Delta }}}_{{FOI}}^{{boldsymbol{ssp}}}) values were scaled and a combined tabular dataset was generated. A random forest regressor model was fitted (RandomForestRegressor(), scikit-learn⁶⁴). For each random forest model, feature importance based on mean impurity reduction was extracted. The partial dependence plots of the top three most important features were developed from average grid values and predicted outcomes using functions from scikit-learn, matplotlib and seaborn packages in python^64,72,73.

Where ({boldsymbol{f}}) are individual features in the dataset (for example, temperature range, annual precipitation, etc.)

Methodology discussion

In our study, we restricted the focus to three NWAs—GTOV, MACV and JUNV—which have been associated with outbreaks of substantial magnitude⁴³. This restriction excludes nine NWAs out of 19 identified mammarenaviruses³⁰. Several of these NWAs have been shown to be emergent, such as Sabia virus in Brazil, Chapare virus in Bolivia, Pirital virus and Tacaribe virus in Venezuela and Guyana. However, only a few cases of these lesser-known NWAs have been reported with confirmed diagnoses or acute deaths investigated⁴³. Of note, rodent species that are reservoirs for one virus may also be reservoirs for other similarly transmitted viruses and are thus non-specific reservoirs. For example, reservoirs of VHF caused by Guanarito virus have been shown to carry viruses from the Bunyaviridae family, as well as other hantaviruses⁷⁴. The issue of putative rodent reservoirs is further complicated by the contentious classification of rodent species in the New World. The continuously expanding species diversity and evolving classification of these rodents (especially around the Andes foothills⁷⁵), combined with different nomenclature and local names, complicates the study and targeting of NWAs reservoir species. We thus applied a species distribution model for rodent reservoirs that was not limited to the extent of historically endemic zones of NWAs, but rather by the geographical extent of reported occurrences of the rodents. A macro-ecological outlook that is reservoir-centric rather than disease outbreak-centric could be more helpful, given the lack of clarity on reservoir species diversity and their classification⁶⁷. Whether multiple reservoirs exist for the same virus, a single reservoir carries multiple viruses, or multiple reservoirs carry multiple viruses, the risk NWA zoonotic spillover heavily depends on the presence, density and habitat patterns of reservoir hosts^67,76.

To investigate the association between changes in climate and changes in the spillover risk represented by FOI, we fitted a random forest (RF) model. The changes in bioclimatic and land use variables (({{boldsymbol{Delta }}}_{{boldsymbol{f}}}^{{boldsymbol{ssp}}})) were multi-collinear given the nature of the SSP scenario estimation of bioclimatic (see ref. ⁶¹ for details) and land use variables. Tree-based ensemble models such as RF are more robust to multicollinearity than classical statistical models, which are normally applied in association studies when iterated multiple times and combined to give a single estimate of permutation feature importance⁷⁷. These ensembles naturally capture non-monotonic effects and higher‑order interactions among climate variables, which are expected in ecological and epidemiological responses to climate change⁷⁸. While our outcomes cannot demonstrate statistically significant associations between changes in features and changes in FOI, the RF model identified feature importances based on permutation importances over 100 model fits each (3 viruses × 2 SSP scenarios × 100 iterations of models = 600 model fits). These feature importance results indicated associated trends in change between climate and FOI.

Similar to the random forest models used to investigate associations between FOI and climate change, we developed and implemented SDMs using tree‑based ensemble algorithms. This method is more robust to multicollinearity among predictors and can accommodate nonlinear responses and interactions without relying on unstable coefficient estimates, compared to classical statistical methods⁷⁹. To address the multicollinearity of bioclimatic and land‑use covariates (details in Supplementary Section C), we combined repeated model refitting on resampled presence–absence data with recursive feature elimination embedded in cross‑validation, thereby selecting a subset of predictors based on out‑of‑sample predictive performance rather than manually defined arbitrary correlation thresholds^80,81. In line with the concerns raised by Dormann et al. that collinearity is pervasive in ecological data and that no single technique fully resolves its effects⁸². Permutation importance per feature was then calculated repeatedly on the RFE‑selected ensembles to quantify each predictor’s contribution under the observed collinearity structure (including standard errors), providing an interpretable ranking of climate and land‑use drivers of projected changes in disease risk (Supplementary Fig. 2.3). By combining these machine-driven techniques, we let out‑of‑sample performance determine which correlated climate predictors were retained, rather than relying on fixed correlation coefficient cut‑offs.