Philippa Kaur

Nadell, C. D., Drescher, K. & Foster, K. R. Spatial structure, cooperation and competition in biofilms. Nat. Rev. Microbiol. 14, 589–600 (2016).Article 
CAS 

Google Scholar 
Palmer, J. D. & Foster, K. R. Bacterial species rarely work together. Science 376, 581–582 (2022).Article 
ADS 
CAS 

Google Scholar 
Pande, S. & Kost, C. Bacterial unculturability and the formation of intercellular metabolic networks. Trends Microbiol. 25, 349–361 (2017).Article 
CAS 

Google Scholar 
Nadell, C. D., Xavier, J. B. & Foster, K. R. The sociobiology of biofilms. FEMS Microbiol. Rev. 33, 206–224 (2009).Article 
CAS 

Google Scholar 
Fritts, R. K., McCully, A. L. & McKinlay, J. B. Extracellular metabolism sets the table for microbial cross-feeding. Microbiol. Mol. Biol. Rev. 85, 135 (2021).Article 

Google Scholar 
D’Souza, G. et al. Ecology and evolution of metabolic cross-feeding interactions in bacteria. Nat. Prod. Rep. 35, 455–488 (2018).Article 

Google Scholar 
Libby, E., Hébert-Dufresne, L., Hosseini, S.-R. & Wagner, A. Syntrophy emerges spontaneously in complex metabolic systems. PLoS Comput. Biol. 15, e1007169 (2019).Article 

Google Scholar 
Staley, J. T. & Konopka, A. Measurement of in situ activities of nonphotosynthetic microorganisms in aquatic and terrestrial habitats. Annu. Rev. Microbiol. 39, 321–346 (1985).Article 
CAS 

Google Scholar 
Zachar, I. Closing the energetics gap. Nat. Ecol. Evol. https://doi.org/10.1038/s41559-022-01839-3 (2022).Article 

Google Scholar 
Zachar, I. & Boza, G. Endosymbiosis before eukaryotes: mitochondrial establishment in protoeukaryotes. Cell. Mol. Life Sci. 77, 3503–3523. https://doi.org/10.1007/s00018-020-03462-6 (2020).Article 
CAS 

Google Scholar 
Zachar, I. & Szathmáry, E. Breath-giving cooperation: critical review of origin of mitochondria hypotheses. Biol. Direct 12, 19. https://doi.org/10.1186/s13062-017-0190-5 (2017).Article 
CAS 

Google Scholar 
Booth, A. & Doolittle, W. F. Eukaryogenesis, how special really?. Proc. Natl. Acad. Sci. 112, 10278–10285 (2015).Article 
ADS 
CAS 

Google Scholar 
Morris, B. E. L., Henneberger, R., Huber, H. & Moissl-Eichinger, C. Microbial syntrophy: Interaction for the common good. FEMS Microbiol. Rev. 37, 384–406 (2013).Article 
CAS 

Google Scholar 
Szathmáry, E. On the propagation of a conceptual error concerning hypercycles and cooperation. J. Syst. Chem. 4, 2208 (2013).Article 

Google Scholar 
Seth, E. C. & Taga, M. E. Nutrient cross-feeding in the microbial world. Front. Microbiol. 5, 350 (2014).Article 

Google Scholar 
Piccardi, P., Vessman, B. & Mitri, S. Toxicity drives facilitation between 4 bacterial species. Proc. Natl. Acad. Sci. 116, 15979–15984 (2019).Article 
ADS 
CAS 

Google Scholar 
Yurtsev, E. A., Conwill, A. & Gore, J. Oscillatory dynamics in a bacterial cross-protection mutualism. Proc. Natl. Acad. Sci. 113, 6236–6241 (2016).Article 
ADS 
CAS 

Google Scholar 
Kehe, J. et al. Positive interactions are common among culturable bacteria. Sci. Adv. 7, 45 (2021).Article 

Google Scholar 
Momeni, B., Xie, L. & Shou, W. Lotka-Volterra pairwise modeling fails to capture diverse pairwise microbial interactions. Elife 6, 25051 (2017).Article 

Google Scholar 
Zengler, K. & Zaramela, L. S. The social network of microorganisms: How auxotrophies shape complex communities. Nat. Rev. Microbiol. 16, 383–390 (2018).Article 
CAS 

Google Scholar 
Koschwanez, J. H., Foster, K. R. & Murray, A. W. Sucrose utilization in budding yeast as a model for the origin of undifferentiated multicellularity. PLoS Biol. 9, e1001122 (2011).Article 
CAS 

Google Scholar 
Ciofu, O., Beveridge, T. J., Kadurugamuwa, J., Walther-Rasmussen, J. & Høiby, N. Chromosomal beta-lactamase is packaged into membrane vesicles and secreted from Pseudomonas aeruginosa. J. Antimicrob. Chemother. 45, 9–13 (2000).Article 
CAS 

Google Scholar 
Xenophontos, C., Harpole, W. S., Küsel, K. & Clark, A. T. Cheating promotes coexistence in a two-species one-substrate culture model. Front. Ecol. Evol. 9, 78006 (2022).Article 

Google Scholar 
West, S. A., Diggle, S. P., Buckling, A., Gardner, A. & Griffin, A. S. The social lives of microbes. Annu. Rev. Ecol. Evol. Syst. 38, 53–77 (2007).Article 

Google Scholar 
Flemming, H.-C. & Wingender, J. The biofilm matrix. Nat. Rev. Microbiol. 8, 623–633 (2010).Article 
CAS 

Google Scholar 
Kümmerli, R. & Brown, S. P. Molecular and regulatory properties of a public good shape the evolution of cooperation. Proc. Natl. Acad. Sci. 107, 18921–18926 (2010).Article 
ADS 

Google Scholar 
Griffin, A. S., West, S. A. & Buckling, A. Cooperation and competition in pathogenic bacteria. Nature 430, 1024–1027 (2004).Article 
ADS 
CAS 

Google Scholar 
Kramer, J., Özkaya, Ö. & Kümmerli, R. Bacterial siderophores in community and host interactions. Nat. Rev. Microbiol. 18, 152–163 (2019).Article 

Google Scholar 
van der Meij, A., Worsley, S. F., Hutchings, M. I. & van Wezel, G. P. Chemical ecology of antibiotic production by actinomycetes. FEMS Microbiol. Rev. 41, 392–416 (2017).Article 

Google Scholar 
Kümmerli, R., Schiessl, K. T., Waldvogel, T., McNeill, K. & Ackermann, M. Habitat structure and the evolution of diffusible siderophores in bacteria. Ecol. Lett. 17, 1536–1544 (2014).Article 

Google Scholar 
Jautzus, T., van Gestel, J. & Kovács, Á. T. Complex extracellular biology drives surface competition in lessigreaterBacillus subtilisless/igreater. Ecol. Lett. 16, 2320–2328. https://doi.org/10.1101/2022.02.28.482363 (2022).Article 
CAS 

Google Scholar 
Sachs, J. L., Mueller, U. G., Wilcox, T. P. & Bull, J. J. The evolution of cooperation. Q. Rev. Biol. 79, 135–160 (2004).Article 

Google Scholar 
Hillesland, K. L. & Stahl, D. A. Rapid evolution of stability and productivity at the origin of a microbial mutualism. Proc. Natl. Acad. Sci. 107, 2124–2129 (2010).Article 
ADS 
CAS 

Google Scholar 
Bruno, J. F., Stachowicz, J. J. & Bertness, M. D. Inclusion of facilitation into ecological theory. Trends Ecol. Evol. 18, 119–125 (2003).Article 

Google Scholar 
Gore, J., Youk, H. & van Oudenaarden, A. Snowdrift game dynamics and facultative cheating in yeast. Nature 459, 253–256 (2009).Article 
ADS 
CAS 

Google Scholar 
Sorg, R. A. et al. Collective resistance in microbial communities by intracellular antibiotic deactivation. PLoS Biol. 14, e2000631 (2016).Article 

Google Scholar 
Karray, F. et al. Extracellular hydrolytic enzymes produced by halophilic bacteria and archaea isolated from hypersaline lake. Mol. Biol. Rep. 45, 1297–1309 (2018).Article 
CAS 

Google Scholar 
Datta, M. S., Sliwerska, E., Gore, J., Polz, M. F. & Cordero, O. X. Microbial interactions lead to rapid micro-scale successions on model marine particles. Nat. Commun. 7, 11965 (2016).Article 
ADS 
CAS 

Google Scholar 
Tarnita, C. E. The ecology and evolution of social behavior in microbes. J. Exp. Biol. 220, 18–24 (2017).Article 

Google Scholar 
Özkaya, Ö., Xavier, K. B., Dionisio, F. & Balbontn, R. Maintenance of microbial cooperation mediated by public goods in single- and multiple-trait scenarios. J. Bacteriol. 199, 22 (2017).Article 

Google Scholar 
Yang, D.-D. et al. Fitness and productivity increase with ecotypic diversity among Escherichia coli strains that coevolved in a simple, constant environment. Appl. Environ. Microbiol. 86, 8 (2020).Article 

Google Scholar 
Pande, S. et al. Fitness and stability of obligate cross-feeding interactions that emerge upon gene loss in bacteria. ISME J. 8, 953–962 (2013).Article 

Google Scholar 
Zhou, K., Qiao, K., Edgar, S. & Stephanopoulos, G. Distributing a metabolic pathway among a microbial consortium enhances production of natural products. Nat. Biotechnol. 33, 377–383 (2015).Article 
CAS 

Google Scholar 
Harcombe, W. R., Chacón, J. M., Adamowicz, E. M., Chubiz, L. M. & Marx, C. J. Evolution of bidirectional costly mutualism from byproduct consumption. Proc. Natl. Acad. Sci. 115, 12000–12004 (2018).Article 
ADS 
CAS 

Google Scholar 
Summers, Z. M. et al. Direct exchange of electrons within aggregates of an evolved syntrophic coculture of anaerobic bacteria. Science 330, 1413–1415 (2010).Article 
ADS 
CAS 

Google Scholar 
Maddamsetti, R., Lenski, R. E. & Barrick, J. E. Adaptation, clonal interference, and frequency-dependent interactions in a long-term evolution experiment with Escherichia coli. Genetics 200, 619–631 (2015).Article 
CAS 

Google Scholar 
Gerrish, P. J. & Lenski, R. E. The fate of competing beneficial mutations in an asexual population. Genetica 102, 127–144 (1998).Article 

Google Scholar 
Popat, R. et al. Quorum-sensing and cheating in bacterial biofilms. Proc. R. Soc. B 279, 4765–4771 (2012).Article 
CAS 

Google Scholar 
Rainey, P. B. & Rainey, K. Evolution of cooperation and conflict in experimental bacterial populations. Nature 425, 72–74 (2003).Article 
ADS 
CAS 

Google Scholar 
Hardin, G. Tragedy of the commons. Science 162, 1243 (1968).Article 
ADS 
CAS 

Google Scholar 
West, S. A., Cooper, G. A., Ghoul, M. B. & Ten Griffin, A. S. recent insights for our understanding of cooperation. Nat. Ecol. Evol. 5, 419–430 (2021).Article 

Google Scholar 
MacArthur, R. Species packing and competitive equilibrium for many species. Theor. Popul. Biol. 1, 1–11 (1970).Article 
CAS 

Google Scholar 
Oliveira, N. M., Niehus, R. & Foster, K. R. Evolutionary limits to cooperation in microbial communities. Proc. Natl. Acad. Sci. 111, 17941–17946 (2014).Article 
ADS 
CAS 

Google Scholar 
Tilman, D. Resource Competition and Community Structure. Monographs in Population Biology, Vol. 17 (Princeton University Press, 1982).
Google Scholar 
Ferenci, T. Trade-off mechanisms shaping the diversity of bacteria. Trends Microbiol. 24, 209–223 (2016).Article 
CAS 

Google Scholar 
Rozen, D. E., Philippe, N., de Visser, J. A., Lenski, R. E. & Schneider, D. Death and cannibalism in a seasonal environment facilitate bacterial coexistence. Ecol. Lett. 12, 34–44 (2009).Article 

Google Scholar 
Brännström, Å., Johansson, J. & von Festenberg, N. The Hitchhiker’s Guide to Adaptive Dynamics. Games 4, 304–328 (2013).Article 
MATH 

Google Scholar 
Ramin, K. I. & Allison, S. D. Bacterial tradeoffs in growth rate and extracellular enzymes. Front. Microbiol. 10, 2956 (2019).Article 

Google Scholar 
Imachi, H. et al. Isolation of an archaeon at the prokaryote–eukaryote interface. Nature 577, 519–525 (2020).Article 
ADS 
CAS 

Google Scholar 
Wintermute, E. H. & Silver, P. A. Emergent cooperation in microbial metabolism. Mol. Syst. Biol. 6, 407 (2010).Article 

Google Scholar 
Libby, E., Kempes, C. & Okie, J. Metabolic compatibility and the rarity of prokaryote endosymbioses. BioRxiv https://doi.org/10.1101/2022.04.14.488272 (2022).Article 

Google Scholar 
Pauli, B., Oña, L., Hermann, M. & Kost, C. Obligate mutualistic cooperation limits evolvability. Nat. Commun. 13, 27630 (2022).Article 

Google Scholar 
Oña, L. & Kost, C. Cooperation increases robustness to ecological disturbance in microbial cross-feeding networks. Ecol. Lett. 25, 1410–1420 (2022).Article 

Google Scholar 
Machado, D. et al. Polarization of microbial communities between competitive and cooperative metabolism. Nat. Ecol. Evol. 5, 195–203 (2021).Article 

Google Scholar 
Mee, M. T., Collins, J. J., Church, G. M. & Wang, H. H. Syntrophic exchange in synthetic microbial communities. Proc. Natl. Acad. Sci. 111, E2149–E2156 (2014).Article 
ADS 
CAS 

Google Scholar 
Goldford, J. E. et al. Emergent simplicity in microbial community assembly. Science 361, 469–474 (2018).Article 
ADS 
CAS 

Google Scholar 
McCutcheon, J. P. The genomics and cell biology of host-beneficial intracellular infections. Annu. Rev. Cell Dev. Biol. 37, 115–142 (2021).Article 
CAS 

Google Scholar 
Sousa, F. L., Neukirchen, S., Allen, J. F., Lane, N. & Martin, W. F. Lokiarchaeon is hydrogen dependent. Nat. Microbiol. 1, 5 (2016).Article 

Google Scholar 
Spang, A. et al. Proposal of the reverse flow model for the origin of the eukaryotic cell based on comparative analyses of Asgard archaeal metabolism. Nat. Microbiol. 4, 1138–1148 (2019).Article 
CAS 

Google Scholar 
Martin, W. & Müller, M. The hydrogen hypothesis for the first eukaryote. Nature 392, 37–41 (1998).Article 
ADS 
CAS 

