Philippa Kaur

1.Polis, G. A. & Strong, D. R. Food web complexity and community dynamics. Am. Nat. 147, 813–846 (1996).Article 

Google Scholar 
2.Lindeman, R. L. The trophic-dynamic aspect of ecology. Ecology 23, 399–417 (1942).Article 

Google Scholar 
3.Coleman, D. C., Andrews, R., Ellis, J. E. & Singh, J. S. Energy flow and partitioning in selected man-managed and natural ecosystems. Agro-Ecosyst. 3, 45–54 (1976).Article 

Google Scholar 
4.Steffan, S. A. et al. Unpacking brown food-webs: Animal trophic identity reflects rampant microbivory. Ecol. Evol. 7, 3532–3541 (2017).PubMed 
PubMed Central 
Article 

Google Scholar 
5.Coleman, D. C. Energetics of detritivory and microbivory in soil in theory and practice. In Food Webs (eds. Polis G. A. & Winemiller K. O.) 39–50 (Springer, 1996).Chapter 

Google Scholar 
6.Moore, J. C. et al. Detritus, trophic dynamics and biodiversity. Ecol. Lett. 7, 584–600 (2004).ADS 
Article 

Google Scholar 
7.Hagen, E. M. et al. A meta-analysis of the effects of detritus on primary producers and consumers in marine, freshwater, and terrestrial ecosystems. Oikos 121, 1507–1515 (2012).Article 

Google Scholar 
8.Danovaro, R., Snelgrove, P. V. R. & Tyler, P. Challenging the paradigms of deep-sea ecology. Trends Ecol. Evol. 29, 465–475 (2014).PubMed 
Article 

Google Scholar 
9.Smith, C. R., De Leo, F. C., Bernardino, A. F., Sweetman, A. K. & Arbizu, P. M. Abyssal food limitation, ecosystem structure and climate change. Trends Ecol. Evol. 23, 518–528 (2008).PubMed 
Article 

Google Scholar 
10.Gage, J. D. & Tyler, P. A. Deep-Sea Biology: A Natural History of Organisms at the Deep-Sea Floor. (Cambridge University Press, 1991).Book 

Google Scholar 
11.De La Rocha, C. L. & Passow, U. Factors influencing the sinking of POC and the efficiency of the biological carbon pump. Deep Res. Part II Top. Stud. Oceanogr. 54, 639–658 (2007).Article 

Google Scholar 
12.Smith, K. L., Ruhl, H. A., Huffard, C. L., Messié, M. & Kahru, M. Episodic organic carbon fluxes from surface ocean to abyssal depths during long-term monitoring in NE Pacific. Proc. Natl. Acad. Sci. USA. 115, 12235–12240 (2018).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
13.Ramirez-Llodra, E. et al. Deep, diverse and definitely different: Unique attributes of the world’s largest ecosystem. Biogeosciences 7, 2851–2899 (2010).ADS 
Article 

Google Scholar 
14.Ruhl, H. A. Abundance and size distribution dynamics of abyssal epibenthic megafauna in the northeast Pacific. Ecology 88, 1250–1262 (2007).PubMed 
Article 
PubMed Central 

Google Scholar 
15.Billett, D. S. M. Deep-sea holothurians. Oceanogr. Mar. Biol. An Annu. Rev. 29, 259–317 (1991).
Google Scholar 
16.Bett, B. J., Malzone, M. G., Narayanaswamy, B. E. & Wigham, B. D. Temporal variability in phytodetritus and megabenthic activity at the seabed in the deep northeast Atlantic. Prog. Oceanogr. 50, 349–368 (2001).ADS 
Article 

Google Scholar 
17.Durden, J. M. et al. Response of deep-sea deposit-feeders to detrital inputs: A comparison of two abyssal time-series sites. Deep. Res. Part II Top. Stud. Oceanogr. 173, 104677 (2020).18.Khripounoff, A. & Sibuet, M. L. nutrition d’echinodermes abyssaux I. Alimentation des holothuries. Mar. Biol. 60, 17–26 (1980).Article 

Google Scholar 
19.Roberts, D., Gebruka, A., Levin, V. & Manship, B. A. D. Feeding and digestive strategies in deposit-feeding holothurians. Oceanogr. Mar. Biol. Annu. Rev. 38, 257–310 (2000).
Google Scholar 
20.FitzGeorge-Balfour, T., Billett, D. S. M., Wolff, G. A., Thompson, A. & Tyler, P. A. Phytopigments as biomarkers of selectivity in abyssal holothurians; interspecific differences in response to a changing food supply. Deep. Res. Part II Top. Stud. Oceanogr. 57, 1418–1428 (2010).21.Miller, R. J., Smith, C. R., Demaster, D. J. & Fornes, W. L. Feeding selectivity and rapid particle processing by deep-sea megafaunal deposit feeders : A 234 Th tracer approach. J. Mar. Res. 58, 653–673 (2000).Article 

Google Scholar 
22.Witte, U. et al. In situ experimental evidence of the fate of a phytodetritus pulse at the abyssal sea floor. Lett. Nat. 424, 763–766 (2003).CAS 
Article 

Google Scholar 
23.Moore, H., Manship, B. & Roberts, D. Gut structure and digestive strategies in three species of abyssal holothurians. in Echinoderm Research 111–119 (Balkema, 1995).24.Deming, J. W. & Colwell, R. R. Barophilic bacteria associated with digestive tracts of abyssal holothurians. Appl. Environ. Microbiol. 44, 1222–1230 (1982).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
25.Amaro, T., Luna, G. M., Danovaro, R., Billett, D. S. M. & Cunha, M. R. High prokaryotic biodiversity associated with gut contents of the holothurian Molpadia musculus from the Nazaré Canyon (NE Atlantic). Deep. Res. Part I Oceanogr. Res. Pap. 63, 82–90 (2012).26.Roberts, D. et al. Sediment distribution, hydrolytic enzyme profiles and bacterial activities in the guts of Oneirophanta mutabilis, Psychropotes longicauda and Pseudostichopus villosus: What do they tell us about digestive strategies of abyssal holothurians?. Prog. Oceanogr. 50, 443–458 (2001).ADS 
Article 

Google Scholar 
27.Sibuet, M., Khripounoff, A., Deming, J., Colwell, R. & Dinet, A. Modification of the gut contents in the digestive tract of abyssal holothurians. In Proceedings of the International Echinoderm Conference. Tampa Bay. (ed Lawrence, J. M.) 421–428 (Balkema, 1982).
Google Scholar 
28.Bradley, C. J. et al. Trophic position estimates of marine teleosts using amino acid compound specific isotopic analysis. Limnol. Oceanogr. Methods 13, 476–493 (2015).Article 

Google Scholar 
29.Ohkouchi, N. et al. Advances in the application of amino acid nitrogen isotopic analysis in ecological and biogeochemical studies. Org. Geochem. 113, 150–174 (2017).CAS 
Article 

Google Scholar 
30.Popp, B. N. et al. Stable isotopes as indicators of ecological change. Terr. Ecol. 1, 173–190 (2007).
Google Scholar 
31.Chikaraishi, Y. et al. Determination of aquatic food-web structure based on compound-specific nitrogen isotopic composition of amino acids. Limnol. Oceanogr. Methods 7, 740–750 (2009).CAS 
Article 

Google Scholar 
32.McClelland, J. W. & Montoya, J. P. Trophic relationships and the nitrogen isotopic composition of amino acids in plankton. Ecology 83, 2173–2180 (2002).Article 

Google Scholar 
33.Steffan, S. A. et al. Microbes are trophic analogs of animals. Proc. Natl. Acad. Sci. U. S. A. 112, 15119–15124 (2015).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
34.Yamaguchi, Y. T. et al. Fractionation of nitrogen isotopes during amino acid metabolism in heterotrophic and chemolithoautotrophic microbes across Eukarya, Bacteria, and Archaea: Effects of nitrogen sources and metabolic pathways. Org. Geochem. 111, 101–112 (2017).CAS 
Article 

Google Scholar 
35.Iken, K., Brey, T., Wand, U., Voigt, J. & Junghans, P. Food web structure of the benthic community at the Porcupine Abyssal Plain (NE Atlantic): A stable isotope analysis. Prog. Oceanogr. 50, 383–405 (2001).ADS 
Article 

Google Scholar 
36.Romero-Romero, S. et al. Seasonal pathways of organic matter within the Avilés submarine canyon: Food web implications. Deep. Res. Part I. Oceanogr. Res. Pap. 117, 1–10 (2016).ADS 
CAS 
Article 

Google Scholar 
37.Reid, W., Wigham, B., McGill, R. & Polunin, N. Elucidating trophic pathways in benthic deep-sea assemblages of the Mid-Atlantic Ridge north and south of the Charlie-Gibbs fracture zone. Mar. Ecol. Prog. Ser. 463, 89–103 (2012).ADS 
Article 

Google Scholar 
38.Drazen, J. et al. Bypassing the abyssal benthic food web: Macrourid diet in the eastern North Pacific inferred from stomach content and stable isotopes analyses. Limnol. Oceanogr. 53, 2644–2654 (2008).ADS 
Article 

Google Scholar 
39.Post, D. Using stable isotopes to estimate trophic position: Models, methods, and assumptions. Ecology 83, 703–718 (2002).Article 

Google Scholar 
40.Witbaard, R., Duineveld, G. C. A., Kok, A., Van Der Weele, J. & Berghuis, E. M. The response of Oneirophanta mutabilis (Holothuroidea) to the seasonal deposition of phytopigments at the Porcupine Abyssal Plain in the Northeast Atlantic. Prog. Oceanogr. 50, 423–441 (2001).ADS 
Article 

Google Scholar 
41.Wigham, B. D., Hudson, I. R., Billett, D. S. M. & Wolff, G. A. Is long-term change in the abyssal Northeast Atlantic driven by qualitative changes in export flux? Evidence from selective feeding in deep-sea holothurians. Prog. Oceanogr. 59, 409–441 (2003).ADS 
Article 

Google Scholar 
42.Hudson, I. R., Wigham, B. D., Billett, D. S. M. & Tyler, P. A. Seasonality and selectivity in the feeding ecology and reproductive biology of deep-sea bathyal holothurians. Prog. Oceanogr. 59, 381–407 (2003).ADS 
Article 

Google Scholar 
43.Lauerman, L. M. L., Smoak, J. M., Shaw, T. J., Moore, W. S. & Smith, K. L. 234Th and 210Pb evidence for rapid ingestion of settling particles by mobile epibenthic megafauna in the abyssal NE Pacific. Limnol. Oceanogr. 42, 589–595 (1997).ADS 
CAS 
Article 

Google Scholar 
44.Roberts, D., Billett, D. S. M., McCartney, G. & Hayes, G. E. Procaryotes on the tentacles of deep-sea holothurians: A novel form of dietary supplementation. Limnol. Oceanogr. 36, 1447–1451 (1991).ADS 
Article 

Google Scholar 
45.Larsen, T., Lee Taylor, D., Leigh, M. B. & O’Brien, D. M. Stable isotope fingerprinting: A novel method for identifying plant, fungal, or bacterial origins of amino acids. Ecology 90, 3526–3535 (2009).46.Plante, C. J., Jumars, P. A. & Baross, J. A. Digestive associations between marine detritivores and bacteria. Annu. Rev. Ecol. Syst. 21, 93–127 (1990).Article 

Google Scholar 
47.Drazen, J. C., Phleger, C. F., Guest, M. A. & Nichols, P. D. Lipid, sterols and fatty acid composition of abyssal holothurians and ophiuroids from the North-East Pacific Ocean: Food web implications. Comp. Biochem. Physiol. – B Biochem. Mol. Biol. 151, 79–87 (2008).48.Amaro, T. et al. Possible links between holothurian lipid compositions and differences in organic matter (OM) supply at the western Pacific abyssal plains. Deep. Res. Part I Oceanogr. Res. Pap. 152, (2019).49.Kharlamenko, V. I., Maiorova, A. S. & Ermolenko, E. V. Fatty acid composition as an indicator of the trophic position of abyssal megabenthic deposit feeders in the Kuril Basin of the Sea of Okhotsk. Deep. Res. Part II Top. Stud. Oceanogr. 154, 374–382 (2018).50.Ginger, M. L. et al. Organic matter assimilation and selective feeding by holothurians in the deep sea: Some observations and comments. Prog. Oceanogr. 50, 407–421 (2001).ADS 
Article 

Google Scholar 
51.McMahon, K. W., Thorrold, S. R., Houghton, L. A. & Berumen, M. L. Tracing carbon flow through coral reef food webs using a compound-specific stable isotope approach. Oecologia 180, 809–821 (2016).ADS 
PubMed 
Article 
PubMed Central 

Google Scholar 
52.Kaufmann, R. S. & Smith, K. L. Activity patterns of mobile epibenthic megafauna at an abyssal site in the eastern North Pacific: Results from a 17-month time-lapse photographic study. Deep. Res. Part I Oceanogr. Res. Pap. 44, 559–579 (1997).53.Kuhnz, L. A., Ruhl, H. A., Huffard, C. L. & Smith, K. L. Rapid changes and long-term cycles in the benthic megafaunal community observed over 24 years in the abyssal northeast Pacific. Prog. Oceanogr. 124, 1–11 (2014).ADS 
Article 

Google Scholar 
54.Vardaro, M. F., Ruhl, H. A. & Smith, K. L. Climate variation, carbon flux, and bioturbation in the abyssal north pacific. Limnol. Oceanogr. 54, 2081–2088 (2009).ADS 
Article 

Google Scholar 
55.Ziegler, A., Mooi, R., Rolet, G. & De Ridder, C. Origin and evolutionary plasticity of the gastric caecum in sea urchins (Echinodermata: Echinoidea). BMC Evol. Biol. 10, 313 (2010).PubMed 
PubMed Central 
Article 

Google Scholar 
56.McCarthy, M. D., Benner, R., Lee, C. & Fogel, M. L. Amino acid nitrogen isotopic fractionation patterns as indicators of heterotrophy in plankton, particulate, and dissolved organic matter. Geochim. Cosmochim. Acta 71, 4727–4744 (2007).ADS 
CAS 
Article 

Google Scholar 
57.Calleja, M. L., Batista, F., Peacock, M., Kudela, R. & McCarthy, M. D. Changes in compound specific δ15N amino acid signatures and d/l ratios in marine dissolved organic matter induced by heterotrophic bacterial reworking. Mar. Chem. 149, 32–44 (2013).CAS 
Article 

Google Scholar 
58.Smith, K. L., Ruhl, H. A., Kaufmann, R. S. & Kahru, M. Tracing abyssal food supply back to upper-ocean processes over a 17-year time series in the northeast Pacific. Limnol. Oceanogr. 53, 2655–2667 (2008).ADS 
Article 

Google Scholar 
59.Huffard, C. L., Kuhnz, L. A., Lemon, L., Sherman, A. D. & Smith, K. L. Demographic indicators of change in a deposit-feeding abyssal holothurian community (Station M, 4000 m). Deep. Res. Part I Oceanogr. Res. Pap. 109, 27–39 (2016).60.Hannides, C. C. S., Popp, B. N., Anela Choy, C. & Drazen, J. C. Midwater zooplankton and suspended particle dynamics in the North Pacific Subtropical Gyre: A stable isotope perspective. Limnol. Oceanogr. 58, 1931–1936 (2013).61.Neto, R. R., Wolff, G. A., Billett, D. S. M., Mackenzie, K. L. & Thompson, A. The influence of changing food supply on the lipid biochemistry of deep-sea holothurians. Deep. Res. Part I Oceanogr. Res. Pap. 53, 516–527 (2006).62.Amaro, T., Witte, H., Herndl, G. J., Cunha, M. R. & Billett, D. S. M. Deep-sea bacterial communities in sediments and guts of deposit-feeding holothurians in Portuguese canyons (NE Atlantic). Deep Res. Part I Oceanogr. Res. Pap. 56, 1834–1843 (2009).ADS 
Article 

