Disease data source
All analyses were conducted at the administrative level, and the exact list of known zoonotic diseases is recorded in the GIDEON database22. GIDEON is currently the most comprehensive and frequently updated infectious disease outbreak database reporting epidemics of human infectious diseases at the global scale and has been widely used in global zoonosis studies42,43 (Last access date, November 9, 2020). The administrative designations used in our analyses were based on the Global Administrative Areas (GADM) database (www.gadm.org, downloaded on November 8, 2020), which includes very detailed boundary data for global countries and major island groups.
Pattern and correlates of zoonosis events worldwide
Number of zoonosis events
GIDEON defines human infectious disease reservoirs as any animal, plant, or substrate supporting the survival and reproduction of infectious agents and promoting transmission to potential susceptible hosts. Its host category therefore includes all human-specific, zoonotic, multihost, and environmental agents. As our main aim was to test the role of established alien animal species in the emergence of zoonotic diseases, we focused on a total of 161 diseases specified in GIDEON’s host designations and definitions as nonhuman zoonotic (n = 115) and multihost (n = 46) diseases (Supplementary Data 1) and excluded diseases with human-specific hosts that do not need animals to persist or be transmitted. The infectious agents of nonhuman zoonotic diseases complete their entire lifecycle in nonhuman hosts but may have the potential to spillover and infect human populations. Infectious agents of multihost diseases can use both human and animal hosts for their development and reproduction. We measured the number of zoonosis events for each jurisdiction according to five host taxonomic groups: mammals, birds, invertebrates, reptiles and amphibians. These zoonoses were mainly caused by bacteria, viruses, parasitic animals and fungi. We excluded zoonoses from the Algae (3 diseases) due to low sample sizes in GIDEON.
Correlates of the number of zoonosis events
Climatic variables
Following a previous study21, we used global environmental stratification (GEnS) as a composite bioclimatic variable generated by stratifying the Earth’s surface into zones with similar climates44. The GEnS database was constructed based on a total of 125 strata across 18 global environmental zones with a spatial resolution of 30 arc seconds (equivalent to approximately 0.86 km2 at the equator). The values in GEnS range from 1 to 18 with a higher value indicating warmer and wetter conditions.
Human population density
We used human population density as one general anthropogenic factor reflecting propagule pressure and human-assisted pathogen movements1,21,45. Human population size data and the land area of each jurisdiction were collected from World Bank Open Data from 2011 to 2020 (available at https://data.worldbank.org/indicator/SP.POP.TOTL, accessed on November 18, 2020). We then calculated the human population density using the human population size divided by the land area.
Native potential host richness and biodiversity loss
Data on the richness of native amphibians, birds, and mammals were derived from the Biodiversity Mapping website (https://biodiversitymapping.org/wordpress/index.php/home/, accessed on August 19, 2020), which were based on studies from Jenkins et al. (2013)’s and Pimm et al. (2014)46,47. The map of reptile diversity is based on an updated database of the global spatial distribution of reptiles48. All diversity maps for each taxon were generated through the calculation of grid-based richness at a spatial resolution of 10 km × 10 km in ArcGIS46. We did not include native invertebrate richness, as global maps for most invertebrate taxa are not yet available. For the loss of native biodiversity, we followed the previous study by first extracting the list of threatened species (NT, EN and VU categories evaluated by the IUCN Red List, access on May 10th, 2021)29, and then calculated the number of threatened species for each taxon distributes in each administrative unite as a proxy of biodiversity loss.
Richness of established alien zoonotic host species
We quantified the richness of established alien animal species from the five main taxonomic groups (mammals, birds, reptiles, amphibians and invertebrates) based on 4,522 establishment events of 795 alien animals in each of 201 jurisdictions according to various databases. Data on 262 established alien reptiles and amphibians were compiled from multiple publications, including Kraus’s compendium49 and other recent updates50. Data on 337 established alien birds after removing all migratory bird species as vagrants were collected from the Global Avian Invasions Atlas (GAVIA)51, which is a comprehensive database of the global distribution of established alien birds. Data on 119 established alien mammals were obtained from the Introduced Mammals of the World database52 and the more recent update53. Data on 77 terrestrial alien invertebrates (66 insects and 11 other groups) across 7 taxa with native and invaded range information were obtained from the Global Invasive Species Database (GISD, http://www.iucngisd.org/gisd/, accessed on July 1, 2020). We calculated the richness of both zoonotic and non-zoonotic alien host species for each order. We first conducted an intensive literature review for each established alien species of each of the four taxa to determine whether they transmit pathogens to humans (Supplementary Data 2). The identification of zoonotic or non-zoonotic host may be influenced by under-sampling in the literature. We therefore incorporated the latest synthesis of human-infecting pathogens in the ‘CLOVER’ dataset to identify zoonotic and non-zoonotic animal hosts54. The CLOVER dataset compiled GMPD255, EID256, HP323 and Shaw57 databases and is currently the most comprehensive dataset on host-pathogen associations. Based on this information, we then categorized each alien species as a ‘zoonotic host’ or ‘non-zoonotic host’. The records of the established alien species were assigned to GADM jurisdictions, and we calculated the richness of the established alien zoonotic and non-zoonotic host species for each taxonomic group within each jurisdiction. In order to increase the statistical power, we conducted subsequent modeling analyses based on four mammalian orders (i.e., Carnivora, Cetartiodactyla, Lagomorpha, and Rodentia), five avian groups (i.e., waterfowl including five orders: Anseriformes, Gruiformes, Pelecaniformes, Phoenicopteriformes and Suliformes; Columbiformes, Galliformes, Passeriformes, Psittaciformes), the order Diptera of the invertebrates, and herpetofauna as a whole, which have established alien populations in at least 50 administrative units.
Climate change
We extracted historical monthly mean temperature and precipitation data recorded between 1901 and 2009 from the University of East Anglia Climate Research Unit (CRU, https://sites.uea.ac.uk/cru/, accessed on November 30, 2020)58. This database provides historical global-scale yearly climatic data with the finest resolution of 0.5° grids. We generated the temperature and precipitation values for all grids in each jurisdiction, calculated the slope of the temperature and precipitation for the time series of the years 1901 to 2009 for each grid and generated the averages based on all grids within each jurisdiction.
Anthropogenic land-use change
We downloaded global land-use data from the Anthromes v2 Dataset (Anthropogenic Biomes version 2, accessed on October 15, 2020) in ESRI GRID format59. We used the 1900 and 2000 data to calculate the temporal changes in land use. By using the reclassify and raster function in ArcGIS, we calculated the percentage of grids in which the land-use type changed to a more anthropogenically influenced type from 1900 to 2000 for each jurisdiction, including 15 scenarios: Wildlands to Seminatural, Wildlands to Rangelands, Wildlands to Croplands, Wildlands to Villages, Wildlands to Dense Settlements, Seminatural to Rangelands, Seminatural to Croplands, Seminatural to Villages, Seminatural to Dense Settlements, Rangelands to Croplands, Rangelands to Villages, Rangelands to Dense Settlements, Croplands to Villages, Croplands to Dense Settlements, and Villages to Dense Settlements.
Sampling effort, reporting bias and incomplete data
A potential issue in quantifying the effects of different predictor variables on the number of zoonosis events is the need to account for the differences in survey effort, reporting bias and incomplete disease data among regions1,21,28. There is a high probability that zoonosis discovery is spatially biased by uneven levels of surveillance across countries, as the global allocation of scientific resources has been focused on rich and developed countries. We thus included the Infectious Disease Vulnerability Index (IDVI), which is a comprehensive metric reflecting the demographic, health care, public health, socioeconomic, and political factors that may have an impact on the capacity of surveillance and detection of infectious diseases in each country60. Second, we followed the methods of a previous study21 to control for reporting biases. We incorporated PubMed citations per disease for each jurisdiction using a Python-based PubCrawler21. In addition, we added the longitude and latitude of the geographic centroid of administrative units to control for spatial autocorrelation as there would be a higher probability of having similar diseases in nearby than distant administrative units61.
Statistical analysis
The number of zoonosis events, native potential host richness, established alien animal richness and human population density were log-transformed to improve linearity. A potential issue in our data analysis is that the numbers of zoonosis events and the numbers of native and alien animal species are strongly influenced by geographical area, as larger countries or regions may host more native or alien animal species and more disease events. We therefore calculated the density of native or alien species richness and the number of zoonosis events using the total number divided by the geographical area of each jurisdiction. Furthermore, the number of zoonosis events may also be influenced by the degree of local disease surveillance. We thus obtained the residuals from a regression correlating zoonosis event density and all disease event density, and used them as the dependent variable for further analyses (Fig. 1). As some of our variables may be expected to be nonlinear, we performed generalized additive mixed model (GAMM) analyses following Mollentze & Streicker 2020’s framework25 to quantify the relationships between different predictor variables and the number of zoonosis events. We started with a full model with zoonosis event density controlling for overall disease surveillance as the response variable and 13 smoothed fixed effects (Fig. 1 and Supplementary Data 4): GEnS, human population density, density of native species richness, biodiversity loss, density of alien zoonotic host richness, density of alien non-zoonotic host richness, climate (temperature and precipitation) change, land-use change, IDVI, PubMed citations, longitude and latitude of geographic centroid of administrative units. The reason why we included the density of alien non-zoonotic host richness as a covariate is because this variable can serve as a positive control for propagule pressure, allowing us to more explicitly test whether zoonotic alien hosts contribute to zoonoses beyond propagule pressure associated with non-zoonotic alien hosts, which cannot directly increase zoonotic diseases. These predictor variables were not highly collinear as their correlation coefficients based on Pearson rank correlation analyses were all <0.65 (Supplementary Fig. 6). Because human history may have a great influence on disease outbreaks, as there may be more human pathogens on continents subject to earlier human settlement45, we followed this literature by including continental identity (i.e., Africa, region of origin and first settlements; most of Europe, Asia, Australia, New Guinea, by approx. 40,000–60,000 BP; North America, by approx. 15,000–25,000 BP; South America, by approx. 1000–5000 years after North America, i.e., 10,000–24,000 BP; and most islands, by approx. 1000–7000 BP) as a random intercept to control for potential pseudoreplication. In addition, to account for the lack of complete independence among disease events caused by the same pathogens or associated with the same introduced host orders (Fig.1), we included pathogen identity and host order as two additional random intercepts. We fitted all models using restricted maximum likelihood method and ranked all candidate models by the Akaike’s Information Criterion (AIC) theoretic approach25,62. Models including all possible combinations of the 13 predictor variables (total 213-1 = 8191 models) were ranked, and the models within 2 AIC unites (i.e., ΔAIC ≤ 2) compared with the top model were considered to be highly supported62. For each model, we computed the standardized estimates of the regression coefficients of the predictor variables with the 95% confidence intervals (CIs), and considered effects statistically significant when the 95% CIs did not overlap zero. In addition to variable significance, we also calculated the proportion of the deviance explained by each predictor variable by comparing the sub-models in the absence of the variable to the full model25. To better compare the coefficients of the different covariates, we standardized each of the predictor variables to a mean of zero and standard deviation of one before it was entered into the model63. Furthermore, we used different levels of thin-plate smoothers with 6, 8, and 10 knots for the fixed-effect variables in GAMMs25. As the results were similar regardless of which knot we used, we present the results from the analyses with 10 knots in the main text but provide the results with 6 and 8 knots in the supporting materials (Supplementary Fig. 3). All analyses were conducted in the gamm4, mgcv, visreg, dplyr, and MuMIn packages in R version 4.0.364 (Supplementary Notes).
To further test whether there are different responses of the zoonosis emergences among alien host groups, we fit an interaction between host order and alien species richness to investigate whether the effect of alien zoonotic host richness varied across taxonomic groups. Finally, to evaluate the potential contribution of alien zoonotic host introductions on historical zoonosis events for each administrative unit, we generated the fitted values along with 95% CIs of the number of zoonosis events in each administrative unit predicted by using the predictors in GAMMs, and subtracted the predicted values excluding zoonotic host introductions from those using all predictor variables.
The observed spatial correlation of zoonosis emergences with alien animal invasions might be problematic because there may be mismatch in the occurrence of alien animal invasions and zoonosis events in time. For instance, zoonosis events at a location might have occurred earlier than the alien animal invasions, which could not have caused the disease despite a strong spatial correlation. We therefore further explored the temporal relationship of alien zoonotic (and non-zoonotic) host introductions and zoonotic diseases over years. To achieve this, we collected the introduction time of each alien zoonotic and non-zoonotic host species in each country or region based on the alien avian introduction database and literatures from alien birds51,63, the Introduced Mammals of the World database52 for alien mammals, and the Global Invasive Species Database (GISD) for alien invertebrates. We obtained the number of new zoonosis events over time from the GIDEON database (Last access date, November 9, 2020), and calculated the number of zoonosis events in each year for different taxa. For each year, the number of new zoonosis events divided by the number of alien host introductions for each country or region was used to evaluate the magnitude of the relationship between alien zoonotic host introductions and the number of new zoonosis events. We then used two approaches to analyze their temporal relationships. Firstly, we conducted a generalized additive modeling analysis where we treated year as the replicate, and zoonotic and non-zoonotic alien introductions as the independent variables, and the number of zoonosis events as the dependent variable, to explore the general trend of alien animal host introductions and the number of zoonosis events along time. Additionally, we applied breakpoint regression analyses using the segmented package in R (Supplementary Notes). In these analyses, the identified the breakpoint reflects the year in which there was a rapid increase in the number of zoonotic or non-zoonotic host species or the number of zoonosis events. For these analyses, we fit the left-horizontal regression and two-slope regression that are widely used in ecological and biogeographical studies65, and applied an AIC-based approach to identify the optimal breakpoint year. We combined data across various host taxa for the temporal analyses because some orders had little data on the timing of alien species establishment.
Reporting summary
Further information on research design is available in the Nature Research Reporting Summary linked to this article.
Source: Ecology - nature.com
