DOI: 10.5281/zenodo.15645100
title: Unveiling the Ecological Drivers of Bat-Borne Viral Epidemics: A Climate-Driven Risk Assessment

description: Bats serve as natural reservoirs for numerous lethal viruses, significantly threatening global public health. Despite the substantial rise in identified bat-hosted viruses, our understanding of their ecological dynamics remains limited, hindering the development of effective preventive strategies for bat-borne viral diseases. By examining 19 underlying variables, this study investigates the complex relationships between 1,350 bat species and 68 viral genera, focusing on potentially infectious (PIVs) and epidemic viruses (PEVs). Phylogenetic and environmental analyses reveal uneven viral diversity across bat families, driven by temperature seasonality, host life-history traits (e.g., longevity, litter size), and migratory behavior. Additionally, socioeconomic factors, including human population density and GDP per capita, amplify viral detectability and spillover risks. While PIV-carrying bats are distributed across diverse biomes, current hotspots for PEV-carrying bats are concentrated in tropical seasonal forests, temperate regions, and areas with high human activity. Climate change projections for 2050 indicate significant shifts in the distribution of bats harboring PEVs, with viruses like SARS-related coronaviruses and Nipah viruses expanding into colder, drier, and higher-latitude regions. Concurrently, high-risk zones for human-bat interactions are expected to grow due to the combined impacts of climate change and urban population expansion. These findings underscore the critical need to integrate ecological, environmental, and socioeconomic data for assessing zoonotic risks. Proactive surveillance and intervention strategies, particularly in dynamic hotspots experiencing rapid environmental and societal changes, will be essential to mitigating the emergence of bat-borne epidemics in an era of global change.
description: Data Collection

We developed a comprehensive dataset of 1,350 bat species, with 357 species carrying various viruses. Host-virus associations were established using the Database of Bat-Associated Viruses (up to October 2024)56 and 180 published studies on bat virus sampling (Appendix 2).

Due to unresolved viral taxonomy, most data were analyzed at the genus level. To maximize data utility, viral diversity was estimated using genus-level classifications. To account for sampling bias, we applied a nonparametric estimator—the Chao2 index (tabula package)—as it is considered a more precise diversity estimator for host–parasite association databases with sampling variation61. To evaluate the potential threat of bat-borne viruses to humans, we classified bats dichotomously according to WHO guidelines using species-level viral data (Supplementary Table 1): Potentially human-infectious viruses (PIVs): viruses with zoonotic potential and confirmed or predicted capacity for limited human-to-human transmission. (b) Potentially human-epidemic viruses (PEVs): viruses with zoonotic potential for sustained human transmission or outbreak initiation. The classification framework adopts a cautious approach because direct bat-to-human transmission events remain rare for many viruses. For instance, while SARS-CoV-2 is phylogenetically linked to bat coronaviruses, its emergence likely involved intermediate hosts or recombination events40.

To further account for sampling biases, we analyzed bat occurrence data from 1990–2023, totaling 2,433,331 records from the Global Biodiversity Information Facility (GBIF)62 and VertNet63 as bat observation counts, and also viral sampling counts, viral detection counts, and prevalence (viral detection counts/ viral sampling counts)19 for each bat species (Appendix 1).

For evaluating the potential influences of environment and life history factors on the relationship between bats and viruses, we collected (1) geographic and climatic drivers, including altitude, wind speed, and 19 bioclimatic data64; (2) socioeconomic drivers, including human population density and GDP per capita65; (3) life history and morphological drivers including generation length, litter size, diet, habitat breadth, body mass66,67, and flight ability indicator--wing aspect ratio68,69. See details provided in Supplementary Table 6.

Data processing

Geographic, climatic, and socioeconomic traits: 19 bioclimatic layers (WorldClim 2.0) were reduced to 6 low-collinearity predictors (Bio4, 10, 11, 15, 18, 19; VIF < 1.92) using vifstep (usdm package). Species-level geographic, climatic, and socioeconomic niches were defined by averaging values across IUCN range polygons (GBIF/ VertNet records), excluding points beyond IUCN boundaries with a 2° buffer.

Life history and morphological traits: Diet was simplified into three categories, including Herbivore (Fruit, Nect., Seed, Plant), Invertivore (Inv.), and Carnivore (Vend., Vect., Vfish). Wing aspect ratio data from EurobaTrait v1.068 and AfrobaT69  datasets were imputed to genus/ family means for species lacking measurements (337 species with original data).

Migration proxy: many bats undertake seasonal migrations spanning hundreds or thousands of kilometers. During these migrations, bats may be exposed to novel pathogens70 or experience stress-induced immunosuppression71, which can alter the composition of their host-associated viruses. However, despite sparse and sporadic observations72,73, the migratory behaviors at the species level have not been well documented. To address this gap, we propose using indirect migration indicators. Specifically, for each bat species, the seasonal range distribution centroids during the seasonal range (breeding: Jun-Aug; wintering: Dec-Feb) distribution (occurrence data) centroids for each species were calculated using gCentroid (rgeos package) as indirect migration indicators.

All continuous variables were normalized to 0–1 scales for comparable effect sizes.