Google Scholar 
63.Silfer, J. A., Engel, M. H., Macko, S. A. & Jumeau, E. J. Stable carbon isotope analysis of amino acid enantiomers by conventional isotope ratio mass spectrometry and combined gas chromatography/isotope ratio mass spectrometry. Anal. Chem. 63, 370–374 (1991).CAS 
Article 

Google Scholar 
64.Jarman, C. L. et al. Diet of the prehistoric population of Rapa Nui (Easter Island, Chile) shows environmental adaptation and resilience. Am. J. Phys. Anthropol. 164, 343–361 (2017).PubMed 
PubMed Central 
Article 

Google Scholar 
65.Dauwe, B., Middelburg, J. J., Herman, P. M. J. & Heip, C. H. R. Linking diagenetic alteration of amino acids and bulk organic matter reactivity. Limnology 44, 1809–1814 (1999).CAS 

Google Scholar 
66.McCarthy, M. D., Benner, R., Lee, C., Hedges, J. I. & Fogel, M. L. Amino acid carbon isotopic fractionation patterns in oceanic dissolved organic matter: An unaltered photoautotrophic source for dissolved organic nitrogen in the ocean?. Mar. Chem. 92, 123–134 (2004).CAS 
Article 

Google Scholar 
67.R: A Language and Environment for Statistical Computing. https://www.R-project.org (R Core Team. R Foundation for Statistical Computing, 2020). More

1.Pereira, H. M. et al. Scenarios for global biodiversity in the 21st century. Science 330, 1496–1501 (2010).ADS 
CAS 
PubMed 
Article 

Google Scholar 
2.Pereira, H. M. et al. Essential biodiversity variables. Science 339, 277–278 (2013).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
3.Stem, C., Margoluis, R., Salafsky, N. & Brown, M. Monitoring and evaluation in conservation: a review of trends and approaches. Conserv. Biol. 19, 295–309 (2005).Article 

Google Scholar 
4.Atwood, T. B. et al. Herbivores at the highest risk of extinction among mammals, birds, and reptiles. Sci. Adv. 6, eabb8458 (2020).ADS 
PubMed 
PubMed Central 
Article 

Google Scholar 
5.Cardillo, M. et al. Human population density and extinction risk in the world’s carnivores. PLoS Biol 2, e197 (2004).PubMed 
PubMed Central 
Article 

Google Scholar 
6.Ripple, W. J. et al. Status and ecological effects of the world’s largest carnivores. Science 343, 1241484 (2014).PubMed 
Article 
CAS 
PubMed Central 

Google Scholar 
7.Wolf, C. & Ripple, W. J. Prey depletion as a threat to the world’s large carnivores. R. Soc. Open Sci. 3, 160252 (2016).ADS 
PubMed 
PubMed Central 
Article 

Google Scholar 
8.Kissling, W. D. et al. Building essential biodiversity variables (EBV s) of species distribution and abundance at a global scale. Biol. Rev. 93, 600–625 (2018).PubMed 
Article 
PubMed Central 

Google Scholar 
9.Nesshöver, C., Livoreil, B., Schindler, S. & Vandewalle, M. Challenges and solutions for networking knowledge holders and better informing decision-making on biodiversity and ecosystem services. Biodivers. Conserv. 25, 1207–1214 (2016).Article 

Google Scholar 
10.Gardner, T. A. et al. The cost-effectiveness of biodiversity surveys in tropical forests. Ecol. Lett. 11, 139–150 (2008).PubMed 
Article 
PubMed Central 

Google Scholar 
11.Field, S. A., Tyre, A. J. & Possingham, H. P. Optimizing allocation of monitoring effort under economic and observational constraints. J. Wildl. Manag. 69, 473–482 (2005).Article 

Google Scholar 
12.Braunisch, V. & Suchant, R. Predicting species distributions based on incomplete survey data: The trade-off between precision and scale. Ecography 33, 826–840 (2010).Article 

Google Scholar 
13.Taberlet, P., Coissac, E., Hajibabaei, M. & Rieseberg, L. H. Environmental DNA. Mol. Ecol. 21, 1789–1793 (2012).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
14.Deiner, K. et al. Long-range PCR allows sequencing of mitochondrial genomes from environmental DNA. Methods Ecol. Evol. 8, 1888–1898 (2017).Article 

Google Scholar 
15.Deiner, K., Walser, J.-C., Mächler, E. & Altermatt, F. Choice of capture and extraction methods affect detection of freshwater biodiversity from environmental DNA. Biol. Conserv. 183, 53–63 (2015).Article 

Google Scholar 
16.Tsuji, S., Takahara, T., Doi, H., Shibata, N. & Yamanaka, H. The detection of aquatic macroorganisms using environmental DNA analysis—A review of methods for collection, extraction, and detection. Environ. DNA 1, 99–108 (2019).Article 

Google Scholar 
17.Sales, N. G. et al. Fishing for mammals: Landscape-level monitoring of terrestrial and semi-aquatic communities using eDNA from riverine systems. J. Appl. Ecol. 57, 707–716 (2020).CAS 
Article 

Google Scholar 
18.Sales, N. G. et al. Assessing the potential of environmental DNA metabarcoding for monitoring Neotropical mammals: a case study in the Amazon and Atlantic Forest, Brazil. Mammal Rev. 50, 221–225 (2020).Article 

Google Scholar 
19.Barnes, M. A. & Turner, C. R. The ecology of environmental DNA and implications for conservation genetics. Conserv. Genet. 17, 1–17 (2016).CAS 
Article 

Google Scholar 
20.Deiner, K. et al. Environmental DNA metabarcoding: Transforming how we survey animal and plant communities. Mol. Ecol. 26, 5872–5895 (2017).PubMed 
Article 
PubMed Central 

Google Scholar 
21.Leempoel, K., Hebert, T. & Hadly, E. A. A comparison of eDNA to camera trapping for assessment of terrestrial mammal diversity. Proc. R. Soc. B Biol. Sci. 287, 20192353 (2020).CAS 
Article 

Google Scholar 
22.Harper, L. R. et al. Environmental DNA (eDNA) metabarcoding of pond water as a tool to survey conservation and management priority mammals. Biol. Conserv. 238, 108225 (2019).Article 

Google Scholar 
23.Rodgers, T. W. & Mock, K. E. Drinking water as a source of environmental DNA for the detection of terrestrial wildlife species. Conserv. Genet. Resour. 7, 693–696 (2015).Article 

Google Scholar 
24.Ushio, M. et al. Environmental DNA enables detection of terrestrial mammals from forest pond water. Mol. Ecol. Resour. 17, e63–e75 (2017).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
25.Williams, K. E., Huyvaert, K. P., Vercauteren, K. C., Davis, A. J. & Piaggio, A. J. Detection and persistence of environmental DNA from an invasive, terrestrial mammal. Ecol. Evol. 8, 688–695 (2018).PubMed 
Article 
PubMed Central 

Google Scholar 
26.Merkes, C. M., McCalla, S. G., Jensen, N. R., Gaikowski, M. P. & Amberg, J. J. Persistence of DNA in carcasses, slime and avian feces may affect interpretation of environmental DNA data. PLoS ONE 9, e113346 (2014).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
27.Pont, D. et al. Environmental DNA reveals quantitative patterns of fish biodiversity in large rivers despite its downstream transportation. Sci. Rep. 8, 1–13 (2018).ADS 
CAS 
Article 

Google Scholar 
28.Zinger, L. et al. Advances and prospects of environmental DNA in neotropical rainforests. Adv. Ecol. Res. 62, 331–373 (2020).Article 

Google Scholar 
29.Withers, P. C., Cooper, C. E., Maloney, S. K., Bozinovic, F. & Cruz-Neto, A. P. Ecological and Environmental Physiology of Mammals Vol. 5 (Oxford University Press, 2016).Book 

Google Scholar 
30.Bicudo, J. E. P., Buttemer, W. A., Chappell, M. A., Pearson, J. T. & Bech, C. Ecological and Environmental Physiology of Birds Vol. 2 (Oxford University Press, 2010).Book 

Google Scholar 
31.Naidoo, R. & Burton, A. C. Relative effects of recreational activities on a temperate terrestrial wildlife assemblage. Conserv. Sci. Pract. 2, e271 (2020).
Google Scholar 
32.Cantera, I. et al. Optimizing environmental DNA sampling effort for fish inventories in tropical streams and rivers. Sci. Rep. 9, 1–11 (2019).CAS 
Article 

Google Scholar 
33.Tarboton, D. G., Bras, R. L. & Rodriguez-Iturbe, I. The fractal nature of river networks. Water Resour. Res. 24, 1317–1322 (1988).ADS 
Article 

Google Scholar 
34.Ishige, T. et al. Tropical-forest mammals as detected by environmental DNA at natural saltlicks in Borneo. Biol. Conserv. 210, 281–285 (2017).Article 

Google Scholar 
35.Joseph, L. N., Field, S. A., Wilcox, C. & Possingham, H. P. Presence–absence versus abundance data for monitoring threatened species. Conserv. Biol. 20, 1679–1687 (2006).PubMed 
Article 
PubMed Central 

Google Scholar 
36.Lacy, R. C. Lessons from 30 years of population viability analysis of wildlife populations. Zoo Biol. 38, 67–77 (2019).PubMed 
Article 

Google Scholar 
37.Gärdenfors, U., Hilton-Taylor, C., Mace, G. M. & Rodríguez, J. P. The application of IUCN Red List criteria at regional levels. Conserv. Biol. 15, 1206–1212 (2001).Article 

Google Scholar 
38.Munro, R., Nielsen, S. E., Price, M., Stenhouse, G. & Boyce, M. S. Seasonal and diel patterns of grizzly bear diet and activity in west-central Alberta. J. Mammal. 87, 1112–1121 (2006).Article 

Google Scholar 
39.Deiner, K. & Altermatt, F. Transport distance of invertebrate environmental DNA in a natural river. PLoS ONE 9, e88786 (2014).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
40.Hunter, M. E. et al. Detection limits of quantitative and digital PCR assays and their influence in presence–absence surveys of environmental DNA. Mol. Ecol. Resour. 17, 221–229 (2017).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
41.Roussel, J.-M., Paillisson, J.-M., Treguier, A. & Petit, E. The downside of eDNA as a survey tool in water bodies. J. Appl. Ecol. 52, 823–826 (2015).CAS 
Article 

Google Scholar 
42.Williams, B. K., Nichols, J. D. & Conroy, M. J. Analysis and Management of Animal Populations (Academic Press, 2002).
Google Scholar 
43.Burton, A. C. et al. Wildlife camera trapping: A review and recommendations for linking surveys to ecological processes. J. Appl. Ecol. 52, 675–685 (2015).Article 

Google Scholar 
44.Morris, W. F. et al. Quantitative Conservation Biology (Sinauer Sunderland, 2002).
Google Scholar 
45.Parks, B. C. South Chilcotin Mountains Park and Big Creek Park Management Plan (2019).46.McLellan, M. L. et al. Divergent population trends following the cessation of legal grizzly bear hunting in southwestern British Columbia, Canada. Biol. Conserv. 233, 247–254 (2019).Article 

Google Scholar 
47.Kays, R. et al. Camera traps as sensor networks for monitoring animal communities. In 2009 IEEE 34th Conference on Local Computer Networks 811–818. https://doi.org/10.1109/LCN.2009.5355046 (IEEE, 2009).48.Hendry, H. & Mann, C. Camelot–intuitive software for camera trap data management. BioRxiv 203216 (2017).49.Valentini, A. et al. Next-generation monitoring of aquatic biodiversity using environmental DNA metabarcoding. Mol. Ecol. 25, 929–942 (2016).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
50.Taberlet, P., Bonin, A., Zinger, L. & Coissac, E. Environmental DNA: For biodiversity research and monitoring (Oxford University Press, Oxford, 2018).Book 

Google Scholar 
51.Boyer, F. et al. obitools: A unix-inspired software package for DNA metabarcoding. Mol. Ecol. Resour. 16, 176–182 (2016).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
52.De Barba, M. et al. DNA metabarcoding multiplexing and validation of data accuracy for diet assessment: Application to omnivorous diet. Mol. Ecol. Resour. 14, 306–323 (2014).PubMed 
Article 
CAS 
PubMed Central 

Google Scholar 
53.Schnell, I. B., Bohmann, K. & Gilbert, M. T. P. Tag jumps illuminated–reducing sequence-to-sample misidentifications in metabarcoding studies. Mol. Ecol. Resour. 15, 1289–1303 (2015).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
54.Wearn, O. & Glover-Kapfer, P. Camera-trapping for conservation: A guide to best-practices. WWF Conserv. Technol. Ser. 1, 2019–2104 (2017).
Google Scholar 
55.R Core Team. R: A language and environment for statistical computing (R Foundation for Statistical Computing, 2020).
Google Scholar 
56.Plummer, M., et al.. JAGS: A program for analysis of Bayesian graphical models using Gibbs sampling. In Proceedings of the 3rd International Workshop on Distributed Statistical Computing vol. 124, 1–10 (Vienna, Austria, 2003).57.Tuomisto, H. A diversity of beta diversities: Straightening up a concept gone awry. Part 1. Defining beta diversity as a function of alpha and gamma diversity. Ecography 33, 2–22 (2010).Article 

Google Scholar 
58.Ahumada, J. A. et al. Community structure and diversity of tropical forest mammals: data from a global camera trap network. Philos. Trans. R. Soc. B Biol. Sci. 366, 2703–2711 (2011).Article 

Google Scholar 
59.Samejima, H., Ong, R., Lagan, P. & Kitayama, K. Camera-trapping rates of mammals and birds in a Bornean tropical rainforest under sustainable forest management. For. Ecol. Manag. 270, 248–256 (2012).Article 

Google Scholar 
60.Parsons, A. W. et al. Do occupancy or detection rates from camera traps reflect deer density?. J. Mammal. 98, 1547–1557 (2017).Article 

Google Scholar 
61.Sollmann, R., Mohamed, A., Samejima, H. & Wilting, A. Risky business or simple solution–Relative abundance indices from camera-trapping. Biol. Conserv. 159, 405–412 (2013).Article 

Google Scholar 
62.Villette, P., Krebs, C. J., Jung, T. S. & Boonstra, R. Can camera trapping provide accurate estimates of small mammal (Myodes rutilus and Peromyscus maniculatus) density in the boreal forest?. J. Mammal. 97, 32–40 (2016).Article 

Google Scholar 
63.Feldhamer, G. A., Thompson, B. C. & Chapman, J. A. Wild Mammals of North America: Biology, Management, and Conservation (The Johns Hopkins University Press, 2003).
Google Scholar 
64.Burnham, K. P. & Anderson, D. R. A Practical Information-Theoretic Approach. Model Selection Multimodel Inference 2nd edn, Vol. 2 (Springer, Berlin, 2002).MATH 

Google Scholar  More