Google Scholar 
López-García, P. & Moreira, D. The Syntrophy hypothesis for the origin of eukaryotes revisited. Nat. Microbiol. 5, 655–667 (2020).Article 

Google Scholar 
Mills, D. B. et al. Eukaryogenesis and oxygen in Earth history. Nat. Ecol. Evol. 6, 520–532 (2022).Article 

Google Scholar 
Liu, Y. et al. Expanded diversity of Asgard archaea and their relationships with eukaryotes. Nature 593, 553–557 (2021).Article 
ADS 
CAS 

Google Scholar 
Zachar, I., Szilágyi, A., Számadó, S. & Szathmáry, E. Farming the mitochondrial ancestor as a model of endosymbiotic establishment by natural selection. Proc. Natl. Acad. Sci. USA. 115, E1504–E1510. https://doi.org/10.1073/pnas.1718707115 (2018).Article 
ADS 
CAS 

Google Scholar 
Cavalier-Smith, T. & Chao, E.E.-Y. Multidomain ribosomal protein trees and the planctobacterial origin of neomura (eukaryotes, archaebacteria). Protoplasma https://doi.org/10.1007/s00709-019-01442-7 (2020).Article 

Google Scholar 
Searcy, D. G. Nutritional syntrophies and consortia as models for the origin of mitochondria. Symb. Mech. Model Syst. 1, 163–183. https://doi.org/10.1007/0-306-48173-1_10 (2002).Article 

Google Scholar 
Müller, N., Timmers, P., Plugge, C. M., Stams, A. J. M. & Schink, B. Syntrophy in methanogenic degradation. Endosymb. Methanog. Archaea 1, 153–192. https://doi.org/10.1007/978-3-319-98836-8_9 (2018).Article 

Google Scholar 
Searcy, D. G. Metabolic integration during the evolutionary origin of mitochondria. Cell Res. 13, 229–238 (2003).Article 
CAS 

Google Scholar 
Flemming, H.-C. & Wuertz, S. Bacteria and archaea on Earth and their abundance in biofilms. Nat. Rev. Microbiol. 17, 247–260 (2019).Article 
CAS 

Google Scholar 
Spang, A. et al. Asgard archaea are the closest prokaryotic relatives of eukaryotes. PLoS Genet. 14, e1007080 (2018).Article 

Google Scholar 
Burns, J. A., Pittis, A. A. & Kim, E. Gene-based predictive models of trophic modes suggest Asgard archaea are not phagocytotic. Nat. Ecol. Evol. 2, 697–704 (2018).Article 

Google Scholar 
Seitz, K. W. et al. Asgard archaea capable of anaerobic hydrocarbon cycling. Nat. Commun. 10, 1 (2019).Article 
ADS 
CAS 

Google Scholar 
Jimenez, P. & Scheuring, I. Density-dependent private benefit leads to bacterial mutualism. Evolution 75, 1619–1635. https://doi.org/10.1111/evo.14241 (2021).Article 

Google Scholar 
Preussger, D., Giri, S., Muhsal, L. K., Oña, L. & Kost, C. Reciprocal fitness feedbacks promote the evolution of mutualistic cooperation. Curr. Biol. 30, 3580-3590.e7 (2020).Article 
CAS 

Google Scholar 
Monaco, H. et al. Spatial-temporal dynamics of a microbial cooperative behavior resistant to cheating. Nat. Commun. 13, 3580 (2022).Article 
ADS 

Google Scholar 
Yanni, D., Márquez-Zacarias, P., Yunker, P. J. & Ratcliff, W. C. Drivers of spatial structure in social microbial communities. Curr. Biol. 29, 545–550 (2019).Article 

Google Scholar  More

Data collectionWe obtained trade data from two different sources: the United States Fish and Wildlife Service (USFWS) Law Enforcement Management Information System (LEMIS)31 and the International Union for Conservation of Nature (IUCN) Red List32. We used the former to obtain data on the live wildlife trade in general and the latter for data on the pet trade specifically. We then matched trade data with our previously compiled global scale datasets of life history traits and introductions in mammals, reptiles and amphibians25,26.We obtained data on the US live wildlife trade from LEMIS by a Freedom of Information Act Request on 12/08/2019. We requested summary data on all US imports and exports of wildlife across all available years (1999-2019) and all trade purposes, including information on species identities and shipment contents (e.g. live individuals, meat, skins, etc.). For each species, we summed the total number of recorded shipments of live individuals (including individuals that died in transit, and live eggs) as a measure of trade frequency. We classified species as in trade if there was at least one shipment of live individuals recorded in the LEMIS database, and as not traded otherwise. The LEMIS dataset is geographically limited to trade by the US, and therefore may not capture the full diversity of species involved in the wildlife trade. For example, the LEMIS database may be missing some species involved in the substantial trade in live wildlife between South–East Asian countries50. However, the US represents one of the most dominant players in the global market for live wildlife16, and by summing both imports and exports we capture demand for species in countries beyond the US to some extent. Supplementary Fig. 2 illustrates the frequency of trade between the US and countries represented in the US LEMIS dataset. LEMIS data should be considered a minimum estimate of the diversity of species involved in the wildlife trade since they mostly record only legal trade (although confiscated shipments are recorded), and shipments are sometimes not identified to the species level16,51,53,53. The LEMIS database also contains some mis-spelled and incorrectly identified species due to human input errors52. To minimise the effect of misidentified shipments on our species level classifications of US trade status, we discarded all LEMIS records that were not identified to the species level (i.e. those identified using genus, common or generic names only), and manually checked the LEMIS data for synonyms and alternate spellings when we could not automatically match any records in LEMIS with species in our life history datasets. Species classified as traded on the basis of a single recorded live shipment in LEMIS are most vulnerable to species level misclassification due to misidentified shipments. The vast majority of traded species have multiple shipments recorded in LEMIS (259/312 [83%] of traded mammals, 265/285 [93%] of traded reptiles and 72/75 [96%] of traded amphibians), reducing the potential impact of shipment level misidentification over the reliability of species level trade classifications. However, to investigate the robustness of our findings to possible errors in species identification in LEMIS, we re-ran our key analyses excluding species classified as traded on the basis of a single live shipment. We found qualitatively the same effects of life history traits on the probability of trade when removing these species as in our full sample (Supplementary Tables 25–27). Despite its limitations, LEMIS is an invaluable resource for identifying broad scale trends in the wildlife trade since few other countries maintain such detailed records, and it is the only large-scale international trade dataset that includes both CITES- and non-CITES-listed species16,41. Including non-CITES listed species in our analyses is important because CITES-listed species represent only a small minority of those in trade14 and are likely to be a biased sample in terms of life history traits, since species vulnerable to extinction typically have slower life histories40.We obtained separate data on the pet trade from the IUCN Red List. The IUCN has assessed the vast majority of mammal, reptile and amphibian species (91%, 79% and 86% respectively54). Here, we classified a species as involved in the pet trade if the IUCN species account included at least one clear description of involvement in the pet trade. Otherwise, we considered a species as not involved in the pet trade. Although LEMIS records the purpose of trade, it uses broad categories (e.g. ‘Commercial’, ‘Personal’, ‘Breeding in captivity’), none of which refers specifically to nor necessarily equates to trade for pets. Therefore, we sought this additional data on the pet trade from the IUCN Red List instead of following the approach of some previous studies which have used LEMIS data as a proxy for the pet trade (e.g. Refs. 15,19). In contrast, the IUCN Red List contains clear textual descriptions of use and trade for many species, allowing us to identify which species are traded specifically for pets32. The IUCN data has further complementary strengths compared with LEMIS in that it is global in scope and includes both legal and illegal trade. We obtained data from the IUCN Red List by manually searching the binomial name of each species in our samples and consulting the ‘Threats’ and ‘Use and Trade’ sections of the species accounts. We classified species as in the pet trade if the information clearly stated this was the case (e.g. “It has been recorded in the pet trade”, “This species appears in the international pet trade”). We discounted descriptions where the information was uncertain (e.g. the species is described as “probably” or “possibly” traded for pets). We did not count as pets those species that the IUCN categorises as used for “Pets/display animals, horticulture” but which are used only for zoos or captive display, such as beluga whales (Delphinapterus leucas). All species described as pets by the IUCN are ‘exotic’, i.e. those without a long history of domestication14, since the IUCN does not list domesticated species.We matched trade data with our previously published global scale datasets on life history traits and introductions25,26. Internationally traded species may or not be released in the wild outside their native range: some may remain in the confines of captivity (e.g. in zoos or kept by private owners). We defined a species as introduced if there was at least one reliable record of its release, by humans, into the wild outside of its native range, either accidentally or intentionally25,26. We included only species with complete data for the same life history traits as used in our prior analyses (mammals: body mass, gestation period, weaning age, neonatal body mass, litter size, litters per year, age at first reproduction and reproductive lifespan; reptiles: body mass, hatchling mass, clutch size, clutches per year, age of sexual maturity, reproductive lifespan and parity; amphibians: snout-vent length, egg size, clutch size, age of sexual maturity and reproductive lifespan) to facilitate direct comparisons with previous results and to allow us to account for covariation between life history traits55. Species with complete life history data represent 7.8%, 3.5% and 1.6% of the total estimated number of species of mammals, reptiles and amphibians respectively56,57,58. These samples are not random as they over-represent orders containing many species of interest and utility to humans (e.g. ungulates, primates, crocodilians) (Supplementary Tables 28–30). However, these biases are unlikely to undermine our results since we examine life history effects on trade and introduction within these samples. Trade and introduction data do not necessarily cover the same time periods: the US dataset covers only the years 1999-present and the IUCN descriptions also typically refer to recent trade. In contrast, our introduction dataset includes both historical and recent introductions25,26. Therefore, the goal of our analyses is not to test causal hypotheses on the direct relationship between trade and introduction but rather to investigate whether the same life history traits predispose species towards both trade and introduction across diverse taxa, locations and circumstances. When combining the datasets and phylogenies59,60,61,62,63, we resolved species name mis-matches by referring to taxonomic information from the IUCN Red List32, the Global Biodiversity Information Facility (GBIF33) and the Integrated Taxonomic Information System (ITIS64). Table 1 summarises final sample sizes and Supplementary Table 1 the degree of overlap between the trade datasets. Most species in the pet trade are also in the general live wildlife trade, but many more species are traded by the US for general purposes than are involved in the pet trade specifically.Finally, we obtained data for a proxy measure of species detectability in order to control for a potential confounding effect on relationships between life history traits and introduction: larger bodied and longer-lived species may be more likely to be recorded by human observers when introduced compared with smaller and shorter-lived species. We obtained data on species occurrence records, geographic range size and population density, assuming that highly detectable species will have a disproportionately large number of recorded observations than expected based on the size of their geographic ranges and average population densities, following similar approaches by e.g. Refs. 65,66. We obtained occurrence records from the Global Biodiversity Information Facility (GBIF33) via the R package rgbif67 selecting only records resulting from human observation. We obtained range sizes (in decimal degrees squared) from the IUCN Red List32 and processed them for analysis using functions from the rgdal package68, excluding areas of uncertain presence (i.e. limiting range to presence code 1, ‘extant’). We obtained population density estimates from the TetraDENSITY database (version 134), a global database of population density estimates for terrestrial vertebrates. Most species in the TetraDENSITY dataset are represented by estimates from multiple different studies (median = 3, range 1–408). We collapsed density estimates to the species level by taking the median value across studies, including all estimates regardless of sampling method to maximise sample size, and converting all units to individuals/km2 to ensure comparability.Statistical analysesTo investigate relationships between life history traits and trade, we run models treating US or pet trade as the outcome variable and life history traits as the predictors. For all analyses, all life history variables were included in the same models to account for covariation among life history traits55. For US trade, where data on trade frequency are available, we run models both in which trade is treated as a binary variable (traded vs. not traded) and as a count variable (frequency of live shipments, including zero values), while for the pet trade, we have no data on trade frequency and so we treat pet trade as a binary variable only. To investigate the effects of life history traits on introduction, we run models in which introduction is the outcome variable and life history traits are the predictors. In introduction models, we only include traded species (running separate models for the set of species in US trade and the set of species in the pet trade). This approach allows us to disentangle effects associated with trade and introduction and thus identify at which stage(s) life history biases emerge. We also run introduction models in which frequency of US trade is included as an additional predictor alongside life history traits, anticipating that highly traded species are more likely to be introduced. Finally, to investigate possible confounding effects of species detectability on relationships between life history traits and introduction, we investigate effects of number of observations, geographic range size and, where sample sizes allowed, population density on the probability of introduction. If highly detectable species are more likely to be recorded as introduced, we expect to find a positive effect of the number of observations (while accounting for geographic range size and population density) on the probability of introduction. If this effect confounds relationships between body mass/lifespan and introduction, the effect of these life history traits on the probability of introduction should disappear when detectability measures are included in the models alongside life history traits. All analyses were conducted using the R statistical programming environment (Version 4.2.069). Plots were coloured using palettes from the viridis package70.To estimate effects of predictor variables, we fit generalized linear mixed models (GLMMs) using Markov chain Monte-Carlo (MCMC) estimation, implemented in the MCMCglmm package35,36. For analyses with binary outcome variables (traded vs. not traded, introduced vs. not introduced) we fit probit models, while for analyses with US trade frequency as the outcome variable we fit hurdle models. Hurdle models estimate two latent variables: the probability that the outcome is zero (on the logit scale), and the probability of the outcome modelled as a Poisson distribution for non-zero values71. This method therefore allows us to estimate effects of life history traits on the probability and frequency of trade in the same model. While the binary component of a hurdle model estimates the probability that outcomes are zero, when reporting results we reverse the sign of coefficients from the binary model for ease of interpretation, so that effects can be interpreted as the probability that the outcome is not zero. Therefore, here predictors with consistent effects on the probability and frequency of trade in hurdle models will have the same sign (so that if, for example, litter size has a positive effect on both the probability and frequency of trade, both coefficients for litter size from the hurdle model will be positive).Datasets comprising biological measures from multiple related species violate the fundamental statistical assumption that observations are independent of one another, since closely related species are more phenotypically similar than expected by chance due to their shared evolutionary history72. To account for the non-independence of species due to shared ancestry, we included a phylogenetic random effect in all models, represented by a variance-covariance (VCV) matrix derived from the phylogeny. The off-diagonal elements of the VCV matrix contain the amount of shared evolutionary history for each pair of species35,37,38 based on the branch lengths of the phylogeny (here proportional to time)59,61,62,63,63. This approach allows us to estimate phylogenetic signal using the heritability (H2) parameter, which measures the proportion of total variance in the latent variable attributable to the phylogeny35,37,38. Heritability is interpreted in the same way as Pagel’s λ in phylogenetic generalized least squares regression35,37,38,72. Specifically, phylogenetic signal is constrained between 0, indicating no phylogenetic effect so that species can be treated as independent, and 1, indicating that similarity between species is directly proportional to their amount of shared evolutionary history35,38,72. As hurdle models estimate two latent variables, for each hurdle model we report two heritability estimates, one for the binary and one for the Poisson component. All continuous independent variables were log-10 transformed due to positively skewed distributions. Although GLMMs do not require normally distributed predictor variables, log-transforming positively skewed life history predictors in phylogenetic comparative analyses allows us to model life history evolution on proportional rather than absolute scales. This is important as it facilitates biologically meaningful comparisons between species across large scales of life history variation73. Further, log-transforming positively skewed predictors helps to meet assumptions of the underlying Brownian motion model of evolutionary change, which assumes that phenotypic change along branches of the phylogeny is normally distributed74.We calculated variance inflation factors (VIFs) using functions from the car R package75 to check for multicollinearity between predictor variables. Where any model reported a variance inflation factor of 5 or above, indicating potentially problematic levels of collinearity76, we re-ran the model removing the variable with the highest VIF iteratively until all the remaining variables had VIFs of More