Environmental, Life-History, and Morphological Predictors Selection

We conducted phylogenetic analyses on 914 bat species (Supplementary Fig. 5) with complete trait data to identify drivers of viral diversity and spillover risk. Virus diversity (accumulated genera per species) was analyzed using phylogenetic linear regression (phylolm), while hosting probability of human-infectious viruses (PIVs) and epidemic-potential viruses (PEVs) were assessed via phylogenetic logistic regression (phyloglm) in the R package phylolm74. Predictors included all variables spanning environmental niches (climate, geographic ranges), socioeconomic factors, life-history, and morphological traits (Supplementary Table 6).

For phylogenetic linear regression (viral diversity): Full model sets were evaluated using automated dredging (dredge) in the R package MuMIn with Pagel’s lambda correlation structure.

For phylogenetic binary logistic regression (PIV/PEV) regressions: (1) AIC ranked single-predictor phylogenetic logistic regressions for each predictor; (2) the lowest-AIC model was iteratively expanded by adding predictors that reduced AIC by ≥2 units and showed significance (p < 0.05). This process continued until no additional models improved upon the base model (AIC < 4 compared to previous or non-significant predictors)75, as shown in process examples in Supplementary Table 7. (3) Research effort covariates (bat observation counts, viral sampling effort, prevalence) were also considered in all models.

Phylogenetic Path Analysis

Using a hypothesis-driven phylogenetic path analysis framework within the R package phylopath76, we assessed the direct and indirect effects of life-history traits and climate on bat-borne virus type and diversity, accounting for phylogenetic influence. The model determination procedure was as follows: (1) a core model structure linking viral diversity, PIV, and PEV probabilities was tested using C-statistic information criterion (CICc, Supplementary Figure. 6); (2) predictors with significant effects in regression analyses were sequentially added to evaluate direct/indirect pathways (Supplementary Table 8); (3) Final models retained only pathways that improved CICc and showed statistical significance (p < 0.05), iterating until no further improvements occurred (Supplementary Tables 8 and 9). The script for these analyses is available in Appendix 3.

Biome Classification of Bat-Virus Associations and Socioeconomic Links

We analyzed the global distribution of 1,307 bat species and their viral status (PIV/PEV carriers vs. non-carriers) using the Whittaker biome classification77. Species were categorized into biomes based on the distribution of centroids.

Socioeconomic associations focused on the overlap of viral host bats and human population based on the risk level of viruses hosted by bats, assessed via bivariate heatmaps using bi_class (bi_class package)78 to compare the distribution of (1) PIV carriers with non-PIV carriers, and (2) PEV carriers with non-PEV-PIV carriers across human population density and GDP per capita. Viral load status was classified into four groups (LL, LH, HL, HH), where "H" (high) and "L" (low) corresponded to values in the top/bottom 33% of pixel-level viral abundance. T-tests were performed to evaluate differences in socioeconomic variables between carrier groups. Details of the analytical script are provided in Appendix 3.

Species Distribution Modeling (SDM) for Risk Projection

(1) Data preparation:

Occurrence records for 42 PEV-hosting bat species (GBIF/VertNet; n >5 per species) were filtered using CoordinateCleaner (CoordinateCleaner package (83))to remove biased locations (e.g., museums). Spatial thinning (0.5° grid resolution) was applied via spThin (spThin package)80 to minimize autocorrelation. Pseudo-absences (1,000 random points/species) and climate predictors (WorldClim 2.0 Bio1–Bio1964; 0.5° resolution) were standardized for model inputs. Our study primarily employs machine learning algorithms, which are significantly less sensitive to collinearity in predictive tasks and can achieve high accuracy without actively removing collinear variables81. In addition, some collinear variables (e.g., temperature and precipitation seasonality) may reflect genuine ecological synergies, and removing such variables could result in the loss of critical environmental drivers (86). So, we chose not to actively remove collinear variables to maximize the retention of data information, which aligns with recent best practices in ecological modeling.

(2) Model implementation:

Multi-algorithm SDMs (SSDM package)83 integrated six methods: GLM, GAM, GBM, CTA, RF, and SVM. Ensemble models were built using 3-fold cross-validation (70% training, 30% testing), retaining only models with AUC > 0.7 and weighting predictions by AUC scores. Future projections (2050) were averaged under SSP2-4.5 and SSP3-7.0 scenarios84. Variable importance and model performance metrics are reported in Supplementary Table 10; Supplementary Fig. 7.

(3) Risk analysis:

Virus-centric range aggregation: Distributions of 11 PEV-hosting bat species were pooled to map potential virus circulation zones (Supplementary Fig. 2). Human exposure hotspots were identified by overlaying PEV richness maps with SSP-based population density projections85, classifying pixels into four risk categories (LL, LH, HL, HH) based on median splits (Supplementary Fig. 4). The category HH (i.e. High PEV richness and High human population density) should be primarily focused due to the heightened public health relevance. Temporal shifts (current vs. 2050) in risk category proportions were quantified to assess emerging public health threats.

All spatial analyses used consistent 0.5° grids; code and parameter details are archived in Appendix 3. All analyses were conducted in R v4.1.2.