Metacommunity model and asymptotic community assemblyWe built a large set of model metacommunities (detailed in full in “Methods”) describing competitive dynamics within a single guild of species across a landscape. Each metacommunity consisted of a set of patches, or local communities, randomly placed in a square arena and linked by a spatial network. The dynamics of each population are governed by three processes: inter- and intraspecific interactions, heterogeneous responses to the environment and dispersal between adjacent patches (Fig. 1). Competition coefficients between species are drawn at random and the population dynamics within each patch are described by a Lotka-Volterra competition model. We control the level of environmental heterogeneity across the network directly by generating an intrinsic growth rate for each species at each patch from a random, spatially correlated distribution. To ensure any turnover is purely autonomous, we keep the environment fixed throughout simulations. Dispersal between neighbouring patches declines exponentially with distance between sites. This formulation allows precise and independent control of key properties of the metacommunity–the number of patches, the characteristic dispersal length and the heterogeneity of the environment.Fig. 1: Elements of the Lotka-Volterra metacommunity model and the emergence of autonomous population dynamics.Environmental heterogeneity, represented by the intrinsic growth rate matrix R, is modelled using a spatially autocorrelated Gaussian random field. A random spatial network, represented by the dispersal matrix D, defines the spatial connectivity of the landscape. The network of species interactions, represented by the competitive overlap matrix A, is modelled by sampling competition coefficients at random (perpendicular bars indicate recipients of a deleterious competitive impact). The resulting dynamics of local population biomasses, given by the colour-coded equation, are numerically simulated. The Hadamard product ‘∘’ represents element-wise matrix multiplication. For large metacommunities, local populations exhibit persistent dynamics despite the absence of external drivers. In the 3D boxes, typical simulated biomass dynamics of dominating species are plotted on linear axes over 2500 unit times. The graphs illustrate the complexity of the autonomous dynamics and the propensity for compositional change (local extinction and colonisation).Full size imageTo populate the model metacommunities, we iteratively introduced species with randomly generated intrinsic growth rates and interspecific interaction coefficients. Between successive regional invasions we simulated the model dynamics, and removed any species whose abundance fell below a threshold across the whole network. Through this assembly process and the eventual onset of ecological structural instability, both average local diversity, the number of species coexisting in a given patch, and regional diversity, the total number of species in the metacommunity, eventually saturate and then fluctuate around an equilibrium value—any introduction of a new species then leads on average to the extinction of one other species (Supplementary Fig. 1). In these intrinsically regulated metacommunities we then studied the phenomenology of autonomous community turnover in the absence of regional invasions or abiotic change.In our metacommunity model, local community dynamics and therefore local limits on species richness depend on a combination of biotic and abiotic filtering (non-uniform responses of species to local conditions)33,34,35 and immigration from adjacent patches, generating so called mass effects in the local community36,37,38. Biotic filtering via interspecific competition is encoded in the interaction coefficients Aij, while abiotic filtering occurs via the spatial variation of intrinsic growth rates Rix. For simplicity, and since predator-prey dynamics are known to generate oscillations39 through mechanisms distinct from those we report here, we restrict our analysis to competitive communities for which all ecological interactions are antagonistic. The off-diagonal elements of the interaction matrix A describe how one species i affects another species j. These are sampled independently from a discrete distribution, such that the interaction strength Aij is set to a constant value in the range 0 to 1 (in most cases 0.5) with fixed probability (connectance, in most cases 0.5) and otherwise set to zero. Intraspecific competition coefficients Aii are set to 1 for all species. This discrete distribution of the interaction terms was chosen for its relative efficiency. In the Supplementary discussion (and Supplementary Fig. 2) we show that outcomes remain unaffected when more complex distributions are modelled. Intrinsic growth rates Rix are sampled from spatially correlated normal distributions with mean 1, autocorrelation length ϕ and variance σ2 (Supplementary Fig. 3).Dispersal is modelled via a spatial connectivity matrix with elements Dxy. The topology of the model metacommunity, expressed through D, is generated by sampling the spatial coordinates of N patches from a uniform distribution ({mathcal{U}}(0,sqrt{N})times {mathcal{U}}(0,sqrt{N})), i.e., an area of size N. Thus, under variation of the number of patches, the inter-patch distances remain fixed on average. Spatial connectivity is defined by linking these patches through a Gabriel graph40, a planar graph generated by an algorithm that, on average, links each local community to four close neighbours41. Avoidance of direct long-distance dispersal and the sparsity of the resulting dispersal matrix permit the use of efficient numerical methods. The exponential dispersal kernel defining Dxy is tuned by the dispersal length ℓ, which is fixed for all species.The dynamics of local population biomasses Bix = Bix(t) are modelled using a system of spatially coupled Lotka-Volterra (LV) equations that, in matrix notation, takes the form23$$frac{d{bf{B}}}{dt}={bf{B}}circ ({bf{R}}-{bf{A}}{bf{B}})+{bf{B}}{bf{D}},$$
(1)
with ∘ denoting element-wise multiplication. Hereafter this formalism is referred to as the Lotka-Volterra Metacommunity Model (LVMCM). Further technical details are provided in Methods and the Supplementary Discussion.In order to numerically probe the impacts of ℓ, ϕ and σ2 on the emergent temporal dynamics, we initially fixed N = 64 and varied each parameter through multiple orders of magnitude (Supplementary Fig. 4). In order to obtain a full characterisation of autonomous turnover in the computationally accessible spatial range (N ≤ 256), we then selected a parameter combination found to generate substantial fluctuations for further analysis. Thereafter we assembled metacommunities of 8–256 patches (Fig. 2a) until regional diversity limits were reached (with tenfold replication) and generated community time series of 104 unit times from which the phenomenology of autonomous turnover could be explored in detail. We found no evidence to suggest that the phenomenology described below depends on this specific parameter combination. While future results may confirm or refute this, autonomous turnover arises over a wide range of parameters (Supplementary Fig. 4) and as such the phenomenon is robust.Fig. 2: Autonomous turnover in model metacommunities.a Typical model metacommunities: a spatial network with N nodes representing local communities (or patches) and edges, channels of dispersal. Patch colour represents the number of clusters in local community state space detected over 104 unit times t using hierarchical clustering of the Bray-Curtis (BC) dissimilarity matrix, Supplementary Fig. 6. b Colour coded matrices of pairwise temporal BC dissimilarity corresponding to the circled patches in (a). Insets represent 102 unit times. For small networks (N = 8) local compositions converge to static fixed points. As metacommunity extent increases, however, persistent dynamics emerge. Initially this autonomous turnover is oscillatory in nature with communities fluctuating between small numbers of states which can be grouped into clusters (16 ≤ N ≤ 32). Intermediate metacommunities (32 ≤ N ≤ 64) manifest “Clementsian” temporal turnover, characterised by sharp transitions in composition, implying species turn over in cohorts. Large metacommunities (N ≥ 128) turn over continuously, implying “Gleasonian” assembly dynamics in which species’ temporal occupancies are independent. c The mean number of local compositional clusters detected for metacommunities of various numbers of patches N (error bars represent standard deviation across all replicated simulations). While the transition from static to dynamic community composition at the local scale is sharp (see text), non-uniform turnover within metacommunities (a) blurs the transition at the regional scale. Aij = 0.5 with probability 0.5, ϕ = 10, σ2 = 0.01, ℓ = 0.5.Full size imageAutonomous turnover in model metacommunitiesFor small (N ≤ 8) metacommunities assembled to regional diversity limits, populations attain equilibria, i.e., converge to fixed points, implying the absence of autonomous turnover23. With increasing metacommunity size N, however, we observe the emergence of persistent population dynamics (Supplementary Fig. 5 and external video) that can produce substantial turnover in local community composition. This autonomous turnover can be represented through Bray-Curtis42 (BC) dissimilarity matrices comparing local community composition through time (Fig. 2b), and quantified by the number of compositional clusters detected in such matrices using hierarchical cluster analysis (Fig. 2a, c).At intermediate spatial scales (Fig. 2, 16 ≤ N ≤ 32) we often find oscillatory dynamics, which can be perfectly periodic or slightly irregular. With increasing oscillation amplitude, these lead to persistent turnover dynamics where local communities repeatedly transition between a small number of distinct compositional clusters (represented in Fig. 2 by stripes of high pairwise BC dissimilarity spanning large temporal ranges). At even larger scales (N ≥ 64) this compositional coherence begins to break down, and for very large metacommunities (N ≥ 128) autonomous dynamics drive continuous acyclic change in community composition. The number of compositional clusters detected over time typically varies within a given metacommunity (Fig. 2a node colour), however we find a clear increase in the average number of compositional clusters, i.e., an increase in turnover, with increasing total metacommunity size (Fig. 2c).Metacommunities in which the boundaries of species ranges along environmental gradients are clumped are termed Clementsian, while those for which range limits are independently distributed are  referred to as Gleasonian43. We consider the block structure of the temporal dissimilarity matrix at intermediate N to represent a form of Clementsian temporal turnover, characterised by sudden significant shifts in community composition. Metacommunity models similar to ours have been found to generate such patterns along spatial gradients44, potentially via an analogous mechanism45. Large, diverse model metacommunities manifest Gleasonian temporal turnover. In such cases, species colonisations and extirpations are largely independent and temporal occupancies predominantly uncorrelated, such that compositional change is continuous, rarely, if ever, reverting to the same state.Mechanistic explanation of autonomous turnoverSurprisingly, the onset and increasing complexity of autonomous turnover as system size N increases (Fig. 2) can be understood as a consequence of local community dynamics alone. To explain this, we first recall relevant theoretical results for isolated LV communities. Then we demonstrate that, in presence of weak propagule pressure, these results imply local community turnover dynamics, controlled by the richness of potential invaders, that closely mirror the dependence on system size seen in full LV metacommunities.Application of methods from statistical mechanics to models of large isolated LV communities with random interactions revealed that such models exhibit qualitatively distinct phases46,47,48. If the number of modelled species, S, interpreted as species pool size, lies below some threshold value determined by the distribution of interaction strengths (Supplementary Fig. 7), these models exhibit a unique linearly stable equilibrium (Unique Fixed Point phase, UFP). Some species may go extinct, but the majority persists48. When pool size S exceeds this threshold, there appear to be no more linearly stable equilibrium configurations. Any community formed by a selection from the S species is either unfeasible (there is no equilibrium with all species present), intrinsically linearly unstable, or invadable by at least one of the excluded species. This has been called the multiple attractor (MA) phase47. However, the implied notion that this part of the phase space is in fact characterised by multiple stable equilibria may be incorrect.Population dynamical models with many species have been shown to easily exhibit attractors called stable heteroclinic networks49, which are characterised by dynamics in which the system bounces around between several unstable equilibria, each corresponding to a different composition of the extant community, implying indefinite, autonomous community turnover (Fig. 3, red line). As these attractors are approached, models exhibit increasingly long intermittent phases of slow dynamics, which, when numerically simulated, can give the impression that the system eventually reaches one of several ‘stable’ equilibria, suggesting that turnover comes to a halt. We demonstrate in the Supplementary discussion that the MA phase of isolated LV models is in fact characterised by such stable heteroclinic networks (Supplementary Figs. 8 and 9). Note, we retain the MA terminology here because the underlying complete heteroclinic networks, interpreted as a directed graph50,51 (Fig. 3, inset), might have multiple components that are mutually unreachable through dynamic transitions52, each representing a different attractor.Fig. 3: Approximate heteroclinic networks underlie autonomous community turnover.The main panel shows two trajectories in the state space of a community of three hypothetical species (population biomasses B1, B2, B3) that are in non-hierarchical competition with each other, such that no species can competitively exclude both others (a “rock-paper-scissors game”17). Without propagule pressure, the system has three unstable equilibrium points (P1, P2, P3) and cycles between these (red curve), coming increasingly close to the equilibria and spending ever more time in the vicinity of each. The corresponding attractor is called a heteroclinic cycle (dashed arrows). Under weak extrinsic propagule pressure (blue curve), the three equilibria and the heteroclinic cycle disappear, yet the system closely tracks the original cycle in state space. Such a cycle can be represented as a graph linking the dynamically connected equilibria (inset). With more interacting species, these graphs can become complex “heteroclinic networks”49,50,51 with trajectories representing complex sequences of species composition during autonomous community turnover.Full size imageIf one now adds to such isolated LV models terms representing weak propagule pressure for all S species (Supplementary Eq. (2)), dynamically equivalent to mass effects occurring in the full metacommunity model (Eq. (1)), then none of the S species can entirely go extinct. The weak influx of biomass drives community states away from the unstable equilibria representing coexistence of subsets of the S species and the heteroclinic network connecting them (blue line in Fig. 3). Typically, system dynamics then still follow trajectories closely tracking the original heteroclinic networks (Fig. 3), but now without requiring boundless time to transition from the vicinity of one equilibrium to the next.The nature and complexity of the resulting population dynamics depend on the size and complexity of the underlying heteroclinic network, and both increase with pool size S. In simulations (Supplementary Fig. 10) we find that, as S increases, LV models with weak propagule pressure pass through the same sequence of states as we documented for LVMCM metacommunities in Fig. 2: equilibria, oscillatory population dynamics, Clementsian and finally Gleasonian temporal turnover.Above we introduced the number of clusters detected in Bray-Curtis dissimilarity matrices of fixed time series length as a means of quantifying the approximate number of equilibria visited during local community turnover. As shown in Fig. 4a, b, this number increases in LV models with S in a manner strikingly similar to its increase in the LVMCM with the number of species present in the ecological neighbourhood of a given patch. Thus, dynamics within a patch are controlled not by N directly but rather by neighbourhood species richness. For a given neighbourhood, species richness depends on the number of connected patches, the total area and therefore total abiotic heterogeneity encompassed, and the connectivity, all of which can vary substantially within a metacommunity of a given size N. As illustrated in Fig. 4b, there is a tendency for neighbourhood richness to be larger in larger metacommunities, leading indirectly to the dependence of metacommunity dynamics on N seen in Fig. 2.Fig. 4: Ecological mass effects drive autonomous turnover.a The number of compositional clusters detected, plotted against the size of the pool of potential invaders S for an isolated LV community using a propagule pressure ϵ of 10−10 and 10−15, fit by a generalised additive model87. For S 1 compositional cluster) occurs at a pools size of around S = 35 species, consistent with the theoretical prediction47 of the transition between the UFP and MA phases (Supplementary Discussion). Close inspection of this threshold reveals an important and hitherto unreported relationship between the transition into the MA phase and local ecological limits set by the onset of ecological structural instability, which is known to regulate species richness in LV systems subject to external invasion pressure23,24: in the Supplementary Discussion we show that the boundary between the UFP and MA phases47 coincides precisely with the onset of structural instability24 (Supplementary Eqs. (3)–(9)).For LVMCM metacommunities, this relationship (demonstrated analytically in the Supplementary Discussion) is numerically confirmed in Fig. 5. During assembly, local species richness increases until it reaches the limit imposed by local structural instability. Further assembly occurs via the “regionalisation” of the biota53—a collapse in average range sizes23 and associated increase in spatial beta diversity—until regional diversity limits are reached23. The emergence of autonomous turnover coincides with the onset of species saturation at the local scale. Autonomous turnover can therefore serve as an indirect indication of intrinsic biodiversity regulation via local structural instability in complex communities.Fig. 5: The emergence of temporal turnover during metacommunity assembly.a Local species richness, defined by reference to source populations only (({overline{alpha }}_{text{src}}), grey) and regional diversity (γ black) for a single metacommunity of N = 32 coupled communities during iterative regional invasion of random species. We quantify local source diversity ({overline{alpha }}_{text{src}}) as the metacommunity average of the number αsrc of non-zero equilibrium populations persisting when immigration is switched off (off-diagonal elements of D set to zero), since this is the component of a local community subject to strict ecological limits to biodiversity. Note the log scale chosen for easy comparison of local and regional species richness. b Increases in regional diversity beyond local limits arise via corresponding increases in spatial turnover (({overline{beta }}_{text{s}}), black). Autonomous temporal turnover (({overline{beta }}_{text{t}}), grey) sets in (crosses a threshold mean Bray Curtis (BC) dissimilarity of 10−2) precisely when average local species richness ({overline{alpha }}_{text{src}}) has reached its limit, reflecting the equivalence of the transition to the MA phase space and the onset of local structural instability. In both panels, the dashed line marks the point at which autonomous temporal turnover was first detected. Aij = 0.3 with probability 0.3, ϕ = 10, σ2 = 0.01, ℓ = 0.5. Both spatial and temporal turnover computed as the mean BC dissimilarity. In each iteration of the assembly model (regional invasion event), 0.1S + 1 species were introduced. Dynamics were simulated for 2 × 104 unit times, with the second 104 unit times analysed for autonomous turnover, and a total of 104 invasions were modelled.Full size imageThus, we have shown that propagule pressure perturbs local communities away from unstable equilibria and drives compositional change. In order to invade, however, species need to be capable of passing through biotic and abiotic filters33,34,35. We would expect, therefore, that turnover would be suppressed in highly heterogeneous or poorly connected environments where mass effects are weak. Indeed, by manipulating the autocorrelation length ϕ and variance σ2 of the abiotic filter represented by the matrix R and the characteristic dispersal length ℓ, we observe a sharp drop-off in temporal turnover in parameter regimes that maximise between-patch community dissimilarity (short environmental correlation or dispersal lengths, Supplementary Fig. 11). Thus, we conclude that it is not species richness or spatial dissimilarity per se that best predict temporal turnover, but the size of the pool of species with positive invasion fitness, i.e., those not repelled by the combined effects of biotic and abiotic filters.The macroecology of autonomous turnoverWe find good correspondence between temporal and spatio-temporal biodiversity patterns emerging in model metacommunities in the absence of external abiotic change and in empirical data (Fig. 6), with quantitative characteristics lying within the ranges observed in natural ecosystems.Fig. 6: Macroecological signatures of autonomous compositional change.A bimodal distribution in temporal occupancy observed in North American birds54 (a) and in simulations (e N = 64, ϕ = 5, σ2 = 0.01, ℓ = 0.5). Intrisically regulated species richness observed in estuarine fish species59 (b) and in simulations (f N = 64, ϕ = 5, σ2 = 0.01, ℓ = 0.5, 1000 unit times t). The decreasing slopes of the STR with increasing sample area12 (c), and the SAR with increasing sample duration12 (d) for various communities and in simulations (g and h N = 256, ϕ = 10, σ2 = 0.01, ℓ = 0.5, spatial window ΔA, temporal windo ΔT). In (c) and (d) we have rescaled the sample area/duration by the smallest/shortest reported value and coloured by community (see original study for details). In (g) and (h) we study the STAR in metacommunities of various size N, represented by colour. Limited spatio-temporal turnover in the smallest metacommunties (blue colours) greatly reduces the exponents of the STAR relative to large metacommunities (red colours). Aij = 0.5 with probability 0.5 in all cases.Full size imageTemporal occupancyThe proportion of time in which species occupy a community tends to have a bi-modal empirical distribution54,55,56 (Fig. 6a). The distribution we found in simulations (Fig. 6e) closely matches the empirical pattern.Community structureTemporal turnover has been posited to play a stabilising role in the maintenance of community structure57,58. In an estuarine fish community59, for example, species richness (Fig. 6b) and the distribution of abundances were remarkably robust despite changes in population biomasses by multiple orders of magnitude. In model metacommunities with autonomous turnover we found, likewise, that local species richness exhibited only small fluctuations around the steady-state mean (Fig. 6f, three random local communities shown) and that the macroscopic structure of the community was largely time invariant (Supplementary Fig. 12). In the light of our results, we propose the absence of temporal change in community properties such as richness or the abundance distribution despite potentially large fluctuations in population abundances59 as indicative of autonomous compositional turnover.The species-time-area-relation, STARThe species-time-relation (STR), typically fit by a power law of the form S ∝ Tw 12,60,61, describes how observed species richness increases with observation time T. The exponent w of the STR has been found to be consistent across taxonomic groups and ecosystems12,13,62, indicative of some general population dynamical mechanism. However, the exponent of the STR decreases with increasing sampling area12, and the exponent of the empirical Species Area Relation (SAR) (S ∝ Az) consistently decreases with increasing sampling duration12 (Fig. 6c, d). We tested for these patterns in a large simulated metacommunity with N = 256 patches by computing the species-time-area-relation (STAR) for nested subdomains and variable temporal sampling windows (see “Methods”). We observed exponents of the nested SAR in the range z = 0.02–0.44 and for the STR a range w = 0.01–0.44 (Supplementary Fig. 13). We also found a clear decrease in the rate of species accumulation in time as a function of sample area and vice-versa (Fig. 6g, h), consistent with the empirical observations. Meta-analyses of these patterns in nature have reported exponents which are remarkably consistent, with z typically in the range 0.1–0.363, and w typically in the range 0.2–0.413, in both cases largely independent of location or taxonomic group13.Thus, the distribution of temporal occupancy, the time invariance of key macroecological structures and the STAR in our model metacommunities match observed patterns. This evidence suggests that such autonomous dynamics cannot be ruled out as an important driver of temporal compositional change in natural ecosystems.Turnover rate in simulated metacommunitiesHow do the turnover rates that we find in our model compare with those observed? Our current analytic understanding of autonomous turnover is insufficient for estimating the rates directly from parameters, but the simulation results provide some indication of the expected order of magnitude, that can be compared with observations. Key for such a comparison is the fact that, because the elements of R are 1 on average, the time required for an isolated single population to reach carrying capacity is ({mathcal{O}}(1)) unit times. Supplementary Fig. 12b suggests that transitions between community states occur at the scale of around 10–50 unit times. This gives a holistic, rule-of-thumb estimate for the expected rate of autonomous turnover, depending on the typical reproductive rates of the guild of interest. In the case of macroinvertebrates, for example, the time required for populations to saturate in population biomass could be of the order of a month or less. By our rule of thumb, this would mean that autonomous community turnover would occur on a timescale of years. In contrast, for slow growing species like trees, where monoculture stands can take decades to reach maximum population biomass, the predicted timescale for autonomous turnover would be on the order of centuries or more. Indeed, macroinvertebrate communities have been observed switching between community configurations with a period of a few years64,65, while the proportional abundance of tree pollen and tree fern spores fluctuates in rain forest bog deposits with a period of the order of 103 years66—suggesting that the predicted autonomous turnover rates are biologically plausible.ConclusionsCurrent understanding of the mechanisms driving temporal turnover in ecological communities is predominantly built upon phenomenological studies of observed patterns2,67,68,69 and is unquestionably incomplete10,59. That temporal turnover can be driven by external forces—e.g., seasonal or long term climate change, direct anthropogenic pressures—is indisputable. A vitally important question is, however, how much empirically observed compositional change is actually due to such forcing. Recent landmark analyses of temporal patterns in biodiversity have detected no systematic change in species richness or structure in natural communities, despite rates of compositional turnover greater than predicted by stochastic null models1,70,71,72. Here we have shown that empirically realistic turnover in model metacommunities can occur via precisely the same mechanism as that responsible for regulating species richness at the local scale. While the processes regulating diversity in natural communities remain insufficiently understood, our theoretical work suggests local structural instability may explain these empirical observations in a unified and parsimonious way. Therefore, we advocate for the application of null models of metacommunity dynamics that account for natural turnover in ecological status assessments and predictions based on ancestral baselines. Future work will involve fitting the model described here to observations by estimating abiotic and biotic parameters from empirical datasets. In the Supplementary Discussion we show how different combinations of parameters lead to different quantitative outcomes (Supplementary Fig. 4), likely representing different types of empirical metacommunities. Understanding where in this parameter space natural systems exist may provide the foundation for a quantitative null model, a baseline expectation of turnover against which observations can be compared.Our simulations revealed a qualitative transition from “small” metacommunities, where autonomous turnover is absent or minimal, to “large” metacommunities with pronounced autonomous turnover (Fig. 2). The precise location of the transition between these cases depends on details such as dispersal traits, the ecological interaction network, and environmental gradients (Supplementary Fig. 4). Taking, for simplicity, regional species richness as a measure of metacommunity size suggests that both ‘small’ and ‘large’ communities in this sense are realised in nature. In our simulations, the smallest metacommunities sustain 10s of species, while the largest have a regional diversity of the order 103, which is not large comparable to the number of tree species in just 0.25 km2 of tropical rainforest (1100–1200 in Borneo and Ecuador73) or of macroinvertebrates in the UK ( >32,00074). Within the ‘small’ category, where autonomous turnover is absent, we would therefore expect to be, e.g., communities of marine mammals or large fish, where just a few species interact over ranges that can extend across entire climatic niches, implying that the effective number of independent “patches” is small and providing few opportunities for colonisation by species from neighbouring communities. Likely to belong to the ‘large’ category are communities of organisms that occur in high diversity with range sizes that are small compared to climatic niches, such as macroinvertebrates. For these, autonomous turnover of local communities can plausibly be expected based on our findings. Empirically distinguishing between these two cases for different guilds will be an important task for the future.For metacommunities of intermediate spatial extent, autonomous turnover is characterised by sharp transitions between cohesive states at the local scale. To date, few empirical analyses have reported such coherence in temporal turnover, perhaps because the taxonomic and temporal resolution required to detect such patterns is not yet widely available. Developments in biomonitoring technologies75 are likely to reveal a variety of previously undetected ecological dynamics, however and by combining high resolution temporal sampling and metagenetic analysis of community composition, a recent study demonstrated cohesive but short-lived community cohorts in coastal plankton76. Such Clementsian temporal turnover may offer a useful signal of autonomous compositional change in real systems.Thus, overcoming previous computational limits to the study of complex metacommunities11,77, we have discovered the existence of two distinct phases of metacommunity ecology—one characterised by weak or absent autonomous turnover, the other by continuous compositional change even in the absence of external drivers. By synthesising a wide range of established ecological theory11,23,24,47,48,49, we have heuristically explained these phases. Our explanation implies that autonomous turnover requires little more than a diverse neighbourhood of potential invaders, a weak immigration pressure, and a complex network of interactions between co-existing species. More

