DOI: 10.5281/zenodo.17598049
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: Supplementary Methods

for “Unveiling the Ecological and Climatic Drivers of Bat-Borne Viral Epidemics: A Climate-Driven Risk Assessment”

 

Yang Wanga,e#*, Liang Lüa,e#, Ghulam Nabia#, Tianlong Caib, Jiaxin Zhanga,  Xu Liua, Gang Songc, Xinglou Yangd, Limin Wanga, Dongming Lia,e*

 

a Hebei Key Laboratory of Animal Physiology, Biochemistry and Molecular Biology, College of Life Sciences, Hebei Normal University, Shijiazhuang 050024, China.

b School of Life Sciences, Westlake University, Hangzhou 310058, China.

c Institute of Zoology, Chinese Academy of Sciences, Beijing 100101, China

d Kunming Institute of Zoology, Chinese Academy of Sciences, Kunming 650023, China

e Hebei Collaborative Innovation Center for Eco-Environment, Hebei Normal University, Shijiazhuang 050024, China

 

#These authors equally contributed to this study.

*Correspondence:Yang Wang, Email: wangy@hebtu.edu.cn; Dongming Li, Email: lidongming@hebtu.edu.cn

 

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) (Chen et al. 2014) and 2070 records in 180 published studies on bat virus sampling (Appendix 1).

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 variation (Teitelbaum et al. 2020). To evaluate the potential threat of bat-borne viruses to humans, we classified bats dichotomously according to WHO guidelines using species-level viral data: 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 events (Wang et al. 2023).

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.org 2023) and VertNet (Bloom 2016) as bat observation counts, and also viral sampling counts, viral detection counts, and prevalence (viral detection counts/ viral sampling counts) (Ruiz-Aravena et al. 2022) 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, and 19 bioclimatic data (Fick and Hijmans 2017); (2) socioeconomic drivers, including human population density and GDP per capita (Kummu et al.); (3) life history and morphological drivers, including generation length, litter size, diet, habitat breadth, body mass (Myhrvold et al. 2015, Cooke et al. 2019), and flight ability indicator--wing aspect ratio (Cosentino et al. 2023, Froidevaux et al. 2023). 

 

Data processing

Geographic, climatic, and socioeconomic traits: 19 bioclimatic layers (WorldClim 2.0) were reduced to 6 low-collinearity predictors (Bio04, 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.0 (Froidevaux et al. 2023) and AfrobaT (Cosentino et al. 2023) 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 pathogens (Luis et al. 2015) or experience stress-induced immunosuppression (Randolph and Barreiro 2018), which can alter the composition of their host-associated viruses. However, despite sparse and sporadic observations (O’Mara et al. 2019, Bach et al. 2022), 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.

 

Predictors Selection

We conducted phylogenetic analyses on 914 bat species (Appendix 2) 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 phylolm (Ho and Ané 2014). Predictors included all variables spanning environmental niches (climate, geographic ranges), socioeconomic factors, life-history, and morphological traits (Appendix 2). 

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) (Burnham and Anderson 2002), as shown in process examples in Table 1. (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 phylopath (van der Bijl 2018), 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, Fig. 1); (2) predictors with significant effects in regression analyses were sequentially added to evaluate direct/indirect pathways (Table 2); (3) Final models retained only pathways that improved CICc and showed statistical significance (p < 0.05), iterating until no further improvements occurred. The script for these analyses is available in Appendix 2.

 

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 classification (Whittaker 1980). 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 (2022)) 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 2.

 

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 (Zizka et al. 2019) to remove biased locations (e.g., museums). Spatial thinning (0.5°grid resolution) was applied via spThin (spThin package) (Aiello-Lamm ens et al. 2015) to minimize autocorrelation. Pseudo-absences (1,000 random points/species) and climate predictors (WorldClim 2.0 Bio1–Bio19 (Fick and Hijmans 2017); 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 variables (Valavi et al.). 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 (Graham 2003). 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) (Schmitt et al. 2017) 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 scenarios (O’Neill et al. 2016). 

(3) Risk analysis:

Virus-centric range aggregation: Distributions of 11 PEV-hosting bat species were pooled to map potential virus circulation zones (Fig. 2). Human exposure hotspots were identified by overlaying PEV richness maps with SSP-based population density projections (Jones and O’Neill 2020), classifying pixels into four risk categories (LL, LH, HL, HH) based on median splits. 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.

 

References 

Aiello-Lammens ME, Boria RA, Radosavljevic A, Vilela B, Anderson RP. 2015. spThin: an R package for spatial thinning of species occurrence records for use in ecological niche models. Ecography 38: 541–545.

Bach P, Voigt CC, Gttsche M, Bach L, Brust V, Hill R, Hüppop O, Lagerveld S, Schmaljohann H, Seebens‐Hoyer A. 2022. Offshore and coastline migration of radio‐tagged Nathusius’ pipistrelles. Conservation Science & Practice 4.

van der Bijl W. 2018. phylopath: Easy phylogenetic path analysis in R. Peerj 6.

Bloom D. 2016. VertNet_Mammalia_Sept2016.