The Campbell et al. model provides an excellent framework to identify species whose extinction leads to community collapse and species whose reintroduction can restore the community (see Fig. 2 for an illustration of these processes). Our first objective, finding the effect of species extinction on the rest of the species in an established community, is achievable using the concept of Logical Domain of Influence (LDOI)41; the LDOI represents the influence of a (set of) fixed node state(s) on the rest of the components in a system. In this section we first present our proposed method to calculate the LDOI for the Boolean threshold functions governing the Campbell et al. model of plant–pollinator community assembly. Then we verify that the simplified logical functions preserve the LDOI and hence can be implemented to further analyze the effect of extinction in plant–pollinator networks. Next, we address one of the main questions that motivated this study: Can stable motif driver set analysis facilitate the identification of keystone species? We discuss the identification of the driver sets of inactive stable motifs and motif groups and present the results of stabilizing these sets to measure the magnitude of the effect of species extinction on the communities. Lastly we discuss possible prevention and mitigation measures based on the knowledge acquired from driver sets of stable motifs and motif groups.Figure 2Illustration of species extinction and restoration in a hypothetical 6-species community. (a) The interaction network (on the left), and the maximal richness community possible for this network (the community with the most established species). Nodes highlighted with green represent established species. (b) The initial extinction of two species, po_1 and po_2 (left) and the community that results after cascading extinctions (right). Nodes highlighted with grey represent extinct species. (c) An intervention to restore pl_2 (left), which induces the restoration of further species, finally leading to a restored community with all the species present (right). The nodes highlighted with teal represent the restored species.Full size imageLDOI in the Boolean threshold modelThe LDOI concept was originally defined on Boolean functions expressed in a disjunctive prime form. Here we extend it to Boolean threshold functions. We implemented it as a breadth first search on the interaction network, as exemplified in Fig. 3. Assume that we want to find the LDOI of a (set of) node(s) (S_0={n_1,dots ,n_N}) and their specific fixed state (Q(S_0)={sigma _{n_1},dots ,sigma _{n_N}}). Starting from the set (S_0), the next set of nodes (S_1) that can acquire a fixed state due to the influence of (Q(S_0)) consists of the nodes that have an incoming edge from the nodes in the set (S_0) in the interaction network. The nodes in set (S_1) are the subject of the first search level. For each node (n_i in S_0) and (n^prime _i in S_1) we assume a “worst case scenario” (i.e., maximal opposition of the effect of (n_i) on (n^prime _i) from other regulators) to find the possible sufficiency relationships between the two. There are five cases:

1.

If (n_i) is a positive regulator of (n^prime _i), then (sigma _{n_i}=1) is a candidate for being sufficient for (sigma _{n^prime _i}=1). We assume that all other positive regulators of (n^prime _i) that have an unknown state (i.e., are not in (Q(S_0))) are inactive and all negative regulators of (n^prime _i) that have an unknown state are active. If (sum _j W_{ij} > 0) under this assumption, then the active state of (n_i) is sufficient to activate (n^prime _i). The virtual node (n^prime _i) that corresponds to (sigma _{n^prime _i}=1) is added to LDOI((Q(S_0))).

2.

If (n_i) is a positive regulator of (n^prime _i), then (sigma _{n_i}=0) is a candidate for being sufficient for (sigma _{n^prime _i}=0). We assume all other positive regulators of (n^prime _i) that have an unknown state are active and all negative regulators of (n^prime _i) that have an unknown state are inactive. If (sum _j W_{ij}le 0) under this assumption, then the inactive state of (n_i) is sufficient to deactivate (n^prime _i). The virtual node (sim n^prime _i) that corresponds to (sigma _{n^prime _i}=0) is added to LDOI((Q(S_0))).

3.

If (n_i) is a negative regulator of (n^prime _i), then (sigma _{n_i}=1) is a candidate for being sufficient for (sigma _{n^prime _i}=0). We assume all positive regulators of (n^prime _i) that have an unknown state are active and all other negative regulators of (n^prime _i) that that have an unknown state are inactive. If (sum _j W_{ij}le 0) under this assumption, then the active state of (n_i) is sufficient to deactivate (n^prime _i). The virtual node (sim n^prime _i) that corresponds to (sigma _{n^prime _i}=0) is added to LDOI((Q(S_0))).

4.

If (n_i) is a negative regulator of (n^prime _i), then (sigma _{n_i}=0) is a candidate for being sufficient for (sigma _{n^prime _i}=1). We assume all positive regulators of (n^prime _i) that have an unknown state are inactive and all other negative regulators of (n^prime _i) that that have an unknown state are active. If (sum _j W_{ij} > 0) under this assumption, then the inactive state of (n_i) is sufficient to activate (n^prime _i). The virtual node (n^prime _i) that corresponds to (sigma _{n^prime _i}=1) is added to the LDOI((Q(S_0))).

5.

If none of the past four sufficiency checks are satisfied, the node (n^prime _i) will be visited again in the next search levels.

The second set of nodes that can be influenced, (S_2), are the nodes that have an incoming edge from the nodes in the set (S_1). The algorithm goes over these nodes in the second search level as described above. This search continues to all the levels of the search algorithm until all nodes are visited (possibly multiple times) and either acquire a fixed state and are added to the LDOI or their state will be left undetermined at the end of the algorithm. In Fig. 3, we illustrate this search to find the LDOI((sim )pl_1). The first search level is (S_1={)po_1, po_3(}); (sim )pl_1 is sufficient to deactivate po_3, but not po_1. As a result, (sim )po_3(in ) LDOI((sim )pl_1). This process continues until all levels are visited and at the end of the algorithm LDOI((sim )pl_1()={sim )po_3, (sim )pl_2, (sim )pl_3, (sim )pl_4, (sim )pl_5, (sim )po_1, (sim )po_2 (}).Figure 3Breadth first search of the interaction network to find the LDOI of a (set of) fixed note state(s) in Boolean threshold functions governing the dynamics of plant–pollinator networks. (a) An interaction network with five plants and 3 pollinators. (b) The breadth first search in the case of starting from the node state (sim )pl_1. The nodes with incoming edges from pl_1 make up (S_1={)po_1, po_3(}). The second sufficiency check is satisfied for node state (sim )po_3, as a result (sim )po_3(in ) LDOI((sim )pl_1). The same process is applied for node po_1, but none of the sufficiency checks are satisfied, so this node will be visited again later. The next level of the search consists of the nodes that have incident edges from (S_1), i.e., (S_2={)pl_2, pl_3, pl_4, pl_5(}). The second sufficiency check is satisfied for all of these nodes and they are all fixed to their inactive state in the LDOI((sim )pl_1). Lastly, we reach (S_3={)po_1, po_2(}). Node po_1 is reached again, and with both its positive regulators fixed to their inactive states the second sufficiency check is satisfied and node po_1 is fixed to its inactive state as well. The same holds for po_2 and hence LDOI((sim )pl_1()={sim )po_3, (sim )pl_2, (sim )pl_3, (sim )pl_4, (sim )pl_5, (sim )po_1, (sim )po_2 (}).Full size imageTo measure the accuracy of the simplification method originally introduced in28, we analyzed logical domains of influence in 6000 networks with 50–70 nodes. These networks are among the largest in our ensembles and have the most complex structures. We randomly selected (sets of) inactive node states, found their LDOIs using the Boolean threshold functions and the simplified Boolean functions, and compared the two resulting LDOIs. We used 8 single node states and 8 combinations of size 2 to 4 for each network. We found that in all cases the LDOI calculated using the simplified Boolean functions matches the LDOI calculated using the Boolean threshold functions.Next, we analyzed (sets of) active node states and their LDOIs in the same ensembles of networks. Similar to the previous analysis, we used 8 single node states and 8 combinations of size 2 to 4 for each network. Our analysis shows that in 77.1% of the cases the LDOI calculated using the simplified Boolean functions matches the LDOI calculated using the Boolean threshold functions. In 22% of the cases the LDOI calculated from the simplified Boolean functions contains the LDOI calculated from the threshold functions, and it also contains extra active node states, overestimating the LDOI by 57.5% on average. These additional members of the LDOI result from the fact that the simplified Boolean functions contain fewer negative regulators than the threshold functions. The guiding principle of the simplification method is that the probability of (H(x)=1) conserves the probability of each node having an active state across all the states it can have. In contrast, the probability of the propagation of the active state is not necessarily preserved and tends to be higher in the simplified Boolean model; thus the LDOI of the active node states is overestimated in some cases.In the rest of the cases (about 1%), the LDOI calculated from the simplified Boolean functions does not fully capture the LDOI calculated from the threshold functions. This again is caused by the sparsification of the negative edges in the simplified Boolean functions. In the threshold functions, the activation of 4 or more negative regulators of a target node combined with one active positive regulator is sufficient to deactivate the target node, i.e., there might be inactive node states in the LDOI of a set of active node states. However, some of these negative regulators drop in the simplified Boolean model and the inactive state of the target node is not necessarily in the LDOI of the set of active node states in the simplified case. This is the rare mechanism by which the simplified model might underestimate the influence of active node states on the rest of the network.In the following section we are interested in analyzing the effect of species extinction on the established community, i.e., we look at the LDOI of (sets of) inactive node states. Observing that the influence of extinction of species is measured correctly in the simplified Boolean models, we conclude that these models can be utilized to further analyze the process of extinction and its ecological implications.Stable motif based identification of species whose loss leads to cascading extinctionsEach stable motif or motif group can have multiple driver sets; stabilization of each driver set leads to the stabilization of the whole motif or motif group. In plant–pollinator interaction networks, the stable motifs either represent a sub-community (when the constituent nodes stabilize in their active states) or the simultaneous extinction of all species in the group (when the constituent nodes stabilize to their inactive states). Stabilization of the nodes in the driver set of an inactive stable motif results in stabilization of all the nodes in the stable motif to their inactive state, i.e., cascading extinction of the constituent species.The knowledge gained from stable motif analysis and the network of functional relationships offers insight into the cascading effect of an extinction that constitutes a driver set of an inactive stable motif. The magnitude of this effect depends on (i) the number of nodes that the inactive stable motif contains and (ii) the number of virtual nodes (including motifs and motif groups) corresponding to inactive species that are logically determined by the stabilization of the inactive stable motif.To investigate the role of stable motifs in the study of species extinction in plant–pollinator networks, we simulated extinctions that drive inactive stable motifs in 6000 networks with the sizes of 50–70 nodes. We considered driver sets of size 1, 2, or 3, and implemented them by fixing the corresponding node(s) to its (their) inactive state. As a point of comparison, we also performed a “control” analysis using the same networks with the same size of initial extinction; however, the candidates of initial extinction are inactive node states that do not drive stable motifs or motif groups. Based on the properties of the drivers of stable motifs, one expects that following the extinction of driver species, cascading extinctions of other species follow, while the same does not necessarily hold for non-driver species. As a result, we expect to observe greater damage to the original community when driver species become extinct.We assume that the “maximal richness community”—the community (attractor) in which the largest number of species managed to establish—is the subject of species extinction. This maximal richness community results from the stabilization of all active stable motifs. All other attractors that have some established species contain a subset of all active stable motifs and thus will contain a subset of the species of the maximal richness community. While for a generic Boolean model with multiple attractors one expects that a perturbed version of the model also has multiple attractors, this specific perturbation of a plant–pollinator model (namely, extinction of species in the maximal richness community) has a single attractor. We prove this by contradiction. Assume there are two separate attractors in the perturbed model, which means that there is at least one node that has opposite states in these two attractors. Note that this bi-stability is the result of the perturbation and not a property of the original system as the maximal richness community (an attractor) is the starting point for the introduced extinction. Specifically, the inactive state of the extinct node has to lead to the stabilization of another node to its active versus inactive states in the two separate attractors. The only case in which the stabilization of an inactive node state can result in the stabilization of an active node state is if there is a negative edge from the former to the latter in the interaction network after simplification. Since the Boolean function in 2 is inhibitor dominant, the negative regulators that remain in the Boolean model must be in their inactive states in the maximal richness attractor. As they are already inactive (extinct), they are not candidates for extinction. The only nodes that are candidates for extinction are the ones that positively regulate other nodes; perturbing the system by fixing these candidates to their inactive states cannot lead to the active state of a target node. In conclusion, bi-stability is not possible.We found the new attractor of the system given the (combination of) inactive node state(s) using the the functions percolate_and_remove_constants() and trap_spaces() from the pyboolnet Python package. We quantify the effects of the initial extinction(s) on the maximal richness attractor by the percentage change in the number of active species, which we call damage percentage. Note that this choice of maximal richness community as the reference and starting point allows us to detect the cascading extinctions following the initial damage.In Fig. 4 the left column plots show the average damage percentage caused by the extinction of 1 (top panel), 2 (middle panel), or 3 (bottom panel) species that represent driver sets of stable motifs and motif groups, while the right column plots illustrate the average damage percentage as a result of the extinction of 1, 2 or 3 species that represent non-driver nodes. Comparing the two columns, one can notice that stabilization of the driver sets of stable motifs and motif groups leads to considerably larger damage to the communities. This is due to the fact that stabilization of driver sets ensures the stabilization of entire inactive stable motifs and motif groups and hence ensures cascading extinctions. Comparing the plots in the left column, we see that the larger the driver sets are, the larger the damage to the community becomes. This is because larger driver sets are more likely to stabilize larger stable motifs and motif groups. This figure illustrates the significance of stable motifs and their driver sets in the study of species extinction in plant–pollinator communities.Figure 4Histogram plots illustrating the average percentage of the damage caused in an established community after the extinction of species. This analysis is performed over 6000 networks with the size of 50–70 nodes. To study the extinction of species we started from the maximal richness community, then we fixed the nodes that correspond to the focal species to the their inactive states. The original extinctions are excluded from the damage percentages. The left column plots show the average damage percentage caused to the maximal richness community by the extinction of a driver set of size 1 (top), 2 (middle), or 3 (bottom) of an inactive stable motif or motif group. For each network, we determined all the relevant driver sets of one stable motif or motif group, we performed the extinction and calculated the resulting damage, then we calculated the average damage percentage over all data points collected for the same network. The right column plots show the average damage percentage caused to the maximal richness community by the extinction of 1 (top), 2 (middle), and 3 (bottom) non-driver, randomly chosen nodes. Each time a randomly selected combination of non-driver nodes were the subject of simultaneous extinction until all combinations are explored and then we calculated the average damage percentage over all data points collected for each network. The number of networks that qualify for each of these 6 categories differ (e.g., some networks have a stable motif with a driver set of size 2 but no stable motif with a driver set of size 3). In the left column 5529, 3212, and 1980 networks and in the right column 5779, 5626, and 5423 networks qualified respectively. The red lines represent the mean value of all the presented data points in each plot.Full size imageIn Fig. 4 left column, the full driver set of one inactive stable motif or motif group was stabilized. However, the species that become extinct might only contain a subset of a driver set of a stable motif or motif group, i.e., they only stabilize a subset of the inactive node states in the stable motif or motif group. We compare the extinction effect caused by the stabilization of a full driver set of four nodes with the effect of the extinction of four nodes that contain a partial driver set in Fig. 5 using the batch of the largest networks in this study, i.e, the batch that contains networks with 30 nodes representing plant species and 40 nodes representing pollinator species. This choice is due to the fact that the existence of stable motifs and motif groups having a driver set of four node states is highly probable in larger networks. As expected, the stabilization of the complete driver set leads to greater damage. Stabilization of the same number of nodes that contain a partial driver set leads to significantly less damage and species loss in the community; the median damage percentage in the case of stabilization of partial driver sets is 22.6% while it is 69.2% in the case of stabilization of the full driver sets. We also note that damage of more than 90% occurs rarely and is only possible when a full driver set is stabilized (see Fig. 5 right plot). This suggests that the motif groups that lead to total extinction tend to have a driver set with more than four nodes; in other words, only the simultaneous extinction of five or more species would lead to total community collapse.Figure 5Histogram plots illustrating the average percentage of the damage caused in an established community after the extinction of species. This analysis is performed over 1000 networks with the size of 70 nodes (30 nodes representing plant species and 40 nodes representing pollinator species). The original extinctions are excluded from the damage percentages. The left plot shows the average damage percentage caused to the maximal richness community by the extinction of 2 species that are a subset of the 4-node driver set of an inactive stable motif or motif group plus 2 randomly selected non-driver species. The right plot shows the damage percentage caused to the maximal richness community by the extinction of 4-node driver sets of the same inactive stable motifs and motif groups. Each time the driver set of one stable motif or motif group was the subject of extinction and we calculated the average damage percentage over all data points collected for each network. 295 networks qualified for this analysis.Full size imageMotif driver set analysis outperforms structural measures in identifying keystone speciesThe literature on ecological networks offers multiple measures that reflect the importance of each species for community stability. One family of such measures is centrality (quantified by the network measures degree centrality and betweenness centrality). Previous studies45,46 have shown that species (nodes) with higher centrality scores are keystone species in ecological communities (i.e., species whose loss would dramatically change or even destroy the community). The nodes with highest in-degree centrality (such as pl_2 in Fig. 6a) represent generalist species that can receive beneficial interactions from multiple sources and survive. The nodes with highest betweenness centrality (such as pl_2 and po_2 in Fig. 6a) represent species that act as connectors and help the community survive. We find that high centrality corresponds to specific patterns in the expanded network: the inactive state of generalist or connector species is often the driver of a cascading extinction. Indeed, stable motif analysis of the expanded network in Fig. 6b confirms that there is an inactive stable motif (highlighted with grey) driven by the minimal set {(sim )pl_2}. The fact that node pl_2 is a stable motif driver means that in the case of the extinction of pl_2 the whole community collapses.To compare the effectiveness of stable motif analysis to the effectiveness of the more studied structural measures to identify keystone species, we performed an analysis similar to the previous section. We compared the magnitude of cascading extinctions in the case of extinction of stable motif driver nodes and of nodes with high values of previously introduced structural importance measures. Specifically, we used node betweenness centrality, node contribution to nestedness47, and mutualistic species rank (MusRank)22 to find crucial species based on their structural properties. For more details on definition and adaptation of these two measures see “Methods”. In this analysis, we used each measure to target species in the simplified Boolean models as follows:

1.

Betweenness centrality: The 10% of species with the highest betweenness centrality are chosen to be candidates for extinction.

2.

Node contribution to nestedness: The species with the most interactions tend to contribute the least to the community nestedness. Targeting them most likely leads to a faster community collapse48. As a result, 10% of species with the lowest contribution to network nestedness are chosen to be candidates for extinction. For more details on this measure, please see “Methods”.

3.

Pollinator MusRank: The pollinator species with the highest MusRank importance are more likely to interact with multiple plants, so the 10% of pollinator species with the highest importance are chosen to be candidates for extinction. For more details on this measure, please see “Methods”.

4.

Plant MusRank: The plant species with the highest MusRank importance are more likely to interact with multiple pollinators, so the 10% of plant species with the highest importance are chosen to be candidates for extinction.

Figure 7 illustrates the results of this analysis in 6000 networks with 50–70 nodes. In each network the 1-node, 2-node, and 3-node driver sets of inactive stable motifs are identified and made extinct. In the same networks 10% of nodes based on betweenness centrality, node contribution to nestedness, and node MusRank score were chosen to be candidates for extinction. To match the “driver set” data, all choices of 1, 2, or 3 nodes in these sets were explored and the damage was averaged over each extinction size for each network. We observe the cascading extinction and calculate the damage percentage relative to the maximal richness attractor. The plot represents the collective data over all initial simultaneous extinction sizes of 1, 2, and 3 species.Comparing the four methods, one notices that the histograms acquired using stable motif driver sets, node betweenness centrality, and node contribution to nestedness are very similar, showing a peak for the 10–20% bin of the damage, and a long tail that reaches a damage percentage of 80–100%. The MusRank score performs less well in identifying the crucial species. Also, the frequency of the higher damage percentages shows that node contribution to nestedness is the closest to the “driver set” method in identifying nodes whose extinction causes the collapse of the whole community, making it the best structural measure out of the three. Nevertheless, the driver set method finds keystone species that cannot be identified via structural measures, as the corresponding damage percentage histogram has the most prominent tail at the right edge of the panel. Indeed, stable motif driver sets identified 82%, 80%, and 546% more species whose extinction leads to 60% or higher damage to the community when compared to betweenness centrality, node nestedness, and node MusRank score based methods respectively.The reason for the higher effectiveness of driver set analysis is illustrated in Fig. 8 in which the MusRank score and node contribution to nestedness are calculated for two example networks. One can see how these two measures might incorrectly identify less vital species. In the left column of Fig. 8, MusRank identifies the node po_2, highlighted with green, as the most important species. However, this node does not have any outgoing edges; its extinction does not lead to any cascading extinction. The inability of the MusRank score to consider the direction of edges causes such misidentification. In the right column, the three nodes highlighted with yellow have the lowest contributions to network nestedness. The expanded network shows that these three nodes together are not able to cause full community collapse, while the three-node driver set of the inactive stable motif can. Since the nestedness definition depends on the number of mutual interactions, it might fail to identify some of the keystone nodes that are necessary to the stability of the community (for more details on node nestedness see “Methods”).Previously it was shown that identifying the stable motifs and their driver sets can successfully steer the system toward a desired attractor or away from unwanted ones37,38,43. Stable motif analysis of the Boolean model offers insight into the dynamical trajectories of the system; hence control strategies can be developed accordingly. In the next section we use stable motif driver sets to suggest control methods and analyze their efficiency.Figure 6Generalist species in the interaction network and the expanded network. (a) A simplified network consisting of 3 plant and 3 pollinator species. pl_2 is a generalist species, i.e., it has two incoming edges indicating that it can survive on either of its sources of pollination, po_1 or po_2. The expanded network in (b) illustrates that the stabilization of the grey stable motif stabilizes all the nodes to their inactive states, and hence causes full community collapse. (sim )pl_2 is the minimal driver set of the grey stable motif, consistent with the strong damage induced by the loss of a generalist species.Full size imageFigure 7Histogram plots illustrating the performance of driver set analysis versus structural measures in identifying keystone species. The analysis was done on 6000 networks with sizes of 50–70 nodes. The starting point is the maximal richness community, i.e., the attractor in which the most species establish. For each network 1, 2, and 3 node(s) were selected and simultaneously fixed to their inactive states. After the cascading damage the new attractor is compared to the maximal richness attractor to calculate the damage percentage. The structural measures—betweenness centrality, node nestedness contribution, and node MusRank score—were calculated for all nodes in each network; the top 10% according to the relevant ordering were candidates to being fixed to their inactive states. The network IDs were matched, i.e., only the networks that had candidate nodes according to all four measures for each extinction size are included in this plot. The total number of data points is 6360. The red solid lines represent the mean and the black dashed lines represent the median over all data points in each plot.Full size imageFigure 8Networks illustrating examples of when structural measures fail to identify keystone species. In both columns simplified networks consisting of 3 plant and 3 pollinator species are presented. The MusRank is calculated for all the nodes in the network in the left column and denoted in the node labels. The expanded network corresponding to this network is shown below. Node contribution to network nestedness is calculated for all the nodes in the network in the right column and denoted in the node labels. Similarly the expanded network that correspond to it is shown below. Note that these two networks have different edges. In the left column MusRank score identifies node po_2, highlighted with green, as the most important, while the expanded network shows that the extinction of po_2 does not cause any further damage to the community, as this node has no outgoing edges. This is due to the fact that MusRank calculation process fails to consider the directed network and replaces all the directed edges with undirected ones. The MusRank score does not identify po_3 as a crucial species; however, virtual node (sim )po_3, outlined with black in the expanded network is a driver of a stable motif that has all other nodes in its LDOI; the extinction of po_3 leads to full community collapse. In the right column, the nodes highlighted with yellow (pl_2, pl_3, and po_2) have the lowest node contribution to nestedness, which predicts that these nodes are likely crucial to the stability of the community. Analyzing the expanded network, one can see that these three nodes together are not able to drive the inactive stable motif highlighted with teal. The minimal driver set for this stable motif, outlined with black, consists of {(sim )po_1, (sim )po_2, (sim )po_3}; together these nodes drive the inactive stable motif and cause full community collapse. The nestedness-based measure was not able to capture the significance of nodes po2 and po_3.Full size imageDamage mitigation measures and strategies for endangered communitiesThere are two substantial questions related to managing the damage induced by species extinction: (1) How can one prevent the damage as much as possible? (2) Once the damage happens, the reintroduction of which species can restore the community and to what extent? In this section we aim to answer these questions in the context of the Campbell et al. model, implementing stable motif based network control. This analysis can inform agricultural and ecological strategies employed to prevent and mitigate damage.Damage preventionOne of the most important questions in ecology is what strategies to use so that we can prevent and avert extinction damage to the community. In this section we analyze how the knowledge from stable motif analysis and driver sets can be implemented to minimize the effect of extinction of keystone species in case of limited resources. Each attractor of the original system can have multiple control sets; stabilizing the node states in each control set ensures that the system reaches that specific attractor. The same information from the attractor control sets can be implemented to prevent the system from converging into unwanted attractors. Zañudo et al. illustrated that by blocking (not allowing to stabilize) the stable motifs that lead to the unwanted attractors, one can decrease the probability (sometimes to zero) that the system arrives in those attractors38. In order to block an attractor, the control sets of that attractor are identified and the negations of the node states in the control sets are externally imposed. This approach eliminates the undesired attractor; however, new attractors might form that are similar to the eliminated attractor. Campbell et al. showed that in order to avoid such new attractors one needs to block the parent motif, which in this case is the largest strongly connected subgraph of the expanded network that contains the inactive virtual nodes44. Here, we investigate how stable motif blocking based attractor control can identify the species whose preservation would offer the highest benefit in avoiding catastrophic damage to the community. This information would aid the development of management strategies in plant–pollinator communities.To avoid all attractors that lead to some degree of species extinction, one needs to block all the driver sets of all inactive stable motifs and motif groups in a given network. Implementing this in 100 randomly selected networks with 25 plant and 25 pollinator nodes, we found that 45.6% of the species in the maximal richness community need to be kept (prevented from extinction) to ensure the lack of cascading extinctions. Given that management resources are usually limited, active monitoring and conservation of almost half of the species in a community seems costly and impractical. Hence, we set a more feasible goal of identifying and blocking the driver set(s) of the largest inactive stable motif or motif group in each network. The same 100 networks containing 50 nodes are the subject of analysis in this section. The reason for performing the analysis in a relatively limited ensemble is that it involves the identification of all driver sets of the largest inactive stable motif or motif group, which is computationally expensive. For each network, the driver set of the largest inactive stable motif or motif group (which corresponds to the extinction of all the species in that group) is identified and blocked (that is, the corresponding species are not allowed to go extinct). Then the same number of species as in the driver set of that stable motif or motif group are selected and stabilized to their inactive state. We considered all combinations of node extinctions outside the blocked subset, calculated the damage percentage relative to the maximal richness community, and then averaged over all data points for each network. As a control, we repeated the analysis without blocking; the size of the initial extinction is the same as in the previous analysis for consistency.Figure 9 shows the result of the analysis described above for 100 networks. The left box and whiskers plot illustrates the damage percentage relative to the maximal richness community when the blocking feature is activated, while the right box and whiskers plot shows the damage percentage relative to the maximal richness community when the blocking is disabled. The average and median damage percentages are 14.96% and 13.04% respectively when the largest inactive stable motif or motif group was blocked and 24.73% and 20.38% when it was not. This (sim )10% difference in the average between the two sets of results, as well as the fewer cases of high-damage outliers in the left plot, demonstrates that by preventing the extinction of species identified by stable motif analysis, one can prevent catastrophic community damage considerably.To estimate the fraction of species that would need to be monitored to prevent their extinction, we compared the size of the maximal richness attractor and the size of the driver set of the largest stable motif. The maximal richness community represents an average of 32% of the original species pool, approximately 15 out of 50 species. The driver sets of the largest stable motifs had an average size of 2.5 node states over all 100 networks, i.e., about 16.6% of the maximal richness community. In ecological terms, given limited resources, the information gained from stable motif driver sets can help direct the conservation efforts toward the keystone species that play a key role in maintaining the rest of the community in a cost-effective manner.Figure 9Box plots comparing the damage communities face if the largest inactive stable motif or motif group is completely blocked, i.e., all the drivers of this inactive stable motif or motif group are prevented from stabilizing versus if the same stable motif or motif group is allowed to stabilize. This analysis was performed over 100 randomly selected networks that contain 25 plant and 25 pollinator nodes. All the driver sets of an inactive stable motif or motif group are identified. From left to right the box and whiskers plots show the average damage percentage relative to the maximal richness community if the largest inactive stable motif is blocked and the same quantity if the largest stable motif or motif group is not blocked respectively. For the left box and whiskers plot, all combinations of inactive node states except the driver sets are considered, and for the right box and whiskers plot all combinations are explored. Due to the computational complexity caused by combinatorial explosion, this analysis was performed over 100 randomly selected 50-node networks.Full size imageRestoration of a group of speciesAlthough human preservation efforts have been directed toward community conservation, there are many industrial activities that lead to ecosystem degradation. Ecologists are interested in developing restoration strategies to be deployed after a stable community is hit by catastrophic damage to recover biodiversity and the ecosystem functions it provides49. Here we propose that stable motif analysis and the driver sets identified from the expanded network can give insight into restoration measures. While we examined the inactive stable motifs in the study of species extinction, here we focus on the active stable motifs as our goal is to restore as much biodiversity as possible.Several network measures have been proposed to identify the species that if re-introduced would restore the community considerably. Two of the most studied algorithms include maximising functional complementarity (or diversity) and maximising functional redundancy50. The first strategy targets the restoration of the species that provide as many functions to the ecosystem as possible; this approach results in a community that has a maximal number of functions provided by different groups of species. Alternatively, maximising the functional redundancy yields a community in which several species perform the same function. While this resultant community might have a limited number of functions, it is robust. Both of these community restoration approaches have been studied extensively (e.g. see21).We hypothesize that restoring the species that constitute driver sets of active stable motifs can help maximise the number of species post-restoration. Since there is evidence that functional diversity correlates with the number of species in the community51, we compare the post-restoration communities identified by stable motif driving with the functional diversity maximisation approach. As discussed in section LDOI in the Boolean threshold model, the Boolean simplification of the threshold functions leads to an overestimation of the LDOI of active node states (compared to the original threshold functions) in some networks. We evaluate the negative effects of this overestimation by checking the effectiveness of the restored species in the original threshold model.The same 6000 networks we examined in the last section were the subject of this analysis. To create an unbiased initial community, we create the damaged communities by eliminating the same number of species from the maximal richness community as the number that will be restored. We identify the inactive stable motif or motif group with the driver set size of 1, 2, or 3 node states that causes the most damage to the maximal richness community. We then eliminate the species corresponding to this driver set to reach the most damaged community for the given size of the initial extinction. This community is the starting point for two analyses. In the stable motif driving approach we stabilized an active stable motif that has a driver set of the same size as the initial extinction to reach a post-restoration community and calculated the percentage of the extinct species that were restored. In the functional diversity maximization based approach we re-introduced the same number of species selected from the to 10% of species in terms of their contribution to functional diversity.To calculate the functional diversity of a community one needs to (1) define and construct a trait matrix, (2) determine the distance (trait dissimilarity) of pairs of species, (3) perform hierarchical clustering based on the distances to create a dendrogram, and (4) calculate the total branch length of the dendrogram, i.e., the sum of the length of all paths51,52. Petchey et al. argued that resource-use traits among plant and pollinator species can be used to classify the organisms into separate functional groups53 and Devoto et al. proposed the use of the adjacency matrix based on the interaction network as the trait matrix21. In this study we do the same and implement the bipartite adjacency matrix to construct the distance matrix.Since the networks of the Campbell et al. model are directed, we modify the algorithm in that we have two separate adjacency matrices, one denoting the edges incoming to plant species and the other denoting the edges incoming to pollinator species. The hierarchical clustering algorithm is then run on each of these matrices separately, resulting in a dendrogram for each adjacency matrix. If extinction occurs in a community, the functional diversity of the survived community can be determined by calculating the total branch length of the subset of the dendrogram that includes only the survived species. The restoration strategy using this method is to re-introduce the nodes whose branches add the most to the total branch length of this subset, i.e., maximise the functional diversity of the survived community54. For more details see “Methods”.In each network, the percentage of the extinct species that were restored was calculated and averaged over all data points for each restoration size and each network. Figure 10 illustrates the results of this investigation. Applied to the simplified Boolean model, the median restoration percentage in the case of active stable motif driver set method (blue plot) is 80%. The functional diversity maximization strategy to restoration (yellow plot) yields a lower median restoration percentage, 73%, as well as a large number of low-restoration outliers. Although one might argue that identifying beneficial species using the functional diversity maximization strategy works well, the higher percentage of the cases of 80–100% restoration in case of the active stable motif driver set analysis indicates that the latter identifies some of the most effective restorative species that are not identified via the former method. As in a minority of cases the simplified Boolean model overestimates the positive impact of the sustained presence of a species (see section LDOI in the Boolean threshold model), we sought to verify the effectiveness of the predicted restoration candidates in the original threshold model. The blue (respectively, yellow) box and whiskers plot on the right represents the restoration percentages of the same species as in the left blue (respectively, yellow) plot when these species are restored in the threshold model. The median of the right blue plot is 70%, while the median of the right yellow is 63%, preserving the advantage of the stable motif driver sets. We conclude that although the simplified Boolean model overestimates the restoration effectiveness of certain driver sets (visible in the fact that the lower whisker of the blue plot on the right goes well below the lower whisker of the blue plot on the left), stable motif driver sets are more effective in both comparisons.Figure 10Box and whiskers plots illustrating the average percentage of the extinct species that are restored following the stable motif driver set restoration strategy (blue) versus the functional diversity based approach (yellow). This analysis is performed over 6000 networks with sizes of 50–70 nodes. Starting from the maximal richness community, for each network one inactive stable motif with a driver set of 1, 2 or 3 nodes was stabilized to reach a new damaged community. This task was performed until the community with the most extinct species was identified. This is the community we set as the starting point for the restoration process using both methods. The pair on the left represents the two methods applied to the simplified Boolean model. For both methods we identified 1, 2, or 3 influential nodes for community restoration and we calculated the percentage of the extinct species that could be restored. The pair on the right represents restoring the same species identified by each method in the previous analysis in the original threshold model. In all analyses the community restoration percentage was averaged over all combinations of the same size, for each network and each method. The IDs of all networks are matched.Full size imageCommunity restoration via attractor controlAs illustrated in section “Restoration of a group of species”, stable motif analysis identifies promising and cost-effective group restoration strategies. In this section we aim to go further and identify interventions that can maximally restore a community. Previous stable motif based network control methods37,38,55 require a search for the smallest set of node states to control the system once the stable motif stabilization trajectories are identified. This smallest set may not contain a node from each stable motif in the sequence. In this work, however, we know that each stable motif or motif group needs to be controlled individually28 because the stabilization of none of the motifs results in the stabilization of another. As a result, the control set of each attractor is the same as the union of the driver sets of all members in the consistent combination corresponding to that attractor.In this section we examined this attractor control method by setting the communities with 70% or more of the species in the maximal richness community as the target, i.e., the attractors that have 70% of the species in the maximal richness community are assumed to be the desired attractors. We then recorded the size of the minimal control set needed to achieve each of these attractors. Note that stabilizing each of these control sets guarantees that the system reaches the corresponding attractor38.For this section, we analyzed 6000 networks that have 50–70 nodes. Figure 11 represents box-and-whiskers plots of the size of the minimal set of species that need to be restored, where the target community sizes are classified into three groups based on the percentage of the species relative to the maximal richness attractor. One can see that in half of the cases, the restoration of either 1 or 2 species manages to restore more than 70% of the maximal richness community. The largest set has 8 species that need to be restored; however, this data point is an outlier. As illustrated, driver set analysis and stable motif based attractor control can efficiently identify the species that play an influential restorative role and suggest management strategies that are effective at the scale of the whole community. To assess the impact of the LDOI inflation on this result, we used the restoration candidates identified by control sets of the attractors of the Boolean model in the threshold functions of a subset of networks. The results of comparing the restoration percentage is shown in Fig. 14. The first quartile, median and third quartile values are 78.26%, 86.6%, and 100% for the simplified Boolean models and 43.78%, 72.41%, and 85.71% for the threshold model.To further compare the results of restoration obtained from the two models we sorted the species in the order of their contribution to community restoration following a catastrophic damage. We randomly selected 100 of the largest (70-node) networks, which have the highest probability of a discrepancy between the threshold functions and the simplified Boolean model. In 72% of the cases the two rankings matched completely, and in the majority of the remaining cases only one species was misplaced in the simplified Boolean model-based ranking. To conclude, there is a significant advantage to the implementation of the simplified Boolean model and the drawback can be addressed by a follow-up checking on the original threshold functions.Figure 11The number of species that need to be restored to save 70% of more of the species in the maximal richness community. In this analysis 6000 networks with 50–70 nodes were the subject. For each networks all the attractors that have 70% or more of the species in the maximal richness attractor are identified and set to be the target attractors. The control set of these attractors are then classified into three groups based on the percentage as illustrated in the figure. From left to right, the box and whiskers represent the size of the control set of attractors that have 70–80%, 80–90%, and 90–100% of the species in the maximal richness attractor respectively.Full size image More