Power spectral density for a general Ornstein–Uhlenbeck processIn the following we develop a method to compute the power spectral density of N-dimensional Ornstein–Uhlenbeck processes,$$frac{d{boldsymbol{xi }}}{dt}={boldsymbol{A}}{boldsymbol{xi }}+{boldsymbol{zeta }}(t),$$
(15)
where ζ(t) is an N-vector of Gaussian white noise with correlations ({mathbb{E}}[{boldsymbol{zeta }}(t){boldsymbol{zeta }}{(t^{prime} )}^{T}]=delta (t-t^{prime} ){boldsymbol{B}}). The matrix A determines the mean behaviour of ξ and is considered to be locally stable, i.e. all eigenvalues of A have negative real part. Using the matrices A and B one can fully determine the power spectral density of fluctuations for the Ornstein-Uhlenbeck process.We are interested in the case that the coefficients Aij and Bij are derived from a complex network of interactions with weights drawn at random, possibly with correlations. This framework encompasses a very general class of models with a wealth of real-world applications including but not limited to the ecological focus we have here. The method we describe exploits the underlying network structure of A and B to deduce a self-consistent scheme of equations whose solution contains information on the power spectral density.We start with the definition of the power spectral density Φ(ω) as the Fourier transform of the covariance ({mathbb{E}}[{boldsymbol{xi }}(t){boldsymbol{xi }}{(t+tau )}^{T}]) at equilibrium,$${mathbf{Phi }}(omega )=int_{-infty }^{infty }{{rm{e}}}^{-{rm{i}}omega tau }{mathbb{E}}[{boldsymbol{xi }}(t){boldsymbol{xi }}(t+tau )]dtau .$$
(16)
From ref. 33 on multivariate Ornstein–Uhlenbeck processes, we know that the power spectral density can also be written in the form of the matrix equation,$${mathbf{Phi }}(omega )={({boldsymbol{A}}-iomega {boldsymbol{I}})}^{-1}{boldsymbol{B}}{({{boldsymbol{A}}}^{T}+iomega {boldsymbol{I}})}^{-1}.$$
(17)
In practice, this equation is difficult to use for large systems as large matrix inversion is analytically intractable and numerical schemes are slow and sometimes unstable. We take an alternative route by recasting Eq. (17) as a complex Gaussian integral reminiscent of problems appearing in the statistical physics of disordered systems. Our approach in the following is to treat ω as a fixed parameter and drop the explicit dependence from our notation. We begin by writing$${mathbf{Phi }}(omega )=frac{| {boldsymbol{A}}-iomega {boldsymbol{I}}{| }^{2}}{{pi }^{N}| {boldsymbol{B}}| }int_{{mathbb{C}}}{e}^{-{{boldsymbol{u}}}^{dagger }{{boldsymbol{Phi }}}^{-1}{boldsymbol{u}}}{boldsymbol{u}}{{boldsymbol{u}}}^{dagger }mathop{prod }limits_{i=1}^{N}d{u}_{i} .$$
(18)
Simplification of the integrand is achieved by unpicking the matrix inversion in the exponent via a Hubbard-Stratonovich transformation46,47. To this end we recast the system in the language of statistical mechanics by introducing N complex-valued ‘spins’ ui and N auxiliary variables vi, with the ‘Hamiltonian’$${mathcal{H}}({boldsymbol{u}},{boldsymbol{v}})=-{{boldsymbol{u}}}^{dagger }({boldsymbol{A}}-{rm{i}}omega ){boldsymbol{v}}+{{boldsymbol{v}}}^{dagger }{({boldsymbol{A}}-{rm{i}}omega )}^{dagger }{boldsymbol{u}}+{{boldsymbol{v}}}^{dagger }{boldsymbol{B}}{boldsymbol{v}} .$$
(19)
Introducing a bracket operator$$langle cdots rangle :=frac{{int}_{{mathbb{C}}}{e}^{-{mathcal{H}}({boldsymbol{u}},{boldsymbol{v}})}(cdots )d{boldsymbol{u}}d{boldsymbol{v}}}{{int}_{{mathbb{C}}}{e}^{-{mathcal{H}}({boldsymbol{u}},{boldsymbol{v}})}d{boldsymbol{u}}d{boldsymbol{v}}} ,$$
(20)
we can obtain succinct expressions for the power spectral density Φ = 〈uu†〉 as well as the resolvent matrix ({boldsymbol{{mathcal{R}}}}={({rm{i}}omega -{boldsymbol{A}})}^{-1}=langle {boldsymbol{u}}{{boldsymbol{v}}}^{dagger }rangle). Thus we may write,$${mathbf{Phi }}=frac{1}{{mathcal{Z}}}{int}_{{mathbb{C}}}{e}^{-{mathcal{H}}({boldsymbol{u}},{boldsymbol{v}})}{boldsymbol{u}}{{boldsymbol{u}}}^{dagger }mathop{prod }limits_{i=1}^{N}d{u}_{i}d{v}_{i} ,$$
(21)
where ({mathcal{Z}}=| {boldsymbol{A}}-iomega {boldsymbol{I}}{| }^{2}/{pi }^{2N}).This construction may seem laborious at first, but it unlocks a powerful collection of statistical mechanics tools, including the ‘cavity method’. Originally, the cavity method has been introduced in order to analyse a model for spin glass systems48,49. Further applications of the method include the analysis of the eigenvalue distribution in sparse matrices50,51,52. We will exploit the network structure in a similar fashion in order to compute the power spectral density.In our analysis, we find that it is convenient to split the Hamiltonian in Eq. (19) into the sum of its local contributions at sites i, ({{mathcal{H}}}_{i}), and contributions from interactions between i and j, ({{mathcal{H}}}_{ij}),$${mathcal{H}}=mathop{sum}limits_{i}{{mathcal{H}}}_{i}+mathop{sum}limits_{i sim j}{{mathcal{H}}}_{ij} .$$
(22)
These terms can be decomposed as ({{mathcal{H}}}_{i}={{boldsymbol{w}}}_{i}^{dagger }{{boldsymbol{chi }}}_{i}{{boldsymbol{w}}}_{i}) and ({{mathcal{H}}}_{ij}={{boldsymbol{w}}}_{i}^{dagger }{{boldsymbol{chi }}}_{ij}{{boldsymbol{w}}}_{j}), where we introduce the compound spins ({{boldsymbol{w}}}_{i}={({u}_{i},{v}_{i})}^{T}) and transfer matrices,$${{boldsymbol{chi }}}_{i} = , left(begin{array}{ll}0&{A}_{ii}+iomega \ -{A}_{ii}+iomega &{B}_{ii}end{array}right) ,\ {{boldsymbol{chi }}}_{ij} = , left(begin{array}{ll}0&{A}_{ji}\ -{A}_{ij}&{B}_{ij}end{array}right) .$$
(23)
Let us focus on the power spectral density of a particular variable ξi, obtained from the diagonal element ϕi = Φii. For this we compute the single-site marginal fi by integrating over all other variables,$${f}_{i}({{boldsymbol{w}}}_{i})=frac{1}{{mathcal{Z}}}{int}_{{mathbb{C}}}{e}^{-{mathcal{H}}}mathop{prod}limits_{jne i}d{{boldsymbol{w}}}_{j}.$$
(24)
Alternatively, ϕi can be obtained as the top left entry of the covariance matrix ({{mathbf{Psi }}}_{i}=langle {{boldsymbol{w}}}_{i}{{boldsymbol{w}}}_{i}^{dagger }rangle). We write the covariance matrix as the integral,$${{mathbf{Psi }}}_{i}={int}_{{mathbb{C}}}{f}_{i}({{boldsymbol{w}}}_{i}){{boldsymbol{w}}}_{i}{{boldsymbol{w}}}_{i}^{dagger }d{{boldsymbol{w}}}_{i} ,$$
(25)
which could also be expressed in terms of a Gaussian integral,$${{mathbf{Psi }}}_{i}=frac{1}{{pi }^{2}| {{mathbf{Psi }}}_{i}| }{int}_{{mathbb{C}}}{e}^{-{{boldsymbol{w}}}_{i}^{dagger }{{mathbf{Psi }}}_{i}^{-1}{{boldsymbol{w}}}_{i}}{{boldsymbol{w}}}_{i}{{boldsymbol{w}}}_{i}^{dagger }d{{boldsymbol{w}}}_{i} .$$
(26)
By comparing Eqs. (25) and (26) we find that$${f}_{i}({{boldsymbol{w}}}_{i})=frac{1}{{pi }^{2}| {{mathbf{Psi }}}_{i}| }{e}^{-{{boldsymbol{w}}}_{i}^{dagger }{{mathbf{Psi }}}_{i}^{-1}{{boldsymbol{w}}}_{i}} .$$
(27)
We now insert Eq. (22) into Eq. (24) and obtain,$${f}_{i}({{boldsymbol{w}}}_{i})=frac{1}{{pi }^{2}| {{mathbf{Psi }}}_{i}| }{e}^{-{{mathcal{H}}}_{i}}{int}_{{mathbb{C}}}mathop{prod}limits_{i sim j}left({e}^{-{{mathcal{H}}}_{ij}-{{mathcal{H}}}_{ji}}{f}_{j}^{(i)}d{{boldsymbol{w}}}_{j}right) ,$$
(28)
where we write ({f}_{j}^{(i)}) for the ‘cavity marginals’,$${f}_{j}^{(i)}({{boldsymbol{w}}}_{j})=frac{1}{{{mathcal{Z}}}^{(i)}}{int}_{{mathbb{C}}}{e}^{-{{mathcal{H}}}^{(i)}}mathop{prod}limits_{kne i,j}d{{boldsymbol{w}}}_{k} .$$
(29)
In essence, the above discussion amounts to organising the 2N integrals in Eq. (21) in a convenient way, with the advantage of providing a simple intuition for the role of the underlying network. The superscript (i) is used to indicate that the quantity corresponds to the cavity network where node i has been removed. We will further use this notation for the ‘cavity covariance matrix’ ({{mathbf{Psi }}}_{jl}^{(i)}) introduced in the following.Next we perform the integration in Eq. (28) and compare to the form in Eq. (27). We thus obtain a recursion formula for the covariance matrix Ψi and the cavity covariance matrices ({{mathbf{Psi }}}_{jl}^{(i)}),$${{mathbf{Psi }}}_{i}={left({{boldsymbol{chi }}}_{i}-mathop{sum}limits_{{{i !sim! j}atop {i! sim! l}}}{{boldsymbol{chi }}}_{ij}{{mathbf{Psi }}}_{jl}^{(i)}{{boldsymbol{chi }}}_{li}right)}^{-1},$$
(30)
where the notation i ~ j indicates that we sum over nodes j connected to node i. Unless there is some specific structure underlying the network, we assume that most real world cases have a ‘tree-like’ structure from the local view point of a single node i. Hence, it is highly unlikely that the nodes j and l are nearby in the cavity network where node i is removed, and thus ({{mathbf{Psi }}}_{jl}^{(i)}) only gives non-zero contributions if j = l. We therefore reduce Eq. (30) and obtain for the covariance matrix,$${{mathbf{Psi }}}_{i}={left({{boldsymbol{chi }}}_{i}-mathop{sum}limits_{i sim j}{{boldsymbol{chi }}}_{ij}{{mathbf{Psi }}}_{j}^{(i)}{{boldsymbol{chi }}}_{ji}right)}^{-1}.$$
(31)
Similarly, the cavity covariance matrix obeys the equation,$${{mathbf{Psi }}}_{j}^{(i)}={left({{boldsymbol{chi }}}_{j}-mathop{sum}limits_{j sim k,kne i}{{boldsymbol{chi }}}_{jk}{{mathbf{Psi }}}_{k}^{(j)}{{boldsymbol{chi }}}_{kj}right)}^{-1}.$$
(32)
Here we use that Ψ(i, j) = Ψ(j) when the nodes i and k are not connected. In other words, removing node j from the cavity network where node i is missing, has the same effect as removing it from the full network. The system in Eq. (31) describes a collection of nonlinear matrix equations that must be solved self-consistently.For networks with high enough connectivity (and to good approximation even with modest connectivity), the removal of a single node does not affect the rest of the network, as its contribution is negligible compared to the full system. Hence the system in Eq. (31) can be reduced to a smaller set of equations approximately satisfied by the matrices Ψi:$${{mathbf{Psi }}}_{i}approx {left({{boldsymbol{chi }}}_{i}-mathop{sum}limits_{i sim j}{{boldsymbol{chi }}}_{ij}{{mathbf{Psi }}}_{j}{{boldsymbol{chi }}}_{ji}right)}^{-1}.$$
(33)
The power spectral density ϕi can be obtained as the top left entry of Ψi.In order to progress further, we now consider specific approximations that help us compute the power spectral density. First, we take a mean-field approach in order to obtain the mean power spectral density for all nodes part of the network; we then use the result for the mean-field in order to compute a close approximation to the local power spectral density of a single node. Later, we adapt the method to partitioned networks where nodes belong to different types of connected groups.Mean fieldFor the following, we assume that all agents in the system behave the same on average. In practice, the terms governed by self-interactions Aii are drawn from the same distribution for all agents. Similarly, the terms including Bii are governed by one distribution. Interaction strengths and connections with other nodes in the network are also sampled equally for all agents (we have explored a large Lotka-Volterra ecosystem as an example of such a network). In the mean-field (MF) formulation we assume that the mean degree and excess degree are approximately equal, and replace all quantities in Eqs. (31) and (32) with their average. Ψi = ΨMF ∀ i. We then obtain the following recursion equation,$${{mathbf{Psi }}}^{{rm{MF}}}={left[{mathbb{E}}[{{boldsymbol{chi }}}_{i}]-{mathbb{E}}left(mathop{sum}limits_{i sim j}{{boldsymbol{chi }}}_{ij}{{mathbf{Psi }}}^{{rm{MF}}}{{boldsymbol{chi }}}_{ji}right)right]}^{-1}.$$
(34)
In order to solve this equation, we parameterise,$${{mathbf{Psi }}}^{{rm{MF}}}=left(begin{array}{ll}phi &r\ -bar{r}&0end{array}right),$$
(35)
where the top left entry ϕ corresponds to the mean power spectral density, and we introduce r as the mean diagonal element of the resolvent matrix ({boldsymbol{{mathcal{R}}}}). Finally by inserting the ansatz of Eq. (35) into Eq. (34) we obtain,$$left(begin{array}{ll}phi &r\ -bar{r}&0end{array}right)^{-1}= , left(begin{array}{ll}0&{mathbb{E}}[{A}_{ii}]+iomega \ -{mathbb{E}}[{A}_{ii}]+iomega &{mathbb{E}}[{B}_{ii}]end{array}right)\ , +cleft(begin{array}{ll}0&bar{r}{mathbb{E}}[{A}_{ij}{A}_{ji}]\ -r{mathbb{E}}[{A}_{ij}{A}_{ji}]&phi {mathbb{E}}[{A}_{ij}^{2}]+(r+bar{r}){mathbb{E}}[{A}_{ij}{B}_{ij}]end{array}right),$$
(36)
where c is the average degree (i.e. number of connections) per node. Moreover, the expectations in the second term are to be taken over connected nodes i ~ j (i.e. non-zero matrix entries).From Eq. (36) above, we obtain the equations,$$frac{phi }{| r{| }^{2}} = , {mathbb{E}}[{B}_{ii}]+cleft(phi {mathbb{E}}[{A}_{ij}^{2}]+2{rm{Re}}(r){mathbb{E}}[{A}_{ij}{B}_{ij}]right),\ frac{bar{r}}{| r{| }^{2}} = , -!{mathbb{E}}[{A}_{ii}]+iomega -cr{mathbb{E}}[{A}_{ij}{A}_{ji}].$$
(37)
We solve the second equation in Eq. (37) for r and write the mean power spectral density in terms of r,$$phi = , | r{| }^{2}frac{{mathbb{E}}[{B}_{ii}]+2c{rm{Re}}(r){mathbb{E}}[{A}_{ij}{B}_{ij}]}{1-c| r{| }^{2}{mathbb{E}}[{A}_{ij}^{2}]},\ r = , frac{1}{2c{mathbb{E}}[{A}_{ij}{A}_{ji}]}left[-{mathbb{E}}[{A}_{ii}]+iomega right.\ , left.-sqrt{{(-{mathbb{E}}[{A}_{ii}]+iomega )}^{2}-4c{mathbb{E}}[{A}_{ij}{A}_{ji}]}right]$$
(38)
This equation informs the first part of the results presented in the main text.Single defect approximationThe single defect approximation (SDA) makes use of the mean-field approximation for the cavity fields, but retains local information about individual nodes. We parameterise similarly to Eq. (35) for a single individual. Moreover, we replace all other quantities with the respective mean-field approximation. Specifically, we obtain$${left(begin{array}{ll}{phi }_{i}^{{rm{SDA}}}&{r}_{i}^{{rm{SDA}}}\ -{bar{r}}_{i}^{{rm{SDA}}}&0end{array}right)}^{-1}= ,left(begin{array}{ll}0&{A}_{ii}+iomega \ -{A}_{ii}+iomega &{B}_{ii}end{array}right)\ , +mathop{sum}limits_{i sim j}left(begin{array}{ll}0&{bar{r}}^{{rm{MF}}}{A}_{ij}{A}_{ji}\ -{r}^{{rm{MF}}}{A}_{ij}{A}_{ji}&{phi }^{{rm{MF}}}{A}_{ij}^{2}+({r}^{{rm{MF}}}+{bar{r}}^{{rm{MF}}}){A}_{ij}{B}_{ij}end{array}right).$$
(39)
We solve this equation for ({phi }_{i}^{{rm{SDA}}},{r}_{i}^{{rm{SDA}}}), which delivers$$frac{{phi }_{i}^{{rm{SDA}}}}{| {r}_{i}^{{rm{SDA}}}{| }^{2}} = , {phi }^{{rm{MF}}}mathop{sum}limits_{i sim j}{A}_{ij}^{2}+2{rm{Re}}({r}^{{rm{MF}}})mathop{sum}limits_{i sim j}{A}_{ij}{B}_{ij}+{B}_{ii} ,\ {r}_{i}^{{rm{SDA}}} = , {left({A}_{ii}+iomega +{bar{r}}^{{rm{MF}}}mathop{sum}limits_{i sim j}{A}_{ij}{A}_{ji}right)}^{-1}.$$
(40)
Partitioned networkPreviously we assumed that all nodes in a network are interchangeable in distribution. However, many real-world applications feature agents with different properties, imposing a high-level structure on the network. We realise this by partitioning nodes into distinct groups that interact with each other (see the section Trophic structure model for a simple example).In order to handle different connected groups we make use of the cavity method as in Eqs. (31) and (32). In particular, we split the sum in the second term on the right-hand side of these equations into contributions from each group in the partitioned network. Let M denote the number of subgroups Vm in a partitioned network then we write,$${{mathbf{Psi }}}_{i}= , {left({{boldsymbol{chi }}}_{i}-mathop{sum }limits_{m}^{M}mathop{sum}limits_{{{i !sim! j}atop {j!in! {V}_{m}}}}{{boldsymbol{chi }}}_{ij}{{mathbf{Psi }}}_{j}^{(i)}{{boldsymbol{chi }}}_{ji}right)}^{-1} ,\ {{mathbf{Psi }}}_{j}^{(i)} = , {left({{boldsymbol{chi }}}_{j}-mathop{sum }limits_{m}^{M}mathop{sum}limits_{{{j !sim! k}atop {k!in !{V}_{m}}}}{{boldsymbol{chi }}}_{jk}{{mathbf{Psi }}}_{k}^{(j)}{{boldsymbol{chi }}}_{kj}right)}^{-1}.$$
(41)
Similar to the previous sections we replace all quantities with a mean-field average ({{mathbf{Psi }}}_{m}^{{rm{MF}}}), but for each group separately. Hence we obtain M equations of the form$${{mathbf{Psi }}}_{i}^{{rm{MF}}}={left[{mathbb{E}}[{{boldsymbol{chi }}}_{i}]-{mathbb{E}}left(mathop{sum }limits_{m}^{M}mathop{sum}limits_{{{i !sim! j}atop {j!in !{V}_{m}}}}{{boldsymbol{chi }}}_{ij}{{mathbf{Psi }}}_{m}^{{rm{MF}}}{{boldsymbol{chi }}}_{ji}right)right]}^{-1}.$$
(42)
In order to compute the mean power spectral density for different groups separately, we use a parameterisation as in Eq. (35) for each group. Therefore we have,$${{mathbf{Psi }}}_{m}^{{rm{MF}}}=left(begin{array}{ll}{phi }_{m}&{r}_{m}\ -{bar{r}}_{m}&0end{array}right),$$
(43)
for all m = 1, …, M. This delivers 2M equations to solve for all rm and ϕm. Numerically this is straightforward, although algebraically long-winded for the general case. However, the equations simplify for special cases. In the section Trophic structure model we demonstrate this method for a bipartite network where a lack of intra-group interactions simplifies the analysis.Large Lotka-Volterra ecosystemModel descriptionFirst, we define the framework for a general Lotka-Volterra ecosystem with N species and a large but finite system size V ≫ 1. Note that this parameter can be interpreted as a scaling factor for the fluctuation amplitude and thus, larger systems exhibit higher stability and quantitative reliability for our analytic results. Let Xi denote the number of individuals and xi = Xi/V the density of species i = 1, …, N. We start from the following set of reactions that define the underlying stochastic dynamics of the system:$$ , {X}_{i},mathop{to }limits^{{b}_{i}},2{X}_{i} ({rm{birth}})\ , 2{X}_{i},mathop{to }limits^{{R}_{ii}},{X}_{i} ({rm{death}})\ , {X}_{i}+{X}_{j},mathop{to }limits^{{R}_{ij}},left{begin{array}{ll}2{X}_{i}+{X}_{j}&({rm{mutualism}}),hfill\ {X}_{i}&({rm{competition}}),\ 2{X}_{i}&({rm{predation}}).hfillend{array}right.$$
(44)
The self-interactions are governed by the birth rate bi  > 0 and density-dependent mortality rate Rii  > 0. Furthermore, we define three interaction types between species i and j, namely mutualism, competition and predation. In the case of mutualistic interactions, both species benefit from each other, whereas competition means that both species have a higher mortality rate, depending on the density of the other species. For predator-prey pairs, one predator species benefits from the death of a prey species. The predator and prey species are chosen randomly, such that species i is equally likely to be a predator or prey of species j.With probability Pc we assign an interaction rate Rij  > 0 to the species pair (i, j), and with probability 1 − Pc there is no interaction between species i and j (i.e. Rij = 0). In other words, each species has on average c = NPc interaction partners. The reaction rates are considered to be i.i.d. random variables drawn from a half-normal distribution (| {mathcal{N}}(0,{sigma }^{2})|), where we write for the mean reaction rate (mu ={mathbb{E}}[{R}_{ij}]=sigma sqrt{2/pi }) and raw second moment ({sigma }^{2}={mathbb{E}}[{R}_{ij}^{2}]). For each interaction pair, the interaction type is chosen such that the proportion of predator-prey pairs is p ∈ [0, 1], and all non-predator-prey interactions are equally distributed between mutualistic and competitive interactions (i.e. the overall proportion of mutualistic/competitive interactions is 1/2(1 − p)). Lastly, we define the symmetry parameter γ = 1 − 2p, where γ = −1 if all interactions are of predator-prey type (p = 1), and similarly γ = +1 if there are no predator-prey interactions (p = 0). In a mixed case where predator-prey and mutualistic/competitive interactions have equal proportion (p = 1/2), we have γ = 0. Later we will see that γ is equivalent to the correlation of signed interaction strengths.In the limit V → ∞, the dynamics of the species density xi obey the ordinary differential equations,$$frac{d{x}_{i}}{dt}={x}_{i}left({b}_{i}+mathop{sum }limits_{j}^{N}{alpha }_{ij}{x}_{j}right),$$
(45)
where αij are the interaction coefficients with ∣αij∣ = ∣αji∣ = Rij. The signs of the interaction coefficients are determined by the type of interaction between species i and j. For mutualistic interactions we have αij = αji  > 0, and αij = αji  More

Fungal culture and insect rearingThe fungus F. verticillioides was isolated from sugarcane plants and cultivated in potato dextrose (PD) medium (Difco, Sparks, NV, USA) at 25 °C with a 12 h photoperiod in climatic chambers. A. nidulans (A4 strain) was used as a control because it is not involved in red rot disease. It was cultivated in minimal medium (MM) [24] and maintained in climatic chambers at 37 °C in the dark.The D. saccharalis was provided by Prof. Dr. José R. P. Parra from the University of São Paulo, Piracicaba. The caterpillars were fed an artificial diet [25] and maintained in a room under controlled conditions (temperature 25 ± 4 °C, relative humidity 60 ± 10% and 14 h of light). Adults were kept in cages covered with white paper sheets, where the eggs were deposited, collected and sanitized with 1% copper sulfate solution daily. Newly hatched caterpillars were transferred to the artificial diet [25].Olfactory preference assayFive days before the experiment, a total of 105 fungal conidia of F. verticillioides or A. nidulans were inoculated in a Falcon tube (15 mL) containing 7 mL of MM. The negative control was sterile MM. Tubes containing fungus-colonized medium and control medium were placed at opposite ends of the Petri dish (15 cm diameter) bottom, lined with moistened filter paper. A group of ten third-instar D. saccharalis caterpillars was released in the central region of the arena. The choice was quantified in the end of the experiment when the caterpillar remained in the Falcon tube to feed. The medium in the tubes represents a food source, once the caterpillars find it, they remain in the chosen tube. The Petri dishes were closed, sealed and kept in a dark room for 5 h at 25 °C; then, the number of caterpillars inside each tube was recorded. The assay was also performed using third-instar Spodoptera frugiperda, to detect specific attractiveness, and with fifth-instar D. saccharalis, to find changes in insect behavior during different immature stages.To confirm insect attraction to fungal volatiles, VOCs collected from F. verticillioides were used to attract D. saccharalis. This assay was performed as described; however, only the control medium was added to the tubes. The hexane solvent was removed from the samples using nitrogen gas and the fungal VOCs were eluted in mineral oil. In addition to the control medium, each tube contained a piece of cotton loaded with either 50 µL of an aerated sample of F. verticillioides VOCs or solvent control (mineral oil). The dishes were placed in the dark for 7 h at 25 °C. All assays were repeated 10 times. Statistical analyses were performed using t-test (p  More

1.Craik, D. J., Daly, N. L., Bond, T. & Waine, C. Plant cyclotides: A unique family of cyclic and knotted proteins that defines the cyclic cystine knot structural motif. J. Mol. Biol. 294, 1327–1336 (1999).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
2.Gran, L. On the effect of a polypeptide isolated from “Kalata-Kalata” (Oldenlandia affinis DC) on the oestrogen dominated uterus. Acta Pharmacol. Toxicol. (Copenh) 33, 400–408 (1973).CAS 
Article 

Google Scholar 
3.Schoepke, T., Hasan Agha, M. I., Kraft, R., Otto, A. & Hiller, K. Haemolytisch aktive Komponenten aus Viola tricolor L. und Viola arvensis murray. Sci. Pharm. 61, 145–153 (1993).CAS 

Google Scholar 
4.Claeson, P., Göransson, U., Johansson, S., Luijendijk, T. & Bohlin, L. Fractionation protocol for the isolation of polypeptides from plant biomass. J. Nat. Prod. 61, 77–81 (1998).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
5.Göransson, U., Luijendijk, T., Johansson, S., Bohlin, L. & Claeson, P. Seven novel macrocyclic polypeptides from Viola arvensis. J. Nat. Prod. 62, 283–286 (1999).PubMed 
Article 
PubMed Central 

Google Scholar 
6.Poth, A. G. et al. Discovery of cyclotides in the Fabaceae plant family provides new insights into the cyclization, evolution, and distribution of circular proteins. ACS Chem. Biol. 6, 345–355 (2011).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
7.Poth, A. G. et al. Cyclotides associate with leaf vasculature and are the products of a novel precursor in Petunia (Solanaceae). J. Biol. Chem. 287, 27033–27046 (2012).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
8.Burman, R. et al. Distribution of circular proteins in plants: Large-scale mapping of cyclotides in the Violaceae. Front. Plant Sci. 6, 20 (2015).ADS 
Article 

Google Scholar 
9.Hernandez, J. F. et al. Squash trypsin inhibitors from Momordica cochinchinensis exhibit an atypical macrocyclic structure. Biochemistry 39, 5722–5730 (2000).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
10.Nguyen, G. K. T. et al. Discovery of linear cyclotides in monocot plant Panicum laxum of Poaceae family provides new insights into evolution and distribution of cyclotides in plants. J. Biol. Chem. 288, 3370–3380 (2013).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
11.Saether, O. et al. Elucidation of the primary and three-dimensional structure of the uterotonic polypeptide kalata B1. Biochemistry 34, 4147–4158 (1995).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
12.Ravipati, A. S. et al. Understanding the diversity and distribution of cyclotides from plants of varied genetic origin. J. Nat. Prod. 80, 1522–1530 (2017).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
13.Gruber, C. W. et al. Distribution and evolution of circular miniproteins in flowering plants. Plant Cell 20, 2471–2483 (2008).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
14.Simonsen, S. M. et al. A continent of plant defense peptide diversity: Cyclotides in Australian Hybanthus (Violaceae). Plant Cell 17, 3176–3189 (2005).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
15.Slazak, B., Jacobsson, E., Kuta, E. & Göransson, U. Exogenous plant hormones and cyclotide expression in Viola uliginosa (Violaceae). Phytochemistry 117, 527–536 (2015).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
16.Lindholm, P. et al. Cyclotides: A novel type of cytotoxic agents. Mol. Cancer Ther. 1, 365–369 (2002).CAS 
PubMed 
PubMed Central 

Google Scholar 
17.Ovesen, R. G. et al. Biomedicine in the environment: Cyclotides constitute potent natural toxins in plants and soil bacteria. Environ. Toxicol. Chem. 30, 1190–1196 (2011).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
18.Pränting, M., Lööv, C., Burman, R., Göransson, U. & Andersson, D. I. The cyclotide cycloviolacin O2 from Viola odorata has potent bactericidal activity against Gram-negative bacteria. J. Antimicrob. Chemother. 65, 1964–1971 (2010).PubMed 
Article 
CAS 
PubMed Central 

Google Scholar 
19.Tam, J. P., Lu, Y. A., Yang, J. L. & Chiu, K. W. An unusual structural motif of antimicrobial peptides containing end-to-end macrocycle and cystine-knot disulfides. Proc. Natl. Acad. Sci. USA 96, 8913–8918 (1999).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
20.Slazak, B. et al. How Does the sweet violet (Viola odorata L.) fight pathogens and pests—cyclotides as a comprehensive plant host defense system. Front. Plant Sci. 9, 20 (2018).Article 

Google Scholar 
21.Colgrave, M. L. et al. Anthelmintic activity of cyclotides: In vitro studies with canine and human hookworms. Acta Trop. 109, 163–166 (2009).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
22.Jennings, C., West, J., Waine, C., Craik, D. & Anderson, M. A. Biosynthesis and insecticidal properties of plant cyclotides: The cyclic knotted proteins from Oldenlandia affinis. Proc. Natl. Acad. Sci. USA. 98, 10614–10619 (2001).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
23.Gilding, E. K. et al. Gene coevolution and regulation lock cyclic plant defence peptides to their targets. New Phytol. 210, 717–730 (2016).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
24.Mylne, J. S., Wang, C. K., van der Weerden, N. L. & Craik, D. J. Cyclotides are a component of the innate defense of Oldenlandia affinis. Biopolymers 94, 635–646 (2010).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
25.Dörnenburg, H. Cyclotide synthesis and supply: From plant to bioprocess. Biopolymers 94, 602–610 (2010).PubMed 
Article 
CAS 
PubMed Central 

Google Scholar 
26.Trabi, M. et al. Variations in cyclotide expression in Viola species. J. Nat. Prod. 67, 806–810 (2004).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
27.Lista de especies silvestres de Canarias (hongos, plantas y animales terrestres). (Consejería de Política Territorial y Medio Ambiente. Gobierno de Canarias., 2001).28.Myers, N., Mittermeier, R. A., Mittermeier, C. G., da Fonseca, G. A. B. & Kent, J. Biodiversity hotspots for conservation priorities. Nature 403, 853–858 (2000).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
29.Gómez, M. V. M., Esquivel, J. L. M., Díaz, J. R. D. & Izquierdo, M. S. Viola guaxarensis (Violaceae): A new Viola from Tenerife, Canary Islands, Spain. Willdenowia 50, 13–21 (2020).Article 

Google Scholar 
30.Rodríguez-Rodríguez, P., De Castro, A. G. F., Seguí, J., Traveset, A. & Sosa, P. A. Alpine species in dynamic insular ecosystems through time: Conservation genetics and niche shift estimates of the endemic and vulnerable Viola cheiranthifolia. Ann. Bot. 123, 505–519 (2019).PubMed 
Article 
PubMed Central 

Google Scholar 
31.Ireland, D. C., Colgrave, M. L. & Craik, D. J. A novel suite of cyclotides from Viola odorata: Sequence variation and the implications for structure, function and stability. Biochem. J. 400, 1–12 (2006).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
32.Burman, R., Gunasekera, S., Strömstedt, A. A. & Göransson, U. Chemistry and biology of cyclotides: Circular plant peptides outside the box. J. Nat. Prod. 77, 724–736 (2014).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
33.Trabi, M. & Craik, D. J. Tissue-specific expression of head-to-tail cyclized miniproteins in Violaceae and structure determination of the root cyclotide Viola hederacea root cyclotide1. Plant Cell 16, 2204–2216 (2004).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
34.Ballard, H. E., Sytsma, K. J. & Kowal, R. R. Shrinking the violets: Phylogenetic relationships of infrageneric groups in Viola (Violaceae) based on internal transcribed spacer DNA sequences. Syst. Bot. 23, 439 (1998).Article 

Google Scholar 
35.Batista, F. & Sosa, P. A. Allozyme diversity in natural populations of Viola palmensis. Webb & Berth (Violaceae) from La Palma (Canary Islands): Implications for conservation genetics. Ann. Bot. 90, 725–733 (2002).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
36.Marcussen, T., Heier, L., Brysting, A. K., Oxelman, B. & Jakobsen, K. S. From gene trees to a dated allopolyploid network: Insights from the angiosperm genus Viola (Violaceae). Syst. Biol. 64, 84–101 (2014).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
37.Marcussen, T., Oxelman, B., Skog, A. & Jakobsen, K. S. Evolution of plant RNA polymerase IV/V genes: Evidence of subneofunctionalization of duplicated NRPD2/NRPE2-like paralogs in Viola (Violaceae). BMC Evol. Biol. 10, 45 (2010).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
38.Gilli, A. Viola anagae Gilli sp. Nov.. Feddes Repert. 89, 595–596 (1979).Article 

Google Scholar 
39.Moreno-Saiz, J. Lista Roja 2008 de la Flora Vascular Española (Dirección General de Medio Natural y Política Forestal, Ministerio de Medio Ambiente, y Medio Rural y Marino, y Sociedad Española de Biología de la Conservación de Plantas, 2008).
Google Scholar 
40.Broussalis, A. M. et al. First cyclotide from Hybanthus (Violaceae). Phytochemistry 58, 47–51 (2001).41.Mulvenna, J. P., Wang, C. & Craik, D. J. CyBase: A database of cyclic protein sequence and structure. Nucleic Acids Res. 34, D192–D194 (2006).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
42.Hellinger, R. et al. Peptidomics of circular cysteine-rich plant peptides—analysis of the diversity of cyclotides from Viola tricolor by transcriptome- and proteome-mining. J. Proteome Res. https://doi.org/10.1021/acs.jproteome.5b00681 (2015).Article 
PubMed 
PubMed Central 

Google Scholar 
43.Slazak, B., Haugmo, T., Badyra, B. & Göransson, U. The life cycle of cyclotides: Biosynthesis and turnover in plant cells. Plant Cell Rep. 39, 1359–1367 (2020).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
44.Colgrave, M. L., Jones, A. & Craik, D. J. Peptide quantification by matrix-assisted laser desorption ionisation time-of-flight mass spectrometry: Investigations of the cyclotide kalata B1 in biological fluids. J. Chromatogr. A 1091, 187–193 (2005).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
45.Marcussen, T. Allozymic variation in the widespread and cultivated Viola odorata (Violaceae) in western Eurasia. Bot. J. Linn. Soc. 151, 563–571 (2006).Article 

Google Scholar 
46.Källback, P., Nilsson, A., Shariatgorji, M. & Andrén, P. E. msIQuant—quantitation software for mass spectrometry imaging enabling fast access, visualization, and analysis of large data sets. Anal. Chem. 88, 4346–4353 (2016).PubMed 
Article 
CAS 
PubMed Central 

Google Scholar 
47.Pohlert, T. PMCMRplus: Calculate Pairwise Multiple Comparisons of Mean Rank Sums Extended.48.Wickham, H. ggplot2: Elegant Graphics for Data Analysis. Media (Springer, 2009). https://doi.org/10.1007/978-0-387-98141-3.Book 
MATH 

Google Scholar 
49.R Development Core Team, R. R A Language and Environment for Statistical Computing, Vol 1 409 (R Foundation for Statistical Computing, 2011).
Google Scholar 
50.Package, T. Package ‘ PMCMRplus ’ R topics documented (2019).51.Kolde, R. pheatmap: Pretty Heatmaps. R package version 1.0.12. (2019). https://cran.r-project.org/package=pheatmap.52.Grabherr, M. G. et al. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat. Biotechnol. 29, 644–652 (2011).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
53.Sigrist, C. J. A. et al. PROSITE: A documented database using patterns and profiles as motif descriptors. Brief. Bioinform. 3, 265–274 (2002).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
54.Rice, P., Longden, I. & Bleasby, A. EMBOSS: The European molecular biology open software suite. Trends Genet. 16, 276–277 (2000).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
55.Sievers, F. et al. Fast, scalable generation of high-quality protein multiple sequence alignments using Clustal Omega. Mol. Syst. Biol. 7, 539 (2011).PubMed 
PubMed Central 
Article 

Google Scholar 
56.Burman, R. et al. Cyclotide proteins and precursors from the genus Gloeospermum: Filling a blank spot in the cyclotide map of Violaceae. Phytochemistry 71, 13–20 (2010).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
57.Levenfors, J. J., Hedman, R., Thaning, C., Gerhardson, B. & Welch, C. J. Broad-spectrum antifungal metabolites produced by the soil bacterium Serratia plymuthica A 153. Soil Biol. Biochem. 36, 677–685 (2004).CAS 
Article 

Google Scholar 
58.Broekaert, W. F., Terras, R. F. G., Cammue, B. P. A. & Vandedeyden, J. An automated quantitative assay for fungal growth inhibition. Most 69, 20 (1990).
Google Scholar 
59.CLSI. M38–A2 reference method for broth dilution antifungal susceptibility testing of filamentous fungi; approved standard—second edition. Clin. Lab. Stand. Inst. 20, 20 (2008).
Google Scholar  More

Data collectionWe extracted hourly air temperature and SH from the North America Land Data Assimilation System project46, a near real-time dataset with a 0.125° × 0.125° grid resolution. We spatially and temporally averaged these data into daily county-level records. SH is the mass of water vapor in a unit mass of moist air (g kg−1). Daily downward UV radiation at the surface, with a wavelength of 0.20–0.44 µm, was extracted from the European Centre for Medium-Range Weather Forecasts ERA5 climate reanalysis47.Other characteristics of each county, including geographic location, population density, demographic structure of the population, socioeconomic factors, proportion of healthcare workers, intensive care unit (ICU) bed capacity, health risk factors, long-term and short-term air pollution, and climate zone were collected from multiple sources. Geographic coordinates, population density, median household income, percent of people older than 60 years, percent Black residents, percent Hispanic residents, percent owner-occupied housing, percent residents aged 25 years and over without a high school diploma, and percent healthcare practitioners or support staff were collected from the U.S. Census Bureau48. Total ICU beds in each county were derived from Kaiser Health News49. The prevalence of smoking and obesity among adults in each county was obtained from the Robert Wood Johnson Foundation’s 2020 County Health Rankings50. We extracted annual PM2.5 concentrations in the U.S. from 2014 to 2018 from the 0.01° × 0.01° grid resolution PM2.5 estimation provided by the Atmospheric Composition Analysis Group51, and calculated average PM2.5 levels during this 5-year period for each county to represent long-term PM2.5 exposure (Supplementary Fig. 5). Short-term air quality data during the study period, including daily mean PM2.5 and daily maximum 8-h O3, were obtained from the United States Environmental Protection Agency52. We categorized study counties into one of five climate zones based on the guide released by U.S. Department of Energy53 (Supplementary Fig. 6).The county-level COVID-19 case and death data were downloaded from the John Hopkins University Coronavirus Resource Center1. The U.S. county-to-county commuting data were available from the U.S. Census Bureau48. Daily numbers of inter-county visitors to points of interest (POI) were provided by SafeGraph54.Data ethicsSafeGraph utilizes data from mobile applications of which users optionally consent to provide their anonymous location data.Estimation of reproduction numberWe estimated the daily reproduction number (Rt) in all 3142 U.S. counties using a dynamic metapopulation model informed by human mobility data31,55. Rt is the mean number of new infections caused by a single infected person, given the public health measures in place, in a population in which everyone is assumed to be susceptible. In the metapopulation model, two types of movement were considered: daily work commuting and random movement. During the daytime, some commuters travel to a county other than their county of residence, where they work and mix with the populations of that county; after work, they return home and mix with individuals in their home, residential county. Apart from regular commuting, a fraction of the population in each county, assumed to be proportional to the number of inter-county commuters, travels for purposes other than work. As the population present in each county is different during daytime and night-time, we modelled the transmission dynamics of COVID-19 separately for these two time periods, each depicted by a set of ordinary differential equations (Supplementary Notes).To account for case underreporting, we explicitly simulated reported and unreported infections, for which separate transmission rates were defined. Recent studies from several countries indicate that asymptomatic cases of COVID-19, which are typically unreported, are less contagious than symptomatic cases56,57,58,59. Studies on the early transmission of SARS-CoV-2 in China18 and the U.S.60 also showed that undocumented infections are less transmissible than documented infections.In order to reflect the spatiotemporal variation of disease transmission rate and reporting, we allowed transmission rates and ascertainment rates to vary across counties and to change over time. The transmission model simulated daily confirmed cases and deaths for each county. To map infections to deaths, we used an age-stratified infection fatality rate (IFR)61 and computed the weekly IFR for each county as a weighted average using state-level age structure of confirmed cases reported by the U.S. Centers for Disease Control and Prevention. We further adjusted for reporting lags using an observational delay model informed by a U.S. line-list COVID-19 data record62.For the period prior to March 15, 2020, we used commuting data from the U.S. census survey to prescribe the inter-county movement in the transmission model48. Starting March 15, the census survey data are no longer representative due to changes in mobility behavior following the implementation of non-pharmaceutical interventions. We, therefore, used estimates of the reduction of inter-county visitors to POI (e.g., restaurants, stores, etc.) from SafeGraph54 to account for the change in inter-county movement on a county-by-county basis. Because there is no direct relationship between population-level mobility patterns and COVID-19 transmission rates63, we did not model local transmission rate as a function of inter-county mobility. Instead, the SafeGraph data were only used to inform the change of population mixing across counties.To infer key epidemiological parameters, we fitted the transmission model to county-level daily cases and deaths reported from March 15, 2020 to December 31, 2020. The estimated reproduction number was computed as follows:$${R}_{t}=beta Dleft[alpha +left(1-alpha right)mu right],$$
(1)
where β is the county-specific transmission rate, μ is the relative transmissibility of unreported infections, α is the county-specific ascertainment rate, and D is the average duration of infectiousness. Note (beta) and (alpha) were defined for each county separately and were allowed to vary over time. Unlike previous studies using effective reproduction number$${R}_{e}=beta Dleft[alpha +left(1-alpha right)mu right]s,$$
(2)
where s is the estimated local population susceptibility, we used reproduction number Rt to exclude the influence of population susceptibility on disease transmission rate.D, (mu), (Z) (the average latency period from infection to contagiousness), and a multiplicative factor adjusting random movement ((theta)) were randomly drawn from the posterior distributions inferred from case data through March 13, 202060: (D=3.56) (3.21–3.83), (mu =0.64) (0.56–0.70), (Z=3.59) (95% CI: 3.28–3.99), and (theta =0.15) (0.12–0.17). (Z) and (theta) are used in ordinary differential equations used to model transmission dynamics (Supplementary Notes).The daily transmission rate (beta) and ascertainment rate (alpha) were estimated sequentially for each county using the ensemble adjustment Kalman filter (EAKF)64. Specifically, parameters ({beta }_{i}) and ({alpha }_{i}) for county (i) were updated each day using incidence and death data. We used the estimates on day (t-1) as the prior parameters on day (t), and then updated the priors to posteriors using the EAKF and observations. The posteriors are the estimated parameter values on day (t). To ensure a smooth parameter estimation, we imposed a (pm 30 %) limit on the daily change of parameters ({beta }_{i}) and ({alpha }_{i}). Other smoothing constraints were tested and the results were similar. To avoid possible inaccurate estimation for counties with few cases, we inferred Rt in the 2669 U.S. counties with at least 400 cumulative confirmed cases as of December 31, 2020 (Supplementary Fig. 7).Statistical analysisAll statistical analyses were conducted with R software (version 3.6.1) using the mgcv and dlnm packages.Association between meteorological factors and R
t
Given the potential non-linear and temporally delayed effects of meteorological factors, a distributed lag non-linear model65 combined with generalized additive mixed models66 was applied to estimate the associations of daily mean temperature, daily mean SH, and daily mean UV radiation with SARS-CoV-2 Rt. To quantify the total contribution, independent effects, and relative importance of meteorological factors (i.e., temperature, SH, and UV radiation), we included all three variables in the same model. To reduce collinearity, we used cross-basis terms rather than the raw variables (Supplementary Tables 5–6). The full model can be expressed as:$$log (E({{{R}}}_{i,j,t}))= alpha +te(s({{rm{latitude}}}_{i}{,{rm{longitude}}}_{i},{rm{k}}=200),s({{rm{time}}}_{t},{rm{k}}=30))+{rm{cb}}.{rm{temperature}}+{rm{cb}}.{rm{SH}}+ {rm{cb}}.{rm{UV}}\ +{beta }_{1}({rm{population}},{rm{density}}_{i})+{beta }_{2}({rm{percent}},{rm{Black}},{rm{residents}}_{i})+{beta }_{3}({rm{percent}},{rm{Hispanic}},{rm{residents}}_{i})\ +{beta }_{4}({rm{percent}},{rm{people}},{rm{older}},{rm{than}},60,{rm{years}}_{i})+{beta }_{5}({rm{median}},{rm{household}},{rm{income}}_{i})\ +{beta }_{6}({rm{percent}},{rm{owner}}-{rm{occupied}},{rm{housing}}_{i})\ +{beta }_{7}({rm{percent}},{rm{residents}},{rm{older}},{rm{than}},25,{rm{years}},{rm{without}},{rm{a}},{rm{high}},{rm{school}},{rm{diploma}}_{i})\ +{beta }_{8}({rm{number}},{rm{of}},{rm{ICU}},{rm{beds}},{rm{per}},10,000,{rm{people}}_{i})+{beta }_{9}({rm{percent}},{rm{healthcare}},{rm{workers}}_{i})\ quad , {beta }_{10}({rm{day}},{rm{when}},100,{rm{cumulative}},{rm{cases}},{rm{per}},100,000,{rm{people}},{rm{was}},{rm{reached}}_{i})+{re}({rm{county}}_{i})+{re}({rm{state}}_{j})$$
(3)
where E(Ri,j,t) refers to the expected Rt in county i, state j, on day t, and α is the intercept. Given the distribution of Rt in our data close to a lognormal distribution (Supplementary Fig. 8), we used log-transformed Rt as the outcome variable, and the Gaussian family in the model. A thin plate spline with a maximum of 200 knots was used to control the coordinates of the centroid of each county; the time trend was controlled by a flexible natural cubic spline over the range of study dates with a maximum of 30 knots; due to the unique pattern of the non-linear time trend of Rt in each county (Supplementary Fig. 4), we constructed tensor product smooths (te) of the splines of geographical coordinates and time, to better control for the temporal and spatial variations (Supplementary Fig. 3).Cb.temperature, cb.SH, and cb.UV are cross-basis terms for the mean air temperature, mean SH and mean UV radiation, respectively. We modeled exposure-response associations (meteorological factors vs. percent change in Rt) using a natural cubic spline with 3 degrees of freedom (df) and modeled the lag-response association using a natural cubic spline with an intercept and 3 df with a maximum lag of 13 days. We adjusted for county-level characteristics, including population density, percent Black residents, percent Hispanic residents, percent people older than 60 years, median household income, percent owner-occupied housing, percent residents older than 25 years without a high school diploma, number of ICU beds per 10,000 people, and percent healthcare workers, given their potential relationship with SARS-CoV-2 transmission67,68,69,70. Day when 100 cumulative cases per 100,000 people was reached in each county was used to approximate local epidemic stage45 (Supplementary Fig. 9). The random effects of state and county were modeled by parametric terms penalized by a ridge penalty (re), to further control for unmeasured state- and county-level confounding. Residual plots were used to diagnose the model (Supplementary Fig. 10). In additional analyses, we included air temperature, SH, and UV radiation in separate models (Supplementary Fig. 2).Based on the estimated exposure-response curves, between the 1st and the 99th percentiles of the distribution of air temperature, SH, and UV radiation, we determined the value of exposure associated with the lowest relative risk of Rt to be the optimum temperature, the optimum SH, or the optimum UV radiation, respectively. The natural cubic spline functions of the exposure-response relationship were then re-centered with the optimum values of meteorological factors as reference values. We report the cumulative relative risk of Rt associated with daily temperature, SH, or UV radiation exposure in the previous two weeks (0– 13 lag days) as the percent changes in Rt when comparing the daily exposure with the optimum reference values (i.e., the cumulative relative risk of Rt equals one and the percent change in Rt equals zero when the temperature, SH, or UV radiation exposure is at its optimum value).Attribution of R
t to meteorological factorsWe used the optimum value of temperature, SH, or UV radiation as the reference value for calculating the fraction of Rt attributable to each meteorological factor; i.e., the attributable fraction (AF). For these calculations, we assumed that the associations of meteorological factors with Rt were consistent across the counties. For each day in each county, based on the cumulative lagged effect (cumulative relative risk) corresponding to the temperature, SH, or UV radiation of that day, we calculated the attributable Rt in the current and next 13 days, using a previously established method71. Specifically, in a given county, the Rt attributable to a meteorological factor (xt) for a given day t was defined as the attributable absolute excess of Rt (AEx,t, the excess reproduction number on day t attributable to the deviation of temperature or SH from the optimum value) and the attributable fraction of Rt (AFx,, the fraction of Rt attributable to the deviation of the meteorological factor from its optimum value), each accumulated over the current and next 13 days. The formulas can be expressed as:$${{AF}}_{x,t}=1-{rm{exp }}left(-mathop{sum }limits_{l=0}^{13}{beta }_{{x}_{t},l}right)$$
(4)
$${{AE}}_{x,t}={{AF}}_{x,t}times mathop{sum }limits_{l=0}^{13}frac{{n}_{t+1}}{13+1},$$
(5)
where nt is the Rt on day t, and ({sum }_{l=0}^{13}{beta }_{{x}_{t},l}) is the overall cumulative log-relative risk for exposure xt on day t obtained by the exposure-response curves re-centered on the optimum values. Then, the total absolute excess of Rt attributable to temperature, SH, or UV radiation in each county was calculated by summing the absolute excesses of all days during the study period, and the attributable fraction was calculated by dividing the total absolute excess of Rt for the county by the sum of the Rt of all days during the study period for the county. The attributable fraction for the 2669 counties combined was calculated in a similar manner at the national level. We derived the 95% eCI for the attributable absolute excess and attributable fraction by 1000 Monte Carlo simulations71. The total fraction of Rt attributable to meteorological factors was the sum of the attributable fraction for temperature, SH, and UV radiation. We also calculated the attributable fractions by month in the study period.Sensitivity analysesWe conducted several sensitivity analyses to test the robustness of our results: (a) the lag dimension was redefined using a natural cubic spline and three equally placed internal knots in the log scale; (b) an alternative four df was used in the cross-basis term for meteorological factors in the exposure-response function; (c) the maximum number of knots was reduced to 25 in the flexible natural cubic spline to control time trend in the tensor product smooths; (d) all demographic and socioeconomic variables were excluded from the model; (e) adjustment for the prevalence of smoking and obesity among adults was included in the model; (f) adjustment for climate zone was included in the model; (g) additional adjustment was made for the average PM2.5 concentration in each county during 2014–201845; (h) additional adjustment was made for daily mean PM2.5, and daily maximum 8-h O3. For daily covariates with available data in only some of the counties or study period, the results of sensitivity analyses were compared to the main model re-run on the same partial dataset.Reporting summaryFurther information on research design is available in the Nature Research Reporting Summary linked to this article. More

1.Crick FHC, Barnett FRSL, Brenner S, Watts-Tobin RJ. General Nature of the Genetic Code for Proteins. Nature. 1961;192:1227–32.CAS 
PubMed 
Article 

Google Scholar 
2.Hershey AD, Chase M. Independent functions of viral protein and nucleic acid in growth of bacteriophage. J Gen Physiol. 1952;36:39–56.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
3.Luria S, Delbrück M. Mutations of Bacteria from Virus Sensitivity to Virus Resistance. Genetics. 1943;28:491–511.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
4.Kortright KE, Chan BK, Koff JL, Turner PE. Phage Therapy: a Renewed Approach to Combat Antibiotic-Resistant Bacteria. Cell Host Microbe. 2019;25:219–32.CAS 
PubMed 
Article 

Google Scholar 
5.Mushegian AR. Are there 10^31 virus particles on Earth, or more, or less? J Bacteriol. 2020;202:e00052–20.PubMed 
PubMed Central 
Article 

Google Scholar 
6.Dennehy JJ. What Can Phages Tell Us about Host-Pathogen Coevolution? Int J Evol Biol. 2012;2012:1–12.Article 

Google Scholar 
7.Jessup CM, Kassen R, Forde SE, Kerr B, Buckling A, Rainey PB, et al. Big questions, small worlds: microbial model systems in ecology. Trends Ecol Evol. 2004;19:189–97.PubMed 
Article 

Google Scholar 
8.Tecon R, Mitri S, Ciccarese D, Or D, Meer JR, van der, Johnson DR. Bridging the Holistic-Reductionist Divide in Microbial Ecology. MSystems. 2019;4:e00265–18.PubMed 
PubMed Central 
Article 

Google Scholar 
9.Bohannan BJM, Lenski RE. Linking genetic change to community evolution: insights from studies of bacteria and bacteriophage. Ecol Lett. 2000;3:362–77.Article 

Google Scholar 
10.Buckling A, Brockhurst MA. Bacteria-Virus Coevolution. In: Orkun S Soyer, editor. Evolutionary Systems Biology. 2012. New York, NY: Springer; 2012. p. 347–70.11.Koskella B, Brockhurst MA. Bacteria-phage coevolution as a driver of ecological and evolutionary processes in microbial communities. FEMS Microbiol Rev. 2014;38:1–16.Article 
CAS 

Google Scholar 
12.De Sordi L, Lourenço M, Debarbieux L. The Battle Within: interactions of Bacteriophages and Bacteria in the Gastrointestinal Tract. Cell Host Microbe. 2019;25:210–8.PubMed 
Article 
CAS 

Google Scholar 
13.Scanlan PD. Bacteria–Bacteriophage Coevolution in the Human Gut: implications for Microbial Diversity and Functionality. Trends Microbiol. 2017;25:614–23.CAS 
PubMed 
Article 

Google Scholar 
14.Breitbart M. Marine viruses: truth or dare. Annu Rev Mar Sci. 2012;4:425–48.Article 

Google Scholar 
15.Pratama AA, van Elsas JD. The ‘neglected’ soil virome–potential role and impact. Trends Microbiol. 2018;26:649–62.CAS 
PubMed 
Article 

Google Scholar 
16.Lourenço M, De Sordi L, Debarbieux L. The diversity of bacterial lifestyles hampers bacteriophage tenacity. Viruses. 2018;10:1–11.Article 
CAS 

Google Scholar 
17.Martiny JBH, Riemann L, Marston MF, Middelboe M. Antagonistic Coevolution of Marine Planktonic Viruses and Their Hosts. Annu Rev Mar Sci. 2014;6:393–414.Article 

Google Scholar 
18.Díaz-Muñoz SL, Koskella B. Bacteria–Phage Interactions in Natural Environments. In: Sariaslani S, Gadd GM, editors. Advances in Applied Microbiology. Cambridge, MA:Academic Press; 2014. p.135–83.19.Avrani S, Schwartz DA, Lindell D. Virus-host swinging party in the oceans. Mob Genet Elem. 2012;2:88–95.Article 

Google Scholar 
20.Winter C, Bouvier T, Weinbauer MG, Thingstad TF. Trade-Offs between Competition and Defense Specialists among Unicellular Planktonic Organisms: the “Killing the Winner” Hypothesis Revisited. Microbiol Mol Biol Rev. 2010;74:42–57.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
21.Hansen MF, Svenningsen SL, Røder HL, Middelboe M, Burmølle M. Big Impact of the Tiny: bacteriophage–bacteria Interactions in Biofilms. Trends Microbiol. 2019;27:739–52.CAS 
PubMed 
Article 

Google Scholar 
22.O’Brien S, Hodgson DJ, Buckling A. The interplay between microevolution and community structure in microbial populations. Curr Opin Biotechnol. 2013;24:821–5.PubMed 
Article 
CAS 

Google Scholar 
23.Brockhurst MA, Koskella B. Experimental coevolution of species interactions. Trends Ecol Evol. 2013;28:367–75.PubMed 
Article 

Google Scholar 
24.Geredew Kifelew L, Mitchell JG, Speck P. Mini-review: efficacy of lytic bacteriophages on multispecies biofilms. Biofouling. 2019;35:472–81.CAS 
PubMed 
Article 

Google Scholar 
25.Miki T, Jacquet S. Complex interactions in the microbial world: Underexplored key links between viruses, bacteria and protozoan grazers in aquatic environments. Aquat Micro Ecol. 2008;51:195–208.Article 

Google Scholar 
26.Johnke J, Cohen Y, de Leeuw M, Kushmaro A, Jurkevitch E, Chatzinotas A. Multiple micro-predators controlling bacterial communities in the environment. Curr Opin Biotechnol. 2014;27:185–90.CAS 
PubMed 
Article 

Google Scholar 
27.Hall AR, Ashby B, Bascompte J, King KC. Measuring Coevolutionary Dynamics in Species-Rich Communities. Trends Ecol Evol. 2020;35:539–50.PubMed 
Article 

Google Scholar 
28.Strauss SY. Ecological and evolutionary responses in complex communities: implications for invasions and eco-evolutionary feedbacks. Oikos. 2014;123:257–66.Article 

Google Scholar 
29.Strauss SY, Irwin RE. Ecological and evolutionary consequences of multispecies plant-animal interactions. Annu Rev Ecol Evol Syst. 2004;35:435–66.Article 

Google Scholar 
30.Inouye B, Stinchcombe JR. Relationships between ecological interaction modifications and diffuse coevolution: similarities, differences, and causal links. Oikos. 2011;95:353–60.Article 

Google Scholar 
31.Barraclough TG. How Do Species Interactions Affect Evolutionary Dynamics Across Whole Communities? Annu Rev Ecol Evol Syst. 2015;46:25–48.Article 

Google Scholar 
32.Bottery MJ, Pitchford JW, Friman V-P. Ecology and evolution of antimicrobial resistance in bacterial communities. ISME J. 2021;15:939–48.PubMed 
Article 

Google Scholar 
33.Gómez P, Bennie J, Gaston KJ, Buckling A. The Impact of Resource Availability on Bacterial Resistance to Phages in Soil. PLoS ONE. 2015;10:e0123752.PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
34.Gorter FA, Scanlan PD, Buckling A. Adaptation to abiotic conditions drives local adaptation in bacteria and viruses coevolving in heterogeneous environments. Biol Lett. 2016;12:20150879.PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
35.Scanlan JG, Hall AR, Scanlan PD. Impact of bile salts on coevolutionary dynamics between the gut bacterium Escherichia coli and its lytic phage PP01. Infect Genet Evol. 2019;73:425–32.CAS 
PubMed 
Article 

Google Scholar 
36.Gómez P, Buckling A. Bacteria-phage antagonistic coevolution in soil. Science. 2011;332:106–9.PubMed 
Article 
CAS 

Google Scholar 
37.Weinbauer MG, Rassoulzadegan F. Are viruses driving microbial diversification and diversity? Environ Microbiol. 2004;6:1–11.PubMed 
Article 

Google Scholar 
38.Johnke J, Baron M, de Leeuw M, Kushmaro A, Jurkevitch E, Harms H, et al. A generalist protist predator enables coexistence in multitrophic predator-prey systems containing a phage and the bacterial predator Bdellovibrio. Front Ecol Evol. 2017;5:1–12.Article 

Google Scholar 
39.R Core Team. R: a Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing; 2020.40.Mumford R, Friman VP. Bacterial competition and quorum-sensing signalling shape the eco-evolutionary outcomes of model in vitro phage therapy. Evol Appl. 2017;10:161–9.CAS 
PubMed 
Article 

Google Scholar 
41.Connell JH. The influence of interspecific competition and other factors on the distribution of the barnacle Chthamalus stellatus. Ecology. 1961;42:710–23.Article 

Google Scholar 
42.Vellend M. Conceptual Synthesis in Community Ecology. Q Rev Biol. 2010;85:183–206.PubMed 
Article 

Google Scholar 
43.Alseth EO, Pursey E, Lujan AM, McLeod I, Rollie C, Westra ER. Bacterial biodiversity drives the evolution of CRISPR-based phage resistance in Pseudomonas aeruginosa. Nature. 2019;574:549–74.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
44.Goldhill DH, Turner PE. The evolution of life history trade-offs in viruses. Curr Opin Virol. 2014;8:79–84.PubMed 
Article 

Google Scholar 
45.Keen EC. Tradeoffs in bacteriophage life histories. Bacteriophage. 2014;4:e28365.PubMed 
PubMed Central 
Article 

Google Scholar 
46.Gómez P, Buckling A. Real-time microbial adaptive diversification in soil. Ecol Lett. 2013;16:650–5.PubMed 
Article 

Google Scholar 
47.Houte S, van, Buckling A, Westra ER. Evolutionary Ecology of Prokaryotic Immune Mechanisms. Microbiol Mol Biol Rev. 2016;80:745–63.PubMed 
PubMed Central 
Article 

Google Scholar 
48.Middelboe M, Hagström A, Blackburn N, Sinn B, Fischer U, Borch NH, et al. Effects of bacteriophages on the population dynamics of four strains of pelagic marine bacteria. Micro Ecol. 2001;42:395–406.CAS 
Article 

Google Scholar 
49.Gómez P, Buckling A. Coevolution with phages does not influence the evolution of bacterial mutation rates in soil. ISME J. 2013;7:2242–4.PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
50.De Sordi L, Khanna V, Debarbieux L. The Gut Microbiota Facilitates Drifts in the Genetic Diversity and Infectivity of Bacterial Viruses. Cell Host Microbe. 2017;22:801–8.e3.CAS 
PubMed 
Article 

Google Scholar 
51.De Sordi L, Lourenço M, Debarbieux L. “I will survive”: A tale of bacteriophage-bacteria coevolution in the gut. Gut Microbes. 2019;10:92–9.CAS 
PubMed 
Article 

Google Scholar 
52.Landsberger M, Gandon S, Meaden S, Chabas H, Buckling A, Westra ER, et al. Anti-CRISPR phages cooperate to overcome CRISPR-Cas immunity. Cell. 2018;174:908–16.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
53.Westra ER, van Houte S, Oyesiku-Blakemore S, Makin B, Broniewski JM, Best A, et al. Parasite exposure drives selective evolution of constitutive versus inducible defense. Curr Biol. 2015;25:1043–9.CAS 
PubMed 
Article 

Google Scholar 
54.Dy RL, Richter C, Salmond GP, Fineran PC. Remarkable mechanisms in microbes to resist phage infections. Annu Rev Virol. 2014;1:307–31.PubMed 
Article 
CAS 

Google Scholar 
55.Rostøl JT, Marraffini L. (Ph)ighting phages: how bacteria resist their parasites. Cell Host Microbe. 2019;25:184–94.PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
56.Burmeister AR, Turner PE. Trading-off and trading-up in the world of bacteria–phage evolution. Curr Biol. 2020;30:R1120–R1124.CAS 
PubMed 
Article 

Google Scholar 
57.Plummer M. JAGS: a program for analysis of Bayesian graphical models using Gibbs sampling. Vienna, Austria: Proc. 3rd Int. Workshop Distrib. Stat. Comput; 2003. p. 1–10.58.Wickham H. ggplot2: elegant Graphics for Data Analysis. Verlag New York: Springer; 2016.59.Wickham H. tidyr: Tidy Messy Data. 2020.60.Plummer M. rjags: Bayesian Graphical Models using MCMC. 2019.61.Wickham H, François R, Henry L, Müller K. dplyr: A Grammar of Data Manipulation. 2020.62.Gandon S, Buckling A, Decaestecker E, Day T. Host-parasite coevolution and patterns of adaptation across time and space. J Evol Biol. 2008;21:1861–6.CAS 
PubMed 
Article 

Google Scholar  More