Burnham KP, Anderson DR. 2002. Model Selection and Multimodel Inference.

Chen L, Liu B, Yang J, Jin Q. 2014. DBatVir: the database of bat-associated viruses. Database: The Journal of Biological Databases and Curation 2014: bau021.

Cooke RSC, Eigenbrod F, Bates AE. 2019. Projected losses of global mammal and bird ecological strategies. Nature Communications.

Cosentino F, Castiello G, Maiorano L. 2023. A dataset on African bats’ functional traits. Scientific Data 10.

Fick SE, Hijmans RJ. 2017. WorldClim 2: new 1-km spatial resolution climate surfaces for global land areas. International Journal of Climatology 37: 4302–4315.

Froidevaux JSP, Toshkova N, Barbaro L, Benítez-López A, Kerbiriou C, Le Viol I, Pacifici M, Santini L, Stawski C, Russo D. 2023. A species-level trait dataset of bats in Europe and beyond. Scientific Data 10.

GBIF.org. 2023. Occurrence data of Chiroptera download from GBIF.  https://doi.org/10.15468/dl.ywmbpq.

Graham MH. 2003. CONFRONTING MULTICOLLINEARITY IN ECOLOGICAL MULTIPLE REGRESSION. Ecology 84.

Ho L si T, Ané C. 2014. A linear-time algorithm for Gaussian and non-Gaussian trait evolution models. Systematic Biology 63: 397–408.

Jones B, O’Neill BC. 2020. Global One-Eighth Degree Population Base Year and Projection Grids Based on the Shared Socioeconomic Pathways, Revision 01.

Kummu M, Taka M, Guillaume JHA. Gridded global datasets for Gross Domestic Product and Human Development Index over 1990–2015. Scientific Data.

Luis AD, O’Shea TJ, Hayman DTS, Wood JLN, Cunningham AA, Gilbert AT, Mills JN, Webb CT. 2015. Network analysis of host-virus communities in bats and rodents reveals determinants of cross-species transmission. Ecology Letters 18: 1153–1162.

Myhrvold NP, Baldridge E, Chan B, Sivam D, Freeman DL, Ernest SKM. 2015. An amniote life‐history database to perform comparative analyses with birds, mammals, and reptiles. Ecology 96.

O’Mara MT, Wikelski M, Kranstauber B, Dechmann DKN. 2019. First three‐dimensional tracks of bat migration reveal large amounts of individual behavioral flexibility. Ecology 100.

O’Neill BC, Tebaldi C, van Vuuren DP, Eyring V, Friedlingstein P, Hurtt G, Knutti R, Kriegler E, Lamarque J-F, Lowe J, Meehl GA, Moss R, Riahi K, Sanderson BM. 2016. The Scenario Model Intercomparison Project (ScenarioMIP) for CMIP6. Geoscientific Model Development 9: 3461–3482.

Prener C, Grossenbacher T, Zehr A. 2022. biscale: Tools and Palettes for Bivariate Thematic Mapping. R package version 1.1.0.9000.

Randolph HE, Barreiro LB. 2018. Holy Immune Tolerance, Batman! Immunity 48: 1074–1076.

Ruiz-Aravena M, McKee C, Gamble A, Lunn T, Morris A, Snedden CE, Yinda CK, Port JR, Buchholz DW, Yeo YY, Faust C, Jax E, Dee L, Jones DN, Kessler MK, Falvo C, Crowley D, Bharti N, Brook CE, Aguilar HC, Peel AJ, Restif O, Schountz T, Parrish CR, Gurley ES, Lloyd-Smith JO, Hudson PJ, Munster VJ, Plowright RK. 2022. Ecology, evolution and spillover of coronaviruses from bats. Nature Reviews Microbiology 20: 299–314.

Schmitt S, Pouteau R, Justeau D, de Boissieu F, Birnbaum P. 2017. ssdm: An r package to predict distribution of species richness and composition based on stacked species distribution models. Methods in Ecology and Evolution 8: 1795–1803.

Teitelbaum CS, Amoroso CR, Huang S, Davies TJ, Nunn CL. 2020. A comparison of diversity estimators applied to a database of host–parasite associations. Ecography.

Valavi R, Guillera-Arroita G, José J. Lahoz-Monfort, Elith J. Predictive performance of presence-only species distribution models: a benchmark study with reproducible code. Ecological Monographs.

Wang J, Pan Y, Yang LF, Yang WH, Lv K, Luo C, Wang J, Kuang G, Wu WC, Gou QY. 2023. Individual bat virome analysis reveals co-infection and spillover among bats and virus zoonotic potential. Nature Communications 14.

Whittaker RH. 1980. Calssification of plant communities. Calssification of plant communities.

Zizka A, Silvestro D, Andermann T, Azevedo J, Duarte Ritter C, Edler D, Farooq H, Herdean A, Ariza M, Scharn R, Svantesson S, Wengström N, Zizka V, Antonelli A. 2019. CoordinateCleaner: Standardized cleaning of occurrence records from biological collection databases. Methods in Ecology and Evolution 10: 744–751.