Good, S. P., Noone, D. & Bowen, G. Hydrologic connectivity constrains partitioning of global terrestrial water fluxes. Science 349, 175–177 (2015).Article 
ADS 
CAS 

Google Scholar 
Humphrey, V. et al. Sensitivity of atmospheric CO2 growth rate to observed changes in terrestrial water storage. Nature 560, 628–631 (2018).Article 
ADS 
CAS 

Google Scholar 
Humphrey, V. et al. Soil moisture–atmosphere feedback dominates land carbon uptake variability. Nature 592, 65–69 (2021).Article 
CAS 

Google Scholar 
Devaraju, N., Bala, G. & Nemani, R. Modelling the influence of land-use changes on biophysical and biochemical interactions at regional and global scales. Plant Cell Environ. 38, 1931–1946 (2015).Article 
CAS 

Google Scholar 
Chu, C. et al. Does climate directly influence NPP globally?. Glob. Change Biol. 22, 12–24 (2016).Article 
ADS 

Google Scholar 
Pan, S. et al. Impacts of climate variability and extremes on global net primary production in the first decade of the 21st century. J. Geogr. Sci. 25, 1027–1044 (2015).Article 

Google Scholar 
Musavi, T. et al. Stand age and species richness dampen interannual variation of ecosystem-level photosynthetic capacity. Nat. Ecol. Evol. 1, 48 (2017).Article 

Google Scholar 
Cheng, J. et al. Vegetation feedback causes delayed ecosystem response to East Asian Summer Monsoon Rainfall during the Holocene. Nat. Commun. 12, 1–9 (2021).ADS 

Google Scholar 
Yu, Y. et al. Observed positive vegetation-rainfall feedbacks in the Sahel dominated by a moisture recycling mechanism. Nat. Commun. 8, 1–9 (2017).Article 
ADS 

Google Scholar 
Betts, R. A., Cox, P. M., Lee, S. E. & Woodward, F. I. Contrasting physiological and structural vegetation feedbacks in climate change simulations. Nature 387, 796–799 (1997).Article 
ADS 
CAS 

Google Scholar 
Forzieri, G. et al. Increased control of vegetation on global terrestrial energy fluxes. Nat. Clim. Chang. 10, 356–362 (2020).Article 
ADS 

Google Scholar 
Ukkola, A. M. et al. Reduced streamflow in water-stressed climates consistent with CO2 effects on vegetation. Nat. Clim. Chang. 6, 75–78 (2016).Article 
ADS 

Google Scholar 
Steffen, W. et al. Trajectories of the earth system in the anthropocene. Proc. Natl. Acad. Sci. USA 115, 8252–8259 (2018).Article 
ADS 
CAS 

Google Scholar 
Bruijnzeel, L. A. Hydrological functions of tropical forests: Not seeing the soil for the trees?. Agric. Ecosyst. Environ. 104, 185–228 (2004).Article 

Google Scholar 
Bierkens, M. F. P. & van den Hurk, B. J. J. M. Groundwater convergence as a possible mechanism for multi-year persistence in rainfall. Geophys. Res. Lett. 34, 2402 (2007).Article 
ADS 

Google Scholar 
Idso, S. B. & Brazel, A. J. Rising atmospheric carbon dioxide concentrations may increase streamflow. Nature 312, 51–53 (1984).Article 
ADS 
CAS 

Google Scholar 
Betts, R. A. et al. Projected increase in continental runoff due to plant responses to increasing carbon dioxide. Nature 448, 1037–1041 (2007).Article 
ADS 
CAS 

Google Scholar 
Swann, A. L. S., Hoffman, F. M., Koven, C. D. & Randerson, J. T. Plant responses to increasing CO2 reduce estimates of climate impacts on drought severity. Proc. Natl. Acad. Sci. USA. 113, 10019–10024 (2016).Article 
ADS 
CAS 

Google Scholar 
Mankin, J. S., Seager, R., Smerdon, J. E., Cook, B. I. & Williams, A. P. Mid-latitude freshwater availability reduced by projected vegetation responses to climate change. Nat. Geosci. 12, 983–988 (2019).Article 
ADS 
CAS 

Google Scholar 
Frank, D. C. et al. Water-use efficiency and transpiration across European forests during the Anthropocene. Nat. Clim. Chang. 5, 579–583 (2015).Article 
ADS 
CAS 

Google Scholar 
Zhang, K. et al. Vegetation greening and climate change promote multidecadal rises of global land evapotranspiration. Sci. Rep. 5, 1–9 (2015).
Google Scholar 
Teuling, A. J., Seneviratne, S. I., Williams, C. & Troch, P. A. Observed timescales of evapotranspiration response to soil moisture. Geophys. Res. Lett. 33, 23 (2006).Article 

Google Scholar 
Teuling, A. J., Uijlenhoet, R., Hupert, F. & Troch, P. A. Impact of plant water uptake strategy on soil moisture and evapotranspiration dynamics during drydown. Geophys. Res. Lett. 33, 3401 (2006).Article 
ADS 

Google Scholar 
Vivoni, E. R. et al. Observed relation between evapotranspiration and soil moisture in the North American monsoon region. Geophys. Res. Lett. 35, 22 (2008).Article 

Google Scholar 
Dirmeyer, P. A., Jin, Y., Csingh, C. & Yan, C. Evolving land-atmosphere interactions over North America from CMIP5 simulations. J. Clim. 26, 7313–7327 (2013).Article 
ADS 

Google Scholar 
Dirmeyer, P. A. et al. Verification of land-atmosphere coupling in forecast models, reanalyses, and land surface models using flux site observations. J. Hydrometeorol. 19, 375–392 (2018).Article 
ADS 

