Philippa Kaur

Detection and validation of high-quality jumbo phage binsDue to the large size of jumbo bacteriophage genomes, it is likely that they are present in multiple distinct contigs in metagenomic datasets and therefore require binning to recover high-quality metagenome-assembled genomes (MAGs) [28]. This has been shown for large DNA viruses that infect eukaryotes, where several recent studies have successfully employed binning approaches to recover viral MAGs [2, 3, 30]. Here, we used the same 1545 high-quality metagenomic assemblies [31] used in a recent study to recover giant viruses of eukaryotes [3], but we modified the bioinformatic pipeline to identify bins of jumbo bacteriophages. These metagenomes were compiled by Parks et al. [31] and included available metagenomes on the NCBI’s Short Read Archive by December 31, 2015 (see Parks et al. [31]). This dataset includes a wide variety of marine metagenomes (n = 469) including many non-Tara metagenomes (n = 165). We focused our benchmarking and distribution analyses on Tara data [29] because of the well-curated metadata and size fractions in this dataset. We first binned the contigs from these assemblies with MetaBat2 [32], which groups contigs together based on similar nucleotide composition and coverage profiles, and focused on bins of at least 200 kilobases in total length. We subsequently identified bins composed of bacteriophage contigs through analysis with VirSorter2 [33], VIBRANT [34], and CheckV [35] (see Methods for details).The occurrence of multiple copies of highly conserved marker genes is typically used to assess the level of contamination present in metagenome-derived genomes of bacteria and archaea [36]. Because bacteriophage lack these marker genes [37], we developed alternative strategies to assess possible contamination in our jumbo phage bins. Firstly, we refined the set of bins by retaining those with no more than 5 contigs that were each at least 5 kilobases in length to reduce the possibility that spurious contigs were put together. Secondly, we assessed the possibility that two strains of smaller phages with similar nucleotide composition may be binned together by aligning the contigs in a bin to each other. Bins that had contigs with high sequence similarity across the majority of their lengths were discarded (Supplementary Fig. 1). Thirdly, we discarded bins if their contigs exhibited aberrant co-abundance profiles in different metagenomes (see Supplementary Methods). To generate these co-abundance profiles, we mapped reads from 225 marine metagenomes provided by Tara Oceans [29] onto the bins. Coverage variation between contigs was benchmarked based on read-mapping results from artificially-fragmented reference genomes present in the samples (See Methods for details). Only bins with coverage variation below our empirically-derived threshold were retained. Using this stringent filtering, we identified 85 bins belonging to jumbo bacteriophages. These bins ranged in length from 202 kbp to 498 kbp, and 31 consisted of a single contig, while 54 consisted of 2–5 contigs (Supplementary Fig. 2).To assess global diversity patterns of jumbo bacteriophages, we combined these jumbo phage bins together with a compiled database of previously-identified jumbo phages that included all complete jumbo Caudovirales genomes on RefSeq (downloaded July 5th, 2020), the INPHARED database [14], a recent survey of cultivated jumbo phages [6], the Al-Shayeb et al. study [4], and marine jumbo phage contigs from metagenomic surveys of GOV 2.0 [26] (60 jumbo phages), ALOHA 2.0 [38] (8 jumbo phages), and one megaphage MAG recovered from datasets of the English Channel [39]. Ultimately, we arrived at a set of 244 jumbo phages, including the 85 bins, that were present in at least one Tara Oceans sample (min. 20% genome covered, see Methods) or deriving from a marine dataset (i.e. ALOHA, GOV 2.0) which we analyzed further in this study and refer to as marine jumbo phages. Statistics on genomic features can be found in Supplementary Dataset 1.Marine jumbo phages belong to distinct groups with diverse infection strategiesBecause bacteriophages lack high-resolution, universal marker genes for classification, such as 16S rRNA in bacteria, phages are often grouped by gene content [40, 41]. Here, we generated a bipartite network that included the 85 bins of jumbo phages with a dataset of available Caudovirales complete genomes in RefSeq (3012 genomes; downloaded July 5th, 2020) and the full set of reference jumbo phages described above. To construct the bipartite network, we compared proteins encoded in all the phage genomes to the VOG database, and each genome was linked to VOG hits that were present (Fig. 1, Supplementary Dataset 2, see Methods for details). To identify groups of phage genomes with similar VOG profiles, we employed a spinglass community detection algorithm [42] to generate genome clusters. Similar methods have recently been used to analyze evolutionary relationships in other dsDNA viruses [41]. The marine jumbo phages of this study clustered into five groups that included both jumbo and non-jumbo phage genomes (Fig. 2a). We refer to these five clusters as Phage Genome Clusters (PGCs): PGC_A, PGC_B, PGC_C, PGC_D, and PGC_E. These PGCs included cultured phages and metagenome-derived jumbo phages found in various environments (i.e. aquatic, engineered) and isolated on a diversity of hosts (i.e. Firmicutes, Proteobacteria, Bacteroidetes) (Fig. 2b, c). Of the marine jumbo phages, 135 belonged to PGC_A, 11 to PGC_B, 90 to PGC_C, 7 to PGC_D, and 1 to PGC_E (Fig. 1b). In addition to this network-based analysis, we also examined phylogenies of the major capsid protein (MCP) and the terminase large subunit (TerL) encoded by the marine jumbo phages and the same reference phage set examined in the network (Fig. 1c, d). With the exception of PGC_A, the marine jumbo phages that belong to the same PGC appeared more closely related to each other than those belonging to different clusters. The polyphyletic placement of jumbo phage PGCs in these marker gene phylogenies is consistent with the view that genome gigantism evolved multiple times, independently within the Caudovirales [6].Fig. 1: Bipartite network and marker gene analyses of jumbo phages.a Network with marine jumbos and references as nodes and edges based on shared VOGs. Marine jumbo phage nodes are colored by PGC as detected with spinglass community detection analysis, other nodes are in gray. Edges and VOG nodes have been omitted to more clearly represent the pattern of phage clustering. Node size corresponds to the natural log of genome length in kilobases. b Barchart of the number of members in each PGC. PGCs with marine jumbo phages are denoted with a star and the number of marine jumbo phages in that PGC. Proportion of marine jumbo phages in that PGC is colored. Phylogenies of TerL (c) and MCP (d) proteins with references and bins. Inner ring and branches are colored by the 5 PGCs that marine jumbo phages belong to. Navy blue circles in the outer ring denote marine jumbo phages.Full size imageFig. 2: Statistics of the Phage Genomes Clusters (PGCs).a Boxplot of genome length in each network cluster (x-axis is PGC number). Star denotes PGC with a marine jumbo phage and the color matches the PGC letters of Fig. 1. b Stacked barplot of the metagenome environment from which each phage derives from in each PGC (x-axis). Reference (yellow) are cultured phages, in black are the bins of jumbo phages from this study. c Stacked barplot of the host phylum of the RefSeq cultured phages in each cluster; metagenomic phages are in gray.Full size imageWe then compared functional content encoded by the marine jumbo phages in the PGCs to identify functional differences that distinguish these groups. PGC_E was excluded from this analysis because this genome cluster contained only a single jumbo phage. Collectively, most genes of the marine jumbo phages could not be assigned a function (mean: 86.60%, std dev: 7.01%; Supplementary Dataset 3), which is common with environmental viruses [43, 44]. Genes with known functions primarily belonged to functional categories related to viral replication machinery, such as information processing and virion structure (Fig. 3a), and these genes drove the variation between the genome clusters of marine jumbo phages (Fig. 3b). A recent comparative genomic analysis of cultivated jumbo phages was able to identify three types of jumbo phages that are defined by different infection strategies and host interactions (referred to as Groups 1–3) [6]. We cross-referenced our PGCs and found that PGCs B, C, and D of this study corresponded to Groups 1, 2, and 3, respectively, suggesting that these genome clusters contain phages with distinct infection and replication strategies. PGC_A corresponded to multiple groups, indicating that this genome cluster contains a particularly broad diversity of phages.Fig. 3: Functional predictions of PGCs.a Functional categories for genes encoded by jumbo phages averaged by PGC. b Heatmap of proportion of genomes in each PGC that contain the listed genes. Listed genes were selected based on containing a known function and having a variance between the PGCs above 0.2. Dendrogram was generated based on hierarchical clustering in pheatmap.Full size imageThe second largest phage cluster with marine jumbo phages, PGC_B, consists of 238 phages (11 (4.6%) marine jumbo phages, including 10 bins generated here), and included cultured phages of Group 1, which is typified by Pseudomonas aeruginosa phage PhiKZ. Supporting this correspondence with Group 1, all marine jumbo phages of PGC_B encoded the same distinct replication and transcription machinery, including a divergent family B DNA polymerase and a multi-subunit RNA polymerase (Fig. 3b, Supplementary Dataset 3). These marine jumbo phages also encoded a PhiKZ internal head protein, and they uniquely encoded shell and tubulin homologs which has recently been found in PhiKZ phages to assist in the formation of a nucleus-like compartment during infection that protects the replicating phage from host defenses [18, 19]. Although we could not confidently predict hosts for the 11 metagenomic marine jumbo phages in this PGC_B (Supplementary Dataset 1), the cultured phages of this genome cluster infect pathogenic bacteria belonging to the phyla Proteobacteria (178 phages) and Firmicutes (6 phages) (Fig. 2c), implicating a potential host range for marine jumbo phages in PGC_B.The next largest phage genome cluster, PGC_C, comprised of 156 phages total (90 marine jumbo phages (57.7%); 4 bins generated from this study) and included reference jumbo phages in Group 2 (31, 19.9%) which are typified by Alphaproteobacteria and Cyanobacteria phages. Likewise, the host range of other cultured phages in PGC_C support the Group 2 correspondence, either infecting Cyanobacteria (139 phages) or Proteobacteria (4 phages) (Fig. 2c). Furthermore, all 3 marine metagenomic phages in PGC_C for which hosts could be predicted were matched to Cyanobacteria hosts (Supplementary Dataset 1). Functional annotations of PGC_C marine jumbo phages revealed nearly all encode a family B DNA polymerase (97.8% of phages) and the photosystem II D2 protein (PF00124, VOG04549) characteristic of cyanophages (90% of phages) (Fig. 3b). This PGC included the reference Prochlorococcus phage P-TIM68 (NC_028955.1), which encodes components of both photosystem I and II as a mechanism to enhance cyclic electron flow during infection [45]. This suggests that an enhanced complement of genes used to manipulate host physiology during infection may be a driver of large genome sizes in this group. Additionally, most of the PGC_C marine jumbo phages encoded lipopolysaccharide biosynthesis proteins (76%), which have been found in cyanophage genomes that may induce a “pseudolysogeny” state, when infected host cells are dormant, by changing the surface of the host cell and preventing additional phage infections [6] (Supplementary Dataset 3). Taken together, most marine jumbo phages of PGC_C likely follow host interactions of jumbo cyanophages, such as potentially manipulating host metabolism by encoding their own photosynthetic genes and potentially inducing a pseudolysogenic state.Finally, phages of PGC_D totaled at 47 phages, of which 7 were marine jumbo phages generated in this study (14.9%). This group included Group 3 jumbo phages (15, 31.9%), which is primarily distinguished by encoding a T7-type DNA polymerase but is not typified by a particular phage type or host range. Supporting this grouping, all marine jumbo phages in this study encoded a T7 DNA polymerase which belongs to family A DNA polymerases (Fig. 3b, Supplementary Dataset 3). Most of the other genes distinctively encoded by the marine jumbo phages in this group included structural genes related to T7 (T7 baseplate, T7 capsid proteins), a eukaryotic DNA topoisomerase I catalytic core (PF01028), and DNA structural modification genes (MmcB-like DNA repair protein, DNA gyrase B). Hosts of metagenomic marine jumbo phages in PGC_D could not be predicted (Supplementary Dataset 1); however, cultured Group 3 jumbo phages in PGC_D all infect Proteobacteria, primarily Enterobacteria and other pathogens.The largest of the phage genome clusters, PGC_A, contained 469 phages, including 135 marine jumbo phages (63 bins from this study). This genome cluster contained the largest jumbo phages, such as Bacillus phage G (498 kb) and the marine megaphage Mar_Mega_1 (656 kb) recently recovered from the English Channel [39]. Unlike other PGCs, PGC_A contained mostly metagenomic phages (401, 85%, Fig. 2b, c). Considering PGC_A contains the largest jumbo phages (Figs. 1b, 2a), the vast genetic diversity in this PGC might explain why few genes were found to distinguish this group. Of the genes unique to PGC_A, only one was present in at least half of the phages (51.9%), which was a Bacterial DNA polymerase III alpha NTPase domain (PF07733). The host ranges of cultured phages from this PGC further reflect the large diversity of this group and included a variety of phyla and genera that can perform complex metabolisms or lifestyles, such as the nitrogen-fixing Cyanobacteria of the Nodularia genus isolated from the Baltic Sea (accessions NC_048756.1 and NC_048757.1) and the Bacteroidetes bacteria Rhodothermus isolated from a hot spring in Iceland (NC_004735.1) [46]. Because this group contains an abundance of metagenome-derived genomes that encode mostly proteins with no VOG annotation (Supplementary Dataset 2), it is possible that it includes several distinct lineages that could not be distinguished using the community detection algorithm of the bipartite network analysis.Relative abundance of jumbo bacteriophages across size fractionsTo explore the distribution of the marine jumbo phages in the ocean, we first examined the size fractions in which the jumbo phages were most prevalent. To remove redundancy for the purposes of read mapping, we examined the 244 jumbo phages at the population-level ( >80% genes shared with >95% average nucleotide identity [24]), corresponding to 142 populations (11 unique to this study, corresponding to 47 bins). We then mapped Tara Oceans metagenomes onto the 142 jumbo phage populations, and 102 of these populations could be detected [min. 20% of genome covered], resulting in 74 populations in PGC_A, 2 in PGC_B, 22 in PGC_C, 3 in PGC_D, and 1 in PGC_E. Out of the 225 Tara Oceans metagenomes examined, 213 (94.6%) contained at least one jumbo phage population (median: 7, Supplementary Dataset 4). Jumbo phages were more frequently detected in samples below 0.22 µm ( More

Cardinale, B. J., Palmer, M. A. & Collins, S. L. Species diversity enhances ecosystem functioning through interspecific facilitation. Nature 415, 426–429 (2002).ADS 
CAS 
PubMed 

Google Scholar 
Gagic, V. et al. Functional identity and diversity of animals predict ecosystem functioning better than species-based indices. Proc. R. Soc. Lond. B Biol. Sci. 282, 20142620 (2015).
Google Scholar 
Oliver, T. H. et al. Biodiversity and resilience of ecosystem functions. Trends Ecol. Evol. 30, 673–684 (2015).PubMed 

Google Scholar 
Oliver, T. H. et al. Declining resilience of ecosystem functions under biodiversity loss. Nat. Commun. 6, 10122 (2015).ADS 
PubMed 

Google Scholar 
Wardle, D. A. Do experiments exploring plant diversity–ecosystem functioning relationships inform how biodiversity loss impacts natural ecosystems?. J. Veg. Sci. 27, 646–653 (2016).
Google Scholar 
Cadotte, M. W., Carscadden, K. & Mirotchnick, N. Beyond species: Functional diversity and the maintenance of ecological processes and services. J. Appl. Ecol. 48, 1079–1087 (2011).
Google Scholar 
Petchey, O. L. & Gaston, K. J. Extinction and the loss of functional diversity. Proc. R. Soc. Lond. B Biol. Sci. 269, 1721–1727 (2002).
Google Scholar 
Rosenfeld, J. S. Functional redundancy in ecology and conservation. Oikos 98, 156–162 (2002).
Google Scholar 
Fonseca, C. R. & Ganade, G. Species functional redundancy, random extinctions and the stability of ecosystems. J. Ecol. 89, 118–125 (2001).
Google Scholar 
Mouillot, D. et al. Functional over-redundancy and high functional vulnerability in global fish faunas on tropical reefs. PNAS 111, 13757–13762 (2014).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Violle, C. et al. Let the concept of trait be functional!. Oikos 116, 882–892 (2007).
Google Scholar 
Laughlin, D. C. Applying trait-based models to achieve functional targets for theory-driven ecological restoration. Ecol. Lett. 17, 771–784 (2014).PubMed 

Google Scholar 
Laughlin, D. C., Strahan, R. T., Huffman, D. W. & Sánchez Meador, A. J. Using trait-based ecology to restore resilient ecosystems: Historical conditions and the future of montane forests in western North America. Restor. Ecol. 25, S135–S146 (2017).
Google Scholar 
Petchey, O. L. & Gaston, K. J. Functional diversity (FD), species richness and community composition. Ecol. Lett. 5, 402–411 (2002).
Google Scholar 
Carmona, C. P., de Bello, F., Mason, N. W. H. & Lepš, J. Traits without borders: Integrating functional diversity across scales. Trends Ecol. Evol. 31, 382–394 (2016).PubMed 

Google Scholar 
Jain, M. et al. The importance of rare species: A trait-based assessment of rare species contributions to functional diversity and possible ecosystem function in tall-grass prairies. Ecol. Evol. 4(104), 112 (2014).
Google Scholar 
Mouillot, D. et al. Rare species support vulnerable functions in high-diversity ecosystems. PLoS Biol. 11, e1001569 (2013).CAS 
PubMed 
PubMed Central 

Google Scholar 
Leitão, R. P. et al. Rare species contribute disproportionately to the functional structure of species assemblages. Proc. R. Soc. B 283, 20160084 (2016).PubMed 
PubMed Central 

Google Scholar 
IUCN/SSC. Guidelines for Reintroductions and Other Conservation Translocations. (IUCN Species Survival Commission, 2013).Bakker, E. S. & Svenning, J.-C. Trophic rewilding: Impact on ecosystems under global change. Philos. Trans. R. Soc. B Biol. Sci. 373, 20170432 (2018).
Google Scholar 
Garrido, P. et al. Experimental rewilding enhances grassland functional composition and pollinator habitat use. J. Appl. Ecol. 56, 946–955 (2019).
Google Scholar 
Svenning, J.-C. et al. Science for a wilder Anthropocene: Synthesis and future directions for trophic rewilding research. Proc. Natl. Acad. Sci. 113, 898–906 (2016).ADS 
CAS 
PubMed 

Google Scholar 
Estes, J. A. et al. Trophic downgrading of planet earth. Science 333, 301–306 (2011).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Ritchie, E. G. et al. Ecosystem restoration with teeth: What role for predators?. Trends Ecol. Evol. 27, 265–271 (2012).PubMed 

Google Scholar 
Chauvenet, A. L. M., Canessa, S. & Ewen, J. G. Setting objectives and defining the success of reintroductions. In Reintroduction of Fish and Wildlife Populations 105–121 (University of California Press, 2016).Ewen, J. G., Soorae, P. S. & Canessa, S. Reintroduction objectives, decisions and outcomes: Global perspectives from the herpetofauna. Anim. Conserv. 17, 74–81 (2014).
Google Scholar 
Kleiman, D. G., Price, M. R. S. & Beck, B. B. Criteria for reintroductions. In Creative Conservation: Interactive Management of Wild and Captive Animals (eds. Olney, P. J. S., Mace, G. M. & Feistner, A. T. C.) 287–303 (Springer Netherlands, 1994). https://doi.org/10.1007/978-94-011-0721-1_14.Hunter, M. L. & Hutchinson, A. The virtues and shortcomings of parochialism: Conserving species that are locally rare, but globally common. Conserv. Biol. 8, 1163–1165 (1994).
Google Scholar 
Brichieri-Colombi, T. A. & Moehrenschlager, A. Alignment of threat, effort, and perceived success in North American conservation translocations. Conserv. Biol. 30, 1159–1172 (2016).PubMed 

Google Scholar 
Thévenin, C., Mouchet, M., Robert, A., Kerbiriou, C. & Sarrazin, F. Reintroductions of birds and mammals involve evolutionarily distinct species at the regional scale. PNAS https://doi.org/10.1073/pnas.1714599115 (2018).Article 
PubMed 
PubMed Central 

Google Scholar 
Seddon, P. J., Soorae, P. S. & Launay, F. Taxonomic bias in reintroduction projects. Anim. Conserv. 8, 51–58 (2005).
Google Scholar 
Thévenin, C., Morin, A., Kerbiriou, C., Sarrazin, F. & Robert, A. Heterogeneity in the allocation of reintroduction efforts among terrestrial mammals in Europe. Biol. Conserv. 241, 108346 (2020).
Google Scholar 
Devictor, V. et al. Spatial mismatch and congruence between taxonomic, phylogenetic and functional diversity: The need for integrative conservation strategies in a changing world. Ecol. Lett. 13, 1030–1040 (2010).PubMed 

Google Scholar 
Crees, J. J., Turvey, S. T., Freeman, R. & Carbone, C. Mammalian tolerance to humans is predicted by body mass: Evidence from long-term archives. Ecology 100, e02783 (2019).PubMed 

Google Scholar 
Sandom, C., Faurby, S., Sandel, B. & Svenning, J.-C. Global late Quaternary megafauna extinctions linked to humans, not climate change. Proc. R. Soc. B Biol. Sci. 281, 20133254 (2014).
Google Scholar 
Wilman, H. et al. EltonTraits 1.0: Species-level foraging attributes of the world’s birds and mammals. Ecology 95, 2027–2027 (2014).
Google Scholar 
Dı́az, S. & Cabido, M. Vive la différence: Plant functional diversity matters to ecosystem processes. Trends Ecol. Evol. 16, 646–655 (2001).
Google Scholar 
Mlambo, M. C. Not all traits are ‘functional’: Insights from taxonomy and biodiversity-ecosystem functioning research. Biodivers. Conserv. 23, 781–790 (2014).
Google Scholar 
van der Plas, F. et al. Plant traits alone are poor predictors of ecosystem properties and long-term ecosystem functioning. Nat. Ecol. Evol. 4, 1602–1611 (2020).PubMed 

Google Scholar 
Lavorel, S. & Garnier, E. Predicting changes in community composition and ecosystem functioning from plant traits: Revisiting the Holy Grail. Funct. Ecol. 16, 545–556 (2002).
Google Scholar 
Luck, G. W., Lavorel, S., McIntyre, S. & Lumb, K. Improving the application of vertebrate trait-based frameworks to the study of ecosystem services. J. Anim. Ecol. 81, 1065–1076 (2012).PubMed 

Google Scholar 
Mouchet, M. et al. Towards a consensus for calculating dendrogram-based functional diversity indices. Oikos 117, 794–800 (2008).
Google Scholar 
Podani, J. & Schmera, D. On dendrogram-based measures of functional diversity. Oikos 115, 179–185 (2006).
Google Scholar 
Maire, E., Grenouillet, G., Brosse, S. & Villéger, S. How many dimensions are needed to accurately assess functional diversity? A pragmatic approach for assessing the quality of functional spaces. Glob. Ecol. Biogeogr. 24, 728–740 (2015).
Google Scholar 
Villéger, S., Maire, E. & Leprieur, F. On the risks of using dendrograms to measure functional diversity and multidimensional spaces to measure phylogenetic diversity: A comment on Sobral et al. (2016). Ecol. Lett. 20, 554–557 (2017).PubMed 

Google Scholar 
Tsirogiannis, C. & Sandel, B. PhyloMeasures: A package for computing phylogenetic biodiversity measures and their statistical moments. Ecography 39, 709–714 (2016).
Google Scholar 
Isaac, N. J., Turvey, S. T., Collen, B., Waterman, C. & Baillie, J. E. Mammals on the EDGE: Conservation priorities based on threat and phylogeny. PLoS ONE 2, e296 (2007).ADS 
PubMed 
PubMed Central 

Google Scholar 
Paradis, E., Claude, J. & Strimmer, K. APE: Analyses of phylogenetics and evolution in R language. Bioinformatics 20, 289–290 (2004).CAS 

Google Scholar 
Hidasi-Neto, J., Loyola, R. & Cianciaruso, M. V. Global and local evolutionary and ecological distinctiveness of terrestrial mammals: Identifying priorities across scales. Divers. Distrib. 21, 548–559 (2015).
Google Scholar 
Dovrat, G., Meron, E., Shachak, M., Golodets, C. & Osem, Y. The relative contributions of functional diversity and functional identity to ecosystem function in water-limited environments. J. Veg. Sci. 30, 427–437 (2019).
Google Scholar 
Funk, J. L. et al. Revisiting the Holy Grail: Using plant functional traits to understand ecological processes. Biol. Rev. 92, 1156–1173 (2017).PubMed 

Google Scholar 
Kuebbing, S. E. & Bradford, M. A. The potential for mass ratio and trait divergence effects to explain idiosyncratic impacts of non-native invasive plants on carbon mineralization of decomposing leaf litter. Funct. Ecol. 33, 1156–1171 (2019).
Google Scholar 
Devictor, V. et al. Defining and measuring ecological specialization. J. Appl. Ecol. 47, 15–25 (2010).
Google Scholar 
Byers, J. E. et al. Using ecosystem engineers to restore ecological systems. Trends Ecol. Evol. 21, 493–500 (2006).PubMed 

Google Scholar 
Jones, C. G., Lawton, J. H. & Shachak, M. Organisms as ecosystem engineers. In Ecosystem Management: Selected Readings (eds. Samson, F. B. & Knopf, F. L.) 130–147 (Springer, 1996). https://doi.org/10.1007/978-1-4612-4018-1_14.Macdonald, D. W. et al. Reintroducing the beaver (Castor fiber) to Scotland: A protocol for identifying and assessing suitable release sites. Anim. Conserv. 3, 125–133 (2000).
Google Scholar 
Wilmers, C. C., Crabtree, R. L., Smith, D. W., Murphy, K. M. & Getz, W. M. Trophic facilitation by introduced top predators: Grey wolf subsidies to scavengers in Yellowstone National Park. J. Anim. Ecol. 72, 909–916 (2003).
Google Scholar 
Dupont, H., Mihoub, J.-B., Bobbé, S. & Sarrazin, F. Modelling carcass disposal practices: Implications for the management of an ecological service provided by vultures. J. Appl. Ecol. 49, 404–411 (2012).
Google Scholar 
Moleon, M. et al. Humans and scavengers: The evolution of interactions and ecosystem services. Bioscience 64, 394–403 (2014).
Google Scholar 
Legras, G., Loiseau, N., Gaertner, J.-C., Poggiale, J.-C. & Gaertner-Mazouni, N. Assessing functional diversity: The influence of the number of the functional traits. Theor. Ecol. 13, 117–126 (2020).
Google Scholar 
Petchey, O. L. & Gaston, K. J. Functional diversity: Back to basics and looking forward. Ecol. Lett. 9, 741–758 (2006).PubMed 

Google Scholar 
Díaz, S. et al. The global spectrum of plant form and function. Nature 529, 167–171 (2016).ADS 
PubMed 

Google Scholar 
Lundgren, E. J. et al. Introduced herbivores restore Late Pleistocene ecological functions. Proc. Natl. Acad. Sci. 117, 7871–7878 (2020).CAS 
PubMed 
PubMed Central 

Google Scholar 
Malhi, Y. et al. Megafauna and ecosystem function from the Pleistocene to the Anthropocene. PNAS 113, 838–846 (2016).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Osborne, P. E. & Seddon, P. J. Selecting suitable habitats for reintroductions: Variation, change and the role of species distribution modelling. Reintrod. Biol. Integr. Sci. Manag. 1, 73–104 (2012).
Google Scholar 
Lipsey, M. K., Child, M. F., Seddon, P. J., Armstrong, D. P. & Maloney, R. F. Combining the fields of reintroduction biology and restoration ecology. Conserv. Biol. 21, 1387–1390 (2007).PubMed 

Google Scholar 
Perino, A. et al. Rewilding complex ecosystems. Science 364, eaav5570 (2019).CAS 
PubMed 

Google Scholar 
Loiseau, N. et al. Global distribution and conservation status of ecologically rare mammal and bird species. Nat. Commun. 11, 5071 (2020).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Cadotte, M. W. & Tucker, C. M. Difficult decisions: Strategies for conservation prioritization when taxonomic, phylogenetic and functional diversity are not spatially congruent. Biol. Conserv. 225, 128–133 (2018).
Google Scholar 
Sarrazin, F. & Barbault, R. Reintroduction: Challenges and lessons for basic ecology. Trends Ecol. Evol. (Amst.) 11, 474–478 (1996).CAS 

Google Scholar  More

Hey, J. On the arbitrary identification of real species. In Speciation and Patterns of Diversity (eds Butlin, R. K. et al.) 15–28 (Cambridge University Press, 2009).
Google Scholar 
Arbogast, B. S., Edwards, S. V., Wakeley, J., Beerli, P. & Slowinski, J. B. Estimating divergence times from molecular data on phylogenetic and population genetic timescales. Annu. Rev. Ecol. Syst. 33, 707–740 (2002).
Google Scholar 
Nielsen, R. & Wakeley, J. Distinguishing migration from isolation: A Markov chain Monte Carlo approach. Genetics 158, 885–896 (2001).CAS 
PubMed 
PubMed Central 

Google Scholar 
Wakeley, J. The effects of subdivision on the genetic divergence of populations and species. Evolution 54, 1092–1101 (2000).CAS 
PubMed 

Google Scholar 
Hey, J. & Nielsen, R. Multilocus methods for estimating population sizes, migration rates and divergence time, with applications to the divergence of Drosophila pseudoobscura and D. persimilis. Genetics 167, 747–760 (2004).CAS 
PubMed 
PubMed Central 

Google Scholar 
Hey, J. Isolation with migration models for more than two populations. Mol. Biol. Evol. 27, 905–920 (2010).CAS 
PubMed 

Google Scholar 
Mailund, T. et al. A new isolation with migration model along complete genomes infers very different divergence processes among closely related great ape species. PLoS Genet. 8, e1003125 (2012).PubMed 
PubMed Central 

Google Scholar 
Igea, J., Aymerich, P., Bannikova, A. A., Gosálbez, J. & Castresana, J. Multilocus species trees and species delimitation in a temporal context: Application to the water shrews of the genus Neomys. BMC Evol. Biol. 15, 209 (2015).PubMed 
PubMed Central 

Google Scholar 
Sánchez-Gracia, A. & Castresana, J. Impact of deep coalescence on the reliability of species tree inference from different types of DNA markers in mammals. PLoS One 7, e30239 (2012).ADS 
PubMed 
PubMed Central 

Google Scholar 
Degnan, J. H. & Rosenberg, N. A. Gene tree discordance, phylogenetic inference and the multispecies coalescent. Trends Ecol. Evol. 24, 332–340 (2009).PubMed 

Google Scholar 
Edwards, S. V. & Beerli, P. Perspective: Gene divergence, population divergence, and the variance in coalescence time in phylogeographic studies. Evolution 54, 1839–1854 (2000).CAS 
PubMed 

Google Scholar 
Peterson, B. K., Weber, J. N., Kay, E. H., Fisher, H. S. & Hoekstra, H. E. Double digest RADseq: An inexpensive method for de novo SNP discovery and genotyping in model and non-model species. PLoS One 7, e37135 (2012).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Andrews, K. R., Good, J. M., Miller, M. R., Luikart, G. & Hohenlohe, P. A. Harnessing the power of RADseq for ecological and evolutionary genomics. Nat. Rev. Genet. 17, 81–92 (2016).CAS 
PubMed 
PubMed Central 

Google Scholar 
Escoda, L., Fernández-González, A. & Castresana, J. Quantitative analysis of connectivity in populations of a semi-aquatic mammal using kinship categories and network assortativity. Mol. Ecol. Resour. 19, 310–326 (2019).PubMed 

Google Scholar 
Bininda-Emonds, O. R. P. Fast genes and slow clades: Comparative rates of molecular evolution in mammals. Evol. Bioinform. 3, 59 (2007).CAS 

Google Scholar 
Welch, J. J., Bininda-Emonds, O. R. P. & Bromham, L. Correlates of substitution rate variation in mammalian protein-coding sequences. BMC Evol. Biol. 8, 53 (2008).PubMed 
PubMed Central 

Google Scholar 
Matassi, G., Sharp, P. M. & Gautier, C. Chromosomal location effects on gene sequence evolution in mammals. Curr. Biol. 9, 786–791 (1999).CAS 
PubMed 

Google Scholar 
Lercher, M. J., Chamary, J. V. & Hurst, L. D. Genomic regionality in rates of evolution is not explained by clustering of genes of comparable expression profile. Genome Res. 14, 1002–1013 (2004).CAS 
PubMed 
PubMed Central 

Google Scholar 
Castresana, J. Genes on human chromosome 19 show extreme divergence from the mouse orthologs and a high GC content. Nucleic Acids Res. 30, 1751–1756 (2002).CAS 
PubMed 
PubMed Central 

Google Scholar 
Benton, M. J., Donoghue, P. C. J. & Asher, R. J. Calibrating and constraining molecular clocks. In The Timetree of Life (eds Hedges, S. B. & Kumar, S.) 35–86 (Oxford University Press, 2009).Bouckaert, R. et al. BEAST 2.5: An advanced software platform for Bayesian evolutionary analysis. PLoS Comput. Biol. 15, e1006650 (2019).CAS 
PubMed 
PubMed Central 

Google Scholar 
Musser, G. G. & Carleton, M. D. Superfamily Muroidea. In Mammal Species of the World. A Taxonomic and Geographic Reference (eds Wilson, D. E. & Reeder, D. M.) 894–1531 (Johns Hopkins University Press, 2005).Pardiñas, U. F. J. et al. Family Cricetidae (True Hamsters, Voles, Lemmings and New World Rats and Mice). In Handbook of the Mammals of the World. Volume 7. Rodents II (eds Wilson, D. E. et al.) 204-279 (Lynx Edicions, 2017).Chevret, P. et al. Genetic structure, ecological versatility, and skull shape differentiation in Arvicola water voles (Rodentia, Cricetidae). J. Zoolog. Syst. Evol. Res. 58, 1323–1334 (2020).
Google Scholar 
Kryštufek, B. et al. Fossorial morphotype does not make a species in water voles. Mammalia 79, 293–303 (2015).
Google Scholar 
Centeno-Cuadros, A., Delibes, M. & Godoy, J. A. Dating the divergence between Southern and European water voles using molecular coalescent-based methods. J. Zool. 279, 404–409 (2009).
Google Scholar 
Castiglia, R. et al. The Italian peninsula hosts a divergent mtDNA lineage of the water vole, Arvicola amphibius s.l., including fossorial and aquatic ecotypes. Hystrix 27, 99–103 (2016).
Google Scholar 
Mahmoudi, A. et al. Evolutionary history of water voles revisited: Confronting a new phylogenetic model from molecular data with the fossil record. Mammalia 84, 171–184 (2020).
Google Scholar 
Cassola, F. Arvicola scherman, Montane Water Vole. The IUCN Red List of Threatened Species e.T136766A115519839 (2016).Somoano, A., Miñarro, M. & Ventura, J. Reproductive potential of a vole pest (Arvicola scherman) in Spanish apple orchards. Spanish J. Agric. Res. 14, e1008 (2016).
Google Scholar 
Somoano, A., Ventura, J. & Miñarro, M. Continuous breeding of fossorial water voles in northwestern Spain: Potential impact on apple orchards. Folia Zool. 66, 37–49 (2017).
Google Scholar 
Ventura, J. & Gosálbez, J. Taxonomic review of Arvicola terrestris (Linnaeus, 1758) (Rodentia, Arvicolidae) in the Iberian Peninsula. Bonn Zool. Beitr. 40, 227–242 (1989).
Google Scholar 
Ventura, J. & Sans-Fuentes, M. A. Geographic variation and divergence in nonmetric cranial traits of Arvicola (Mammalia, Rodentia) in southwestern Europe. Z. Säugetierkunde 62, 99–107 (1997).
Google Scholar 
Gómez, A. & Lunt, D. H. Refugia within refugia: Patterns of phylogeographic concordance in the Iberian Peninsula. In Phylogeography of Southern European Refugia (eds S. Weiss & N. Ferrand) 155–188 (Springer, 2007).Batsaikhan, N. et al. Arvicola amphibius, European Water Vole. The IUCN Red List of Threatened Species e.T2149A197271401 (2016).Cuenca-Bescós, G., Agustí, J., Lira, J., Rubio, M. M. & Rofes, J. A new species of water vole from the early Pleistocene of Southern Europe. Acta Palaeontol. Pol. 55, 565–580 (2010).
Google Scholar 
Cubo, J., Ventura, J. & Casinos, A. A heterochronic interpretation of the origin of digging adaptations in the northern water vole, Arvicola terrestris (Rodentia: Arvicolidae). Biol. J. Linn. Soc. 87, 381–391 (2006).
Google Scholar 
Catchen, J. M., Hohenlohe, P. A., Bassham, S., Amores, A. & Cresko, W. A. Stacks: An analysis tool set for population genomics. Mol. Ecol. 22, 3124–3140 (2013).PubMed 
PubMed Central 

Google Scholar 
Pritchard, J. K., Stephens, M. & Donnelly, P. Inference of population structure using multilocus genotype data. Genetics 155, 945–959 (2000).CAS 
PubMed 
PubMed Central 

Google Scholar 
Altschul, S. F. et al. Gapped BLAST and PSI-BLAST: A new generation of protein database search programs. Nucleic Acids Res. 25, 3389–3402 (1997).CAS 
PubMed 
PubMed Central 

Google Scholar 
Yates, A. D. et al. Ensembl 2020. Nucleic Acids Res. 48, D682–D688 (2020).CAS 

Google Scholar 
Aghová, T. et al. Fossils know it best: Using a new set of fossil calibrations to improve the temporal phylogenetic framework of murid rodents (Rodentia: Muridae). Mol. Phylogenet. Evol. 128, 98–111 (2018).PubMed 

Google Scholar 
Hey, J. et al. Phylogeny estimation by integration over isolation with migration models. Mol. Biol. Evol. 35, 2805–2818 (2018).CAS 
PubMed 
PubMed Central 

Google Scholar 
Phifer-Rixey, M., Harr, B. & Hey, J. Further resolution of the house mouse (Mus musculus) phylogeny by integration over isolation-with-migration histories. BMC Evol. Biol. 20, 120 (2020).PubMed 
PubMed Central 

Google Scholar 
Hey, J. The divergence of chimpanzee species and subspecies as revealed in multipopulation isolation-with-migration analyses. Mol. Biol. Evol. 27, 921–933 (2010).CAS 
PubMed 

Google Scholar 
Kumar, S. & Subramanian, S. Mutation rates in mammalian genomes. Proc. Natl. Acad. Sci. U.S.A. 99, 803–808 (2002).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Uchimura, A. et al. Germline mutation rates and the long-term phenotypic effects of mutation accumulation in wild-type laboratory mice and mutator mice. Genome Res. 25, 1125–1134 (2015).CAS 
PubMed 
PubMed Central 

Google Scholar 
Milholland, B. et al. Differences between germline and somatic mutation rates in humans and mice. Nat. Commun. 8, 15183 (2017).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Wright, B. R. et al. A demonstration of conservation genomics for threatened species management. Mol. Ecol. Resour. 20, 1526–1541 (2020).PubMed 

Google Scholar 
Escoda, L. & Castresana, J. The genome of the Pyrenean desman and the effects of bottlenecks and inbreeding on the genomic landscape of an endangered species. Evol. Appl. 14, 1898–1913 (2021).PubMed 
PubMed Central 

Google Scholar 
Arnold, B., Corbett-Detig, R. B., Hartl, D. & Bomblies, K. RADseq underestimates diversity and introduces genealogical biases due to nonrandom haplotype sampling. Mol. Ecol. 22, 3179–3190 (2013).CAS 
PubMed 

Google Scholar 
Cariou, M., Duret, L. & Charlat, S. How and how much does RAD-seq bias genetic diversity estimates?. BMC Evol. Biol. 16, 240 (2016).PubMed 
PubMed Central 

Google Scholar 
Campbell, C. R. et al. Pedigree-based and phylogenetic methods support surprising patterns of mutation rate and spectrum in the gray mouse lemur. Heredity 127, 233–244 (2021).CAS 
PubMed 

Google Scholar 
Scornavacca, C. et al. Orthomam v10: Scaling-up orthologous coding sequence and exon alignments with more than one hundred mammalian genomes. Mol. Biol. Evol. 36, 861–862 (2019).CAS 
PubMed 
PubMed Central 

Google Scholar 
Willis, S. C., Hollenbeck, C. M., Puritz, J. B., Gold, J. R. & Portnoy, D. S. Haplotyping RAD loci: An efficient method to filter paralogs and account for physical linkage. Mol. Ecol. Resour. 17, 955–965 (2017).CAS 
PubMed 

Google Scholar 
O’Leary, S. J., Puritz, J. B., Willis, S. C., Hollenbeck, C. M. & Portnoy, D. S. These aren’t the loci you’e looking for: Principles of effective SNP filtering for molecular ecologists. Mol. Ecol. 27, 3193–3206 (2018).
Google Scholar 
Dahl-Jensen, D. et al. Eemian interglacial reconstructed from a Greenland folded ice core. Nature 493, 489–494 (2013).ADS 
CAS 

Google Scholar 
Clark, P. U. et al. The last glacial maximum. Science 325, 710–714 (2009).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Pinho, C. & Hey, J. Divergence with gene flow: Models and data. Annu. Rev. Ecol. Evol. Syst. 41, 215–230 (2010).
Google Scholar 
Balmori-de la Puente, A. et al. Size increase without genetic divergence in the Eurasian water shrew Neomys fodiens. Sci. Rep. 9, 17375 (2019).ADS 
PubMed 
PubMed Central 

Google Scholar 
Langmead, B. & Salzberg, S. L. Fast gapped-read alignment with Bowtie 2. Nat. Methods 9, 357–359 (2012).CAS 
PubMed 
PubMed Central 

Google Scholar 
Foll, M. & Gaggiotti, O. A genome-scan method to identify selected loci appropriate for both dominant and codominant markers: A Bayesian perspective. Genetics 180, 977–993 (2008).PubMed 
PubMed Central 

Google Scholar 
Freedman, A. H. et al. Genome sequencing highlights the dynamic early history of dogs. PLoS Genet. 10, e1004016 (2014).PubMed 
PubMed Central 

Google Scholar 
Felsenstein, J. PHYLIP-phylogeny inference package (version 3.4). Cladistics 5, 164–166 (1989).
Google Scholar 
Zheng, X. et al. A high-performance computing toolset for relatedness and principal component analysis of SNP data. Bioinformatics 28, 3326–3328 (2012).CAS 
PubMed 
PubMed Central 

Google Scholar 
Evanno, G., Regnaut, S. & Goudet, J. Detecting the number of clusters of individuals using the software STRUCTURE: A simulation study. Mol. Ecol. 14, 2611–2620 (2005).CAS 
PubMed 

Google Scholar 
Jakobsson, M. & Rosenberg, N. A. CLUMPP: A cluster matching and permutation program for dealing with label switching and multimodality in analysis of population structure. Bioinformatics 23, 1801–1806 (2007).CAS 
PubMed 

Google Scholar 
Goudet, J. HIERFSTAT, a package for R to compute and test hierarchical F-statistics. Mol. Ecol. Notes 5, 184–186 (2005).
Google Scholar 
Katoh, K. & Standley, D. M. MAFFT multiple sequence alignment software version 7: Improvements in performance and usability. Mol. Biol. Evol. 30, 772–780 (2013).CAS 
PubMed 
PubMed Central 

Google Scholar 
Castresana, J. Selection of conserved blocks from multiple alignments for their use in phylogenetic analysis. Mol. Biol. Evol. 17, 540–552 (2000).CAS 
PubMed 

Google Scholar 
Stamatakis, A. RAxML version 8: A tool for phylogenetic analysis and post-analysis of large phylogenies. Bioinformatics 30, 1312–1313 (2014).CAS 
PubMed 
PubMed Central 

Google Scholar 
Brown, R. P. & Yang, Z. Rate variation and estimation of divergence times using strict and relaxed clocks. BMC Evol. Biol. 11, 271 (2011).PubMed 
PubMed Central 

Google Scholar 
Rambaut, A., Drummond, A. J., Xie, D., Baele, G. & Suchard, M. A. Posterior Summarization in Bayesian Phylogenetics Using Tracer 1.7. Syst. Biol. 67, 901–904 (2018).CAS 
PubMed 
PubMed Central 

Google Scholar 
Hey, J. & Wang, K. The effect of undetected recombination on genealogy sampling and inference under an isolation-with-migration model. Mol. Ecol. Resour. 18, 489 (2019).
Google Scholar 
QGIS_Development_Team. QGIS Geographic Information System. Open Source Geospatial Foundation Project. http://qgis.osgeo.org (2021).IUCN. Arvicola scherman. The IUCN Red List of Threatened Species. Version 6.2. https://www.iucnredlist.org. Downloaded on 04 September 2019. (2019).IUCN. Arvicola amphibius. The IUCN Red List of Threatened Species. Version 6.2. https://www.iucnredlist.org. Downloaded on 10 July 2019. (2019). More

Following the work in36, we construct an epidemiological model which tracks the disease dynamics and population of two species of hosts following the introduction of a pathogen. The native host (hereafter simply referred to as “type 1”) is vulnerable to the disease, but due to being well adapted to the native habitat has high fecundity when uninfected. The invasive host (hereafter referred to as “type 2”), has coevolved defenses to the pathogen that increase both its tolerance of and resistance to the disease, but is not inherently as well-adapted to the habitat in the absence of infection (i.e., its intrinsic rate of growth in the new habitat is lower than that of the native).Our initial conditions correspond to a population of uninfected type 1 hosts with a small number of both uninfected and infected type 2 hosts, representing an invasion by a novel competitor carrying a novel pathogen into the type 1 population. We consider a vector-borne pathogen, and make the simplifying assumption that there is an already abundant competent vector species in the habitat. (For this initial formulation, we considered a scenario of mosquito-borne infections in birds, such as avian malaria37 or West Nile virus38, to motivate concrete choices.)The model couples two biological dynamics: the daily vector-borne spread of the disease among hosts, and a yearly host breeding cycle. We simulate in discrete time-steps that represent days using an SIR model taking into account the interactions between the disease, the two species of host, and the vectors. The model also includes a passive death rate for hosts of vectors, which increases for hosts while infected. While the vectors are assumed to breed daily, the hosts reproduce as part of an assumed annual breeding season, every (t_c) time-steps (typically equal to 365). These dynamics were informed by considering an annually breeding bird population in a tropical environment, however, they are not meant to reflect the realism of any one biological system. They are chosen here merely to allow a clean interpretation of modeled scenarios. Future models should explore the impact of greater variety in the dynamics of possible vector and host reproductive patterns.Epidemiological modelThe model tracks eight variables corresponding to combinations of host species and vectors with their infection status. Hosts may be of type 1 or 2, and are either susceptible to the disease ((S_1, S_2)), currently infected ((I_1, I_2)), or recovered ((R_1, R_2)). We assume that recovery is complete and recovered individuals suffer no residual effects from their infection aside from a lifelong immunity to becoming reinfected. (We later set the recovery rate for host type 1 to 0, so (R_1 = 0) at all times, but leave it defined for the sake of generality.) For simplicity, we model using only one stage of infection in which individuals are both infectious and symptomatic. The model also tracks the status of the vector population, which may either be susceptible ((S_v)) or infected ((I_v)). We assume that vectors do not recover from the disease, but also suffer no negative effects from being infected, acting only as carriers.For convenience of notation, we denote the total number of hosts$$begin{aligned} H = S_1 + I_1 + R_1 + S_2 + I_2 + R_2 end{aligned}$$and the relative frequencies of infection within their respective population$$begin{aligned} F_1 = frac{I_1}{H}, F_2 = frac{I_2}{H},F_v = frac{I_v}{S_v+I_v} end{aligned}$$which allows some equations to be written more compactly. Table 1 shows a summary of these variables.Table 1 Variables.Full size tableThe model also has several constant parameters that affect the dynamics. (beta _j) determines the probability that hosts of type j become infected when bitten by a single infected vector. We typically set (beta _1 > beta _2), making type 2 hosts less likely to become infected.Likewise, (delta _j) determines the probability that a vector becomes infected when biting an infected host of type j.(b_j) determines the bite rate for vectors on host type j. We assume that each vector bites the same number of hosts per day, so each vector’s probability of becoming infected depends only on the frequency of infection among hosts, while each host will be bitten more if there are more vectors.(gamma _j) determines the proportion of infected hosts of type j that recover from the disease each day. We typically set (gamma _1 = 0 < gamma _2), meaning infected hosts of type 1 do not recover, while infected type 2 recover after an average of (1/gamma _2) days.(mu _{j-}) determines the daily death rate for uninfected hosts of type j and (mu _{j+}) determines the death rate for infected host of type j. We typically set (mu _{1-} = mu _{2-}< mu _{2+} < mu _{1+}), meaning uninfected hosts have the same death rate regardless of type, infected type 2 have a higher death rate than uninfected hosts, and infected type 1 have the highest. (Both susceptible and recovered hosts are considered to be uninfected.) Table 2 shows a summary of parameters related to the SIR dynamics.Equation 1 shows continuous ordinary differential equations approximating the dynamics. Note that the actual model instantiates these in discrete time-steps using the forward Euler method with (h = 1).$$ begin{aligned}&frac{dS_1}{dt} = - S_1 beta _1 b_1 I_v /H - S_1 mu _{1-} \&frac{dI_1}{dt} = S_1 beta _1 b_1 I_v /H - gamma _1 I_1 - I_1 mu _{1+} \&frac{dR_1}{dt} = I_1 gamma _1 - R_1 mu _{1-} \&frac{dS_2}{dt} = -S_2 beta _2 b_2 I_v /H - S_2 mu _{2-} \&frac{dI_2}{dt} = S_2 beta _2 b_2 I_v /H - I_2 gamma _2 - I_2 mu _{2+} \&frac{dR_2}{dt} = I_2 gamma _2 - R_2 mu _{2-}\&frac{dS_v}{dt} = alpha _v H -S_v delta _1 b_1 F_1 -S_v delta _2 b_2 F_2 -S_v mu _v\&frac{dI_v}{dt} = S_v delta _1 b_1 F_1 + S_v delta _2 b_2 F_2 - I_v mu _v\ end{aligned} $$ (1) Table 2 Parameters for SIR dynamics.Full size tableFollowing a standard SIR model, susceptible hosts can become infected, and infected hosts become recovered, but each equation also contains a negative term corresponding to deaths. Thus, the total population of hosts is strictly decreasing in this time-frame. We assume that the vectors breed on a much shorter timescale than hosts, so we include a term for their births here, while host births are implemented by a yearly breeding event. We assume no vertical disease transmission, so all new vectors begin in the susceptible category. We assume that the daily birthrate for each vector increases with access to hosts, and decreases with competition among other vectors for hosts and breeding sites, so we set it equal to (frac{alpha _v H}{S_v + I_v}), where (alpha _v) is a constant scaling factor. Since the birthrate for each vector contains the total number of vectors in its denominator, the total number of vector births in the population will simply be (alpha _v H).A population with a larger number of hosts will be able to sustain a larger number of vectors. For a population with a constant number of hosts, the equilibrium vector population will be proportional to the number hosts: aH where (a = frac{alpha _v}{mu _v}) is the equilibrium vector density (number of vectors per host). For instance if (a = 2), then in equilibrium there will be twice as many vectors as hosts. Given a fixed number of hosts, the population of vectors will asymptotically approach the equilibrium value. In practice the total number of hosts is constantly changing, so the population of vectors will chase after this moving equilibrium, though for our standard parameters (alpha _v) and (mu _v) are sufficiently large such that this will occur on a short timescale, and the population of vectors remains close to the current equilibrium value.Breeding eventTable 3 shows a summary of parameters related to the breeding event. Every (t_c) days (typically 365), a breeding event occurs according to the following process.Table 3 Parameters for breeding event.Full size tableLet$$begin{aligned}&Delta S_1 = t_c alpha _{1-}(S_1+R_1)+t_calpha _{1+} I_1 \&Delta S_2 = t_c alpha _{2-}(S_2+R_2)+t_calpha _{2+} I_2 \ end{aligned}$$be the number of new host offspring of each type born this generation. In order to maintain consistency of temporal units among the parameters, each birthrate parameter is multiplied by (t_c). Let H be the current total number of hosts. Let$$begin{aligned} c = {left{ begin{array}{ll} 0 &{} hbox {if } H ge kappa \ 1 &{} hbox {if } H + Delta S_1 + Delta S_2 le kappa \ frac{kappa -H}{Delta S_1 + Delta S_2} &{} hbox {otherwise} \ end{array}right. } end{aligned}$$be the proportion of offspring that survive to adulthood. (None, if the population is already above carrying capacity. All, if the difference between the reproducing population size and the carrying capacity exceeds the new births. If the population is approaching carrying capacity, juvenile mortality scales proportionally so that the population will hit carrying capacity but not exceed it.)Then$$begin{aligned}&S_1 + c Delta S_1 rightarrow S_1 \&S_2 + c Delta S_2 rightarrow S_2 \ end{aligned}$$We assume there is no vertical disease transmission, so all new hosts begin in the susceptible category. We assume that the host population is iteroparous, such that the new offspring and the existing adult population both carry over to the next generation. If the new population would exceed the carrying capacity, we assume the limited space or supplies reduces the number of successful offspring so that the population exactly reaches the carry capacity by reduction in juvenile survival rather than population-wide competition that could also reduce the adult population.The carrying capacity is therefore what drives the interspecific host competition. Because births of both species are summed and then normalized by the total number of births, the higher the birthrate of one host, the larger a fraction of the available space it will capture during the breeding event. Similarly, the lower the death-rate of a host, the less space it frees up for the next breeding event. Even if one host species would be able to sustain a stable population on its own, the presence of a more fit competitor can lead to the extinction of the less fit type by driving its effective birth rate down.Immune-reproductive trade-offs and boundary conditionsWe assume that host type 1 is evolutionarily stable in the absence of the disease; an uninfected monoculture population below the carrying capacity will have at least as many births as deaths each cycle. In a continuous version of this model where births and deaths happened simultaneously, this might be defined by (alpha _{1-} ge mu _{1-}) . However in our model, the population spends many days decreasing due to deaths before the next breeding event occurs. The population exponentially decays throughout the cycle, and then jumps up during the breeding event. The number of new host births is proportional to the number of hosts at the start of the breeding event, which will be the lowest value of any other time during the cycle. Thus, the birth rate needs to be high enough that the surviving hosts can compensate despite their diminished numbers. Taking this into account, we get the condition$$begin{aligned}&alpha _{1-} ge frac{1-(1- mu _{1-})^{t_c}}{(1-mu _{1-})^{t_c}} \ end{aligned}$$Which is a higher bound on (alpha _{1-}) than the simpler one above, but will be close to it if (mu _{1-}) and (t_c) are small.To implement the scenario in which type 2 has increased resistance and tolerance to the disease at the expense of overall fecundity, we implement the following boundary conditions:$$begin{aligned}&beta _1 > beta _2 \&0 = gamma _1< gamma _2 \&mu _{1-} = mu _{2-}< mu _{2+} < mu _{1+} \&alpha _{1-} > alpha _{2-} > alpha _{2+} > alpha _{1+} end{aligned}$$Type 2 hosts are less likely to contract the disease, and are able to recover from it, while type 1 lack the immunological strength to eradicate it completely. Additionally, while both types of host are weakened by the disease, type 2 suffer fewer negative effects. However, this stronger immune response comes at the cost of reducing their birth rate when compared to healthy type 1 hosts.Due to the heterogeneous population, there is ambiguity in defining (R_0) for the disease. The two types of host have different transmission rates and durations of infection, and will therefore be responsible for different amounts of disease spread. To resolve this, we define several related values. Let (R_0^j) be the (R_0) of the disease in a homogeneous population of type j hosts: the average number of hosts infected (indirectly, through vectors) from a single infected host in a population consisting entirely of type j hosts.$$begin{aligned}&R_0^1 = frac{delta _1 beta _1 a b_1^2}{mu _v mu _{1+}} \&R_0^2 = frac{delta _2 beta _2 a b_2^2}{mu _v (mu _{2+}+gamma _2)} end{aligned}$$We simplify the equation for (R_0^1) since (gamma _1 = 0). We define w to be the frequency of host type 1: (w := (S_1 + I_1)/H). Then (R_0) for the vectors is$$begin{aligned} R_0^v = R_0^1 w + R_0^2 (1-w) end{aligned}$$which will also be the effective (R_0) of the disease for the hosts in the mixed population.For simplicity of results, we restrict to the case where type 1 is more infectious overall than type 2, in particular (R_0^1 > R_0^2). This allows us to avoid edge cases in simulation outcomes which are beyond the scope of this paper. We intend to lift this restriction and study these outcomes in future work.NoteAlthough usual epidemiological model formulations can rely on the value 1 as the boundary condition for (R_0) to determine the epidemic potential of an outbreak, in this case we are calculating effective (R_0) in a dynamic host population, such that the decrease in disease spread due to saturation from recovered hosts and already infected hosts increases the actual thresholds. More accurate criteria require a technical and somewhat cumbersome analysis, which we leave for a future paper. More

Tusting, L. S. et al. Mapping changes in housing in sub-Saharan Africa from 2000 to 2015. Nature 568, 391–394 (2019).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Lozano, R. et al. Measuring progress from 1990 to 2017 and projecting attainment to 2030 of the health-related Sustainable Development Goals for 195 countries and territories: A systematic analysis for the Global Burden of Disease Study 2017. Lancet 392, 2091–2138 (2018).
Google Scholar 
Kassebaum, N. J. et al. Global, regional, and national levels and causes of maternal mortality during 1990–2013: A systematic analysis for the Global Burden of Disease Study 2013. Lancet 384, 980–1004 (2014).PubMed 
PubMed Central 

Google Scholar 
Dhiman, S. Are malaria elimination efforts on right track? An analysis of gains achieved and challenges ahead. Infect. Dis. Poverty 8, 1–19 (2019).
Google Scholar 
WHO. World Malaria Report 2018 (WHO, 2018).
Google Scholar 
Janko, M. M. et al. The links between agriculture, Anopheles mosquitoes, and malaria risk in children younger than 5 years in the Democratic Republic of the Congo: A population-based, cross-sectional, spatial study. Lancet Planet. Health 2, e74–e82 (2018).PubMed 
PubMed Central 

Google Scholar 
Gerland, P. et al. World population stabilization unlikely this century. Science 346, 234–237 (2014).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Jayne, T. S., Chamberlin, J. & Headey, D. D. Land pressures, the evolution of farming systems, and development strategies in Africa: A synthesis. Food Policy 48, 1–17 (2014).
Google Scholar 
Ordway, E. M., Asner, G. P. & Lambin, E. F. Deforestation risk due to commodity crop expansion in sub-Saharan Africa. Environ. Res. Lett. 12, 044015 (2017).ADS 

Google Scholar 
Chaves, L. S. M. et al. Global consumption and international trade in deforestation-associated commodities could influence malaria risk. Nat. Commun. 11, 1–10 (2020).
Google Scholar 
Adenle, A. A., Azadi, H. & Manning, L. The era of sustainable agricultural development in Africa: Understanding the benefits and constraints. Food Rev. Int. 34, 411–433 (2018).
Google Scholar 
Ijumba, J. N. & Lindsay, S. W. Impact of irrigation on malaria in Africa: Paddies paradox. Med. Vet. Entomol. 15, 1–11 (2001).CAS 
PubMed 

Google Scholar 
Warra, A. A. & Prasad, M. N. V. African perspective of chemical usage in agriculture and horticulture—their impact on human health and environment. In Agrochemicals, Detection Treatment and Remediation 401–436 (Elsevier, 2020).
Google Scholar 
Fornace, K. M., Diaz, A. V., Lines, J. & Drakeley, C. J. Achieving global malaria eradication in changing landscapes. Malar. J. 20, 1–14 (2021).
Google Scholar 
Burkett-Cadena, N. D. & Vittor, A. Y. Deforestation and vector-borne disease: Forest conversion favors important mosquito vectors of human pathogens. Basic Appl. Ecol. 26, 101–110 (2018).PubMed 

Google Scholar 
Lindblade, K. A., Walker, E. D., Onapa, A. W., Katungu, J. & Wilson, M. L. Land use change alters malaria transmission parameters by modifying temperature in a highland area of Uganda. Trop. Med. Int. Heal. 5, 263–274 (2000).CAS 

Google Scholar 
Yasuoka, J. & Levins, R. Impact of deforestation and agricultural development on anopheline ecology and malaria epidemiology. Am. J. Trop. Med. Hyg. 76, 450–460 (2007).PubMed 

Google Scholar 
Guerra, C. A., Snow, R. W. & Hay, S. I. A global assessment of closed forests, deforestation and malaria risk. Ann. Trop. Med. Parasitol. 100, 189–204 (2006).CAS 
PubMed 

Google Scholar 
Laporta, G. Z., de Prado, P. I. K. L., Kraenkel, R. A., Coutinho, R. M. & Sallum, M. A. M. Biodiversity can help prevent malaria outbreaks in tropical forests. PLoS Negl. Trop. Dis. 7, e2139 (2013).PubMed 
PubMed Central 

Google Scholar 
Baeza, A., Santos-Vega, M., Dobson, A. P. & Pascual, M. The rise and fall of malaria under land-use change in frontier regions. Nat. Ecol. Evol. 1, 0108 (2017).
Google Scholar 
Patz, J. A., Graczyk, T. K., Geller, N. & Vittor, A. Y. Effects of environmental change on emerging parasitic diseases. Int. J. Parasitol. https://doi.org/10.1016/S0020-7519(00)00141-7 (2000).Article 
PubMed 

Google Scholar 
Sogoba, N. et al. Spatial analysis of malaria transmission parameters in the rice cultivation area of Office du Niger, Mali. Am. J. Trop. Med. Hyg. 76, 1009–1015 (2007).PubMed 

Google Scholar 
Mwangangi, J. M. et al. Anopheles larval abundance and diversity in three rice agro-village complexes Mwea irrigation scheme, central Kenya. Malar. J. 9, 1–10 (2010).
Google Scholar 
Diuk-Wasser, M. A. et al. Patterns of irrigated rice growth and malaria vector breeding in Mali using multi-temporal ERS-2 synthetic aperture radar. Int. J. Remote Sens. 27, 535–548 (2006).CAS 
PubMed 
PubMed Central 

Google Scholar 
Briët, O. J. T., Dossou-Yovo, J., Akodo, E., Van De Giesen, N. & Teuscher, T. M. The relationship between Anopheles gambiae density and rice cultivation in the savannah zone and forest zone of Côte d’Ivoire. Trop. Med. Int. Heal. 8, 439–448 (2003).
Google Scholar 
Klinkenberg, E., McCall, P. J., Wilson, M. D., Amerasinghe, F. P. & Donnelly, M. J. Impact of urban agriculture on malaria vectors in Accra, Ghana. Malar. J. 7, 1–9 (2008).
Google Scholar 
Keiser, J. et al. Effect of irrigation and large dams on the burden of malaria on a global and regional scale. Am. J. Trop. Med. Hyg. 72, 392–406 (2005).PubMed 

Google Scholar 
Kyei-Baafour, E. et al. Impact of an irrigation dam on the transmission and diversity of Plasmodium falciparum in a seasonal malaria transmission area of Northern Ghana. J. Trop. Med. 2020, 1–8 (2020).
Google Scholar 
Kibret, S. Time to revisit how dams are affecting malaria transmission. Lancet Planet. Heal. 2, e378–e379 (2018).
Google Scholar 
Kibret, S., Lautze, J., McCartney, M., Nhamo, L. & Yan, G. Malaria around large dams in Africa: Effect of environmental and transmission endemicity factors. Malar. J. 18, 1–12 (2019).
Google Scholar 
Kibret, S., Wilson, G. G., Ryder, D., Tekie, H. & Petros, B. Malaria impact of large dams at different eco-epidemiological settings in Ethiopia. Trop. Med. Health 45, 1–14 (2017).
Google Scholar 
Keiser, J., Singer, B. H. & Utzinger, J. Reducing the burden of malaria in different eco-epidemiological settings with environmental management: A systematic review. Lancet Infect. Dis. 5, 695–708 (2005).PubMed 

Google Scholar 
Ijumba, J. N., Shenton, F. C., Clarke, S. E., Mosha, F. W. & Lindsay, S. W. Irrigated crop production is associated with less malaria than traditional agricultural practices in Tanzania. Trans. R. Soc. Trop. Med. Hyg. 96, 476–480 (2002).CAS 
PubMed 

Google Scholar 
Ijumba, J. N., Mosha, F. W. & Lindsay, S. W. Malaria transmission risk variations derived from different agricultural practices in an irrigated area of northern Tanzania. Med. Vet. Entomol. 16, 28–38 (2002).CAS 
PubMed 

Google Scholar 
Klinkenberg, E. et al. Malaria and irrigated crops, Accra, Ghana. Emerg. Infect. Dis. 11, 1290–1293 (2005).PubMed 
PubMed Central 

Google Scholar 
USAID. The DHS Program—DHS Methodology. https://www.dhsprogram.com/What-We-Do/Survey-Types/DHS-Methodology.cfm (1984).Siraj, A. S. et al. Altitudinal changes in malaria incidence in highlands of Ethiopia and Colombia. Science 343, 1154–1158 (2014).ADS 
CAS 
PubMed 

Google Scholar 
Pascual, M., Ahumada, J. A., Chaves, L. F., Rodó, X. & Bouma, M. Malaria resurgence in the East African highlands: Temperature trends revisited. Proc. Natl. Acad. Sci. U. S. A. 103, 5829–5834 (2006).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Lindblade, K. A. et al. Sustainability of reductions in malaria transmission and infant mortality in Western Kenya with use of insecticide-treated bednets 4 to 6 years of follow-up. J. Am. Med. Assoc. 291, 2571–2580 (2004).CAS 

Google Scholar 
Mordecai, E. A. et al. Optimal temperature for malaria transmission is dramatically lower than previously predicted. Ecol. Lett. 16, 22–30 (2013).PubMed 

Google Scholar 
Weiss, D. J. et al. Re-examining environmental correlates of Plasmodium falciparum Malaria endemicity: A data-intensive variable selection approach. Malar. J. 14, 1–18 (2015).
Google Scholar 
Bauhoff, S. & Busch, J. Does deforestation increase malaria prevalence? Evidence from satellite data and health surveys. World Dev. 127, 104734 (2020).
Google Scholar 
Austin, K. F., Bellinger, M. O. & Rana, P. Anthropogenic forest loss and malaria prevalence: A comparative examination of the causes and disease consequences of deforestation in developing nations. AIMS Environ. Sci. 4, 217–231 (2017).
Google Scholar 
Kabaria, C. W., Gilbert, M., Noor, A. M., Snow, R. W. & Linard, C. The impact of urbanization and population density on childhood Plasmodium falciparum parasite prevalence rates in Africa. Malar. J. 16, 1–10 (2017).
Google Scholar 
Herrera, D. et al. Upstream watershed condition predicts rural children’s health across 35 developing countries. Nat. Commun. 8, 1–8 (2017).CAS 

Google Scholar 
Van Ittersum, M. K. et al. Can sub-Saharan Africa feed itself?. Proc. Natl. Acad. Sci. U. S. A. 113, 14964–14969 (2016).PubMed 
PubMed Central 

Google Scholar 
Ickowitz, A. Shifting cultivation and deforestation in tropical Africa: Critical reflections. Dev. Change 37, 599–626 (2006).
Google Scholar 
Kar, N. P., Kumar, A., Singh, O. P., Carlton, J. M. & Nanda, N. A review of malaria transmission dynamics in forest ecosystems. Parasit. Vectors 7, 1–12 (2014).
Google Scholar 
Phalan, B. et al. Crop expansion and conservation priorities in tropical countries. PLoS ONE 8, e51759 (2013).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Afrane, Y. A. et al. Does irrigated urban agriculture influence the transmission of malaria in the city of Kumasi, Ghana?. Acta Trop. 89, 125–134 (2004).PubMed 

Google Scholar 
De Silva, P. M. & Marshall, J. M. Factors contributing to urban malaria transmission in sub-saharan Africa: A systematic review. J. Trop. Med. 2012, 1–10 (2012).
Google Scholar 
Kibret, S., Wilson, G. G., Tekie, H. & Petros, B. Increased malaria transmission around irrigation schemes in Ethiopia and the potential of canal water management for malaria vector control. Malar. J. 13, 1–12 (2014).
Google Scholar 
Dongus, S. et al. Urban agriculture and Anopheles habitats in Dar es Salaam, Tanzania. Geospat. Health 3, 189–210 (2009).PubMed 

Google Scholar 
Bhatt, S. et al. The effect of malaria control on Plasmodium falciparum in Africa between 2000 and 2015. Nature 526, 207–211 (2015).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Georganos, S. et al. Modelling the wealth index of demographic and health surveys within cities using very high-resolution remotely sensed information. Remote Sens. 11, 2543 (2019).ADS 

Google Scholar 
Pascual, M. & Baeza, A. What happens when forests fall?. Elife 10, e67863 (2021).PubMed 
PubMed Central 

Google Scholar 
MacDonald, A. J. & Mordecai, E. A. Amazon deforestation drives malaria transmission, and malaria burden reduces forest clearing. Proc. Natl. Acad. Sci. U. S. A. 116, 22212–22218 (2019).CAS 
PubMed 
PubMed Central 

Google Scholar 
Valle, D. & Clark, J. Conservation efforts may increase malaria burden in the Brazilian Amazon. PLoS ONE 8, e57519 (2013).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Luis, A. D., Kuenzi, A. J. & Mills, J. N. Species diversity concurrently dilutes and amplifies transmission in a zoonotic host–pathogen system through competing mechanisms. Proc. Natl. Acad. Sci. U. S. A. 115, 7979–7984 (2018).CAS 
PubMed 
PubMed Central 

Google Scholar 
Miller, E. & Huppert, A. The effects of host diversity on vector-borne disease: The conditions under which diversity will amplify or dilute the disease risk. PLoS ONE https://doi.org/10.1371/journal.pone.0080279 (2013).Article 
PubMed 
PubMed Central 

Google Scholar 
Yamana, T. K. & Eltahir, E. A. B. Incorporating the Effects of Humidity in a Mechanistic Model of Anopheles Gambiae Mosquito Population Dynamics in the Sahel Region of Africa. http://www.parasitesandvectors.com/content/6/1/235. https://doi.org/10.1186/1756-3305-6-235 (2013).Wielgosz, B., Kato, E. & Ringler, C. Agro-ecology, household economics and malaria in Uganda: Empirical correlations between agricultural and health outcomes. Malar. J. 13, 1–11 (2014).
Google Scholar 
Asale, A., Duchateau, L., Devleesschauwer, B., Huisman, G. & Yewhalaw, D. Zooprophylaxis as a control strategy for malaria caused by the vector Anopheles arabiensis (Diptera: Culicidae): A systematic review. Infect. Dis. Poverty 6, 1–14 (2017).
Google Scholar 
Halliday, F., Rohr, J. & Laine, A.-L. Biodiversity loss underlies the dilution effect of biodiversity. Ecol. Lett. https://doi.org/10.1101/2020.04.20.050377 (2020).Article 
PubMed 
PubMed Central 

Google Scholar 
Pienkowski, T., Dickens, B. L., Sun, H. & Carrasco, L. R. Empirical evidence of the public health benefits of tropical forest conservation in Cambodia: A generalised linear mixed-effects model analysis. Lancet Planet. Health 1, e180–e187 (2017).PubMed 

Google Scholar 
Shah, H. A., Huxley, P., Elmes, J. & Murray, K. A. Agricultural land-uses consistently exacerbate infectious disease risks in Southeast Asia. Nat. Commun. 10, 4299 (2019).ADS 
PubMed 
PubMed Central 

Google Scholar 
Parham, P. E. & Michael, E. Modeling the effects of weather and climate change on malaria transmission. Environ. Health Perspect. 118, 620–626 (2010).PubMed 

Google Scholar 
Krefis, A. C. et al. Modeling the relationship between precipitation and malaria incidence in children from a holoendemic area in Ghana. Am. J. Trop. Med. Hyg. 84, 285–291 (2011).PubMed 
PubMed Central 

Google Scholar 
Drakeley, C. J. et al. Altitude-Dependent and -Independent Variations in Plasmodium falciparum Prevalence in Northeastern Tanzania. J. Infect. Dis. 191, 1589–1598 (2005).PubMed 

Google Scholar 
Masuda, K. Length of maternal schooling and children’s risk of malaria infection: Evidence from a natural experiment in Uganda. BMJ Glob. Health 5, 4–11 (2020).
Google Scholar 
Ma, C. et al. Is maternal education a social vaccine for childhood malaria infection? A cross-sectional study from war-torn Democratic Republic of Congo. Pathog. Glob. Health 111, 98–106 (2017).PubMed 
PubMed Central 

Google Scholar 
Njau, J. D., Stephenson, R., Menon, M. P., Kachur, S. P. & McFarland, D. A. Investigating the important correlates of maternal education and childhood malaria infections. Am. J. Trop. Med. Hyg. 91, 509–519 (2014).PubMed 
PubMed Central 

Google Scholar 
Degarege, A., Fennie, K., Degarege, D., Chennupati, S. & Madhivanan, P. Improving socioeconomic status may reduce the burden of malaria in sub Saharan Africa: A systematic review and meta-analysis. PLoS ONE 14, 1–26 (2019).
Google Scholar 
Sonko, S. T. et al. Does socio-economic status explain the differentials in malaria parasite prevalence? Evidence from the Gambia. Malar. J. 13, 1–12 (2014).
Google Scholar 
Tusting, L. S. et al. Housing improvements and malaria risk in Sub-Saharan Africa: A multi-country analysis of survey data. PLoS Med. 14, 1–15 (2017).
Google Scholar 
Yang, D. et al. Drinking water and sanitation conditions are associated with the risk of malaria among children under five years old in sub-Saharan Africa: A logistic regression model analysis of national survey data. J. Adv. Res. 21, 1–13 (2020).PubMed 

Google Scholar 
Hay, S. I., Guerra, C. A., Tatem, A. J., Atkinson, P. M. & Snow, R. W. Urbanization, malaria transmission and disease burden in Africa. Nat. Rev. Microbiol. 3, 81–90 (2011).
Google Scholar 
Murray, C. J. L. et al. Global malaria mortality between 1980 and 2010: A systematic analysis. Lancet 379, 413–431 (2012).PubMed 

Google Scholar 
Nankabirwa, J. et al. Malaria in school-age children in Africa: An increasingly important challenge. Trop. Med. Int. Heal. 19, 1294–1309 (2014).
Google Scholar 
Okiro, E. A. et al. Age patterns of severe paediatric malaria and their relationship to Plasmodium falciparum transmission intensity. Malar. J. 8, 1–11 (2009).
Google Scholar 
Fullman, N., Burstein, R., Lim, S. S., Medlin, C. & Gakidou, E. Nets, spray or both? the effectiveness of insecticide-treated nets and indoor residual spraying in reducing malaria morbidity and child mortality in sub-Saharan Africa. Malar. J. 12, 1 (2013).
Google Scholar 
Agusto, F. B. et al. The impact of bed-net use on malaria prevalence. J. Theor. Biol. 320, 58–65 (2013).ADS 
PubMed 
MATH 

Google Scholar 
Hughes, R. A., Heron, J., Sterne, J. A. C. & Tilling, K. Accounting for missing data in statistical analyses: Multiple imputation is not always the answer. Int. J. Epidemiol. 48, 1294–1304 (2019).PubMed 
PubMed Central 

Google Scholar 
Beck-Johnson, L. M. et al. The effect of temperature on Anopheles mosquito population dynamics and the potential for malaria transmission. PLoS ONE 8, e79276 (2013).ADS 
PubMed 
PubMed Central 

Google Scholar 
Hien, D. F. D. S. et al. Plant-mediated effects on mosquito capacity to transmit human malaria. PLoS Pathog. 12, 1–17 (2016).
Google Scholar 
Donnelly, B., Berrang-Ford, L., Ross, N. A. & Michel, P. A systematic, realist review of zooprophylaxis for malaria control. Malar. J. 14, 1–16 (2015).
Google Scholar 
Curtis, P. G., Slay, C. M., Harris, N. L., Tyukavina, A. & Hansen, M. C. Classifying drivers of global forest loss. Science 361, 1108–1111 (2018).ADS 
CAS 
PubMed 

Google Scholar 
Townes, L. R., Mwandama, D., Mathanga, D. P. & Wilson, M. L. Elevated dry-season malaria prevalence associated with fine-scale spatial patterns of environmental risk: A case-control study of children in rural Malawi. Malar. J. 12, 1 (2013).
Google Scholar 
Brock, P. M. et al. Predictive analysis across spatial scales links zoonotic malaria to deforestation. Proc. R. Soc. B 286, 20182351 (2019).PubMed 
PubMed Central 

Google Scholar 
Zhang, Z., Manjourides, J., Cohen, T., Hu, Y. & Jiang, Q. Spatial measurement errors in the field of spatial epidemiology. Int. J. Health Geogr. 15, 1–12 (2016).
Google Scholar 
Rockström, J. et al. Managing water in rainfed agriculture: The need for a paradigm shift. Agric. Water Manag. 97, 543–550 (2010).
Google Scholar 
Rockström, J., Barron, J. & Fox, P. Water productivity in rain-fed agriculture: Challenges and opportunities for smallholder farmers in drought-prone tropical agroecosystems. Water Product. Agric. limits Oppor. Improv. 85199, 145–162. https://doi.org/10.1079/9780851996691.0145 (2009).Article 

Google Scholar 
Steinmann, P., Keiser, J., Bos, R., Tanner, M. & Utzinger, J. Schistosomiasis and water resources development: Systematic review, meta-analysis, and estimates of people at risk. Lancet. Infect. Dis. 6, 411–425 (2006).PubMed 

Google Scholar 
Halstead, N. T. et al. Agrochemicals increase risk of human schistosomiasis by supporting higher densities of intermediate hosts. Nat. Commun. 9, 837 (2018).ADS 
PubMed 
PubMed Central 

Google Scholar 
Rohr, J. R. et al. Emerging human infectious diseases and the links to global food production. Nat. Sustain. 2, 445–456 (2019).PubMed 
PubMed Central 

Google Scholar 
Sokolow, S. H. et al. Reduced transmission of human schistosomiasis after restoration of a native river prawn that preys on the snail intermediate host. Proc. Natl. Acad. Sci. U. S. A. 112, 9650–9655 (2015).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Rasolofoson, R. A., Hanauer, M. M., Pappinen, A., Fisher, B. & Ricketts, T. H. Impacts of forests on children’s diet in rural areas across 27 developing countries. Sci. Adv. 4, 1–10 (2018).
Google Scholar 
Doxsey-Whitfield, E. et al. Taking advantage of the improved availability of census data: A first look at the gridded population of the world, version 4. Pap. Appl. Geogr. 1, 226–234 (2015).
Google Scholar 
Hansen, M. C. et al. High-resolution global maps of 21st-century forest cover change. Science 850, 850–854 (2013).ADS 

Google Scholar 
Hijmans, R. J., Cameron, S. E., Parra, J. L., Jones, P. G. & Jarvis, A. Very high resolution interpolated climate surfaces for global land areas. Int. J. Climatol. 25, 1965–1978 (2005).
Google Scholar 
Hollister, M. J. Package ‘elevatr’ Title Access Elevation Data from Various APIs. (2018).Bontemps, S. et al. Consistent global land cover maps for climate modelling communities: Current achievements of the ESA’s land cover CCI. Proc. ESA Living Planet Symp. 13, 9–13 (2013).
Google Scholar 
Mahende, C. et al. Performance of rapid diagnostic test, blood-film microscopy and PCR for the diagnosis of malaria infection among febrile children from Korogwe District, Tanzania. Malar. J. 15, 1–7 (2016).
Google Scholar 
Stauffer, W. M. et al. Diagnostic performance of rapid diagnostic tests versus blood smears for malaria in US clinical practice. Clin. Infect. Dis. 49, 908–913 (2009).PubMed 

Google Scholar 
Yankson, R., Anto, E. A. & Chipeta, M. G. Geostatistical analysis and mapping of malaria risk in children under 5 using point-referenced prevalence data in Ghana. Malar. J. 18, 1–12 (2019).
Google Scholar 
Gatton, M. L. et al. Impact of Plasmodium falciparum gene deletions on malaria rapid diagnostic test performance. Malar. J. 19, 1–11 (2020).
Google Scholar 
Austin, K. F. Export agriculture is feeding malaria: A cross-national examination of the environmental and social causes of malaria prevalence. Popul. Environ. 35, 133–158 (2013).
Google Scholar 
Tyukavina, A. et al. Aboveground carbon loss in natural and managed tropical forests from 2000 to 2012. Environ. Res. Lett. 10, 074002 (2015).ADS 

Google Scholar 
Ayele, D. G., Zewotir, T. T. & Mwambi, H. G. Prevalence and risk factors of malaria in Ethiopia. Malar. J. 11, 1 (2012).
Google Scholar 
Acheson, E. S. & Kerr, J. T. Nets versus spraying: A spatial modelling approach reveals indoor residual spraying targets Anopheles mosquito habitats better than mosquito nets in Tanzania. PLoS ONE 13, 1–19 (2018).
Google Scholar 
Siraj, A. S. et al. Temperature and population density determine reservoir regions of seasonal persistence in highland malaria. Proc. R. Soc. B 282, 20151383 (2015).PubMed 
PubMed Central 

Google Scholar 
Ishengoma, D. S. et al. Trends of Plasmodium falciparum prevalence in two communities of Muheza district North-eastern Tanzania: Correlation between parasite prevalence, malaria interventions and rainfall in the context of re-emergence of malaria after two decades of progressive. Malar. J. 17, 1–10 (2018).
Google Scholar 
Weiss, D. J. et al. Air temperature suitability for Plasmodium falciparum malaria transmission in Africa 2000–2012: A high-resolution spatiotemporal prediction. Malar. J. 13, 1–11 (2014).
Google Scholar 
Watts, A. G. et al. Elevation as a proxy for mosquito-borne zika virus transmission in the Americas. PLoS ONE 12, 1–16 (2017).
Google Scholar 
Shah, H. A., Dritsaki, M., Pink, J. & Petrou, S. Psychometric properties of Patient Reported Outcome Measures (PROMs) in patients diagnosed with Acute Respiratory Distress Syndrome (ARDS). Health Qual. Life Outcomes 14, 15 (2016).PubMed 
PubMed Central 

Google Scholar 
Eneanya, O. A. et al. Environmental suitability for lymphatic filariasis in Nigeria. Parasites Vectors 11, 1–13 (2018).
Google Scholar 
Craney, T. A. & Surles, J. G. Model-dependent variance inflation factor cutoff values. Qual. Eng. 14, 391–403 (2002).
Google Scholar 
Anderson, D. & Burnham, K. Model Selection and Multimodel Inference (Springer, 2002).MATH 

Google Scholar 
Giam, X. & Olden, J. D. Quantifying variable importance in a multimodel inference framework. Methods Ecol. Evol. 7, 388–397 (2016).
Google Scholar 
Guo, G. & Zhao, H. Multilevel modeling for binary data. Annu. Rev. Sociol. 26, 441–462 (2000).
Google Scholar 
Li, B., Lingsma, H. F., Steyerberg, E. W. & Lesaffre, E. Logistic random effects regression models: A comparison of statistical packages for binary and ordinal outcomes. BMC Med. Res. Methodol. 11, 1–11 (2011).
Google Scholar 
Symonds, M. R. E. & Moussalli, A. A brief guide to model selection, multimodel inference and model averaging in behavioural ecology using Akaike’s information criterion. Behav. Ecol. Sociobiol. 65, 13–21 (2011).
Google Scholar 
Burnham, K. P., Anderson, D. R. & Huyvaert, K. P. AIC model selection and multimodel inference in behavioral ecology: Some background, observations, and comparisons. Behav. Ecol. Sociobiol. 65, 23–35 (2011).
Google Scholar 
Shmueli, G. To explain or to predict?. Stat. Sci. 25, 289–310 (2010).MathSciNet 
MATH 

Google Scholar 
Hamlet, A. et al. The seasonal influence of climate and environment on yellow fever transmission across Africa. PLoS Negl. Trop. Dis. 12, 1–17 (2018).
Google Scholar 
Foley, J. A. et al. Global consequences of land use. Science 309, 570–574 (2005).ADS 
CAS 
PubMed 

Google Scholar 
Norton, E. C., Dowd, B. E. & Maciejewski, M. L. Marginal effects: Quantifying the effect of changes in risk factors in logistic regression models. JAMA 320, 84–85 (2018).PubMed 

Google Scholar 
RStudio Team. R Studio: Integrated Development for R (2015).Bates, D., Mächler, M., Bolker, B. & Walker, S. Fitting linear mixed-effects models using lme4. J. Stat. Softw. 67, 1–48 (2015).
Google Scholar  More

Since the bone discontinuities noted in the three vertebrae analyzed show no clear sign of bone overgrowth, it is pivotal to rule out the possibility that we are dealing with preservation damages before proposing an accurate diagnosis for the lesions. The close-up view examination of the abnormalities shows that their edges have clear signs of smoothing and rounding (Fig. 1), which represent important evidence of osteoblastic activity18,19. Additionally, the similar color of the cortical damage and normal bone can be used as secondary evidence to rule out post-mortem processes as a possible origin of the alterations, since recent destructive processes are lighter than the rest of the bone19. Therefore, as taphonomic processes can be ruled out, the pointed evidence strongly suggests that the discontinuities observed are of pathological origin. More specifically, these breaks found in all three vertebrae are indicative of bone fracture.Based on fracture analysis criteria applied here20, which consider the location and morphological pattern of the fractures, we classified the fractures noted in all vertebrae as traumas belonging to Type A (vertebral body compression), Group A2 (split fractures), and subgroup A2.1 (sagittal split fracture). This diagnosis implies that the traumatic episode was likely caused by a compressive force on the vertebral column, which split the vertebral bodies in the sagittal plane. This type of injury is considered stable—i.e., the fracture does not have a tendency to displace after reduction—and neurological deficit is uncommon20,22,23. Although stable traumas cause only moderate pain, without generating significant movement limitations20, the Eremotherium individual here analyzed died with unhealed bones, as there is no evidence of callus formation.The absence of other skeletal signs that point to the presence of another type of disease concomitantly to the fractures allows us to reject the possibility that they have been generated as a result of a pre-existing disease (e.g., infection, neoplasm). We also consider that the vertebral injuries were not caused by repetitive force (stress fractures) because this type of injury is commonly characterized as a nondisplaced line or crack in the bone, called hairline fracture3. Those refer to situations where the broken bone fragments are not visibly out of alignment and exhibit very little relative displacement21. Although the Eremotherium vertebrae fractures’ can be described as nondisplaced, they also have a noticeable gap between their edges that is mostly narrow with wider parts in the middle, something found in split fractures20 but that is not characteristic of hairline fractures. Lastly, the subgroup C1.2.1 (rotational sagittal split fracture) might be a source of confusion due to similar morphological pattern with subgroup A2.1 (sagittal split fracture). However, in subgroup C1.2.1 there are compressive and rotational forces acting simultaneously, producing total separation into two parts20, which clearly did not occur in the vertebrae analyzed here.In humans, compression fractures are most commonly caused by osteoporosis, although infection, neoplasm and trauma can also be etiological factors23,24,25. However, as aforementioned, the absence of other pathological skeletal marks is an important characteristic to take note as it serves to disregard the possibility of the fractures’ genesis to be secondary to another pathology. As such, in this case, osteoporosis, infection and neoplasm are unlikely etiologies. On the other hand, a compression fracture in a healthy individual is commonly generated after a severe traumatic event such as a fall from great height23,26. This scenario seems to better explain the origin of the vertebral fractures in the case of the Eremotherium ground sloth herein studied.The three fractured vertebrae were recovered in the Toca das Onças site (Fig. 2), a small cave considered as one of the richest paleontological sites of the Brazilian Quaternary15. Two complete skeletons of Eremotherium laurillardi and fragments belonging to at least thirteen other individuals, together with several other bones assigned to different smaller species are known to this cave14. It comprises of a single dry chamber that can only be entered through vertical entrances approximately 4.5 m high (Figs. 2b–d and 3). Two different hypotheses concerning the depositional process of Toca da Onças were previously proposed: (1) the animals climbed down into the cave in search of water14; or (2) due to the vertical character of the cave entrance, it could have functioned as a natural trap where animals accidentally fell into the cave15.Figure 2Location map of the Toca das Onças site and images of the cave. (a) Detail of the location, (b) cave entrance area view, (c) view from inside the cave, (d) Cave entrance detail. Scale bars 10 m in (b) and 5 m in (c). This figure was generated by Adobe Photoshop CS6 software (https://www.adobe.com/br/products/photoshop.html).Full size imageFigure 3Schematic representation of the Toca das Onças site. (a) Ground plan of the cave illustrating its morphology and dimension, (b) Cross-section illustrating the abyss-shaped entrance.Full size imageThe first hypothesis would indicate that the animal fell into the cave during an attempt to climb down. However, there is no report in the literature indicating that Eremotherium laurillardi could have been a climbing animal. In addition, the vertical morphology of the cave entrance would be a limiting factor for climbing behavior (see Fig. 3).Therefore, based on the type of fracture (compression sagittal split fracture) observed in the three vertebrae of Eremotherium as well as the inferred origin mechanism (fall from a great height), the presence of the individual here analyzed in the fossil accumulation of Toca das Onças is more likely explained by the second hypothesis. This idea is not particularly new as ‘entrapment due to fall’ has been described as a fossil accumulation mode to several other caves worldwide (e.g.,27,28). However, the use of bones fractures as an indicator of fossil accumulation mode is an interesting novelty. Of course, a detailed taphonomic investigation in the Toca das Onças still needs to be conducted in order to accurately interpret the formation of this important Quaternary fossil accumulation from Brazil.In sum, we suggest that the animal accidentally fell into the cave, fractured at least three sequential vertebrae (12th, 13th thoracic vertebrae and 1st lumbar vertebra) after the impact on the ground, survived for a while, but succumbed trapped inside the cave without food and water (Fig. 4). Other animals found in the cave, but without signs of bone fracture, may have fallen and not fractured their bones or not survived after the fall, especially the smaller ones. Finally, the proposal of falls to explain the unusual record of giant ground sloth fossils preserving much of its skeleton in caves, as reported for Toca das Onças site, contrasts with the better-documented pattern of skeletal accumulation via hydraulic action.Figure 4Artistic reconstruction of the suggested fall of the individual Eremotherium laurillardi into the cave. Artwork by Júlia d’Oliveira.Full size image More

Our compiled dataset consisted of 1674 records of marine mammal records after removing duplicate reports. It included 660 reports of sightings, 59 reports of induced mortalities or hunting records, 240 reports of incidental mortalities, 632 unique stranding records (live / dead), and 83 records which could not be categorised because of incomplete information.SightingsA total of 660 opportunistic sightings (number of individuals, ni = 3299) were recorded throughout the Indian coastline between 1748 and 2017 (Fig. 1a, 2a, 3a). Sighting data on the east coast (species = 18, ni = 1105) was mostly restricted to Odisha and Tamil Nadu (representing 97% of total east coast sightings). On the west coast (ni = 1297), Maharashtra (ni = 549), Gujarat (ni = 248) and Karnataka (ni = 307) contributed to highest sighting records (representing 85% of total west coast sightings). Sightings from the islands also contributed to 24.85% of the dataset (Andaman & Nicobar Islands = 24.37%, Lakshadweep = 0.48%). Highest incidence of sightings was for DFP (ni = 1894) followed by dugongs (ni = 959), BW (ni = 58) and SBW (ni = 17).Figure 1Marine mammal records obtained from data compiled between years 1748 – 2017 along the east coast, west coast and the islands of India for the groups i.e., baleen whales (BW), dolphins and finless porpoise (DFP), sperm and beaked whales (SBW) and dugongs, given as color-coded stacked bars where (a) sighting records—records where live animals were sighted (b) induced mortalities—records where animals were reported hunted or killed or were driven ashore, (c) incidental mortalities—records where animals were found dead after entanglement in fishing nets or being struck by vessels and (d) stranding records—records where dead or live animals were found washed ashore, or floating near shore or stranded alive and were attempted for rescue.Full size imageFigure 2Marine mammal records obtained every year from the data compiled between years 1748–2017 along Indian coastline given as cumulative numbers for each group i.e., baleen whales (BW), dolphins and finless porpoise (DFP), sperm and beaked whales (SBW) and dugongs, as color-coded stacked bars, where (a) sighting records—records where live animals were sighted (b) induced mortalities—records where animals were reported hunted or killed or were driven ashore, (c) incidental mortalities—records where animals were found dead after entanglement in fishing nets or being struck by vessels and (d) stranding records—records where dead or live animals were found washed ashore, or floating near shore or stranded alive and were attempted for rescue.Full size imageFigure 3Bubble plots showing distribution of marine mammal records obtained from data compiled between years 1748–2017 along the Indian coastline for each group i.e., baleen whales (BW), dolphins and finless porpoise (DFP), sperm and beaked whales (SBW) and dugongs, as color-coded stacked bars, where (a) sighting—records where live animals were sighted (b) induced mortalities—records where animals were reported hunted or killed or were driven ashore, (c) incidental mortalities—records where animals were found dead after entanglement in fishing nets or being struck by vessels and (d) strandings—records where dead or live animals were found washed ashore, or floating near shore or stranded alive and were attempted for rescue. Size of the bubble indicates number of individuals. These maps were created using ArcGIS 10.5 (https://desktop.arcgis.com/en/arcmap/10.3/map/working-with-layers/about-symbolizing-layers-to-represent-quantity.htm).Full size imageInduced mortalitiesA total of 59 incidences (ni = 102) were recorded of marine mammals being hunted/ captured between the years 1748–2017 (Fig. 1b, 2b, 3b). The total number of animals hunted/ captured deliberately is similar along east coast (ni = 33), west coast (ni = 29) and islands (ni = 36). Out of all marine mammal species, 90% of the animals hunted at the east coast were dugong D. dugon (ni = 30, all from Tamil Nadu). On the west coast, records of hunting incidences of finless porpoise Neophocaena phocaenoides were highest (79% of total records on west coast, Goa ni = 17, Kerala ni = 4, Karnataka and Maharashtra ni = 1). In the islands (i.e., Andaman and Nicobar Islands), 94% of the hunting records were of dugongs (ni = 34).Incidental mortalitiesA total of 240 net entanglements (ni = 1356) were reported along the Indian coast between the years 1748 and 2017 (Fig. 1c, 2c, 3c). Similar counts of individuals entangled along east (ni = 670) and west coast (ni = 654) were obtained with low reporting from the islands (ni = 26). Fourteen species were reported entangled from both east and west coast with only 4 species recorded from the islands. D. dugon was found to be most frequently entangled along the east coast (63 incidences, ni = 594, contributing to 56% of the total numbers on east coast), followed by Tursiops sp. (11 incidences, ni = 14, 9% of the east coast dataset). On the west coast, Tursiops sp. was the most frequently entangled (18 incidences, ni = 117, contributing to 18% of the west coast dataset), followed by N. phocaenoides (17 incidences, ni = 34, contributing to 17% of the dataset). The total number of DFP being entangled from west coast (ni = 623) were higher than east coast (ni = 68). More dugong individuals were entangled along east coast (i.e., from Tamil Nadu; ni = 594) as compared to the west coast (i.e., Gujarat; ni = 3) and Islands (i.e., Andaman and Nicobar; ni = 19). D. dugon was the most frequently entangled species in the islands (19 incidences, ni = 19, contributing to 79% of the total numbers in islands dataset) followed by false killer whale Pseudorca crassidens (3 incidences, ni = 5, contributing to 12% of the islands dataset). Very few BW or SBW (11 incidences, ni = 11) were recorded accidently entangled throughout the Indian coastline.StrandingsMarine mammals stranding reports consisted of 91.93% dead (ni = 581) and 8.07% live strandings (ni = 51) (Figs. 1d, 2d, 3d). Considering mass strandings as strandings with ni  > 2 (excluding mother and calf;33,34), 8.5% of all reports were mass strandings (21 strandings, ni = 1054). Most of the records did not have information about the sex of the stranded animal (83%), the age class (88%) or the state of decomposition of the carcass (53%). Highest strandings were reported of dugongs (strandings = 190, ni = 228), followed by BW (strandings = 178, ni =  = 190), DFP (strandings = 157, ni =  = 552) and SBW (strandings = 47, individuals = 48). There were 54 incidences (ni = 54, 9% of total stranding data) where the animal was not identified reliably to include in either of the groups.Species composition and frequencies of strandings were different on east coast, west coast and in the islands (Fig. 1, Table 1). Twenty-two species were reported as stranded on the east coast with D. dugon as the most frequently stranded species (83 incidences, ni = 107, ~ 29% of all records), followed by Indo-Pacific humpback dolphin Sousa chinensis, (31 incidences, ni = 108, ~ 10% of all records). On the west coast, out of 20 species reported as stranded, Balaenoptera musculus was most frequent (28 incidences, ni = 29, ~ 12% of all records) followed by N. phocaenoides (23 incidences, ni = 39, ~ 10% of all records). In the islands, 13 species were reported as stranded, D. dugon (93 incidences, ni = 102, contributing to 77% of the total animals found on the islands) followed by strandings of sperm whale Physeter macrocephalus (8 incidences, ni = 8, contributing to 6% of the data; Table 1).

a. Baleen whales

Table 1 Number of stranding events reported for marine mammals between 1748–2017 in India from the east coast, the west coast and Lakshadweep and Andaman & Nicobar archipelagos.Full size tableA total of 178 BW strandings (ni = 190) were reported. Most species were unidentified (east coast ni= 27, west coast ni = 58, islands ni = 4; i.e., 47% of the data). Identified strandings comprised of 6 species (see Table 1), some of which were later found to be misidentification (no confirmed evidence for common Minke Whale Balaenoptera acutorostrata, Sei Whale Balaenoptera borealis and Fin Whale Balaenoptera physalus from Indian waters; MMRCNI, 2018). Higher number of strandings occurred on the west coast (ni = 126), as compared to east coast (ni = 60). The east and west coast reported all six species of BW, whereas only three species stranded on the islands. B. borealis (misidentified) was the most stranded species across the east coast (12 incidences, ni = 12, contributing to 11% of the data) whereas blue whale Balaenoptera musculus was the most frequent across the west coast (28 incidences, ni = 29, contributing to 11% of the data). Baleen whale strandings were rare in the islands (4 incidences, ni = 4).Forty-seven SBW strandings (ni = 48) were reported along the Indian coast. More SBW stranded on the east coast (ni = 23) as compared to the west coast (ni = 13) and the islands (ni = 12). P. macrocephalus was most frequently reported (70% of all SBW records, east coast ni = 20, west coast ni = 6, islands ni = 8).There were 157 strandings (ni =552) of DFP belonging to 14 species. Twenty-one of these events were mass strandings (ni  > 2). The largest mass stranding event (ni = 147) occurred of short-finned pilot whale Globicephala macrorhynchus along the west coast (Tamil Nadu). Higher number of DFP strandings were recorded from east coast (ni = 418) as compared to west coast (ni = 83) and the islands (ni = 51; Table 1). East coast received a higher diversity of stranded DFP (number of species = 11) as compared to west coasts (number of species = 9) and the islands (number of species = 3). S. chinensis was the most frequently stranded species along the east coast (31 incidences, ni = 108, contributing to 33% of the data) whereas N. phocaenoides was the most frequent along the west coast (23 incidences, ni = 39, contributing to 37% of the data; Table 1).

d. Dugongs

The current distribution of dugongs in India is in the shallow coastal waters of Gujarat, Tamil Nadu and Andaman & Nicobar Islands37,38. There are 190 stranding events recorded between the years 1893 and 2017. The highest number of stranded dugongs were recorded from Tamil Nadu (ni = 107) closely followed by Andaman and Nicobar Islands (ni = 102) and few records from Gujarat (ni = 19).Temporal stranding patternsOur analysis of temporal trends for the last 42 years (1975–2017) showed that the mean number of strandings along the Indian coast was 11.25 ± SE 1.39 / year. The number of stranding reports show an increasing trend for two decades after 1975, dropping between 1995 and 2004. We observed a distinct rise in strandings post 2005 (18.23 ± SE 2.98 / year) with the highest reports from 2015–17 (27.66 ± SE 8.51/year) (Fig. 4).

a. Baleen whales

Figure 4A beanplot of decadal trends in marine mammal stranding in India from data compiled between years 1975–2017. Data prior to 1975 was discontinuous over the years to be considered for decadal trends. The data for last decade considered here includes only two years (2015–17) where increased reporting is evident. The bold horizontal lines indicate the mean number of strandings in each decade whereas the smaller horizontal lines indicate stranding numbers recorded for each year within the decade.Full size imageOn the west coast, mean stranding rate throughout the years (1975–2017) was 0.0010 ± SE 0.0014 strandings/km, and a steady rise was observed in rate of reported strandings after 2010. A seasonal trend was observed as well, with a peak in the month of September (sr = 0.0061 ± SE 0.0016 strandings/km), i.e., towards the end of monsoon season, and lowest strandings were recorded in the month of June (sr = 0.0016 ± SE 0.006 strandings/ km) (Fig. 5).Figure 5Temporal patterns (annual and monthly stranding rates / 100 km of coastline) in strandings of marine mammal records obtained from data compiled between years 1975–2017 along east and west coast of India for each group where (a) annual stranding rate and (b) monthly stranding rate for baleen whales (BW); (c) annual stranding rate and (d) monthly stranding rate for dolphins and finless porpoise (DFP); (e) annual stranding rate and (f) monthly stranding rate for sperm and beaked whales (SBW) and (g) annual stranding rate and (h) monthly stranding rate for dugongs.Full size imageThe mean stranding rate of BW on the east coast through 1975–2017 was 0.0013 ± SE 0.0017 strandings/km, but no specific trends were observed according to years or seasons. Stranding rates of BW did not differ between east and west coast (Mann–Whitney U test, U = 390, U standardized = -0.025, p value  > 0.05).The stranding rates of SBW differed significantly along both the coasts (Mann Whitney U test, U = 192, U standardized = 0.0, p value  More

Present fire-climate classificationTo identify the different regions of the planet with suitable climatic conditions for fire activity, we compare the global distribution of climate indicators based on temperature and precipitation, with satellite-derived GFED4 burned area data21 (Fig. 1). Starting from four general climates (Tr-tropical, Ar-arid, Te-temperate and Bo-boreal) based on the Köppen–Geiger climate classification main categories22, we create four fire-prone classes using climate thresholds to define the patterns observed in Fig. 1. Each category is characterised by the element that boosts fire activity during the FS: low precipitation, high temperatures or a combination of both. The classification is made by contrasting the probability distribution of the climatic variables at data points associated with high fire activity vs. points with low fire activity within the main Köppen-Geiger categories (see Threshold Selection in Methods section for a detailed explanation).Fig. 1: Burned area observations and climate drivers.a 1996–2016 maximum annual burned area (BAmax) and monthly burned area time series for selected regions. b Average monthly precipitation percentage from the annual total for the fire season (PPFS). c Average monthly temperature anomaly from the annual mean for the fire season (TAFS).Full size imageThe environmental conditions associated with fire occurrence emerge more clearly in this comparison, yielding the different threshold sets in Table 1 that determine the fire-prone months at any location (the selection method is detailed in the Methods section). We define those years with at least 1-month meeting the thresholds, as fire-prone years (FPY). Depending on the number of FPY at each location, the categories of Table 1 are sub-divided into recurrent (r), occasional (o) and infrequent (i) (see Methods). The average number of fire-prone months during the FPY is defined as the potential FS length (PFSL), i.e., the season with climatic characteristics prone to fire activity.Table 1 Fire classification defining criteria.Full size tableFigure 2a depicts the global map of the burned areas classified according to the selected thresholds (Table 1). Savanna fires are responsible for the largest proportion of burned area on the global scale21. The FS in these areas is longer than in other climates (see Supplementary Fig. 1) and, despite savanna fires being also dependent on ignition patterns and human policies and practices, the FS is tied to a pronounced seasonal cycle of precipitation23,24,25, with fire occurring mainly during the dry part of the cycle. Because of this, the Tropical – dry season fire class (Tr-ds) coincides with the distribution of the tropical savanna climate. In Fig. 2, boreal fires are represented as hot season fires (Bo-hs) due to the large positive temperature anomaly existing in those locations during the FS (Fig. 1c). In fact, temperature variations explain much of the variability in boreal burned area26,27. Temperate fires are classified as dry and hot season (Te-dhs) because they affect regions where the dry season coincides with the warm season (Fig. 1b, c). Here, high temperatures and precipitation seasonality determine fire activity and inter-annual burned area variability, e.g., in Western North America28,29,30,31 and Southern Europe32,33. Fire activity in arid regions occurs during warm months, but the relation with precipitation is more complex. The FS is associated with a hot season in cooler (MAT  27.5 °C), the FS starts right at the beginning of the dry season (e.g., the Sahel, Supplementary Fig. 12) while where MATs are more moderate, between 18.5 and 27.5 °C, it takes longer to develop (e.g., Central Australia and the Kalahari desert, Supplementary Figs. 12 and 13). Due to the dependency between fires and the existence of fuel in arid climates, we named this class Arid fuel limited (Ar-fl). A more in-depth discussion about the definition of this fire-climate class can be found in the section entitled Threshold selection for each climate of the Supplementary Information.Fig. 2: Fire-prone region classification.a With observed burned area data as a reference: not classified (NC, white) and misclassified (C, black) areas with BAmax = 0 ha, unclassified (NC, grey) and classified (Tr-ds, Ar-fl, Te-dhs and Bo-hs) areas with BAmax  > 0 ha. Each class is subdivided into three subcategories depending on the recurrence of the fire-prone conditions: recurrent (r), occasional (o) and infrequent (i). b Present (1996–2016) fire-prone climatic regions. c Future (2070–2099) fire-prone climatic regions with shaded grey representing a  0 ha) or fireless (BA = 0 ha). This reveals a two-way relation between fires and climate: fires take place under specific climatic conditions, and most places with these climatic conditions are indeed fire-prone, which supports our earlier hypothesis. Fire activity is controlled by weather, resources to burn and ignitions, as represented through the fire regime triangle12,20. On broad temporal scales and large spatial scales, temperature and precipitation have an important impact on fire because these climate variables influence vegetation type and the abundance, composition, moisture content, and structure of fuels34. Although ignitions may be driving fires to a greater extent than temperature or precipitation at specific locations or events35, they do not seem to limit fire activity at coarse spatial and temporal resolutions, implying that where fuels are sufficient and atmospheric conditions are conducive to combustion, the potential for ignition exists, either by lightning or human causes13,20. For all these reasons, we can identify specific climates that are prone to fires.The areas classified as fire-prone in Fig. 2b comprise 99.26% of the observed global mean annual burned area in Supplementary Fig. 2. This percentage is above 85% for all four general climates (Supplementary Fig. 20). The percentage of land area with non-zero burned area data classified as fire-prone is 91.22%. Considering for each location only the obtained FPY, the percentage of the observed burned area classified is 90.36%. Furthermore, the PFS obtained in the fire-climate classification (Fig. 3b) also correlates well with the timing of observed fire incidence, as globally 87.91% of the observed mean burned area occurs during the identified months of PFS at classified fire-prone locations.Fig. 3: Potential fire season.a Future minus present potential fire season length (PFSL) difference in months (ΔPFSL). b Present potential fire season. c Future potential fire season.Full size imageUnclassified regions (in grey in Fig. 2a) correspond for the most part to those with the least burned area or those where agricultural practices modify the climatic seasonality of fires. In addition, as the classification is conceived from a climatic point of view, locations with fire activity associated with specific meteorological conditions that are not appreciable at the monthly temporal resolution, are probably not well identified. For example, a week of extremely high temperatures could be almost unnoticeable in the monthly mean temperature, but not in fire activity. Similarly, months with the same total precipitation may have different fire activity if the precipitation falls concentrated in a few days or is distributed throughout the month. Furthermore, the short temporal sampling period of the burned area data could also be influencing our results. Locations with long fire cycles may not be well represented in the data.Future fire-climate classificationA future fire-climate classification map is derived by applying the thresholds obtained in the present fire-climate classification to future climatology variables from multiple coupled model intercomparison project phase 5 (CMIP5) global circulation model (GCM) outputs, considering the RCP8.5 scenario (the worst-case climate change scenario of the CMIP5). Two contrasting approaches can be taken for analysing future fire activity, one that considers quick vegetation adaptation to the new climatic conditions, and another that does not. These two approaches clearly diverge in the boreal regions, where the biome (mainly taiga) is strongly conditioned by the low temperatures and where future temperature changes at the end of the 21st century will have a greater amplitude. It is expected that the boreal forest of these areas will not be immediately replaced by a temperate mixed forest where the average annual temperature exceeds the range of values typical of the taiga biome. Terrestrial vegetation compositional and structural change could occur during the 21st century where vegetation disturbance is accelerated or amplified by human activity, but equilibrium states may not be reached until the 22nd century or beyond36.Based on the assumption that during the future period (2070–2099) the vegetation will not be fully adapted to the new climatic conditions, and since the present Köppen–Geiger climate classification (on which we base our Tr, Ar, Te and Bo categories) closely corresponds to the different existent biomes22, we analyse only the projected changes in the specific fire-climate classification variables, maintaining the general division of Tropical, Arid, Temperate and Boreal regions as is in present climate conditions. The future fire-climate classification is shown in Fig. 2c.We note that we determine future fire activity from relationships of the latter with the present climate; however, these relationships might not be stationary. Our approach does not contemplate possible future changes in precipitation frequency if they are not noticeable in monthly precipitation amounts. Areas with the rising incidence of extreme precipitation events due to global warming37 may experience an increase in monthly precipitation but a decrease in rainy days, which may lead us to consider the conditions there less favourable for fire activity than they actually will be.Future changes in global fire activityModelled future fire-prone regions experience significant variations with respect to the present (Fig. 2b, c). Due to global warming, the Bo-hs fire class pertaining to boreal forests would spread over a larger area, reaching most of Northern Scandinavia and undergoing a southward and northward expansion in Canada, Alaska and Russia. This category may experience a percentual expansion of 47.0% according to our results. This expansion is more accentuated for the combination of the highest recurrence subcategories Bo-hs-r and Bo-hs-o, reaching a value of 111.5%.The conjunction of Te-dhs-r and Te-dhs-o fire classes of midlatitudes also undergoes a considerable expansion of 24.5% in the area (Fig. 2b, c). The most remarkable changes are expected in Southern China and Southern Europe. A large part of Europe transitions from an infrequent fire category to a more frequent fire category with Csa and Csb Mediterranean climates38.The Tr-ds fire classes with frequent fire-prone conditions in the Tropics presents fewer spatial changes (Fig. 2b, c), with a spatial contraction of 6.3%. The most important differences are found in South America. Some of the climate model results considered here indicate also that some parts of the Eastern Amazon rainforest will move from a non-fire class to Tr-ds fire class, as other studies have suggested39.The Arid fire-prone classes Ar-fl-r and Ar-fl-o would increase its area by 5.0%. Projected changes in the extent of this class are very sensitive to changes in annual precipitation, conducive to vegetation and fuel reduction or increment, thus there is significant uncertainty in the proximity of desert regions (Fig. 2c).Clearer conclusions can be drawn from the FPY and PFSL calculation (Figs. 3 and 4). The number of months meeting the set of conditions in Table 1 yields the estimated PFSL (Fig. 3b), and the number of years with at least 1-month meeting the thresholds, the FPY. In the boreal regions, we obtain a general lengthening of the PFS. The PFS of these areas is conditioned by temperature, so the amplified warming of Artic zones40 is expected to make the FS longer. Notwithstanding, in certain parts of Eastern Asia, the intense warming is counterbalanced by an increase of the precipitation in certain warm months (see Supplementary Figs. 21 and 22), leading to a slight shortening of our estimated PFS. There is evidence, however, that temperature increases may lead to drier fuels in the future despite the precipitation increase, thus augmenting fire risk, as some investigations have shown for Canada41. Our results agree in general with several other studies that have previously pointed towards an increase of the FSL in boreal areas1,17,42, even when some suggest a more pronounced lengthening in more northerly latitudes1,17. In terms of the frequency of years with fire-prone conditions, the conclusions are even clearer. A general increase of the FPY is observed, especially for northerly latitudes, where the differences reach values of more than +4 years per decade (Fig. 4a). This possible increase in fire activity in boreal areas may result in significant peatland combustion and a release of the large quantities of soil carbon that they store into the atmosphere43. These greenhouse gas emissions may create a positive feedback loop, leading to a further increase in temperature, which in turn will enhance boreal wildfire incidence and more peatland burning.Fig. 4: Fire-prone years.a Future minus a present number of years with at least one month classified as fire-prone per decade (ΔFPY). b Present fire-prone years per decade. c Future fire-prone years per decade.Full size imageThe Te-dhs fire class, corresponding to temperate climates, would also experience a general lengthening of the PFS (Fig. 3). A future precipitation decline may be especially significant in Southern Europe (Supplementary Fig. 21), associated with an increased anticyclonic circulation yielding more stable conditions44, while the temperature rise would be quite homogeneous among all Te-dhs fire-climate class areas. The FS drought intensification around the Mediterranean, together with the general warming (Supplementary Fig. 21), would lead to a lengthening of the PFS of around 2 months (Fig. 3a), but summer months could also experience this precipitation decline (Supplementary Fig. 22), meaning that the FS would be more severe. The Western US, which has already experienced over the last decades the lengthening of the FS45 and the increase of large fires46 and extreme wildfire weather47,48 due to climate change, may also experience an FS lengthening by the end of the 21st century. Some authors18,48,49,50 have studied projected fire future changes from other points of view (occurrence of very large fires, wildfire potential, etc.), finding also a general increase of fire severity by the end of the century in some of these Te-dhs fire regions. The interannual recurrence of fire-prone conditions will significantly increase in countries like France, Italy or Eastern China (Fig. 4a).The PFSL of the Tropical Tr-ds fire-climate class presents slight differences between present and future values (Fig. 3). Some areas of the Northern African savanna may experience a shortening of the PFS, while Southern Africa shows a lengthening. A dipole pattern of wetting in tropical Eastern Africa and drying in Southern Africa51 could be the reason for these future changes. There is a contrasting influence of ENSO in present African fire patterns52, which suggests that the future pattern of precipitation variations in Central Africa may be associated with ENSO future changes under climate change conditions53. Although the quantification of ENSO changes in a warmer climate is still an issue that continues to be investigated, an expansion and strengthening of ENSO teleconnections is confirmed by some authors53,54,55. The general increase in precipitation along all seasons in western equatorial Africa would lead to a significant decrease in the recurrence of interannual fire-prone conditions (Fig. 4a).Our results show that fire-prone areas in Temperate and especially Boreal climates are projected to undergo the most significant expansion and lengthening of the potential FS at the end of the XXI century driven by rising temperatures. In the Tropics, little change is expected in these respects. Notwithstanding, global warming is likely to make fire risk more severe mostly everywhere, and in particular in some regions such as Mediterranean Europe and the Eastern Amazon, where an important decrease in precipitation is also predicted during the PFS. More favourable fire conditions will potentially increment fire activity and burned areas in many places. In others, especially in the Tropics, increasing suppression efforts and a cease to agricultural and pastoral practices like vegetation clearing by fire, replaced by more intensive farming, could counteract the impact of a warmer climate. A reduction of these human-caused fires in the Tropics could bring global burned area down2, despite rising trends elsewhere, given the vast contribution of Tropical fires to the burned areas at the global scale (Fig. 1). More