Google Scholar 
Friedlingstein, P. et al. Positive feedback between future climate change and the carbon cycle. Geophys. Res. Lett. 28, 1543–1546 (2001).Article 
ADS 
CAS 

Google Scholar 
Arora, K. et al. Carbon-concentration and carbon-climate feedbacks in CMIP6 models and their comparison to CMIP5 models. Biogeosciences 17, 4173–4222 (2020).Article 
ADS 
CAS 

Google Scholar 
Song, X., Wang, D. Y., Li, F. & Zeng, X. D. Evaluating the performance of CMIP6 Earth system models in simulating global vegetation structure and distribution. Adv. Clim. Chang. Res. 12, 584–595 (2021).Article 

Google Scholar 
Levine, P. A., Randerson, J. T., Swenson, S. C. & Lawrence, D. M. Evaluating the strength of the land-atmosphere moisture feedback in Earth system models using satellite observations. Hydrol. Earth Syst. Sci. 20, 4837–4856 (2016).Article 
ADS 

Google Scholar 
Wei, N. et al. Evolution of uncertainty in terrestrial carbon storage in earth system models from CMIP5 to CMIP6. J. Clim. 35, 5483–5499 (2022).Article 
ADS 

Google Scholar 
Smith, N. G. et al. Toward a better integration of biological data from precipitation manipulation experiments into Earth system models. Rev. Geophys. 52, 412–434 (2014).Article 
ADS 

Google Scholar 
Yuan, K., Zhu, Q., Riley, W. J., Li, F. & Wu, H. Understanding and reducing the uncertainties of land surface energy flux partitioning within CMIP6 land models. Agric. For. Meteorol. 319, 108920 (2022).Article 
ADS 

Google Scholar 
Baker, J. C. A. et al. An assessment of land-atmosphere interactions over south america using satellites, reanalysis, and two global climate models. J. Hydrometeorol. 22, 905–922 (2021).Article 
ADS 

Google Scholar 
Mooley, D. A. & Parthasarathy, B. Fluctuations in All-India summer monsoon rainfall during 1871?1978. Clim. Change 6, 287–301 (1984).Article 
ADS 

Google Scholar 
Guhathakurta, P. & Rajeevan, M. Trends in the rainfall pattern over India. Int. J. Climatol. 28, 1453–1469 (2008).Article 

Google Scholar 
Sarkar, S. & Kafatos, M. Interannual variability of vegetation over the Indian sub-continent and its relation to the different meteorological parameters. Remote Sens. Environ. 90, 268–280 (2004).Article 
ADS 

Google Scholar 
Roy, P. S. et al. New vegetation type map of India prepared using satellite remote sensing: Comparison with global vegetation maps and utilities. Int. J. Appl. Earth Obs. Geoinf. 39, 142–159 (2015).ADS 

Google Scholar 
Koster, R. D. et al. Regions of strong coupling between soil moisture and precipitation. Science 305, 1138–1140 (2004).Article 
ADS 
CAS 

Google Scholar 
Paul, S. et al. Weakening of Indian summer monsoon rainfall due to changes in land use land cover. Sci. Rep. 6, 1–10 (2016).Article 

Google Scholar 
Pathak, A., Ghosh, S., Kumar, P. & Murtugudde, R. Role of oceanic and terrestrial atmospheric moisture sources in intraseasonal variability of indian summer monsoon rainfall. Sci. Rep. 7, 12729 (2017).Article 
ADS 

Google Scholar 
Pradhan, R., Singh, N. & Singh, R. P. Onset of summer monsoon in Northeast India is preceded by enhanced transpiration. Sci. Rep. 9, 1–11 (2019).Article 

Google Scholar 
Bonan, G. B. Forests and climate change: Forcings, feedbacks, and the climate benefits of forests. Science 320, 1444–1449 (2008).Article 
ADS 
CAS 

Google Scholar 
Le Quéré, C. et al. Global carbon budget 2017. Earth Syst. Sci. Data 10, 405–448 (2018).Article 
ADS 

Google Scholar 
Chen, C. et al. China and India lead in greening of the world through land-use management. Nat. Sustain. 2, 122–129 (2019).Article 

Google Scholar 
Green, J. K. et al. Regionally strong feedbacks between the atmosphere and terrestrial biosphere. Nat. Geosci. 10, 410–414 (2017).Article 
ADS 
CAS 

Google Scholar 
Pathak, A. et al. Role of oceanic and land moisture sources and transport in the seasonal and interannual variability of summer monsoon in India. J. Clim. 30, 1839–1859 (2017).Article 
ADS 

Google Scholar 
Myers, N., Mittermeler, R. A., Mittermeler, C. G., Da Fonseca, G. A. B. & Kent, J. Biodiversity hotspots for conservation priorities. Nature 403, 853–858 (2000).Article 
ADS 
CAS 

Google Scholar 
Venkateswarlu, B. & Prasad, J. V. N. Carrying capacity of Indian agriculture: issues related to rainfed agriculture. Curr. Sci. 102, 6 (2012).
Google Scholar 
Pai, D. S. et al. Development of a new high spatial resolution (0.25° × 0.25°) long period (1901–2010) daily gridded rainfall data set over India and its comparison with existing data sets over the region. Mausam 65, 1–18 (2014).Article 

Google Scholar 
Rodríguez-Fernández, N. J. et al. Long term global surface soil moisture fields using an SMOS-trained neural network applied to AMSR-E data. Remote Sens. 8, 959 (2016).Article 
ADS 

Google Scholar 
Martens, B. et al. GLEAM v3: Satellite-based land evaporation and root-zone soil moisture. Geosci. Model Dev. 10, 1903–1925 (2017).Article 
ADS 

Google Scholar 
Miralles, D. G. et al. Global land-surface evaporation estimated from satellite-based observations. Hydrol. Earth Syst. Sci. 15, 453–469 (2011).Article 
ADS 

Google Scholar 
Doelling, D. R. et al. Geostationary enhanced temporal interpolation for CERES flux products. J. Atmos. Ocean. Technol. 30, 1072–1090 (2013).Article 
ADS 

Google Scholar 
Doelling, D. R. et al. Advances in geostationary-derived longwave fluxes for the CERES synoptic (SYN1deg) product. J. Atmos. Ocean. Technol. 33, 503–521 (2016).Article 
ADS 

Google Scholar 
Running, S. W., Mu, Q. & Zhao, M. MOD17A2H MODIS/Terra Gross Primary Productivity 8-Day L4 Global 500m SIN Grid V006. NASA EOSDIS Land Processes DAAC. (2015). https://doi.org/10.5067/MODIS/MOD17A2H.006. Accessed 22 May 2021.Pathak, A., Ghosh, S. & Kumar, P. Precipitation recycling in the Indian subcontinent during summer monsoon. J. Hydrometeorol. 15, 2050 (2014).Article 
ADS 

Google Scholar 
Paul, S., Ghosh, S., Rajendran, K. & Murtugudde, R. Moisture supply from the western ghats forests to water deficit east coast of India. Geophys. Res. Lett. 45, 4337–4344 (2018).Article 
ADS 

Google Scholar 
Sebastian, D. E. et al. Multi-scale association between vegetation growth and climate in India: A wavelet analysis approach. Remote Sens. 11, 2073 (2019).Article 

Google Scholar 
Tabari, H. & Hosseinzadeh Talaee, P. Sensitivity of evapotranspiration to climatic change in different climates. Glob. Planet. Change 115, 16–23 (2014).Article 
ADS 

Google Scholar 
Roy, A., Das, S. K., Tripathi, A. K., Singh, N. U. & Barman, H. K. Biodiversity in North East India and their conservation. Progress. Agric. 15, 182 (2015).Article 

Google Scholar 
Schimel, D., Stephens, B. B. & Fisher, J. B. Effect of increasing CO2 on the terrestrial carbon cycle. Proc. Natl. Acad. Sci. USA. 112, 436–441 (2015).Article 
ADS 
CAS 

Google Scholar 
Verma, A., Chandel, V. & Ghosh, S. Climate drivers of the variations of vegetation productivity in India. Environ. Res. Lett. 17, 084023 (2022).Article 
ADS 

Google Scholar 
Dimri, A. P. et al. Western disturbances: A review. Rev. Geophys. 53, 225–246 (2015).Article 
ADS 

Google Scholar 
Joseph, J., Scheidegger, J. M., Jackson, C. R., Barik, B. & Ghosh, S. Is flood to drip irrigation a solution to groundwater depletion in the Indo-Gangetic plain?. Environ. Res. Lett. 17, 104002 (2022).Article 
ADS 

Google Scholar 
Sahastrabuddhe, R., Ghosh, S., Saha, A. & Murtugudde, R. A minimalistic seasonal prediction model for Indian monsoon based on spatial patterns of rainfall anomalies. Clim. Dyn. 52, 3661–3681 (2019).Article 

Google Scholar 
Liang, X., Lettenmaier, D. P., Wood, E. F. & Burges, S. J. A simple hydrologically based model of land surface water and energy fluxes for general circulation models. J. Geophys. Res. 99, 14415 (1994).Article 
ADS 

Google Scholar 
Friedl, M. A. & Sulla-Menashe, D. MCD12Q1 MODIS/Terra+Aqua Land Cover Type Yearly L3 Global 500m SIN Grid V006. NASA EOSDIS Land Processes DAAC. (2019). https://doi.org/10.5067/MODIS/MCD12Q1.006. Accessed 22 May 2021.Myneni, R., Knyazikhin, Y. & Park, T. MODIS/Terra Leaf Area Index/FPAR 8-Day L4 Global 500m SIN Grid V061. NASA EOSDIS Land Processes DAAC. (2021) https://doi.org/10.5067/MODIS/MOD15A2H.061. Accessed 22 May 2021.Schaaf, C. & Wang, Z. MCD43A3 MODIS/Terra+Aqua BRDF/Albedo Daily L3 Global – 500m V006. NASA EOSDIS Land Processes DAAC. https://doi.org/10.5067/MODIS/MCD43A3.006. (2015). https://www.umb.edu/spectralmass/terra_aqua_modis/v006. Accessed 22 May 2021.Didan, K., Barreto Munoz, A., Solano, R. & Huete, A. MODIS Vegetation Index User’s Guide (MOD13 Series).Liu, S.-J., Zhang, J.-H., Tian, G.-H. & Cai, D.-X. Detection Fractional Vegetation Cover Changes Using MODIS Data. in 2008 Congress on Image and Signal Processing 707–710 (IEEE, 2008). https://doi.org/10.1109/CISP.2008.46. More

The Global Assessment Report on Biodiversity and Ecosystem Services: Summary for Policy-Makers (IPBES, 2019)DeFries, R. & Nagendra, H. Ecosystem management as a wicked problem. Science 356, 265–270 (2017).Article 
CAS 

Google Scholar 
Turkelboom, F. et al. When we cannot have it all: ecosystem services trade-offs in the context of spatial planning. Ecosyst. Serv. 29, 566–578 (2018).Article 

Google Scholar 
Lee, H. & Lautenbach, S. A quantitative review of relationships between ecosystem services. Ecol. Indic. 66, 340–351 (2016).Article 

Google Scholar 
Bennett, E. M., Peterson, G. D. & Gordon, L. J. Understanding relationships among multiple ecosystem services. Ecol. Lett. 12, 1394–1404 (2009).Article 

Google Scholar 
Goldstein, J. H. et al. Integrating ecosystem-service tradeoffs into land-use decisions. Proc. Natl Acad. Sci. USA 109, 7565–7570 (2012).Article 
CAS 

Google Scholar 
Vallet, A., Locatelli, B. & Pramova, E. Ecosystem Services and Social Equity: Who Controls, Who Benefits and Who Loses? (CIFOR, 2020); https://doi.org/10.17528/cifor/007849Neyret, M. et al. Assessing the impact of grassland management on landscape multifunctionality. Ecosyst. Serv. 52, 101366 (2021).Linders, T. E. W. et al. Stakeholder priorities determine the impact of an alien tree invasion on ecosystem multifunctionality. People Nat. 3, 658–672 (2021).Article 

Google Scholar 
Herzig, A., Ausseil, A.-G. & Dymond, J. in Ecosystem Services in New Zealand—Conditions and Trends (ed. Dymond, J. R.) 511–523 (Manaaki Whenua Press, 2014).Chan, K. M. A., Shaw, M. R., Cameron, D. R., Underwood, E. C. & Daily, G. C. Conservation planning for ecosystem services. PLoS Biol. 4, e379 (2006).Article 

Google Scholar 
Pennington, D. N. et al. Cost-effective land use planning: optimizing land use and land management patterns to maximize social benefits. Ecol. Econ. 139, 75–90 (2017).Article 

Google Scholar 
Hölting, L. et al. Including stakeholders’ perspectives on ecosystem services in multifunctionality assessments. Ecosyst. People 16, 354–368 (2020).Article 

Google Scholar 
Plieninger, T. et al. Exploring futures of ecosystem services in cultural landscapes through participatory scenario development in the Swabian Alb, Germany. Ecol. Soc. 18, 39 (2013).Article 

Google Scholar 
Tasser, E., Schirpke, U., Zoderer, B. M. & Tappeiner, U. Towards an integrative assessment of land-use type values from the perspective of ecosystem services. Ecosyst. Serv. 42, 101082 (2020).Article 

Google Scholar 
Sayer, J. et al. Ten principles for a landscape approach to reconciling agriculture, conservation, and other competing land uses. Proc. Natl Acad. Sci. USA 110, 8349–8356 (2013).Article 
CAS 

Google Scholar 
Vallet, A. et al. Linking equity, power, and stakeholders: roles in relation to ecosystem services. Ecol. Soc. 24, 14 (2019).Article 

Google Scholar 
Allan, E. et al. Land use intensification alters ecosystem multifunctionality via loss of biodiversity and changes to functional composition. Ecol. Lett. 18, 834–843 (2015).Article 

Google Scholar 
Hector, A. & Bagchi, R. Biodiversity and ecosystem multifunctionality. Nature 448, 188–190 (2007).Article 
CAS 

Google Scholar 
Manning, P. et al. Redefining ecosystem multifunctionality. Nat. Ecol. Evol. 2, 427–436 (2018).Article 

Google Scholar 
Raudsepp-Hearne, C., Peterson, G. D. & Bennett, E. M. Ecosystem service bundles for analyzing tradeoffs in diverse landscapes. Proc. Natl Acad. Sci. USA 107, 5242–5247 (2010).Article 
CAS 

Google Scholar 
Daniel, T. C. et al. Contributions of cultural services to the ecosystem services agenda. Proc. Natl Acad. Sci. USA 109, 8812–8819 (2012).Article 
CAS 

Google Scholar 
Gunton, R. M. et al. Beyond ecosystem services: valuing the invaluable. Trends Ecol. Evol. 32, 249–257 (2017).Article 

Google Scholar 
Peter, S., Le Provost, G., Mehring, M., Müller, T. & Manning, P. Cultural worldviews consistently explain bundles of ecosystem service prioritisation across rural Germany. People Nat. 4, 218–230 (2022).Article 

Google Scholar 
Haines-Young, R. & Potschin, M. in Ecosystem Ecology (eds Raffaelli, D. G. & Frid, C. L. J.) 110–139 (Cambridge Univ. Press, 2010).Fischer, M. et al. Implementing large-scale and long-term functional biodiversity research: the Biodiversity Exploratories. Basic Appl. Ecol. 11, 473–485 (2010).Article 

Google Scholar 
Wilson, E. O. Half-Earth: Our Planet’s Fight for Life (Norton, 2017).Bastin, J.-F. et al. The global tree restoration potential. Science 365, 76–79 (2019).Article 
CAS 

Google Scholar 
Clapp, J. & Moseley, W. G. This food crisis is different: COVID-19 and the fragility of the neoliberal food security order. J. Peasant Stud. 47, 1393–1417 (2020).Article 

Google Scholar 
Kirwan, J. & Maye, D. Food security framings within the UK and the integration of local food systems. J. Rural Stud. 29, 91–100 (2013).Article 

Google Scholar 
Ellis, E. C. To conserve nature in the Anthropocene, half Earth is not nearly enough. One Earth 1, 163–167 (2019).Article 

Google Scholar 
Boetzl, F. A. et al. A multitaxa assessment of the effectiveness of agri-environmental schemes for biodiversity management. Proc. Natl Acad. Sci. USA 118, e2016038118 (2021).Tyllianakis, E. & Martin-Ortega, J. Agri-environmental schemes for biodiversity and environmental protection: how we are not yet ‘hitting the right keys’. Land Use Policy 109, 105620 (2021).Article 

Google Scholar 
Arroyo-Rodríguez, V. et al. Designing optimal human-modified landscapes for forest biodiversity conservation. Ecol. Lett. 23, 1404–1420 (2020).Article 

Google Scholar 
Gilroy, J. J. et al. Cheap carbon and biodiversity co-benefits from forest regeneration in a hotspot of endemism. Nat. Clim. Change 4, 503–507 (2014).Article 

Google Scholar 
Lindenmayer, D. B. et al. Avoiding bio-perversity from carbon sequestration solutions: avoiding bio-perversity in carbon markets. Conserv. Lett. 5, 28–36 (2012).Article 

Google Scholar 
Stoll-Kleemann, S. & O’Riordan, T. in The Encyclopedia of the Anthropocene Vol. 3 (eds DellaSala, D. A. & Goldstein, M. I.) 347–353 (Elsevier, 2018).Schaich, H., Bieling, C. & Plieninger, T. Linking ecosystem services with cultural landscape research. GAIA 19, 269–277 (2010).Article 

Google Scholar 
O’Connor, L. M. J. et al. Balancing conservation priorities for nature and for people in Europe. Science 372, 856–860 (2021).Article 

Google Scholar 
Büscher, B. et al. Half-Earth or Whole Earth? Radical ideas for conservation, and their implications. Oryx 51, 407–410 (2017).Article 

Google Scholar 
van der Plas, F. et al. Towards the development of general rules describing landscape heterogeneity–multifunctionality relationships. J. Appl. Ecol. 56, 168–179 (2019).Article 

Google Scholar 
Almeida, I., Rösch, C. & Saha, S. Converting monospecific into mixed forests: stakeholders’ views on ecosystem services in the Black Forest Region. Ecol. Soc. 26, 28 (2021).Meyer, M. A. & Früh-Müller, A. Patterns and drivers of recent agricultural land-use change in southern Germany. Land Use Policy 99, 104959 (2020).Article 

Google Scholar 
Kastner, T. et al. Global agricultural trade and land system sustainability: implications for ecosystem carbon storage, biodiversity, and human nutrition. One Earth 4, 1425–1443 (2021).Rasmussen, L. V. et al. Social–ecological outcomes of agricultural intensification. Nat. Sustain. 1, 275–282 (2018).Article 

Google Scholar 
Lindborg, R. et al. How spatial scale shapes the generation and management of multiple ecosystem services. Ecosphere 8, e01741 (2017).Article 

Google Scholar 
Duarte, G. T., Santos, P. M., Cornelissen, T. G., Ribeiro, M. C. & Paglia, A. P. The effects of landscape patterns on ecosystem services: meta-analyses of landscape services. Landsc. Ecol. 33, 1247–1257 (2018).Article 

Google Scholar 
Le Provost, G. et al. The supply of multiple ecosystem services requires biodiversity across spatial scales. Nat. Ecol. Evol. https://doi.org/10.1038/s41559-022-01918-5 (2022).Martin, D. A. et al. Land-use trajectories for sustainable land system transformations: identifying leverage points in a global biodiversity hotspot. Proc. Natl Acad. Sci. USA 119, e2107747119 (2022).Article 
CAS 

Google Scholar 
Seabloom, E. W., Borer, E. T. & Tilman, D. Grassland ecosystem recovery after soil disturbance depends on nutrient supply rate. Ecol. Lett. 23, 1756–1765 (2020).Article 

Google Scholar 
Messinger, J. & Winterbottom, B. African forest landscape restoration initiative (AFR100): restoring 100 million hectares of degraded and deforested land in Africa. Nat. Faune 30, 14–17 (2016).
Google Scholar 
Whittingham, M. J. The future of agri-environment schemes: biodiversity gains and ecosystem service delivery? J. Appl. Ecol. 48, 509–513 (2011).Article 

Google Scholar 
Le Clec’h, S. et al. Assessment of spatial variability of multiple ecosystem services in grasslands of different intensities. J. Environ. Manage. 251, 109372 (2019).Article 

Google Scholar 
Forschungsethische Grundsätze und Prüfverfahren in den Sozial‐ und Wirtschaftswissenschaften Output 9, Berufungsperiode 5 (German Data Forum, 2017).Strukturdaten Reutlingen—Statistisches Bundesamt (Bundeswahlleiter, 2020); https://www.bundeswahlleiter.de/europawahlen/2019/strukturdaten/bund-99/land-8/kreis-8415.htmlStrukturdaten Uckermark—Statistisches Bundesamt (Bundeswahlleiter, 2020); https://www.bundeswahlleiter.de/europawahlen/2019/strukturdaten/bund-99/land-12/kreis-12073.htmlStrukturdaten Unstrut-Hainich-Kreis—Statistisches Bundesamt (Bundeswahlleiter, 2020); https://www.bundeswahlleiter.de/europawahlen/2019/strukturdaten/bund-99/land-16/kreis-16064.htmlBlüthgen, N. et al. A quantitative index of land-use intensity in grasslands: integrating mowing, grazing and fertilization. Basic Appl. Ecol. 13, 207–220 (2012).Article 

Google Scholar 
Ostrowski, A., Lorenzen, K., Petzold, E. & Schindler, S. Land use intensity index (LUI) calculation tool of the Biodiversity Exploratories project for grassland survey data from three different regions in Germany since 2006, BEXIS 2 module. Zenodo https://doi.org/10.5281/zenodo.3865579 (2020).Schall, P. et al. The impact of even‐aged and uneven‐aged forest management on regional biodiversity of multiple taxa in European beech forests. J. Ecol. 55, 267–278 (2018).Statistisches Jahrbuch über Ernährung, Landwirtschaft und Forsten der Bundesrepublik Deutschland Vol. 63 (Bundesministerium für Ernährung und Landwirtschaft, 2019).Simons, N. K. & Weisser, W. W. Agricultural intensification without biodiversity loss is possible in grassland landscapes. Nat. Ecol. Evol. https://doi.org/10.1038/s41559-017-0227-2 (2017).Zinke, O. Heupreise steigen: Futter für die Bauern knapp und teuer. Agrarheute https://www.agrarheute.com/markt/futtermittel/heupreise-steigen-futter-fuer-bauern-knapp-teuer-571946 (2020).Bois de Chez Nous (Lignum, 2021); https://www.lignum.ch/files/images/Downloads_francais/Shop/20010_Bois_de_chez_nous.pdfGerman Timber Company—Internationaler Holzhandel (German Timber Company, 2021); https://www.germantimber.company/. Accessed 2021-11-24Holzeinschlag nach Holzartengruppen, Holzsorten, ausgewählten Besitzarten (Statistisches Bundesamt, 2022); https://www.destatis.de/DE/Themen/Branchen-Unternehmen/Landwirtschaft-Forstwirtschaft-Fischerei/Wald-Holz/Tabellen/holzeinschlag-deutschland.htmlJahresjagdstrecke Bundesrepublik Deutschland, 2019–2020 (Deutsche Jagdverband, 2020); https://www.jagdverband.de/sites/default/files/2021-01/2021-01_Infografik_Jahresjagdstrecke_Bundesrepublik_Deutschland_2019_2020.jpgHeinze, E. et al. Habitat use of large ungulates in northeastern Germany in relation to forest management. For. Ecol. Manage. 261, 288–296 (2011).Article 

Google Scholar 
Conant, R. T., Cerri, C. E. P., Osborne, B. B. & Paustian, K. Grassland management impacts on soil carbon stocks: a new synthesis. Ecol. Appl. 27, 662–668 (2017).Article 

Google Scholar 
Hermes, J., Albert, C. & von Haaren, C. Mapping and Assessing Local Recreation as a Cultural Ecosystem Service in Germany. UVP-Report https://doi.org/10.17442/uvp-report.034.08 (2020).Hermes, J., Albert, C. & von Haaren, C. Assessing the aesthetic quality of landscapes in Germany. Ecosyst. Serv. 31, 296–307 (2018).Article 

Google Scholar 
Ehrhart, S. & Schraml, U. Perception and evaluation of natural forest dynamics. Allg. Forst Jagdztg. 185, 166–183 (2014).
Google Scholar 
Villanueva-Rivera, L. J. & Pijanowski, B. C. soundecology: Soundscape ecology. R package version 1.3.3 (2018).Meyer, S., Wesche, K., Krause, B. & Leuschner, C. Dramatic losses of specialist arable plants in central Germany since the 1950s/60s—a cross-regional analysis. Divers. Distrib. 19, 1175–1187 (2013).Article 

Google Scholar 
Sasaki, K., Hotes, S., Kadoya, T., Yoshioka, A. & Wolters, V. Landscape associations of farmland bird diversity in Germany and Japan. Glob. Ecol. Conserv. 21, e00891 (2020).Article 

Google Scholar 
Peña, L., Casado-Arzuaga, I. & Onaindia, M. Mapping recreation supply and demand using an ecological and a social evaluation approach. Ecosyst. Serv. 13, 108–118 (2015).Article 

Google Scholar 
Schägner, J. P., Brander, L., Paracchini, M.-L., Hartje, V. & Maes, J. Mapping recreational ecosystem services and its values across Europe: a combination of GIS and meta-analysis. In European Association of Environmental and Resource Economists 22nd Annual Conference (2016).R Core Team. R: A Language and Environment for Statistical Computing v.4.2.1 (R Foundation for Statistical Computing, 2022).Rust Programming Language https://www.rust-lang.org/ v 1.44Le Provost, G. et al. Contrasting responses of above- and belowground diversity to multiple components of land-use intensity. Nat. Commun. 12, 3918 (2021).Article 

Google Scholar 
Gini, C. On the measurement of concentration and variability of characters (English translation from Italian by Fulvio de Santis in 2005). Metron 63, 1–38 (1914). More

Plant Ecology and Nature Conservation Group, Wageningen University, Wageningen, the NetherlandsFons van der Plas & Liesje MommerSystematic Botany and Functional Biodiversity, Life Science, Leipzig University, Leipzig, GermanyThomas Schröder-Georgi, Alexandra Weigelt, Kathryn Barry & Christian WirthGerman Centre for Integrative Biodiversity Research Halle-Jena-Leipzig, Leipzig, GermanyAlexandra Weigelt, Kathryn Barry, Adriana Alzate, Nico Eisenhauer, Anke Hildebrandt, Christiane Roscher & Christian WirthTerrestrial Ecology Research Group, School of Life Sciences Weihenstephan, Technical University of Munich, Munich, GermanySebastian Meyer & Wolfgang WeisserAquaculture and Fisheries Group, Wageningen University and Research Centre, Wageningen, the NetherlandsAdriana AlzateAgroécologie, AgroSup Dijon, Institut National de la Recherche Agronomique, Université de Bourgogne, Université de Bourgogne Franche-Comté, Dijon, FranceRomain L. BarnardEidgenössische Technische Hochschule Zürich, Zurich, SwitzerlandNina BuchmannDepartment of Experimental Plant Ecology, Institute for Water and Wetland Research, Radboud University Nijmegen, Nijmegen, the NetherlandsHans de KroonInstitute of Ecology and Evolution, University Jena, Jena, GermanyAnne Ebeling & Winfried VoigtInstitute of Biology, Leipzig University, Leipzig, GermanyNico EisenhauerHumboldt-Universität zu Berlin, Berlin, GermanyChristof EngelsInstitute of Plant Sciences, University of Bern, Bern, SwitzerlandMarkus FischerMax Planck Institute for Biogeochemistry, Jena, GermanyGerd Gleixner, Ernst-Detlef Schulze & Christian WirthHelmholtz Centre for Environmental Research, Leipzig, GermanyAnke HildebrandtFriedrich Schiller University Jena, Jena, GermanyAnke HildebrandtGeoecology, University of Tübingen, Tübingen, GermanyEva Koller-France & Yvonne OelmannInstitute of Geography and Geoecology, Karlsruhe Institute of Technology, Karlsruhe, GermanySophia Leimer & Wolfgang WilckeEcotron Européen de Montpellier, Centre National de la Recherche Scientifique, Montferrier-sur-Lez, FranceAlexandru MilcuCentre d’Ecologie Fonctionnelle et Evolutive, Unité Mixte de Recherche 5175 (Centre National de la Recherche Scientifique-Université de Montpellier-Université Paul-Valéry Montpellier-Ecole Pratique des Hautes Etudes), Montpellier, FranceAlexandru MilcuDepartment of Evolutionary Biology and Environmental Studies, University of Zurich, Zurich, SwitzerlandPascal A. NiklausUFZ, Helmholtz Centre for Environmental Research, Department Physiological Diversity, Leipzig, GermanyChristiane RoscherInstitute of Landscape Ecology, University of Münster, Münster, GermanyChristoph ScherberCentre for Biodiversity Monitoring, Zoological Research Museum Alexander Koenig, Bonn, GermanyChristoph ScherberGeobotany, Faculty of Biology, University of Freiburg, Freiburg, GermanyMichael Scherer-LorenzenCentre of Biodiversity and Sustainable Land Use, University of Göttingen, Göttingen, GermanyStefan ScheuJ.F. Blumenbach Institute of Zoology and Anthropology, Animal Ecology, University of Göttingen, Göttingen, GermanyStefan ScheuDepartment of Geography, University of Zurich, Zurich, SwitzerlandBernhard SchmidInstitute of Ecology, College of Urban and Environmental Sciences, Peking University, Beijing, ChinaBernhard SchmidLeuphana University Lüneburg, Institute of Ecology, Lüneburg, GermanyVicky TempertonAgroecology, Department of Crop Sciences, University of Göttingen, Göttingen, GermanyTeja TscharntkeF.v.d.P. wrote the initial draft of the manuscript. T.S.-G., A.W., K.B., S.M., A.A., R.L.B., N.B., H.d.K., A.E., N.E., C.E., M.F., G.G., A.H., E.K.-F., S.L., A.M., L.M., P.A.N., Y.O., C.R., C.S., M.S.-L., S.S., B.S., E.-D.S., V.T., T.T., W.V., W. Weisser, W. Wilcke and C.W. helped edit the manuscript. More

West-Eberhard, M. J. Developmental Plasticity and Evolution (Oxford University Press, 2003).Book 

Google Scholar 
Piersma, T. & Van Gils, J. A. The Flexible Phenotype: A Body-Centred Integration of Ecology, Physiology, and Behaviour (Oxford University Press, 2011).
Google Scholar 
Piersma, T. & Drent, J. Phenotypic flexibility and the evolution of organismal design. Trends Ecol. Evol. 18, 228–233 (2003).Article 

Google Scholar 
Tabari, H. Climate change impact on flood and extreme precipitation increases with water availability. Sci. Rep. 10, 1–10 (2020).
Google Scholar 
Beck, H. E. et al. Present and future Köppen–Geiger climate classification maps at 1-km resolution. Sci. Data 5, 180214 (2018).Article 

Google Scholar 
Bogert, C. M. Thermoregulation in reptiles, a factor in evolution. Evolution 3, 195–211 (1949).Article 
CAS 

Google Scholar 
Rensch, B. Das Prinzip geographischer Rassenkreise und das Problem der Artbildung (Gebrueder Borntraeger, 1929).
Google Scholar 
Clusella Trullas, S., van Wyk, J. H. & Spotila, J. R. Thermal melanism in ectotherms. J. Therm. Biol. 32, 235–245 (2007).Article 

Google Scholar 
Delhey, K. A review of Gloger’s rule, an ecogeographical rule of colour: Definitions, interpretations and evidence. Biol. Rev. 94, 1294–1316 (2019).
Google Scholar 
Stuart-Fox, D., Newton, E. & Clusella-Trullas, S. Thermal consequences of colour and near-infrared reflectance. Philos. Trans. R. Soc. B: Biol. Sci. 372, 20160345 (2017).Article 

Google Scholar 
Friedman, N. R. & Remês, V. Ecogeographical gradients in plumage coloration among Australasian songbird clades. Glob. Ecol. Biogeogr. 26, 261–274 (2017).Article 

Google Scholar 
Delhey, K. Darker where cold and wet: Australian birds follow their own version of Gloger’s rule. Ecography 41, 673–683 (2018).Article 

Google Scholar 
Galván, I., Rodríguez-Martínez, S. & Carrascal, L. M. Dark pigmentation limits thermal niche position in birds. Funct. Ecol. 32, 1531–1540 (2018).Article 

Google Scholar 
Medina, I. et al. Reflection of near-infrared light confers thermal protection in birds. Nat. Commun 9, 3610 (2018).Article 
ADS 

Google Scholar 
Aldrich, J. W. & James, F. C. Ecogeographic variation in the American Robin (Turdus migratorius). Auk 108, 230–249 (1991).
Google Scholar 
Morales, H. E. et al. Neutral and selective drivers of colour evolution in a widespread Australian passerine. J. Biogeogr. 44, 522–536 (2017).Article 

Google Scholar 
Griffith, S. C., Owens, I. P. & Burke, T. Environmental determination of a sexually selected trait. Nature 400, 358–360 (1999).Article 
ADS 
CAS 

Google Scholar 
Fargallo, J. A., Laaksonen, T., Korpimäki, E. & Wakamatsu, K. A melanin-based trait reflects environmental growth conditions of nestling male Eurasian kestrels. Evol. Ecol. 21, 157–171 (2007).Article 

Google Scholar 
Fargallo, J. A., Martínez, F., Wakamatsu, K., Serrano, D. & Blanco, G. Sex-dependent expression and fitness consequences of sunlight derived color phenotypes. Am. Nat. 191, 726–743 (2018).Article 

Google Scholar 
Beebe, W. Geographic variation in birds, with especial reference to the effects of humidity. Zoologica 1, 3–41 (1907).
Google Scholar 
Bieber, H. Fellverdunklung beim hauskaninchen nach kälteeinwirkung. Zeitschrift für Säugetierkunde 38, 33–38 (1972).
Google Scholar 
Johnston, R. F. & Selander, R. K. House sparrows: Rapid evolution of races in North America. Science 144, 548–550 (1964).Article 
ADS 
CAS 

Google Scholar 
Galván, I., Wakamatsu, K. & Alonso-Álvarez, C. Black bib size is associated with feather content of pheomelanin in male house sparrows. Pigment Cell Melanoma Res. 27, 1159–1161 (2014).Article 

Google Scholar 
Endler, J. A. On the measurement and classification of colour in studies of animal colour patterns. Biol. J. Linn. Soc. 41, 315–352 (1990).Article 

Google Scholar 
Montgomerie, R. Analyzing colors. In Bird Colouration I. Mechanisms and Measurements (eds Hill, E. G. & McGraw, K. J.) (Harvard University Press, 2006).
Google Scholar 
McGraw, K. J., Dale, J. & Mackillop, E. A. Social environment during molt and the expression of melanin-based plumage pigmentation in male house sparrows (Passer domesticus). Behav. Ecol. Sociobiol. 53, 116–122 (2003).Article 

Google Scholar 
Lessells, C. M. & Boag, P. T. Unrepeatable repeatabilities a common mistake. Auk 104, 116–121 (1987).Article 

Google Scholar 
Anderson, T. R. Biology of the Ubiquitous House Sparrow (Oxford University Press, 2006).Book 

Google Scholar 
Gelman, A. & Hill, J. Data Analysis Using Regression and Multilevel/Hierarchical Models (Cambridge University Press, 2006).Book 

Google Scholar 
Nakagawa, S., Ockendon, N., Gillespie, D. O., Hatchwell, B. J. & Burke, T. Assessing the function of house sparrows’ bib size using a flexible meta-analysis method. Behav. Ecol. 18, 831–840 (2007).Article 

Google Scholar 
Hill, G. E. & McGraw, K. J. Bird Coloration, Volume I: Mechanisms and Measurements (Harvard University Press, 2006).Book 

Google Scholar 
D’Alba, L. & Shawkey, M. D. Melanosomes: Biogenesis, properties, and evolution of an ancient organelle. Physiol. Rev. 99, 1–19 (2018).Article 

Google Scholar 
Delhey, K., Burger, C., Fiedler, W. & Peters, A. Seasonal changes in colour: A comparison of structural, melanin- and carotenoid-based plumage colours. PLoS ONE 5, e11582 (2010).Article 
ADS 

Google Scholar 
Galván, I., Mousseau, T. A. & Møller, A. P. Bird population declines due to radiation exposure at Chernobyl are stronger in species with pheomelanin-based coloration. Oecologia 165, 827–835 (2011).Article 
ADS 

Google Scholar 
Meunier, J., Pinto, S. F., Burri, R. & Roulin, A. Eumelanin-based coloration and fitness parameters in birds: A meta-analysis. Behav. Ecol. Sociobiol. 65, 559–567 (2011).Article 

Google Scholar 
Roulin, A., Almasi, B., Meichtry-Stier, K. S. & Jenni, L. Eumelanin- and pheomelanin-based colour advertise resistance to oxidative stress in opposite ways. J. Evol. Biol. 24, 2241–2247 (2011).Article 
CAS 

Google Scholar 
Gasparini, J. et al. Strength and cost of an induced immune response are associated with a heritable melanin-based colour trait in female tawny owls. J. Anim. Ecol. 78, 608–616 (2009).Article 

Google Scholar 
Fargallo, J. A. et al. Sex-specific phenotypic integration: Endocrine profiles, coloration, and behavior in fledgling boobies. Behav. Ecol. 25, 76–87 (2013).Article 

Google Scholar 
Wittkopp, P. J. & Beldade, P. Development and evolution of insect pigmentation: Genetic mechanisms and the potential consequences of pleiotropy. Semin. Cell Dev. Biol. 20, 65–71 (2009).Article 
CAS 

Google Scholar 
Hubbard, J. K., Uy, J. A. C., Hauber, M. E., Hoekstra, H. E. & Safran, R. J. Vertebrate pigmentation: From underlying genes to adaptive function. Trends Genet. 26, 231–239 (2010).Article 
CAS 

Google Scholar 
McKinnon, J. S. & Pierotti, M. E. Colour polymorphism and correlated characters: Genetic mechanisms and evolution. Mol. Ecol. 19, 5101–5125 (2010).Article 

Google Scholar 
Poston, J. P., Hasselquist, D., Stewart, I. R. & Westneat, D. F. Dietary amino acids influence plumage traits and immune responses of male house sparrows, Passer domesticus, but not as expected. Anim. Behav. 70, 1171–1181 (2005).Article 

Google Scholar 
McGraw, K. J. Dietary mineral content influences the expression of melanin-based ornamental coloration. Behav. Ecol. 18, 137–142 (2007).Article 

Google Scholar 
Fargallo, J. A., Martínez-Padilla, J., Toledano-Díaz, A., Santiago-Moreno, J. & Dávila, J. A. Sex and testosterone effects on growth, immunity and melanin coloration of nestling Eurasian kestrels. J. Anim. Ecol. 76, 201–209 (2007).Article 

Google Scholar 
Fitze, P. S. & Richner, H. Differential effects of a parasite on ornamental structures based on melanins and carotenoids. Behav. Ecol. 13, 401–407 (2002).Article 

Google Scholar 
Roulin, A., Altwegg, R., Jensen, H., Steinsland, I. & Schaub, M. Sex-dependent selection on an autosomal melanic female ornament promotes the evolution of sex ratio bias. Ecol. Lett. 13, 616–626 (2010).Article 

Google Scholar 
Sharma, A. Effect of ambient humidity on UV/visible photodegradation of melanin thin films. Photochem. Photobiol. 86, 852–855 (2010).Article 
CAS 

Google Scholar 
Burtt, E. H. The adaptiveness of animal colors. Bioscience 31, 723–729 (1981).Article 

Google Scholar 
Heppner, F. The metabolic significance of differential absorption of radiant energy by black and white birds. Condor 72, 50–59 (1970).Article 

Google Scholar 
Clusella-Trullas, S., Terblanche, J. S., Blackburn, T. M. & Chown, S. L. Testing the thermal melanism hypothesis: A macrophysiological approach. Funct. Ecol. 22, 232–238 (2008).Article 

Google Scholar 
Zink, R. M. & Remsen, J. V. Evolutionary processes and patterns of geographic variation in birds. Curr. Ornithol. 4, 1–69 (1986).
Google Scholar 
Burtt, E. H. & Ichida, J. M. Gloger’s rule, feather-degrading bacteria, and color variation among song sparrows. Condor 106, 681–686 (2004).Article 

Google Scholar 
Ruiz-De-Castaneda, R., Burtt, E. H. Jr., Gonzalez-Braojos, S. & Moreno, J. Bacterial degradability of an intrafeather unmelanized ornament: A role for feather-degrading bacteria in sexual selection?. Biol. J. Linn. Soc. 105, 409–419 (2012).Article 

Google Scholar 
Goldstein, G. et al. Bacterial degradation of black and white feathers. Auk 121, 656–659 (2004).Article 

Google Scholar 
Ducrest, A. L., Keller, L. & Roulin, A. Pleiotropy in the melanocortin system, coloration and behavioural syndromes. Trends Ecol. Evol. 23, 502–510 (2008).Article 

Google Scholar 
Kim, S. Y., Fargallo, J. A., Vergara, P. & Martínez-Padilla, J. Multivariate heredity of melanin-based coloration, body mass and immunity. Heredity 111, 139–146 (2013).Article 
CAS 

Google Scholar 
Horrocks, N. P. C. et al. Environmental proxies of antigen exposure explain variation in immune investment better than indices of pace of life. Oecologia 177, 281–290 (2015).Article 
ADS 

Google Scholar 
McLean, N., Van Der Jeugd, H. P. & van de Pol, M. High intra-specific variation in avian body condition responses to climate limits generalisation across species. PLoS ONE 13, e0192401 (2018).Article 

Google Scholar 
Gardner, J. L. et al. Spatial variation in avian bill size is associated with humidity in summer among Australian passerines. Clim. Change Responses 3, 11 (2016).Article 

Google Scholar 
Gerson, A. R. et al. Flight at low ambient humidity increases protein catabolism in migratory birds. Science 333, 1434–1436 (2011).Article 
ADS 
CAS 

Google Scholar  More

In 2008, I was investigating the methane bubbling up on the beaches and in shallow waters of Mocha Island, off the coast of central Chile. I became intrigued by how microorganisms could thrive in methane-rich areas and changed my research focus from marine biology to extreme environments. I wanted to understand how methane acts as a source of energy and carbon for microbes.Since then, I have explored a number of bizarre environments. In 2010, I went in a submarine down to 200 metres in the Black Sea, one of the world’s largest anoxic water bodies. There, I found mats of filamentous bacteria that survive on sulfur compounds.In 2017, I studied the microbes in Canada’s tailing ponds, artificial lakes of water, sand and clay waste that are left behind after petroleum extraction. And I sampled the microorganisms living in 100 °C Antarctic hot springs in 2022.I came home to Chile in 2018 and began collaborating with an international team researching the geomicrobiology of thermal features, including hot springs, geysers and volcanoes. After travelling with the group to Argentina’s active volcanic region, I got funding to explore the microbial communities that exist beneath hydrothermal vents in southern Chile, where the oceanic crust is subducting beneath the continental plate.In this image, I am in the Atacama Desert in South America, the driest non-polar desert on the planet. I am measuring 80–100 °C steam released from a fumarole containing yellow sulfur, which crystallizes at its opening as the vapour cools. I also sampled sub-surface microbes that are flushed out with the fluids. We’ll sequence their DNA to assess the microbial communities and their biological interactions.My goal is to learn more about subsurface microbes in extreme environments. I want to understand how microbial forces shaped the planet and how these communities might shift in the future with climate change. More