Ecology

1.Haley, N. J. & Hoover, E. A. Chronic wasting disease of cervids: current knowledge and future perspectives. Annu. Rev. Anim. Biosci. 3, 305–325 (2015).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
2.Hannaoui, S., Schatzl, H. M. & Gilch, S. Chronic wasting disease: emerging prions and their potential risk. PLoS Pathog. 13, e1006619 (2017).PubMed 
Article 
CAS 
PubMed Central 

Google Scholar 
3.United States Geological Survey. Expanding Distribution of Chronic Wasting Disease. https://www.usgs.gov/centers/nwhc/science/expanding-distribution-chronic-wasting-disease?qt-science_center_objects=0#qt-science_center_objects (2021).4.Benestad, S. L., Mitchell, G., Simmons, M., Ytrehus, B. & Vikøren, T. First case of chronic wasting disease in Europe in a Norwegian free-ranging reindeer. Vet. Res. 47, 1–7 (2016).Article 

Google Scholar 
5.Gough, K. C. & Maddison, B. C. Prion transmission: prion excretion and occurrence in the environment. Prion 4, 275–282 (2010).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
6.Donaldson, D. S., Sehgal, A., Rios, D., Williams, I. R. & Mabbott, N. A. Increased abundance of m cells in the gut epithelium dramatically enhances oral prion disease susceptibility. PLoS Pathog. 12, e1006075 (2016).PubMed 
Article 
CAS 
PubMed Central 

Google Scholar 
7.Press, C. M., Heggebø, R. & Espenes, A. Involvement of gut-associated lymphoid tissue of ruminants in the spread of transmissible spongiform encephalopathies. Adv. Drug Deliv. Rev. 56, 885–899 (2004).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
8.Corr, S. C., Gahan, C. C. G. M. & Hill, C. M-cells: origin, morphology and role in mucosal immunity and microbial pathogenesis. FEMS Immunol. Med. Microbiol. 52, 2–12 (2008).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
9.Mabbott, N. A., Donaldson, D. S., Ohno, H., Williams, I. R. & Mahajan, A. Microfold (M) cells: important immunosurveillance posts in the intestinal epithelium. Mucosal Immunol. 6, 666–677 (2013).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
10.Maignien, T., Lasmézas, C. I., Beringue, V., Dormont, D. & Deslys, J. P. Pathogenesis of the oral route of infection of mice with scrapie and bovine spongiform encephalopathy agents. J. Gen. Virol. 80(Pt 11), 3035–3042 (1999).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
11.Bennett, K. M. et al. Induction of colonic m cells during intestinal inflammation. Am. J. Pathol. 186, 1166–1179 (2016).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
12.Donaldson, D. S. & Mabbott, N. A. The influence of the commensal and pathogenic gut microbiota on prion disease pathogenesis. J. Gen. Virol. 97, 1725–1738 (2016).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
13.Terahara, K. et al. Comprehensive gene expression profiling of peyer’s patch m cells, villous m-like cells, and intestinal epithelial cells. J. Immunol. 180, 7840–7846 (2008).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
14.Sigurdson, C. J. et al. Bacterial colitis increases susceptibility to oral prion disease. J. Infect. Dis. 199, 243–252 (2009).PubMed 
Article 
PubMed Central 

Google Scholar 
15.Tahoun, A. et al. Salmonella transforms follicle-associated epithelial cells into m cells to promote intestinal invasion. Cell Host Microb. 12, 645–656 (2012).CAS 
Article 

Google Scholar 
16.Kamada, N., Chen, G. Y., Inohara, N. & Núñez, G. Control of pathogens and pathobionts by the gut microbiota. Nat. Immunol. 14, 685–690 (2013).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
17.Round, J. L. & Mazmanian, S. K. The gut microbiota shapes intestinal immune responses during health and disease. Nat. Rev. Immunol. 9, 313–323 (2009).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
18.Ogbonnaya, E. S. et al. Adult hippocampal neurogenesis is regulated by the microbiome. Biol. Psychiatr. 78, e7-9 (2015).Article 

Google Scholar 
19.Heijtz, R. D. et al. Normal gut microbiota modulates brain development and behavior. Proc. Natl. Acad. Sci. U. S. A. 108, 3047–3052 (2011).ADS 
CAS 
Article 

Google Scholar 
20.Erny, D. et al. Host microbiota constantly control maturation and function of microglia in the CNS. Nat. Neurosci. 18, 965–977 (2015).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
21.Fung, T. C., Olson, C. A. & Hsiao, E. Y. Interactions between the microbiota, immune and nervous systems in health and disease. Nat. Neurosci. 20, 145–155 (2017).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
22.Chu, Y. & Kordower, J. H. The prion hypothesis of Parkinson’s disease. Curr. Neurol. Neurosci. Rep. 15, 28 (2015).PubMed 
Article 
CAS 
PubMed Central 

Google Scholar 
23.Goedert, M. Alzheimer’s and Parkinson’s diseases: the prion concept in relation to assembled Aβ, tau, and α-synuclein. Science 349, 1255555 (2015).PubMed 
Article 
CAS 
PubMed Central 

Google Scholar 
24.Herva, M. E. & Spillantini, M. G. Parkinson’s disease as a member of prion-like disorders. Virus Res. 207, 38–46 (2015).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
25.Tan, J. M. M., Wong, E. S. P. & Lim, K.-L. Protein misfolding and aggregation in Parkinson’s disease. Antioxid. Redox Signal. 11, 2119–2134 (2009).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
26.Weickenmeier, J., Jucker, M., Goriely, A. & Kuhl, E. A physics-based model explains the prion-like features of neurodegeneration in Alzheimer’s disease, Parkinson’s disease, and amyotrophic lateral sclerosis. J. Mech. Phys. Solids 124, 264–281 (2019).ADS 
Article 

Google Scholar 
27.Olanow, C. W. & Brundin, P. Parkinson’s disease and alpha synuclein: is Parkinson’s disease a prion-like disorder?. Mov. Disord. 28, 31–40 (2013).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
28.D’Argenio, V. & Sarnataro, D. Microbiome influence in the pathogenesis of prion and Alzheimer’s diseases. Int. J. Mol. Sci. 20, 4704 (2019).PubMed Central 
Article 
CAS 

Google Scholar 
29.Kang, D.-W. et al. Microbiota transfer therapy alters gut ecosystem and improves gastrointestinal and autism symptoms: an open-label study. Microbiome 5, 10 (2017).PubMed 
Article 
PubMed Central 

Google Scholar 
30.Rowin, J., Xia, Y., Jung, B. & Sun, J. Gut inflammation and dysbiosis in human motor neuron disease. Physiol. Rep. 5, e13443 (2017).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
31.Sampson, T. R. et al. Gut microbiota regulate motor deficits and neuroinflammation in a model of Parkinson’s disease. Cell 167, 1469-1480.e12 (2016).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
32.Vogt, N. M. et al. Gut microbiome alterations in Alzheimer’s disease. Sci. Rep. 7, 13537 (2017).ADS 
PubMed 
Article 
CAS 

Google Scholar 
33.Lev, M., Raine, C. S. & Levenson, S. M. Enhanced survival of germfree mice after infection with irradiated scrapie brain. Experientia 27, 1358–1359 (1971).CAS 
PubMed 
Article 

Google Scholar 
34.Wade, W. F., Dees, C., German, T. L. & Marsh, R. F. Effect of bacterial flora and mouse genotype (euthymic or athymic) on scrapie pathogenesis. J. Leukoc. Biol. 40, 525–532 (1986).CAS 
PubMed 
Article 

Google Scholar 
35.Bradford, B. M., Tetlow, L. & Mabbott, N. A. Prion disease pathogenesis in the absence of the commensal microbiota. J. Gen. Virol. 98, 1943–1952 (2017).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
36.Guan, Y. et al. Comparison of the gut microbiota composition between wild and captive sika deer (Cervus nippon hortulorum) from feces by high-throughput sequencing. AMB Express 7, 212 (2017).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
37.USDA APHIS|Cervids: Chronic Wasting Disease. https://www.aphis.usda.gov/aphis/ourfocus/animalhealth/animal-disease-information/cervid/cervids-cwd/cervid-cwd (2020).38.Keane, D. P. et al. Chronic wasting disease in a Wisconsin white-tailed deer farm. J. Vet. Diagn. Invest. 20, 698–703 (2008).PubMed 
Article 
PubMed Central 

Google Scholar 
39.Ethanol Precipitation Protocol—MRC Holland Technical Support. https://support.mrcholland.com/kb/articles/ethanol-precipitation-protocol.40.Apprill, A. & Parada, A. E. 16S Illumina amplicon protocol: Earth microbiome project. http://press.igsb.anl.gov/earthmicrobiome/protocols-and-standards/16s/.41.Boylen, E. et al. Reproducible, interactive, scalable and extensible microbiome data science using QIIME 2. Nat. Biotechnol. 37, 852–857 (2019).Article 
CAS 

Google Scholar 
42.Yilmaz, P. et al. The SILVA and ‘all-species living tree project (LTP)’ taxonomic frameworks. Nucl. Acids Res. https://doi.org/10.1093/nar/gkt1209 (2014).Article 
PubMed 
PubMed Central 

Google Scholar 
43.Quast, C. et al. The SILVA ribosomal RNA gene database project: Improved data processing and web-based tools. Nucl. Acids Res. https://doi.org/10.1093/nar/gks1219 (2013).Article 
PubMed 
PubMed Central 

Google Scholar 
44.Callahan, B. J. et al. DADA2: High resolution sample inference from Illumina amplicon data. Nat. Methods 13, 581–583 (2016).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
45.Turnbaugh, P. J. et al. An obesity-associated gut microbiome with increased capacity for energy harvest. Nature 444, 1027–1031 (2006).ADS 
PubMed 
PubMed Central 
Article 

Google Scholar 
46.Min, B. R., Gurung, N., Shange, R. & Solaiman, S. Potential role of rumen microbiota in altering average daily gain and feed efficiency in meat goats fed simple and mixed pastures using bacterial tag-encoded FLX amplicon pyrosequencing. J. Anim. Sci. 97, 3523–3534 (2019).PubMed 
PubMed Central 
Article 

Google Scholar 
47.Clayton, J. B. et al. Associations between nutrition, gut microbiome, and health in a novel nonhuman primate model. Sci. Rep. 8, 11159 (2018).ADS 
PubMed 
Article 
CAS 
PubMed Central 

Google Scholar 
48.Moeller, A. H. et al. Social behavior shapes the chimpanzee pan-microbiome. Sci. Adv. 2, e1500997 (2016).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
49.Goodrich, J. K. et al. Human genetics shape the gut microbiome. Cell 159, 789–799 (2014).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
50.Loo, W. T., García-Loor, J., Dudaniec, R. Y., Kleindorfer, S. & Cavanaugh, C. M. Host phylogeny, diet, and habitat differentiate the gut microbiomes of Darwin’s finches on Santa Cruz Island. Sci. Rep. 9, 18781 (2019).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
51.David, L. A. et al. Diet rapidly and reproducibly alters the human gut microbiome. Nature 505, 559–563 (2014).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
52.Delgado, M. L. et al. Intestinal microbial community dynamics of white-tailed deer (Odocoileus virginianus) in an agroecosystem. Microb. Ecol. 74, 496–506 (2017).PubMed 
Article 
PubMed Central 

Google Scholar 
53.Rothschild, D. et al. Environment dominates over host genetics in shaping human gut microbiota. Nature 555, 210–215 (2018).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
54.Rogers, L. L., Mooty, J. J. & Dawson, D. Foods of White-Tailed Deer in the Upper Great Lakes Region: A Review (North Central Forest Experiment Station, Forest Service, U.S. Dept. of Agriculture, 1981).Book 

Google Scholar 
55.Gibson, K. M. et al. Gut microbiome differences between wild and captive black rhinoceros—implications for rhino health. Sci. Rep. 9, 7570 (2019).ADS 
PubMed 
Article 
CAS 
PubMed Central 

Google Scholar 
56.Guo, W. et al. Comparative study of gut microbiota in wild and captive giant pandas (Ailuropoda melanoleuca). Genes (Basel) 10, 827 (2019).CAS 
Article 

Google Scholar 
57.Hale, V. L. et al. Gut microbiota in wild and captive Guizhou snub-nosed monkeys, Rhinopithecus brelichi. Am. J. Primatol. 81, e22989 (2019).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
58.Prabhu, V. R., Wasimuddin, W., Kamalakkannan, R., Arjun, M. S. & Nagarajan, M. Consequences of domestication on gut microbiome: A comparative study between wild gaur and domestic mithun. Front. Microbiol. 11, 133 (2020).PubMed 
Article 
PubMed Central 

Google Scholar 
59.Amato, K. R. et al. The gut microbiota appears to compensate for seasonal diet variation in the wild black howler monkey (Alouatta pigra). Microb. Ecol. 69, 434–443 (2015).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
60.Clayton, J. B. et al. The gut microbiome of nonhuman primates: Lessons in ecology and evolution. Am. J. Primatol. 80, e22867 (2018).PubMed 
Article 
PubMed Central 

Google Scholar 
61.Khafipour, E. et al. Effects of grain feeding on microbiota in the digestive tract of cattle. Anim. Front. 6, 13–19 (2016).Article 

Google Scholar 
62.Liu, C. et al. Dynamic alterations in yak rumen bacteria community and metabolome characteristics in response to feed type. Front. Microbiol. 10, 1116 (2019).PubMed 
Article 
PubMed Central 

Google Scholar 
63.Wu, G. D. et al. Linking long-term dietary patterns with gut microbial enterotypes. Science 334, 105–108 (2011).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
64.Beier, P. Sex differences in quality of white-tailed deer diets. J. Mammal. 68, 323–329 (1987).Article 

Google Scholar 
65.Walker, A. W. et al. Dominant and diet-responsive groups of bacteria within the human colonic microbiota. ISME J. 5, 220–230 (2011).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
66.Yang, X. et al. Seasonal breeding leads to changes for gut microbiota diversity in the wild ground squirrel (Spermophilus dauricus). https://www.researchsquare.com/article/rs-96089/v1 (2020). https://doi.org/10.21203/rs.3.rs-96089/v1.67.Antwis, R. E., Edwards, K. L., Unwin, B., Walker, S. L. & Shultz, S. Rare gut microbiota associated with breeding success, hormone metabolites and ovarian cycle phase in the critically endangered eastern black rhino. Microbiome 7, 27 (2019).PubMed 
Article 
PubMed Central 

Google Scholar 
68.Gordon, I. R. Controlled Reproduction in Horses, Deer, and Camelids (Cab International, 1997).
Google Scholar 
69.Samuel, M. D. & Storm, D. J. Chronic wasting disease in white-tailed deer: Infection, mortality, and implications for heterogeneous transmission. Ecology 97, 3195–3205 (2016).PubMed 
Article 
PubMed Central 

Google Scholar 
70.Miller, M. W., Hobbs, N. T. & Tavener, S. J. Dynamics of prion disease transmission in mule deer. Ecol. Appl. 16, 2208–2214 (2006).PubMed 
Article 
PubMed Central 

Google Scholar 
71.Gandy, K. A. O., Zhang, J., Nagarkatti, P. & Nagarkatti, M. The role of gut microbiota in shaping the relapse-remitting and chronic-progressive forms of multiple sclerosis in mouse models. Sci. Rep. 9, 1–17 (2019).CAS 
Article 

Google Scholar 
72.Ebringer, A., Rashid, T., Wilson, C., Boden, R. & Thompson, E. A possible link between multiple sclerosis and Creutzfeldt-Jakob disease based on clinical, genetic, pathological and immunological evidence involving Acinetobacter bacteria. Med. Hypotheses 64, 487–494 (2005).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
73.Ricci, S. et al. Impact of supplemental winter feeding on ruminal microbiota of roe deer Capreolus capreolus. wbio 2019, 1–11 (2019).Article 

Google Scholar 
74.Sun, C.-H., Liu, H.-Y., Liu, B., Yuan, B.-D. & Lu, C.-H. Analysis of the gut microbiome of wild and captive Père David’s deer. Front. Microbiol. 10, 2331 (2019).PubMed 
Article 

Google Scholar 
75.Barichella, M. et al. Unraveling gut microbiota in Parkinson’s disease and atypical parkinsonism. Mov. Disord. 34, 396–405 (2019).PubMed 
Article 

Google Scholar 
76.Pietrucci, D. et al. Dysbiosis of gut microbiota in a selected population of Parkinson’s patients. Parkinsonism Relat. Disord. 65, 124–130 (2019).PubMed 
Article 

Google Scholar 
77.Radisavljevic, N., Cirstea, M. & Brett Finlay, B. Bottoms up: The role of gut microbiota in brain health. Environ. Microbiol https://doi.org/10.1111/1462-2920.14506 (2018).Article 
PubMed 

Google Scholar 
78.Geirnaert, A. et al. Butyrate-producing bacteria supplemented in vitro to Crohn’s disease patient microbiota increased butyrate production and enhanced intestinal epithelial barrier integrity. Sci. Rep. 7, 11450 (2017).ADS 
PubMed 
Article 
CAS 

Google Scholar 
79.Vacca, M. et al. The controversial role of human gut Lachnospiraceae. Microorganisms 8, 573 (2020).CAS 
PubMed Central 
Article 
PubMed 

Google Scholar 
80.Zeng, H., Ishaq, S. L., Zhao, F.-Q. & Wright, A.-D.G. Colonic inflammation accompanies an increase of β-catenin signaling and Lachnospiraceae/Streptococcaceae bacteria in the hind gut of high-fat diet-fed mice. J. Nutr. Biochem. 35, 30–36 (2016).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
81.Zhang, T., Li, Q., Cheng, L., Buch, H. & Zhang, F. Akkermansia muciniphila is a promising probiotic. Microb. Biotechnol. 12, 1109–1125 (2019).PubMed 
PubMed Central 
Article 

Google Scholar 
82.Zhou, K. Strategies to promote abundance of Akkermansia muciniphila, an emerging probiotics in the gut, evidence from dietary intervention studies. J. Funct. Foods 33, 194–201 (2017).PubMed 
PubMed Central 
Article 

Google Scholar 
83.Ou, Z. et al. Protective effects of Akkermansia muciniphila on cognitive deficits and amyloid pathology in a mouse model of Alzheimer’s disease. Nutr. Diabetes 10, 1–10 (2020).Article 
CAS 

Google Scholar 
84.Cirstea, M., Radisavljevic, N. & Finlay, B. B. Good bug, bad bug: Breaking through microbial stereotypes. Cell Host Microb. 23, 10–13 (2018).CAS 
Article 

Google Scholar 
85.Berer, K. et al. Gut microbiota from multiple sclerosis patients enables spontaneous autoimmune encephalomyelitis in mice. PNAS 114, 10719–10724 (2017).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
86.Cekanaviciute, E. et al. Gut bacteria from multiple sclerosis patients modulate human T cells and exacerbate symptoms in mouse models. PNAS 114, 10713–10718 (2017).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
87.Heintz-Buschart, A. et al. The nasal and gut microbiome in Parkinson’s disease and idiopathic rapid eye movement sleep behavior disorder. Mov. Disord. 33, 88–98 (2018).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
88.Hill-Burns, E. M. et al. Parkinson’s disease and PD medications have distinct signatures of the gut microbiome. Mov. Disord. 32, 739–749 (2017).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
89.Belzer, C. & de Vos, W. M. Microbes inside—from diversity to function: The case of Akkermansia. ISME J. 6, 1449–1458 (2012).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
90.Donaldson, D. S., Pollock, J., Vohra, P., Stevens, M. P. & Mabbott, N. A. Microbial stimulation reverses the age-related decline in M cells in aged mice. iScience 23, 101147 (2020).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
91.Ganesh, B. P., Klopfleisch, R., Loh, G. & Blaut, M. Commensal Akkermansia muciniphila exacerbates gut inflammation in Salmonella Typhimurium-infected gnotobiotic mice. PLoS ONE 8, e74963 (2013).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
92.Donaldson, D. S. et al. M cell-depletion blocks oral prion disease pathogenesis. Mucosal Immunol. 5, 216–225 (2012).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
93.Derrien, M., Belzer, C. & de Vos, W. M. Akkermansia muciniphila and its role in regulating host functions. Microb. Pathog. 106, 171–181 (2017).PubMed 
Article 
PubMed Central 

Google Scholar 
94.Nagalingam, N. A., Kao, J. Y. & Young, V. B. Microbial ecology of the murine gut associated with the development of DSS-colitis. Inflamm. Bowel Dis. 17, 917–926 (2011).PubMed 
Article 
PubMed Central 

Google Scholar 
95.Haley, N. J. et al. Estimating relative CWD susceptibility and disease progression in farmed white-tailed deer with rare PRNP alleles. PLoS ONE 14, e0224342 (2019).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
96.Huttenhower, C. et al. Structure, function and diversity of the healthy human microbiome. Nature 486, 207–214 (2012).ADS 
CAS 
Article 

Google Scholar 
97.Wagner Mackenzie, B., Waite, D. W. & Taylor, M. W. Evaluating variation in human gut microbiota profiles due to DNA extraction method and inter-subject differences. Front. Microbiol. 6, 130 (2015).PubMed 
Article 
PubMed Central 

Google Scholar 
98.Kawada, Y., Naito, Y., Andoh, A., Ozeki, M. & Inoue, R. Effect of storage and DNA extraction method on 16S rRNA-profiled fecal microbiota in Japanese adults. J. Clin. Biochem. Nutr. 64, 106–111 (2019).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
99.Sinha, R. et al. Collecting fecal samples for microbiome analyses in epidemiology studies. Cancer Epidemiol. Biomark. Prev. 25, 407–416 (2016).Article 

Google Scholar  More

1.deMenocal, P. B. & Stringer, C. Climate and the peopling of the world. Nature 538, 49–50 (2016).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
2.Pisor, A. C. & Jones, J. H. Human adaptation to climate change: an introduction to the special issue. Am. J. Hum. Biol. n/a, e23530 (2020).3.Rick, T. C. & Sandweiss, D. H. Archaeology, climate, and global change in the Age of Humans. PNAS 117, 8250–8253 (2020).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
4.Shennan, S. & Sear, R. Archaeology, demography and life history theory together can help us explain past and present population patterns. Philos. Trans. R. Soc. B: Biol. Sci. 376, 20190711 (2021).CAS 
Article 

Google Scholar 
5.Boivin, N. & Crowther, A. Mobilizing the past to shape a better Anthropocene. Nat. Ecol. Evolution 5, 273–284 (2021).Article 

Google Scholar 
6.Bocquet‐Appel, J. Recent Advances in Paleodemography (Springer, Dordrecht, 2008).7.Chamberlain, A. T. Demography in Archaeology (Cambridge University Press, 2006).8.Drennan, R. D., Berrey, C. A. & Peterson, C. E. Regional Settlement Demography in Archaeology (Eliot Werner Publications, 2015).9.Kintigh, K. W. et al. Grand challenges for archaeology. PNAS 111, 879–880 (2014).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
10.Bocquet‐Appel, J. Paleoanthropological traces of a neolithic demographic transition. Curr. Anthropol. 43, 637–650 (2002).Article 

Google Scholar 
11.Crema, E. R. & Kobayashi, K. A multi-proxy inference of Jōmon population dynamics using bayesian phase models, residential data, and summed probability distribution of 14C dates. J. Archaeol. Sci. 117, 105136 (2020).Article 

Google Scholar 
12.Schmidt, I. et al. Approaching prehistoric demography: proxies, scales and scope of the Cologne Protocol in European contexts. Philos. Trans. R. Soc. B: Biol. Sci. 376, 20190714 (2021).Article 

Google Scholar 
13.Schiffels, S. & Durbin, R. Inferring human population size and separation history from multiple genome sequences. Nat. Genet. 46, 919–925 (2014).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
14.White, A. J. et al. An evaluation of fecal stanols as indicators of population change at Cahokia, Illinois. J. Archaeol. Sci. 93, 129–134 (2018).CAS 
Article 

Google Scholar 
15.Shennan, S. et al. Regional population collapse followed initial agriculture booms in mid-Holocene Europe. Nat. Commun. 4, 1–8 (2013).ADS 
Article 
CAS 

Google Scholar 
16.Timpson, A. et al. Reconstructing regional population fluctuations in the European Neolithic using radiocarbon dates: a new case-study using an improved method. J. Archaeol. Sci. 52, 549–557 (2014).Article 

Google Scholar 
17.Crema, E. R., Habu, J., Kobayashi, K. & Madella, M. Summed probability distribution of 14 C dates suggests regional divergences in the population dynamics of the jomon period in Eastern Japan. PLoS ONE 11, e0154809 (2016).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
18.Crema, E. R., Bevan, A. & Shennan, S. Spatio-temporal approaches to archaeological radiocarbon dates. J. Archaeol. Sci. 87, 1–9 (2017).CAS 
Article 

Google Scholar 
19.Chaput, M. A. & Gajewski, K. Radiocarbon dates as estimates of ancient human population size. Anthropocene 15, 3–12 (2016).Article 

Google Scholar 
20.Carleton, W. C. Evaluating Bayesian Radiocarbon‐dated Event Count (REC) models for the study of long‐term human and environmental processes. Journal of Quaternary Science 36, 110–123 (2021).21.Brown, W. A. The past and future of growth rate estimation in demographic temporal frequency analysis: Biodemographic interpretability and the ascendance of dynamic growth models. J. Archaeol. Sci. 80, 96–108 (2017).Article 

Google Scholar 
22.Carleton, W. C. & Groucutt, H. S. Sum things are not what they seem: problems with point-wise interpretations and quantitative analyses of proxies based on aggregated radiocarbon dates. Holocene 0959683620981700. https://doi.org/10.1177/0959683620981700 (2020).23.Crema, E. R. & Bevan, A. Inference from large sets of radiocarbon dates: software and methods. Radiocarbon 63, 23–39 (2021).Article 

Google Scholar 
24.Williams, A. N. The use of summed radiocarbon probability distributions in archaeology: a review of methods. J. Archaeol. Sci. 39, 578–589 (2012).Article 

Google Scholar 
25.Ward, I. & Larcombe, P. Sedimentary unknowns constrain the current use of frequency analysis of radiocarbon data sets in forming regional models of demographic change. Geoarchaeology 36, 546–570 (2021).Article 

Google Scholar 
26.de Souza, J. G. & Riris, P. Delayed demographic transition following the adoption of cultivated plants in the eastern La Plata Basin and Atlantic coast, South America. J. Archaeol. Sci. 125, 105293 (2021).Article 

Google Scholar 
27.Fernández-López de Pablo, J. et al. Palaeodemographic modelling supports a population bottleneck during the Pleistocene-Holocene transition in Iberia. Nat. Commun. 10, 1872 (2019).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
28.Goldberg, A., Mychajliw, A. M. & Hadly, E. A. Post-invasion demography of prehistoric humans in South America. Nature 532, 232–235 (2016).ADS 
CAS 
PubMed 
Article 

Google Scholar 
29.Lima, M. et al. Ecology of the collapse of Rapa Nui society. Proc. R. Soc. B: Biol. Sci. 287, 20200662 (2020).CAS 
Article 

Google Scholar 
30.Prates, L., Politis, G. G. & Perez, S. I. Rapid radiation of humans in South America after the last glacial maximum: a radiocarbon-based study. PLoS ONE 15, e0236023 (2020).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
31.Riris, P. Dates as data revisited: a statistical examination of the Peruvian preceramic radiocarbon record. J. Archaeol. Sci. 97, 67–76 (2018).Article 

Google Scholar 
32.Riris, P. & Arroyo-Kalin, M. Widespread population decline in South America correlates with mid-Holocene climate change. Sci. Rep. 9, 6850 (2019).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
33.Crema, E. R. & Shoda, S. A Bayesian approach for fitting and comparing demographic growth models of radiocarbon dates: A case study on the Jomon-Yayoi transition in Kyushu (Japan). PLOS ONE 16, e0251695 (2021).34.Timpson, A., Barberena, R., Thomas, M. G., Méndez, C. & Manning, K. Directly modelling population dynamics in the South American Arid Diagonal using 14 C dates. Philos. Trans. R. Soc. B: Biol. Sci. 376, 20190723 (2021).CAS 
Article 

Google Scholar 
35.Bernabeu Aubán, J., García Puchol, O., Barton, M., McClure, S. & Pardo Gordó, S. Radiocarbon dates, climatic events, and social dynamics during the Early Neolithic in Mediterranean Iberia. Quat. Int. 403, 201–210 (2016).Article 

Google Scholar 
36.Bevan, A. et al. Holocene fluctuations in human population demonstrate repeated links to food production and climate. Proc. Natl Acad. Sci. USA 114, E10524–E10531 (2017).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
37.Bird, D. et al. A first empirical analysis of population stability in North America using radiocarbon records. Holocene 30, 1345–1359 (2020).ADS 
Article 

Google Scholar 
38.Capuzzo, G., Zanon, M., Corso, M. D., Kirleis, W. & Barceló, J. A. Highly diverse Bronze Age population dynamics in Central-Southern Europe and their response to regional climatic patterns. PLoS ONE 13, e0200709 (2018).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
39.de Souza, J. G. et al. Climate change and cultural resilience in late pre-Columbian Amazonia. Nat. Ecol. Evol. 3, 1007–1017 (2019).PubMed 
Article 

Google Scholar 
40.Jørgensen, E. K. The palaeodemographic and environmental dynamics of prehistoric Arctic Norway: an overview of human-climate covariation. Quat. Int. 549, 36–51 (2020).Article 

Google Scholar 
41.Kelly, R. L., Surovell, T. A., Shuman, B. N. & Smith, G. M. A continuous climatic impact on Holocene human population in the Rocky Mountains. PNAS 110, 443–447 (2013).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
42.Roberts, N. et al. Human responses and non-responses to climatic variations during the last Glacial-Interglacial transition in the eastern Mediterranean. Quat. Sci. Rev. 184, 47–67 (2018).ADS 
Article 

Google Scholar 
43.Wang, C., Lu, H., Zhang, J., Gu, Z. & He, K. Prehistoric demographic fluctuations in China inferred from radiocarbon data and their linkage with climate change over the past 50,000 years. Quat. Sci. Rev. 98, 45–59 (2014).ADS 
Article 

Google Scholar 
44.Warden, L. et al. Climate induced human demographic and cultural change in northern Europe during the mid-Holocene. Sci. Rep. 7, 15251 (2017).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
45.Stewart, M., Carleton, W. C. & Groucutt, H. S. Climate change, not human population growth, correlates with Late Quaternary megafauna declines in North America. Nat. Commun. 12, 965 (2021).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
46.Weninger, B., Clare, L., Jöris, O., Jung, R. & Edinborough, K. Quantum theory of radiocarbon calibration. World Archaeol. 47, 543–566 (2015).Article 

Google Scholar 
47.Weninger, B. & Edinborough, K. Bayesian 14C-rationality, Heisenberg uncertainty, and Fourier Transform: the beauty of radiocarbon calibration. Doc. Praehist. 47, 536–559 (2020).Article 

Google Scholar 
48.Tavaré, S., Balding, D. J., Griffiths, R. C. & Donnelly, P. Inferring coalescence times from DNA sequence data. Genetics 145, 505–518 (1997).PubMed 
PubMed Central 
Article 

Google Scholar 
49.Beaumont, M. A., Zhang, W. & Balding, D. J. Approximate Bayesian computation in population. Genet. Genet. 162, 2025–2035 (2002).
Google Scholar 
50.Carrignon, S., Brughmans, T. & Romanowska, I. Tableware trade in the Roman East: exploring cultural and economic transmission with agent-based modelling and approximate Bayesian computation. PLoS ONE 15, e0240414 (2020).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
51.Crema, E. R., Edinborough, K., Kerig, T. & Shennan, S. J. An approximate Bayesian computation approach for inferring patterns of cultural evolutionary change. J. Archaeol. Sci. 50, 160–170 (2014).Article 

Google Scholar 
52.Crema, E. R., Kandler, A. & Shennan, S. Revealing patterns of cultural transmission from frequency data: equilibrium and non-equilibrium assumptions. Sci. Rep. 6, 39122 (2016).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
53.Rubio-Campillo, X. Model selection in historical research using approximate Bayesian computation. PLoS ONE 11, e0146491 (2016).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
54.Tsutaya, T., Shimomi, A., Fujisawa, S., Katayama, K. & Yoneda, M. Isotopic evidence of breastfeeding and weaning practices in a hunter–gatherer population during the Late/Final Jomon period in eastern Japan. J. Archaeol. Sci. 76, 70–78 (2016).CAS 
Article 

Google Scholar 
55.Porčić, M. & Nikolić, M. The Approximate Bayesian Computation approach to reconstructing population dynamics and size from settlement data: demography of the Mesolithic-Neolithic transition at Lepenski Vir. Archaeol. Anthropol. Sci. 1–18. https://doi.org/10.1007/s12520-014-0223-2 (2015).56.Porčić, M., Blagojević, T., Pendić, J. & Stefanović, S. The Neolithic Demographic Transition in the Central Balkans: population dynamics reconstruction based on new radiocarbon evidence. Philos. Trans. R. Soc. B: Biol. Sci. 376, 20190712 (2021).Article 
CAS 

Google Scholar 
57.DiNapoli, R. J., Rieth, T. M., Lipo, C. P. & Hunt, T. L. A model-based approach to the tempo of “collapse”: The case of Rapa Nui (Easter Island). J. Archaeol. Sci. 116, 105094 (2020).Article 

Google Scholar 
58.Hunt, T. L. & Lipo, C. The Archaeology of Rapa Nui (Easter Island). in The Oxford Handbook of Prehistoric Oceania (eds. Cochrane, E. E. & Hunt, T. L.) 416–449 (Oxford University Press, 2018).59.Kirch, P. V. The Evolution of Polynesian Chiefdoms (Cambridge University Press, 1984).60.Ponting, C. A Green History of the World: The Environment and the Collapse of Great Civilizations. (St. Martin’s Press, 1991).61.Boersema, J. J. The Survival of Easter Island: Dwindling Resources and Cultural Resilience (Cambridge University Press, 2015).62.Boersema, J. J. An earthly paradise? Easter Island (Rapa Nui) as seen by the eighteenth-century European explorers. in Cultural and Environmental Change on Rapa Nui (eds. Haoa Cardinali, S. et al.) 157–178 (Routledge, 2018).63.Boersema, J. J. & Huele, R. Pondering the population numbers of Easter Island’s Past. in Easter Island and the Pacific: Cultural and Environmental Dynamics. In Proc 9th International Conference on Easter Island and the Pacific, Held in the Ethnological Museum, Berlin, Germany (eds. Vogt, B. et al.) 83–92 (Rapa Nui Press, 2019).64.Puleston, C. O. et al. Rain, sun, soil, and sweat: a consideration of population limits on Rapa Nui (Easter Island) before European Contact. Front. Ecol. Evol. 5, 1–14 (2017).Article 

Google Scholar 
65.Lipo, C. P., DiNapoli, R. J. & Hunt, T. L. Commentary: rain, sun, soil, and sweat: a consideration of population limits on Rapa Nui (Easter Island) before European Contact. Front. Ecol. Evol. 25, 1–3 (2018).66.Hunt, T. L. Rethinking Easter Island’s ecological catastrophe. J. Archaeol. Sci. 34, 485–502 (2007).Article 

Google Scholar 
67.Rull, V. The deforestation of Easter Island. Biol. Rev. 95, 124–141 (2020).Article 

Google Scholar 
68.Brandt, G. & Merico, A. The slow demise of Easter Island: insights from a modeling investigation. Front. Ecol. Evol. 13, 1–12 (2015).69.Diamond, J. Collapse: How Societies Choose to Fail or Succeed (Viking, 2005).70.Bahn, P. & Flenley, J. Easter Island, Earth Island: the Enigmas of Rapa Nui (Rowman & Littlefield, 2017).71.Rull, V. Natural and anthropogenic drivers of cultural change on Easter Island: review and new insights. Quat. Sci. Rev. 150, 31–41 (2016).ADS 
Article 

Google Scholar 
72.Cañellas-Boltà, N. et al. Vegetation changes and human settlement of Easter Island during the last millennia: a multiproxy study of the Lake Raraku sediments. Quat. Sci. Rev. 72, 36–48 (2013).ADS 
Article 

Google Scholar 
73.Rull, V. Drought, freshwater availability and cultural resilience on Easter Island (SE Pacific) during the Little Ice Age. Holocene. https://doi.org/10.1177/0959683619895587 (2020).74.Yan, H. et al. A record of the Southern Oscillation Index for the past 2,000 years from precipitation proxies. Nat. Geosci. 4, 611–614 (2011).ADS 
CAS 
Article 

Google Scholar 
75.Mulrooney, M. A. An island-wide assessment of the chronology of settlement and land use on Rapa Nui (Easter Island) based on radiocarbon data. J. Archaeol. Sci. 40, 4377–4399 (2013).Article 

Google Scholar 
76.Stevenson, C. M. et al. Variation in Rapa Nui (Easter Island) land use indicates production and population peaks prior to European contact. Proc. Natl Acad. Sci. USA 112, 1025–1030 (2015).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
77.DiNapoli, R. J., Lipo, C. P. & Hunt, T. L. Revisiting warfare, monument destruction, and the ‘Huri Moai’ phase in Rapa Nui (Easter Island) culture history. Journal of Pacific Archaeology 12, 1–24 (2021).78.Hogg, A. G. et al. SHCal13 Southern Hemisphere calibration, 0–50,000 years cal BP. Radiocarbon 55, 1889–1903 (2013).CAS 
Article 

Google Scholar 
79.Hogg, A. G. et al. SHCal20 Southern Hemisphere calibration, 0–55,000 Years cal BP. Radiocarbon 62, 759–778 (2020).CAS 
Article 

Google Scholar 
80.Bork, H.-R., Mieth, A. & Tschochner, B. Nothing but stones? A review of the extent and technical efforts of prehistoric stone mulching on Rapa Nui. Rapa Nui J. 18, 10–14 (2004).
Google Scholar 
81.Ladefoged, T. N. et al. Soil nutrient analysis of Rapa Nui gardening. Archaeol. Ocean. 45, 80–85 (2010).Article 

Google Scholar 
82.Ladefoged, T. N., Flaws, A. & Stevenson, C. M. The distribution of rock gardens on Rapa Nui (Easter Island) as determined from satellite imagery. J. Archaeol. Sci. 40, 1203–1212 (2013).Article 

Google Scholar 
83.Mieth, A. & Bork, H. R. History, origin and extent of soil erosion on Easter Island (Rapa Nui). Catena 63, 244–260 (2005).Article 

Google Scholar 
84.Stevenson, C. M., Jackson, T. L., Mieth, A., Bork, H.-R. & Ladefoged, T. N. Prehistoric and early historic agriculture at Maunga Orito, Easter Island (Rapa Nui), Chile. Antiquity 80, 919–936 (2006).Article 

Google Scholar 
85.Wozniak, J. A. Subsistence strategies on Rapa Nui (Easter Island): prehistoric gardening practices on Rapa Nui and how they relate to current farming practices. in Cultural and Environmental Change on Rapa Nui (eds. Haoa-Cardinali, S. et al.) 87–112 (Routledge, 2018).86.Tromp, M. & Dudgeon, J. V. Differentiating dietary and non-dietary microfossils extracted from human dental calculus: the importance of sweet potato to ancient diet on Rapa Nui. J. Archaeol. Sci. 54, 54–63 (2015).Article 

Google Scholar 
87.Brosnan, T., Becker, M. W. & Lipo, C. P. Coastal groundwater discharge and the ancient inhabitants of Rapa Nui (Easter Island), Chile. Hydrogeol. J. 27, 519–534 (2019).ADS 
CAS 
Article 

Google Scholar 
88.DiNapoli, R. J. et al. Rapa Nui (Easter Island) monument (ahu) locations explained by freshwater sources. PLoS ONE 14, e0210409 (2019).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
89.Hixon, S., DiNapoli, R. J., Lipo, C. P. & Hunt, T. L. The ethnohistory of freshwater use on Rapa Nui (Easter Island, Chile). J. Polynesian Soc. 128, 163–189 (2019).Article 

Google Scholar 
90.Brown, A. A. & Crema, E. R. Māori population growth in pre-contact New Zealand: regional population dynamics inferred from summed probability distributions of radiocarbon dates. J. Isl. Coast. Archaeol. 0, 1–19 (2019).
Google Scholar 
91.McFadden, C., Walter, R., Buckley, H. & Oxenham, M. F. Temporal trends in the Colonisation of the Pacific: Palaeodemographic Insights. J. World Prehist. https://doi.org/10.1007/s10963-021-09152-w (2021).Article 

Google Scholar 
92.Kirch, P. V. & Rallu, J.-L. The Growth and Collapse of Pacific Island Societies: Archaeological and Demographic Perspectives. (University of Hawai’i Press, 2007).93.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 
94.Sherwood, S. C. et al. New excavations in Easter Island’s statue quarry: Soil fertility, site formation and chronology. J. Archaeological Sci. 111, 104994 (2019).Article 

Google Scholar 
95.Simpson, D. F. Jr. & Dussubieux, L. A collapsed narrative? Geochemistry and spatial distribution of basalt quarries and fine–grained artifacts reveal communal use of stone on Rapa Nui (Easter Island). J. Archaeol. Sci.: Rep. 18, 370–385 (2018).
Google Scholar 
96.Bevan, A. & Crema, E. R. Modifiable reporting unit problems and time series of long-term human activity. Philos. Trans. R. Soc. B: Biol. Sci. 376, 20190726 (2021).CAS 
Article 

Google Scholar 
97.Davies, B., Holdaway, S. J. & Fanning, P. C. Modelling the palimpsest: an exploratory agent-based model of surface archaeological deposit formation in a fluvial arid Australian landscape. Holocene 26, 450–463 (2016).ADS 
Article 

Google Scholar 
98.Commendador, A. S., Dudgeon, J. V., Fuller, B. T. & Finney, B. P. Radiocarbon dating human skeletal material on Rapa Nui: evaluating the effect of uncertainty in marine-derived carbon. Radiocarbon 56, 277–294 (2014).CAS 
Article 

Google Scholar 
99.Stevenson, C. M., Williams, C., Carpenter, E., Hunt, C. S. & Novak, S. W. Architecturally modified caves on Rapa Nui: post-European contact ritual spaces? Rapa Nui J. 32, 1–36 (2019).Article 

Google Scholar 
100.Heaton, T. J. et al. Marine20—the marine radiocarbon age calibration curve (0–55,000 cal BP). Radiocarbon 62, 779–820 (2020).CAS 
Article 

Google Scholar 
101.Beck, J. W., Hewitt, L., Burr, G. S., Loret, J. & Hochstetter, F. T. Mata ki te rangi: eyes towards the heavens. in Easter Island: Scientific Exploration Into the World’s Environmental Problems in Microcosm (eds. Loret, J. & Tanacredi, J. T.) 93–112 (Kluwer Academic/Plenum Publishers, 2003).102.Burr, G. S. et al. Modern and Pleistocene reservoir ages inferred from South Pacific corals. Radiocarbon 51, 319–335 (2009).CAS 
Article 

Google Scholar 
103.DiNapoli, R. J. et al. Marine reservoir corrections for the Caribbean demonstrate high intra- and inter-island variability in local reservoir offsets. Quat. Geochronol. 61, 101126 (2021).Article 

Google Scholar 
104.Surovell, T. A., Byrd Finley, J., Smith, G. M., Brantingham, P. J. & Kelly, R. Correcting temporal frequency distributions for taphonomic bias. J. Archaeol. Sci. 36, 1715–1724 (2009).Article 

Google Scholar 
105.Sunnåker, M. et al. Approximate Bayesian computation. PLoS Comput. Biol. 9, e1002803 (2013).MathSciNet 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
106.Toni, T., Welch, D., Strelkowa, N., Ipsen, A. & Stumpf, M. P. H. Approximate Bayesian computation scheme for parameter inference and model selection in dynamical systems. J. R. Soc. Interface 6, 187–202 (2009).PubMed 
Article 

Google Scholar 
107.Rick, J. W. Dates as data: an examination of the peruvian preceramic radiocarbon record. Am. Antiquity 52, 55–73 (1987).Article 

Google Scholar 
108.R Core Team. R: a Language and Environment for Statistical Computing. (R Foundation for Statistical Computing, 2020). More

The modelConsider a unit-mass population of agents who repeatedly interact in pairs to play a symmetric stage game. The set of strategies available to each agent is finite and denoted by (S equiv {1, ldots , n}). A population state is a vector (x in X equiv {x in {mathbb{R}}^n_+: sum _{i in S} x_i = 1}), with (x_i) the fraction of the population playing strategy (i in S). Payoffs are described by a function (F: S times S rightarrow {mathbb{R}}), where F(i, j) is the payoff received by an agent playing strategy i when the opponent plays strategy j. As a shorthand, we refer to an undirected pair of individuals, one playing i and the other playing j, as an ij pair. The set of all possible undirected pairs is denoted by (mathscr {P}).The interaction structure is modeled as a function (p : X times mathscr {P} rightarrow left[ 0, 1/2 right] ) subject to (sum _{ij in mathscr {P}} p_{ij}(x)=1/2) (since the mass of pairs is half the mass of agents), with (p_{ij}(x)) indicating the mass of ij pairs formed in state x. Note that the mass of ij pairs can never exceed (min {x_i,x_j}), that is, (p_{ij}(x) le min {x_i,x_j}) for all x. We assume that p is continuous in X, and that (p_{ij}(x) > 0) if and only if (x_i > 0) and (x_j > 0 )—meaning that the probability of an ij pair being formed is strictly positive if and only if strategies i and j are played by someone. In the case of uniform random matching, (p_{ii} = x_i^2/2) and (p_{ij} = x_i x_j) for any i and (j ne i).The revision protocol is modeled as a function (phi : X times S times S rightarrow [-1,1]), where (phi _{ij}(x) in [-1,1]) is the probability that an ij pair will turn into an ii pair minus the probability that it will turn into a jj pair, conditional on the population state being x and an ij pair being formed. We assume that (phi ) is continuous in X. We note that by construction (phi _{ij}=-phi _{ji}) for all (i,j in S), and hence (phi _{ii}=0) for all (i in S). Our main assumption on the revision protocol is the following, which is met, among others, by pairwise proportional imitative and imitate-if-better rules22.
Assumption 1

For every (x in X), (phi _{ij}(x) > 0) if (F(i,j) > F(j,i)).
In what follows we consider a dynamical system in continuous time with state space X, characterized by the following equation of motion.

Definition 1

(Pairwise interact-and-imitate dynamics—PIID) For every (x in X) and every (i in S):$$begin{aligned} dot{x}_i = sum _{j in S} p_{ij}(x) phi _{ij}(x). end{aligned}$$
(1)

Main findingsGlobal asymptotic convergenceIn any purely imitative dynamics, if (x_i(t)=0), then (x_i(t^{prime})=0) for every (t^{prime} > t). This implies that we cannot hope for global asymptotic convergence in a strict sense. Thus, to assess convergence towards a certain state x in a meaningful way, we restrict our attention to those states where all strategies that have positive frequency in x have positive frequency as well. We denote by (X_x) the set of states whose support contains the support of x.

Definition 2

(Supremacy) Strategy (iin S) is supreme if (F(i,j) >F(j,i)) for every (j in S setminus {i}).
We note that under PIID, the concept of supremacy is closely related to that of asymmetry33,34, in that (F(i,j) > F(j,i)) implies that agents can only switch from strategy j to strategy i.

Proposition 1

If (i in S) is a supreme strategy, then state (x^* equiv left{ x in X : x_i = 1 right} ) is globally asymptotically stable for the dynamical system with state space (X_{x^*}) and PIID as equation of motion.
Relation to replicator dynamicsTo further characterize the dynamics induced by the pairwise interact-and-imitate protocol, we make two additional assumptions. First, matching is uniformly random, meaning that everyone in the population has the same probability of interacting with everyone else; formally, (p_{ii} = x_i^2/2) and (p_{ij} = x_i x_j) for all i and (j ne i). Second, the probability that an agent has to imitate the opponent is proportional to the difference in their payoffs if the opponent’s payoff exceeds her own, and is zero otherwise. As a consequence, (phi _{ij} = F(i,j) – F(j,i)) up to a proportionality factor. Let

(F left( i, x right) :=sum _j x_j F left( i, j right) ),

(F left( x, i right) :=sum _j x_j F left( j, i right) ), and

( F left( x, x right) :=sum _i sum _j x_i x_j F left( i, j right) ).

Under these assumptions, at any point in time, the motion of (x_i) is described by:$$begin{aligned} dot{x}_i&= sum _{j ne i} x_j x_i left[ F left( i, j right) – F left( j, i right) right] = x_i sum _{j} x_j left[ F left( i, j right) – F left( j, i right) right] nonumber \&= x_i left[ F left( i, x right) – F left( x, i right) right] , end{aligned}$$
(2)
which is a modified replicator equation. According to (2), for every strategy i chosen by one or more agents in the population, the rate of growth of the fraction of i-players, (dot{x}_i / x_i), equals the difference between the expected payoff from playing i in state x and the average payoff received by those who are matched against an agent playing i. In contrast, under standard replicator dynamics35, the fraction of agents playing i varies depending on the excess payoff of i with respect to the current average payoff in the whole population, i.e., (dot{x}_i = x_i left[ F left( i, x right) – F left( x, x right) right] ).A noteworthy feature of replicator dynamics is that they are always payoff monotone: for any (i,j in S), the proportions of agents playing i and j grow at rates that are ordered in the same way as the expected payoffs from the two strategies36. In the case of PIID, this result fails.

Proposition 2

Pairwise-Interact-and-Imitate dynamics need not satisfy payoff monotonicity.
To verify this, it is sufficient to consider any symmetric (2 times 2) game where (F left( i, j right) > F left( j, i right) ) but (F left( j, x right) > F left( i, x right) ) for some (x in X), meaning that i is the supreme strategy but j yields a higher expected payoff in state x. See Fig. 1 for an example where, in the case of uniform random matching, the above inequalities hold for any x; if strategies are updated according to the interact-and-imitate protocol, then this game only admits switches from i to j, therefore violating payoff monotonicity. Proposition 2 can have important consequences, including the survival of pure strategies that are strictly dominated.Survival of strictly dominated strategiesAn recurring topic in evolutionary game theory is to what extent does support exist for the idea that strictly dominated strategies will not be played. It has been shown that if strategy i does not survive the iterated elimination of pure strategies strictly dominated by other pure strategies, then the fraction of the population playing i will converge to zero in all payoff monotone dynamics37,38. This result does not hold in our case, as PIID is not payoff monotone.More precisely, under PIID, a strictly dominated strategy may be supreme and, therefore, not only survive but even end up being adopted by the whole population. This suggests that from an evolutionary perspective, support for the elimination of dominated strategies may be weaker than is often thought. Our result contributes to the literature on the conditions under which evolutionary dynamics fail to eliminate strictly dominated strategies in some games, examining a case which has not yet been studied39.To see that a strictly dominated strategy may be supreme, consider the simple example shown in Fig. 1. Here each agent has a strictly dominant strategy to play A; however, since the payoff from playing B against A exceeds that from playing A against B, strategy B is supreme. Thus, by Proposition 1, the population state in which all agents choose B is globally asymptotically stable.Figure 1A game where the supreme strategy is strictly dominated.Full size imageFigure 1 can also be used to comment on the relation between a supreme strategy and an evolutionary stable strategy, which is a widely used concept in evolutionary game theory40,41. Indeed, while B is the supreme strategy, A is the unique evolutionary stable strategy because it is strictly dominant. However, if F(B, A) were reduced below 2, holding everything else constant, then B would become both supreme and evolutionary stable. We therefore conclude that no particular relation holds between evolutionary stability and supremacy: neither one property implies the other, nor are they incompatible.ApplicationsHaving obtained general results for the class of finite symmetric games, we now restrict the discussion to the evolution of behavior in social dilemmas. We show that if the conditions of Proposition 1 are met, then inefficient conventions emerge in the Prisoner’s Dilemma, Stag Hunt, Minimum Effort, and Hawk–Dove games. Furthermore, this result holds both without and with the assumption that agents interact assortatively.Ineffectiveness of assortmentConsider the (2 times 2) game represented in Fig. 2. If (c > a > d > b), then mutual cooperation is Pareto superior to mutual defection but agents have a dominant strategy to defect. The resulting stage game is the Prisoner’s Dilemma, whose unique Nash equilibrium is (B, B). Moreover, since (F (B,A) > F(A,B)), B is the supreme strategy and the population state in which all agents defect is globally asymptotically stable.We stress that defection emerges in the long run for every matching rule satisfying our assumptions, and therefore also in the case of assortative interactions. Assortment reflects the tendency of similar people to clump together, and can play an important role in the evolution of cooperation42,43,44,45. Intuitively, when agents meet assortatively, the risk of cooperating in a social dilemma may be offset by a higher probability of playing against other cooperators. However, under PIID, this is not the case: the decision whether to adopt a strategy or not is independent of expected payoffs, and like-with-like interactions have no effect except to reduce the frequency of switches from A to B.Figure 2A (2 times 2) stage game.Full size imageEmergence of the maximin conventionIf (a > c > b), (a > d) and (d > b), then the game in Fig. 2 becomes a Stag Hunt game, which contrasts risky cooperation and safe individualism. The payoffs are such that both (left( A, Aright) ) and (left( B, Bright) ) are strict Nash equilibria, that (left( A, Aright) ) is Pareto superior to (left( B, Bright) ), and that B is the maximin strategy, i.e., the strategy which maximizes the minimum payoff an agent could possibly receive. We also assume that (a + c ne c + d), so that one of A and B is risk dominant46. If (a + b > c + d), then A (Stag) is both payoff and risk dominant. When the opposite inequality holds, the risk dominant strategy is B (Hare).Since (F (B,A) > F(A,B)), B is supreme independently of whether or not it is risk dominant to cooperate. This can result in large inefficiencies because, in the long run, the process will converge to the state in which all agents play the riskless strategy regardless of how rewarding social coordination is. As in the case of the Prisoner’s Dilemma, this holds for all matching rules satisfying our assumptions.Evolution of effort exertionIn a minimum effort game, agents simultaneously choose a strategy i, usually interpreted as a costly effort level, from a finite subset S of ({mathbb{R}}). An agent’s payoff depends on her own effort and on the minimum effort in the pair:$$begin{aligned} F left( i, j right) = alpha min left{ i, j right} – beta i , end{aligned}$$where (beta > 0) and (alpha > beta ) are the cost and benefit of effort, respectively. From a strategic viewpoint, this game can be seen as an extension of the Stag Hunt to cases where there are more than two actions. The best response to a choice of j by the opponent is to choose j as well, and coordinating on any common effort level gives a Nash equilibrium. Nash outcomes can be Pareto-ranked, with the highest-effort equilibrium being the best possible outcome for all agents. Thus, choosing a high i is rationalizable and potentially rewarding but may also result in a waste of effort.Under PIID, any (i > j) implies (phi _{ij} < 0) by Assumption 1, meaning that agents will tend to imitate the opponent when the opponent’s effort is lower than their own. The supreme strategy is therefore to exert as little effort as possible, and the population state in which all agents choose the minimum effort level is the unique globally asymptotically stable state.Emergence of aggressive behaviorConsider again the payoff matrix shown in Fig. 2. If (c > a > b > d), then the stage game is a Hawk–Dove game, which is often used to model the evolution of aggressive and sharing behaviors. Interactions can be framed as disputes over a contested resource. When two Doves (who play A) meet, they share the resource equally, whereas two Hawks (who play B) engage in a fight and suffer a cost. Moreover, when a Dove meets a Hawk, the latter takes the entire prize. Again we have that (F (A,B) < F(B,A)), implying that B is the supreme strategy and that the state where all agents play Hawk is the sole asymptotically stable state.The inefficiency that characterizes the (B, B) equilibrium in the Hawk–Dove game arises from the cost that Hawks impose on one another. This can be viewed as stemming from the fact that neither agent owns the resource prior to the interaction or cares about property. A way to overcome this problem may be to introduce a strategy associated with respect for ownership rights, the Bourgeois, who behaves as a Dove or Hawk depending on whether or not the opponent owns the resource41. If we make the standard assumption that each member of a pair has a probability of 1/2 to be an owner, then in all interactions where a Bourgeois is involved there is a 50 percent chance that she will behave hawkishly (i.e., fight for control over the resource) and a 50 percent chance that she will act as a Dove.Let R and C denote the agent chosen as row and column player, respectively, and let (omega _R) and (omega _C) be the states of the world in which R and C owns the resource. The payoffs of the resulting Hawk–Dove–Bourgeois game are shown in Fig. 3. If agents behave as expected payoff maximizers, then All Bourgeois can be singled out as the unique asymptotically stable state. Under PIID, this is not so; depending on who owns the resource, an agent playing C against an opponent playing B may either fight or avoid conflict and let the opponent have the prize. It is easy to see that (F left( C, B mid omega _R right) = F left( B,C mid omega _C right) = d), meaning that the payoff from playing C against B, conditional on owning the resource, equals the payoff from playing B against C conditional on not being an owner. In contrast, the payoff from playing C against B, conditional on not owning the resource, is always worse than that of the opponent, i.e., (F left( C, B mid omega _C right) = b < c = F left( B, C mid omega _R right) ). Thus, in every state of the world, B (Hawk) yields a payoff that is greater or equal to that from C (Bourgeois). Moreover, since (F left( B,A right) > F left( A, B right) ) in both states of the world, strategy B is weakly supreme by Definition 4, and play unfolds as an escalation of hawkishness and fights.Figure 3The Hawk–Dove–Bourgeois game.Full size image More

1.Igwe, J. C. & Onyegbado, C. C. A review of palm oil mill effluent (pome) water treatment. Glob. J. Environ. Res. 1, 54–62 (2007).
Google Scholar 
2.World Wild Fund (WWF). Overview WWF Statement on the 2020 Palm Oil Buyers Scorecard. https://www.worldwildlife.org/industries/palm-oil (2020). Accessed 22 Feb 2021.3.CNUCED. Huile de palme. New York. https://www.surunctad.org/commodities (2016). Accessed 10 Jan 2020.4.Hassan, M. A., Njeshu, G., Raji, A., Zhengwuvi, L. & Salisu, J. Small-Scale Palm Oil Processing in West and Central Africa: Development and Challenges. J. Appl. Sci. Environ. Sust. 2, 102–114 (2016).
Google Scholar 
5.Bala, J. D., Lalung, J., Al-Gheethi, A. A. S., Kaizar, H. & Ismail, N. Reduction of organic load and biodegradation of palm oil mill effluent by aerobic indigenous mixed microbial consortium isolated from palm oil mill effluent (POME). Water Conserv. Sci. Eng. 3, 139. https://doi.org/10.1007/s41101-018-0043-9 (2018).Article 

Google Scholar 
6.Nwoko, O. C., Ogunyemi, S. & Nkwocha, E. E. Effect of pre-treatment of palm oil mill effluent (POME) and cassava mill effluent (CME) on the growth of tomato (Lycopersicum esculentum). J. Appl. Sci. Environ. 14, 67. https://doi.org/10.4314/JASEM.V14I1.56493 (2010).Article 

Google Scholar 
7.Singh, G., Huan, L. K., Leng, T. & Kow D. L. Oil Palm and the Environment: A Malaysian Perspective. (Kuala Lumpur,
Malaysia, Malaysian Oil Palm Growers’ Council, 1999).8.Poku, K. Small-Scale Palm Oil Processing in Africa. Fao Agricultural Services Bulletin 148. http://www.fao.org/3/Y4355E/y4355e00.htm (2002) (ISSN 1010-1365). Accessed 22 Feb 2021.9.Ibekwe, A. M., Grieve, C. M. & Lyon, S. R. Characterization of microbial communities and composition in constructed dairy wetland wastewater effluent. Appl. Environ. Microbiol. 69, 5060. https://doi.org/10.1128/AEM.69.9.5060-5069.2003 (2003).CAS 
Article 
PubMed 

Google Scholar 
10.Sharuddin, S. S. et al. Bacterial community shift revealed Chromatiaceae and Alcaligenaceae as potential bioindicators in the receiving river due to palm oil mill effluent final discharge. Ecol. Indic. 82, 526–529. https://doi.org/10.1016/j.ecolind.2017.07.038 (2017).CAS 
Article 

Google Scholar 
11.CIAPOL. Arrêté N°011264/MINEEF/CIAPOL/SDIIC du 04 Nov.2008 portant réglementation des rejets et emissions des installations classées pour la protection de l’environnement, 11 (2008).
12.Soleimaninanadegani, M. & Manshad, S. Enhancement of biodegradation of palm oil mill effluents by local isolated microorganisms. Int. Sch. Res. Notices. 2014, Article ID 727049. https://doi.org/10.1155/2014/727049 (2014).Article 

Google Scholar 
13.Nwachukwu, J. N., Njoku, U. O., Agu, C. V., Okonkwo, C. C. & Obidiegwu, C. J. Impact of palm oil mill effluent (POME) contamination on soil enzyme activities and physicochemical properties. Res. J. Environ. Toxicol. 12, 34–41. https://doi.org/10.3923/rjet.2018.34.41 (2018).CAS 
Article 

Google Scholar 
14.Hii, K. L., Yeap, S. P. & Mashitah, M. D. Cellulase production from palm oil mill effluent in Malaysia: Economical and technical perspectives. Eng. Life Sci. 12, 7–28. https://doi.org/10.1002/elsc.201000228 (2012).CAS 
Article 

Google Scholar 
15.Ma, Q. et al. Identification of the microbial community composition and structure of coal-mine wastewater treatment plants. Microbiol. Res. 175, 1–5. https://doi.org/10.1016/j.micres.2014.12.013 (2015).CAS 
Article 
PubMed 

Google Scholar 
16.Wang, X., Hu, M., Xia, Y., Wen, X. & Kun, D. K. Pyrosequencing analysis of bacterial diversity in 14 wastewater treatment systems in China. Bioresour. Technol. 78, 7042–7047. https://doi.org/10.1128/AEM.01617-12 (2012).CAS 
Article 

Google Scholar 
17.Wang, Z. et al. Abundance and diversity of bacterial nitrifiers and denitrifiers and their functional genes in tannery wastewater treatment plants revealed by high-throughput sequencing. PLoS One 9, e113603. https://doi.org/10.1371/journal.pone.0113603 (2014).ADS 
CAS 
Article 
PubMed 

Google Scholar 
18.Caporaso, J. et al. Ultra-high-throughput microbial community analysis on the Illumina HiSeq and MiSeq platforms. ISME J. 6, 1621–1624. https://doi.org/10.1038/ismej.2012.8 (2012).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
19.Rana, S., Singh, L., Wahid, Z. & Liu, H. A recent overview of palm oil mill effluent management via bioreactor configurations. Curr. Pollut. Rep. 3, 254–267. https://doi.org/10.1007/s40726-017-0068-2 (2017).CAS 
Article 

Google Scholar 
20.Vuono, D. C. et al. Disturbance and temporal partitioning of the activated sludge metacommunity. ISME J. 9, 425–435. https://doi.org/10.1038/ismej.2014.139 (2015).CAS 
Article 
PubMed 

Google Scholar 
21.Jang, H. M., Kim, J. H., Ha, J. H. & Park, J. M. Bacterial and methanogenic archaeal communities during the single-stage anaerobic digestion of high-strength food wastewater. Bioresour. Technol. 165, 174–182. https://doi.org/10.1016/j.biortech.2014.02.028 (2014).CAS 
Article 
PubMed 

Google Scholar 
22.Mohd-Nor, D. et al. Dynamics of microbial populations responsible for biodegradation during the full-scale treatment of palm oil mill effluent. Microbes Environ. 34, 121. https://doi.org/10.1264/jsme2.ME18104 (2019).Article 
PubMed 
PubMed Central 

Google Scholar 
23.Sun, Z. et al. Identification and characterization of the dominant lactic acid bacteria from kurut: The naturally fermented yak milk in Qinghai, China. J. Gen. Appl. Microbiol. 56, 1–10. https://doi.org/10.2323/jgam.56.1 (2010).Article 
PubMed 

Google Scholar 
24.Webster, N. S. & Taylor, M. W. Marine sponges and their microbial symbionts: Love and other relationships. Environ. Microbiol. 14, 335–346 (2012).CAS 
Article 

Google Scholar 
25.Morrow, K. M., Fiore, C. L. & Lesser, M. P. Environmental drivers of microbial community shifts in the giant barrel sponge, Xestospongia muta, over a shallow to mesophotic depth gradient. Environ. Microbiol. 18, 2025–2038. https://doi.org/10.1111/1462-2920.13226 (2016).CAS 
Article 
PubMed 

Google Scholar 
26.Parman, A., Isa, M. N. M., Farah, F. B., Noorbatcha, B. A. & Salleh, H. M. Comparative metagenomics analysis of palm oil mill effluent (pome) using three different bioinformatics pipelines. IIUM Eng. J. 20, 1–11. https://doi.org/10.31436/iiumej.v20i1.909 (2019).Article 

Google Scholar 
27.Mwaikono, K. S. et al. High-throughput sequencing of 16S rRNa gene reveals substantial bacterial diversity on the municipal dumpsite. BMC Microbiol. 16, 145. https://doi.org/10.1186/s12866-016-0758-8 (2016).Article 
PubMed 

Google Scholar 
28.Silva-Bedoya, L. M., Sánchez-Pinzón, M. S., Cadavid-Restrepo, G. E. & Moreno-Herrera, C. X. Bacterial community analysis of an industrial wastewater treatment plant in Colombia with screening for lipid-degrading microorganisms. Microbiol. Res. 192, 313. https://doi.org/10.1016/j.micres.2016.08.006 (2016).CAS 
Article 
PubMed 

Google Scholar 
29.Lam, M. K. & Lee, K. T. Renewable and sustainable bioenergies production from palm oil mill effluent (POME): Win–win strategies toward better environmental protection. Biotechnol. Adv. 29, 124–141. https://doi.org/10.1016/j.biotechadv.2010.10.001 (2011).CAS 
Article 
PubMed 

Google Scholar 
30.Baharuddin, A. S., Wakisaka, M., Shirai, A.-A.Y.S., Abdul, R. & Hassan, M. A. Co-composting of empty fruit bunches and partially treated palm oil mill effluents in pilot scale. Int. J. Agric. Res. 4, 69–78. https://doi.org/10.3923/ijar.2009.69.78 (2009).CAS 
Article 

Google Scholar 
31.Morikawa-Sakura, M. S. et al. Application of Lactobacillus plantarum ATCC 8014 for wastewater treatment in fisheries industry processing. Jpn. J. Water Treat. Biol. 49, 1–10. https://doi.org/10.2521/jswtb.49.1 (2013).Article 

Google Scholar 
32.Ren, Z., You, W., Wu, S., Poetsch, A. & Xu, C. Secretomic analyses of Ruminiclostridium papyrosolvens reveal its enzymatic basis for lignocellulose degradation. Biotechnol. Biofuels 12, 183. https://doi.org/10.1186/s13068-019-1522-8 (2019).CAS 
Article 
PubMed 

Google Scholar 
33.Lee, J. Z., Logan, A., Terry, S. & Spear, J. R. Microbial response to single-cell protein production and brewery wastewater treatment. Microb. Biotechnol. 8, 65. https://doi.org/10.1111/1751-7915.12128 (2015).CAS 
Article 
PubMed 

Google Scholar 
34.Ye, L. & Zhang, T. Bacterial communities in different sections of a municipal wastewater treatment plant revealed by 16S rDNA 454 pyrosequencing. Appl. Microbiol. Biotechnol. 97, 2681–2690 (2013).CAS 
Article 

Google Scholar 
35.Stubbs, S., Mao, L., Waddington, R. J. & Embery, G. Hydrolytic and depolymerising enzyme activity of Prevotella intermedia and Prevotella nigrescens. Oral Dis. 2, 272. https://doi.org/10.1111/j.1601-0825.1996.tb00237.x (1996).CAS 
Article 
PubMed 

Google Scholar 
36.Komagata, K., Iino, T. & Yamada, Y. The family Acetobacteraceae. In The Prokaryotes (eds Rosenberg, E. et al.) 3–78 (Springer, 2014).Chapter 

Google Scholar 
37.Pires, J. F., Cardoso, L. S., Schwan, R. F. & Silva, C. F. Diversity of microbiota found in coffee processing wastewater treatment plant. World J. Microbiol. Biotechnol. 33, 211. https://doi.org/10.1007/s11274-017-2372-9 (2017).Article 
PubMed 

Google Scholar 
38.Song, Z. Q., Wang, F. P. & Zhi, X. Y. Bacterial and archaeal diversities in Yunnan and Tibetan hot springs, China. Environ. Microbiol. 15, 1160–1175 (2013).CAS 
Article 

Google Scholar 
39.Li, J., Liu, R., Tao, Y. & Li, G. Archaea in wastewater treatment: Current research and emerging technology. Archaea 2018, 1. https://doi.org/10.1155/2018/6973294 (2018).CAS 
Article 

Google Scholar 
40.Khan, M. A., Khan, S. T. & Sequeira, M. C. Comparative analysis of bacterial and archaeal population structure by illumina sequencing of 16S rRNA genes in three municipal anaerobic sludge digesters. Res. Sq. https://doi.org/10.21203/rs.3.rs-60183/v1 (2020).Article 

Google Scholar 
41.Mladenovska, Z., Dabrowski, S. & Ahring, B. K. Anaerobic digestion of manure and mixture of manure with lipids: Biogas reactor performance and microbial community analysis. Water Sci. Technol. 48, 271–278 (2013).Article 

Google Scholar 
42.Gerardi, M. H. Wastewater Bacteria (Wiley, 2006).Book 

Google Scholar 
43.Parada, A. E., Needham, D. M. & Fuhrman, J. A. Every base matters: Assessing small subunit rRNA primers for marine microbiomes with mock communities, time series and global field samples. Environ. Microbiol. 18, 1403–1414. https://doi.org/10.1111/1462-2920.13023 (2016).CAS 
Article 
PubMed 

Google Scholar 
44.Andrews, S. FastQC: a quality control tool for high throughput sequence data (Online). https://www.bioinformatics.babraham.ac.uk/projects/fastqc (2010). Accessed 15 Sept 2019.45.R Development Core Team. R: A Language and Environment for Statistical Computing. (R Foundation for Statistical Computing, 2019). Accessed 8 Jan 2020.46.Callahan, B. J. et al. DADA2: High-resolution sample inference from illumina amplicon data. Nat. Methods 13, 581. https://doi.org/10.1038/nmeth.3869 (2016).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
47.Quast, C. et al. The SILVA ribosomal RNA gene database project: Improved data processing and web-based tools. Nucleic Acids Res. 41, 590. https://doi.org/10.1093/nar/gks1219 (2012).CAS 
Article 

Google Scholar 
48.Paradis, E., Julien, C. & Korbinian, S. APE: Analyses of phylogenetics and evolution in R language. Bioinformatics 20, 289. https://doi.org/10.1093/bioinformatics/btg412 (2004).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
49.McMurdie, P. J. & Holmes, S. Phyloseq: An R package for reproducible interactive analysis and graphics of microbiome census data. PLoS One 8, 61217. https://doi.org/10.1371/journal.pone.0061217 (2013).ADS 
CAS 
Article 

Google Scholar 
50.Oksanen, J. F. et al. vegan: Community Ecology Package. R package version 2.4-0. https://CRAN.R-project.org/package=vegan (2018). Accessed 8 Jan 2020.51.Lahti, L. & Sudarshan, S. Tools for microbiome analysis in R. Version 1.10.0. https://www.microbiome.github.com/microbiome (2017). Accessed 8 Jan 2020.52.Kenkel, N. C. & Orloci, L. Applying metric and nonmetric multidimensional scaling to ecological studies: Some new results. Ecology 67, 919. https://doi.org/10.2307/1939814 (1986).Article 

Google Scholar 
53.Benjamini, Y. & Hochberg, Y. Controlling the false discovery rate: A practical and powerful approach to multiple testing. J. R. Stat. Soc. 57, 289–300 (1995).MathSciNet 
MATH 

Google Scholar  More

1.Fargione, J. E. et al. Natural climate solutions for the United States. Sci. Adv. 4, eaat1869 (2018).ADS 
PubMed 
PubMed Central 
Article 

Google Scholar 
2.Cardinale, B. J. et al. Biodiversity loss and its impact on humanity. Nature 486, 59–67 (2012).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
3.Conti, G. & Díaz, S. Plant functional diversity and carbon storage – an empirical test in semi-arid forest ecosystems. J. Ecol. 101, 18–28 (2013).CAS 
Article 

Google Scholar 
4.Mensah, S., Veldtman, R., Assogbadjo, A. E., Glèlè Kakaï, R. & Seifert, T. Tree species diversity promotes aboveground carbon storage through functional diversity and functional dominance. Ecol. Evol. 6, 7546–7557 (2016).PubMed 
Article 
PubMed Central 

Google Scholar 
5.Islam, M., Dey, A. & Rahman, M. Effect of Tree Diversity on Soil Organic Carbon Content in the Homegarden Agroforestry System of North-Eastern Bangladesh. Small-scale 14, 91–101 (2015).Article 

Google Scholar 
6.Poorter, L. et al. Diversity enhances carbon storage in tropical forests. Glob. Ecol. Biogeogr. 24, 1314–1328 (2015).Article 

Google Scholar 
7.Rahman, M. M., Kabir, M. E., Jahir Uddin Akon, A. S. M. & Ando, K. High carbon stocks in roadside plantations under participatory management in Bangladesh. Glob. Ecol. Conserv. 3, 412–423 (2015).Article 

Google Scholar 
8.McKee, K. L. Interspecific Variation in Growth, Biomass Partitioning, and Defensive Characteristics of Neotropical Mangrove Seedlings: response To Light and Nutrient Availability. Am. J. Bot. 82, 299–307 (1995).Article 

Google Scholar 
9.Kauffman, J. B. et al. Total ecosystem carbon stocks of mangroves across broad global environmental and physical gradients. Ecol. Monogr. 90, 1–18 (2020).Article 

Google Scholar 
10.Tinh, P. H. et al. A comparison of soil carbon stocks of intact and restored mangrove forests in Northern Vietnam. Forests 11, 1–10 (2020).
Google Scholar 
11.Saintilan, N. Above- and below-ground biomasses of two species of mangrove on the Hawkesbury River stuary, New South Wales. Mar. Freshw. Res. 48, 147–152 (1997).CAS 
Article 

Google Scholar 
12.Tamooh, F. et al. Below-ground root yield and distribution in natural and replanted mangrove forests at Gazi bay, Kenya. Ecol. Manag. 256, 1290–1297 (2008).Article 

Google Scholar 
13.MacKenzie, R. A. et al. Sedimentation and belowground carbon accumulation rates in mangrove forests that differ in diversity and land use: a tale of two mangroves. Wetl. Ecol. Manag. 24, 245–261 (2016).Article 

Google Scholar 
14.Atwood, T. B. et al. Global patterns in mangrove soil carbon stocks and losses. Nat. Clim. Chang. 7, 523–528 (2017).ADS 
CAS 
Article 

Google Scholar 
15.Lovelock, C. E. & Duarte, C. M. Dimensions of blue carbon and emerging perspectives. Biol. Lett. 15, 1–5 (2019).Article 

Google Scholar 
16.Donato, D. C. et al. Mangroves among the most carbon-rich forests in the tropics. Nat. Geosci. 4, 293–297 (2011).ADS 
CAS 
Article 

Google Scholar 
17.Adame, M. F. et al. Root biomass and production of mangroves surrounding a karstic oligotrophic coastal lagoon. Ecol. Manag. 256, 1290–1297 (2014).
Google Scholar 
18.Sharma, S. et al. The impacts of degradation, deforestation and restoration on mangrove ecosystem carbon stocks across Cambodia. Sci. Total Environ. 706, 135416 (2020).19.Ruiz-benito, P. et al. Diversity increases carbon storage and tree productivity in Spanish forests. 1–12 (2013) https://doi.org/10.1111/geb.12126.20.Mace, G. Biodiversity Policy Challenges. (2009) https://doi.org/10.1126/science.1180935.21.Tilman, D. et al. The influence of functional diversity and composition on ecosystem processes. Sci. (80-.) 277, 1300–1302 (1997).CAS 
Article 

Google Scholar 
22.Cavanaugh, K. C. et al. Carbon storage in tropical forests correlates with taxonomic diversity and functional dominance on a global scale: Biodiversity and aboveground carbon storage. Glob. Ecol. Biogeogr. 23, 563–573 (2014).Article 

Google Scholar 
23.Ruiz-Jaen, M. C. & Potvin, C. Can we predict carbon stocks in tropical ecosystems from tree diversity? Comparing species and functional diversity in a plantation and a natural forest. N. Phytol. 189, 978–987 (2011).Article 

Google Scholar 
24.Ali, A., Chen, H. Y. H., You, W.-H. & Yan, E.-R. Multiple abiotic and biotic drivers of aboveground biomass shift with forest stratum. Ecol. Manag. 436, 1–10 (2019).Article 

Google Scholar 
25.Chave, J. et al. Towards a worldwide wood economics spectrum. Ecol. Lett. 12, 351–366 (2009).PubMed 
Article 

Google Scholar 
26.Franck, J. & Jérôme, C. Inferring the parameters of the neutral theory of biodiversity using phylogenetic information and implications for tropical forests. Ecol. Lett. 12, 239–248 (2009).Article 

Google Scholar 
27.Poorter, L., Bongers, L. & Bongers, F. Architecture of 54 moist-forest tree species: traits, trade-offs, and functional groups. Ecology 87, 1289–1301 (2006).PubMed 
Article 

Google Scholar 
28.Fatoyinbo, T. E., Simard, M., Washington-Allen, R. A. & Shugart, H. H. Landscape-scale extent, height, biomass, and carbon estimation of Mozambique’s mangrove forests with Landsat ETM+ and Shuttle Radar Topography Mission elevation data. J. Geophys. Res. Biogeosci. 113, G02S06 (2008).29.Díaz, S. et al. The global spectrum of plant form and function. Nature 529, 167–171 (2016).ADS 
PubMed 
Article 
CAS 

Google Scholar 
30.Simard, M. et al. Mangrove canopy height globally related to precipitation, temperature and cyclone frequency. Nat. Geosci. 12, 40–45 (2019).ADS 
CAS 
Article 

Google Scholar 
31.Díaz, S. et al. The plant traits that drive ecosystems: evidence from three continents. J. Veg. Sci. 15, 295–304 (2004).Article 

Google Scholar 
32.Lasky, J. R., Uriarte, M., Boukili, V. K. & Chazdon, R. L. Trait-mediated assembly processes predict successional changes in community diversity of tropical forests. Proc. Natl Acad. Sci. U. S. A. 111, 5616–5621 (2014).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
33.Grime, J. P. Benefits of plant diversity to ecosystems: immediate, filter and founder effects. J. Ecol. 86, 902–910 (1998). vol.Article 

Google Scholar 
34.Pradisty, N. A., Amir, A. A. & Zimmer, M. Plant species- and stage-specific differences in microbial decay of mangrove leaf litter: the older the better? Oecologia (2021) https://doi.org/10.1007/s00442-021-04865-3.35.Hossain, M. et al. Nutrient Dynamics Associated with Leaching and Microbial Decomposition of Four Abundant Mangrove Species Leaf Litter of the Sundarbans, Bangladesh. Wetlands 34, 439–448 (2014).Article 

Google Scholar 
36.Chanda, A. et al. Mangrove associates versus true mangroves: a comparative analysis of leaf litter decomposition in Sundarban. Wetl. Ecol. Manag. 24, 293–315 (2016).Article 

Google Scholar 
37.Alongi, D. M. Global Significance of Mangrove Blue Carbon in Climate Change Mitigation. Sci 2, 67 (2020).Article 

Google Scholar 
38.Lovelock, C. E. & Reef, R. Variable Impacts of Climate Change on Blue. Carbon One Earth 3, 195–211 (2020).Article 

Google Scholar 
39.Alongi, D. Impact of Global Change on Nutrient Dynamics in Mangrove Forests. Forests 9, 596 (2018).Article 

Google Scholar 
40.Rahman, M. M. & Rahaman, M. M. Impacts of Farakka barrage on hydrological flow of Ganges river and environment in Bangladesh. Sustain. Water Resour. Manag. 1–14 (2017) https://doi.org/10.1007/s40899-017-0163-y.41.Gilman, E. L., Ellison, J., Duke, N. C. & Field, C. Threats to mangroves from climate change and adaptation options: a review. Aquat. Bot. 89, 237–250 (2008).Article 

Google Scholar 
42.Kirui, B., Kairo, J., Skov, M., Mencuccini, M. & Huxham, M. Effects of species richness, identity and environmental variables on growth in planted mangroves in Kenya. Mar. Ecol. Prog. Ser. 465, 1–10 (2012).ADS 
Article 

Google Scholar 
43.Ball, M. C. Ecophysiology of mangroves. Trees 2, 129–142 (1988). vol.Article 

Google Scholar 
44.Grace, J. B. et al. Integrative modelling reveals mechanisms linking productivity and plant species richness. Nature 529, 390–393 (2016).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
45.Rahman, M. M., Khan, M. N. I., Hoque, A. K. F. & Ahmed, I. Carbon stock in the Sundarbans mangrove forest: spatial variations in vegetation types and salinity zones. Wetl. Ecol. Manag. 23, 269–283 (2015).Article 

Google Scholar 
46.Chave, J. et al. Tree allometry and improved estimation of carbon stocks and balance in tropical forests. Oecologia 145, 87–99 (2005).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
47.Komiyama, A., Poungparn, S. & Kato, S. Common allometric equations for estimating the tree weight of mangroves. J. Trop. Ecol. 21, 471–477 (2005).Article 

Google Scholar 
48.Hossain, M., Siddique, M. R. H., Saha, S. & Abdullah, S. M. R. Allometric models for biomass, nutrients and carbon stock in Excoecaria agallocha of the Sundarbans, Bangladesh. Wetl. Ecol. Manag. 23, 765–774 (2015).CAS 
Article 

Google Scholar 
49.Loreau, M., Hector, A., The, C. J. & Funct, D. Partitioning Selection and Complementarity in Biodiversity Experiments Partitioning selection and complementarity in biodiversity experiments. (2001) https://doi.org/10.1038/35083573.50.Phelps, J., Webb, E. L. & Adams, W. M. Biodiversity co-benefits of policies to reduce forest-carbon emissions. Nat. Clim. Chang. 2, 497–503 (2012).ADS 
Article 

Google Scholar 
51.Zimmer, M. Ecosystem design: when mangrove ecology meets human needs. Coast. Res. Libr. 25, 367–376 (2018).Article 

Google Scholar 
52.Rahman, M. S., Sass-Klaassen, U., Zuidema, P. A., Chowdhury, M. Q. & Beeckman, H. Salinity drives growth dynamics of the mangrove tree Sonneratia apetala Buch. -Ham. in the Sundarbans, Bangladesh. Dendrochronologia 62, 125711 (2020).Article 

Google Scholar 
53.Suwa, R., Deshar, R. & Hagihara, A. Forest structure of a subtropical mangrove along a river inferred from potential tree height and biomass. Aquat. Bot. 91, 99–104 (2009).Article 

Google Scholar 
54.Sparks, D. L. et al. Total Carbon, Organic Carbon, and Organic Matter. in SSSA Book Series (Soil Science Society of America, American Society of Agronomy, 1996).55.Yakub, M., Omar Ali, M. & Bhattacharjee, D. K. Strength properties of some Bangladesh timber species. (Govt. of the People’s Republic of Bangladesh, Forest Research Institute, 1972).56.Nandy (Datta), P. & Ghose, M. Photosynthesis and water-use efficiency of some mangroves from Sundarbans, India. J. Plant Biol. 44, 213–219 (2001).Article 

Google Scholar 
57.Shannon, C. E. A mathematical theory of communication. Bell Syst. Tech. J. 27, 379–423 (1948).MathSciNet 
MATH 
Article 

Google Scholar 
58.Laliberté, E. & Legendre, P. A distance-based framework for measuring functional diversity from multiple traits. Ecology 91, 299–305 (2010).Article 

Google Scholar 
59.Oliveira, A. D. et al. Does functional trait diversity predict above-ground biomass and productivity of tropical forests? Testing three alternative hypotheses. 191–201 (2015) https://doi.org/10.1111/1365-2745.12346.60.Demján, P. & Dreslerová, D. Modelling distribution of archaeological settlement evidence based on heterogeneous spatial and temporal data. J. Archaeol. Sci. 69, 100–109 (2016).Article 

Google Scholar 
61.Hossain, G. M. & Bhuiyan, M. A. H. Spatial and temporal variations of organic matter contents and potential sediment nutrient index in the Sundarbans mangrove forest, Bangladesh. KSCE J. Civ. Eng. 20, 163–174 (2016).Article 

Google Scholar 
62.Ggraham, M. H. Confronting Multicollinearity in Ecological Multiple Regression. Ecology 84, 2809–2815 (2003).63.Rosseel, Y. Lavaan:anRpackageforstructuralequationmodeling and more. Version 0.5-12 (BETA). J. Stat. Softw. 48, 1–36 (2012).Article 

Google Scholar 
64.Grace, J. B. & Bollen, K. A. Interpreting the Results from Multiple Regression and Structural Equation Models. Bull. Ecol. Soc. Am. 86, 283–295 (2005).Article 

Google Scholar  More

Network reconstruction from observational dataThe typical field study of plant–pollinator interactions involves recording instances of potential pollinators (such as insects) visiting plants within a prescribed observation area and over a prescribed period of time. We will refer to these records as visitation data. Network ecologists analyze visitation data by constructing networks of plant and pollinator species, where a connection between two species indicates that a plant-pollinator interaction exists between them.However, the meaning of edges in ecological networks is not always clear31. One popular way to transform visitation data into networks is to connect two species when they interact “enough”—say when a pollinator species is seen on the reproductive organ of a plant species a specified number of times—but in this case the precise meaning of an edge will depend on the details of the data collection and the choices made in the analysis. How many visits do we take as evidence of a plant–pollinator interaction? A single visit is probably not enough—it might well be an error or misobservation. Is two enough, or ten, or a hundred? And what about observations that were missed entirely? Other methods of analysis transform the data in different ways, for instance encoding them as weighted networks, possibly with some statistical processing along the way32. Even in this case, however, the edges still just count numbers of visits (perhaps transformed in some way), so the resulting networks are effectively histograms in disguise, recording only potential interactions rather than true biological connections.A more principled approach to network construction begins with a clear definition of what relationship (or relationships) a network’s edges encode33. We argue that network ecology often calls for a network of preferred interactions. In the context of plant-pollinator networks the edges of such a network indicate that pollinators preferentially visit certain plant species and they encode a variety of mechanisms that constrain species interactions, such as temporal or spatial uncoupling (i.e., species that do not co-occur in either time or space), constraints due to trait mismatches (e.g., proboscis size very different from corolla size), and physiological-biochemical constraints that prevent the interactions (e.g., chemical barriers). (One can regard preferred interactions as being the opposite of the “forbidden links” described in refs. 34,35,36). Preferred interactions are arguably the relevant ones for instance when analyzing the reaction of a network to abrupt changes: when one removes a plant species from a system, for example, the pollinators that prefer it will have to modify their behavior7,37,38. The interactions we consider are binary—either a species prefers another species or it doesn’t—so the network does not encode varying strengths of interaction.While the data gathered in a typical field study are certainly reflective of preferred interactions, they are, for many reasons, not perfect measurements of networks of preferred interactions13,17. First, there may be observational errors. While the observers performing the work are usually highly trained individuals, they may nonetheless make mistakes. They may confuse one species for another, which is particularly easy to do for small-bodied insects, or smaller species may be overlooked altogether. Observers may make correct observations but record them wrongly. And there will be statistical fluctuations in the number of visits of an insect species to a plant species over any finite time. For rare interactions there may even be no visits at all if we are unlucky. The insects themselves may also appear to make “mistakes” by visiting plants that they typically do not pollinate. These and other factors mean that the record of observed visits is an inherently untrustworthy guide to the true structure of the network of preferred interactions. Here we develop a statistical method for making estimates of network structure despite these limitations of the data.Model of plant–pollinator dataConsider a typical plant–pollinator study in which some number np of plant species, labeled by i = 1…np, and some number na of animal pollinator species, labeled by j = 1…na, are under observation for a set amount of time, producing a record of observed visits such that Mij is the number of times plant species i is visited by pollinator species j. Collectively the Mij can be regarded as a data matrix M with np rows and na columns. This is the input to our calculation.The unknown quantity, the thing we would like to understand, is the network of plant–pollinator interactions. We can think of this network as composed of two sets of nodes, one representing plant species and the other pollinator species, with connections or edges joining each pollinator to the plants it pollinates. In the language of network science this is a bipartite network, meaning that edges run only between nodes of unlike kinds—plants and pollinators—and never between two plants or two pollinators. Such a network can be represented by a second matrix B, called the incidence matrix, with the same size as the data matrix and elements Bij = 1 if plant i is preferentially visited by pollinator j and 0 otherwise.The question we would like to answer is this: What is the structure of the network, represented by B, given the data M? It is not straightforward to answer this question directly, but it is relatively easy to answer the reverse question. If we imagine that we know B, then we can say what the probability is that we make a specific set of observations M. And if we can do this then the methods of Bayesian inference allow us to invert the calculation and compute B from a knowledge of M and hence achieve our goal. The procedure is as follows.Consider a specific plant-pollinator species pair i, j. How many times do we expect to see j visit i if there is, or is not, a preferred interaction between i and j? The answer will depend on several factors. First, and most obviously, we expect the number of visits to be higher if j is in fact a pollinator of i. That is, we expect Mij to be larger if Bij = 1 than if Bij = 0. Second, we expect there to be more visits if there is greater sampling effort—for instance if the period of observation is longer or if the land area over which observations take place is larger15,16,26,27. Third, we expect to see more visits for more abundant plant and pollinator species than for less abundant ones, as demonstrated by several studies28,30. And fourth, as discussed above, we expect there to be some random variation in the number of visits, driven by fluctuations in individual behavior and the environment. These are the primary features that we incorporate into our model. It is possible to add others to handle specific situations (see ref. 39 and the Methods), but we focus on these four here.We translate these factors into a mathematical model of plant–pollinator interaction as follows. The random variations in the numbers of visits will follow a Poisson distribution for each plant–pollinator pair i, j, parameterized by a single number, the distribution mean μij, provided only that measurements are made sufficiently far apart to be independent (which under normal conditions they will be). We expect μij to depend on the factors discussed above and we introduce additional parameters to represent this dependence. First we introduce a parameter r to represent the change in the average number of visits when two species are connected (Bij = 1), versus when they are not (Bij = 0). We write the factor by which the number of visits is increased as 1 + r with r ≥ 0, so that r = 0 implies no increase and successively larger values of r give us larger increases. Second, we represent the effect of sampling effort by an overall constant C that multiplies the mean μij. The same constant is used for all i and j, since the same sampling effort is devoted to all plant–pollinator pairs. Third, we assume that the number of visits is proportional to the abundance of the relevant plant and pollinator species: twice as many pollinators of species j, for instance, will mean twice as many visits by that species, and similarly for the abundance of the plant species13. Thus the number of visits will be proportional to σiτj, for some parameters σi and τj representing the abundances of plant i and pollinator j, respectively, in suitable units (which we will determine shortly).Putting everything together, the mean number of observed visits to plant i by pollinator j is$${mu }_{ij}=C{sigma }_{i}{tau }_{j}(1+r{B}_{ij}),$$
(1)
and the probability of observing exactly Mij visits is drawn from a Poisson distribution with this mean:$$P({M}_{ij}| {mu }_{ij})=frac{{mu }_{ij}^{{M}_{ij}}}{{M}_{ij}!} {e}^{-{mu }_{ij}}.$$
(2)
This equation gives us the probability distribution of a single element Mij of the data matrix. Then, combining Eqs. (1) and (2), the data likelihood—the probability of the complete data matrix M—is given by the product over all species thus:$$P({boldsymbol{M}}| {boldsymbol{B}},theta )=mathop{prod}limits_{i,j}frac{{left[C{sigma }_{i}{tau }_{j}(1+r{B}_{ij})right]}^{{M}_{ij}}}{{M}_{ij}!} {e}^{-C{sigma }_{i}{tau }_{j}(1+r{B}_{ij})},$$
(3)
where θ is a shorthand collectively denoting all the parameters of the model: C, r, σ and τ. Our model is thus effectively a model of an entire network, rather than single interactions, in contrast with other recent approaches to the modeling of network data reliability17,18,32.There are two important details to note about this model. First, the definition in Eq. (1) does not completely determine C, σ, and τ because we can increase (or decrease) any of these parameters by a constant factor without changing the resulting value of μij if we simultaneously decrease (or increase) one or both of the others. In the language of statistics we say that the parameters are not “identifiable.” We can rectify this problem by fixing the normalization of the parameters in any convenient fashion. Here we do this by stipulating that σi and τj sum to one, thus:$$mathop{sum }limits_{i=1}^{{n}_{p}}{sigma }_{i}=mathop{sum }limits_{j=1}^{{n}_{a}}{tau }_{j}=1.$$
(4)
In effect, this makes σi and τj measures of relative abundance, quantifying the fraction of individual organisms that belong to each species, rather than the total number. (This definition differs from traditional estimates of pollinator abundance that define the abundance of a pollinator species in terms of its number of observed visits.) Second, there may be other species-level effects on the observed number of visits in addition to abundance, such as the propensity for observers to overlook small-bodied pollinators. There is, at least within the data used in this paper, no way to tell these effects from true variation in abundance—no way to tell for example if there are truly fewer individuals of a species or if they are just hard to see and hence less often observed. As a result, the abundance parameters in our model actually capture a combination of effects on observation frequency. This does not affect the accuracy of the model, which works just as well either way, but it does mean that we have to be cautious about interpreting the values of the parameters in terms of actual abundance. This point is discussed further in the applications below.Bayesian reconstructionThe likelihood of Eq. (3) tells us the probability of the data M given the network B and parameters θ. What we actually want to know is the probability of the network and parameters given the data, which we can calculate by applying Bayes’ rule in the form$$P({boldsymbol{B}},theta | {boldsymbol{M}})=frac{P({boldsymbol{M}}| {boldsymbol{B}},theta )P({boldsymbol{B}}| theta )P(theta )}{P({boldsymbol{M}})}.$$
(5)
This is the posterior probability that the network has structure B and parameter values θ given the observations that were made. There are three important parts to the expression: the likelihood P(M∣B, θ), the prior probability of the network P(B∣θ), and the prior probability of the parameters P(θ). The denominator P(M) we can ignore because it depends on the data alone and will be constant (and hence irrelevant for our calculations) once M is determined by the observations.Of the three non-constant parts, the first, the likelihood, we have already discussed—it is given by Eq. (3). For the prior on the network P(B∣θ) we make the conservative assumption—in the absence of any knowledge to the contrary—that all edges in the network are a priori equally likely. If we denote the probability of an edge by ρ, then the prior probability on the entire network is$$P({boldsymbol{B}}| theta )=mathop{prod}limits_{i,j}{(1-rho )}^{1-{B}_{ij}}{rho }^{{B}_{ij}}.$$
(6)
We consider ρ an additional parameter which is to be inferred from the data and which we will henceforth include, along with our other parameters, in the set θ.To complete Eq. (5), we also need to choose a prior P(θ) over the parameters. We expect there to be some limit on the value of r, which we impose using a minimally informative prior with finite mean (this distribution turns out to be the exponential distribution). For the remaining parameters we use uniform priors. With these choices, we then have everything we need to compute the posterior probability, Eq. (5).Once we have the posterior probability there are a number of things we can do with it. The simplest is just to maximize it with respect to the unknown quantities B and θ to find the most likely structure for the network and the most likely parameter values, given the data. This, however, misses an opportunity for more detailed inference and can moreover give misleading results. In most cases there will be more than one value of B and θ with high probability under Eq. (5): there may be a unique maximum of the probability, a most likely value, but there are often many other values that have nearly as high probability and offer plausible network structures competitive with the most likely one. To get the most complete picture of the structure of the network we should consider all these plausible structures.For example, if all plausible structures are similar to one another in their overall shape then we can be quite confident that this shape is reflective of the true preferred interactions between plant and pollinator species. If plausible structures are widely varying, however, then we have many different candidates for the true structure and our certainty about that structure is correspondingly lower. In other words, by considering the complete set of plausible structures we can not only make an estimate of the network structure but also say how confident we are in that estimate, in effect putting “error bars” on the network.How do we specify these errors bars in practice? One way is to place posterior probabilities on individual edges in the network. For example, when considering the edge connecting plant i and pollinator j, we would not ask “Is there an edge?” but rather “What is the probability that there is an edge?” Within the formulation outlined above, this probability is given by the average$$P({B}_{ij}=1| {boldsymbol{M}})=mathop{sum}limits_{{boldsymbol{B}}}int {B}_{ij}P({boldsymbol{B}},theta | {boldsymbol{M}})dtheta ,$$
(7)
where the sum runs over all possible incidence matrices and the integral over all parameter values. More generally we can compute the average of any function f(B, θ) of the matrix B and/or the parameters θ thus:$$leftlangle f({boldsymbol{B}},theta )rightrangle =mathop{sum}limits_{{boldsymbol{B}}}int f({boldsymbol{B}},theta ) P({boldsymbol{B}},theta | {boldsymbol{M}})dtheta .$$
(8)
Functions of the matrix and functions of the parameters can both be interesting—the matrix tells us about the structure of the network but the parameters, as we will see, can also reveal important information.Computing averages of the form (8) is unfortunately not an easy task. A closed-form expression appears out of reach and the brute-force approach of performing the sums and integrals numerically over all possible networks and parameters is computationally intractable in all but the most trivial of cases. The sum over B alone involves ({2}^{{n}_{p}{n}_{a}}) terms, which is normally a very large number.Instead therefore we use an efficient Monte Carlo sampling technique to approximate the answers. We generate a sample of network/parameter pairs (B1, θ1), …, (Bn, θn), where each pair appears with probability proportional to the posterior distribution of Eq. (5). Then we approximate the average of f(B, θ) as$$leftlangle f({boldsymbol{B}},theta )rightrangle simeq frac{1}{n}mathop{sum }limits_{i=1}^{n}f({{boldsymbol{B}}}_{i},{theta }_{i}).$$
(9)
Under very general conditions, this estimate will converge to the true value of the average asymptotically as the number of Monte Carlo samples n becomes large. Full details of the computations are given in Materials and Methods, and an extensive simulation study of the model is presented in Supplementary Note 1.Checking the modelInherent in the discussion so far is the assumption that the data can be well represented by our model. In other words, we are assuming there is at least one choice of the network B and parameters θ such that the model will generate data similar to what we see in the field. This assumption could be violated if our model is a poor one, but there is nothing in the method described above that would tell us so. To be fully confident in our results we need to be able not only to infer the network structure, but also to check whether that structure is a good match to the data. The Bayesian toolbox comes with a natural procedure for doing this. Given a set of high-probability values of B and θ generated by the method, we can use them in Eq. (3) to compute the likelihood P(M∣B, θ) of a data set M and then sample possible data sets from this probability distribution, in effect recreating data as they would appear if the model were in fact correct. We can then compare these data to the original field data to see if they are similar: if they are then our model has done a good job of capturing the structure in the data.In the parlance of Bayesian statistics this approach is known as a posterior–predictive assessment40. It amounts to calculating the probability$$P({widetilde{M}}_{ij}| {boldsymbol{M}})=mathop{sum}limits_{{boldsymbol{B}}}int P({widetilde{M}}_{ij}| {boldsymbol{B}},theta )P({boldsymbol{B}},theta | {boldsymbol{M}})dtheta$$
(10)
that pollinator species j makes ({widetilde{M}}_{ij}) visits to plant species i in artificial data sets generated by the model, averaged over many sets of values of B and θ. We can then use this probability to calculate the average value of ({widetilde{M}}_{ij}) thus:$$langle {widetilde{M}}_{ij}rangle =mathop{sum}limits_{{widetilde{M}}_{ij}}{widetilde{M}}_{ij} P({widetilde{M}}_{ij}| {boldsymbol{M}}).$$
(11)
The averages for all plant–pollinator pairs can be thought of as the elements of a matrix (langle widetilde{{boldsymbol{M}}}rangle), which we can then compare to the actual data matrix M, or alternatively we can calculate a residue ({boldsymbol{M}}-langle widetilde{{boldsymbol{M}}}rangle). If (langle widetilde{{boldsymbol{M}}}rangle) and M are approximately equal, or equivalently if the residue is small, then we consider the model a good one.To quantify the level of agreement between the fit and the data we can also compute the discrepancy40 between the artificial data and M as$${X}^{2}=mathop{sum}limits_{ij}frac{{({M}_{ij}-{langle widetilde{M}_{ij}rangle })}^{2}}{{langle widetilde{M}_{ij}rangle }}.$$
(12)
Under the hypothesis that the model is correct, X2 follows a chi-squared distribution with np × na degrees of freedom40. A good fit between model and data is signified by a value of X2 that is much smaller than its expectation value of np × na. Note that the calculation of (P({widetilde{M}}_{ij}| {boldsymbol{M}})) in Eq. (10) is of the same form as the one in Eq. (8), with (f({boldsymbol{B}},theta )=P({widetilde{M}}_{ij}| {boldsymbol{B}},theta )), which means we can calculate (P({widetilde{M}}_{ij}| {boldsymbol{M}})) in the same way we calculate other average quantities, using Monte Carlo sampling and Eq. (9).Application to visitation data setsWell-sampled dataTo demonstrate how the method works in practice, we first consider a large data set of plant–pollinator interactions gathered by Kaiser-Bunbury and collaborators41 at a set of study sites on the island of Mahé in the Seychelles. The data describe the interactions of plant and pollinator species observed over a period of eight months across eight different sites on the island. The data also include measurements of floral abundances for all observation periods and all sites. Our method for inferring network structure does not make use of the abundance measurements, but we discuss them briefly at the end of this section.The study by Kaiser-Bunbury et al. focused particularly on the role of exotic plant species in the ecosystem and on whether restoring a site by removing exotic species would significantly impact the resilience and function of the plant–pollinator network. To help address these questions, half of the sites in the study were restored in this way while the rest were left unrestored as a control group.As an illustration of our method we apply it to data from one of the restored sites, as observed over the course of a single month in December 2012 (the smallest time interval for which data were available). We pick the site named “Trois-Frères” because it is relatively small but also well sampled. Our calculation then proceeds as shown in Fig. 1. There were 8 plant and 21 pollinator species observed at the site during the month, giving us an 8 × 21 data matrix M as shown in Fig. 1a. (Following common convention, the plots of matrices in this paper are drawn with rows and columns ordered by decreasing numbers of observed interactions, so that the largest elements of the data matrix—the darkest squares—are in the top and left of the plot.)Fig. 1: Illustration of the method of this paper applied to data from the study of Kaiser-Bunbury et al.41.a We start with a data matrix M that records the number of interactions between each plant species and pollinator species. Species pairs that are never observed to interact (Mij = 0) are shown in white. b We then draw 2000 samples from the distribution of Eq. (5), four of which are shown in the figure. Each sample consists of a binary incidence matrix B, values for the relative abundances σ and τ (shown as the orange and blue bar plots, respectively), and values for the parameters C, r, and ρ (not shown). c We combine the samples using Eqs. (7)–(9) to give an estimate of the probability of each edge in the network and the complete parameter set θ. For the data set studied here our estimates of the expected values of the parameters C, r, and ρ are 〈C〉 = 20.2, 〈r〉 = 45.9, and 〈ρ〉 = 0.244.Full size imageNow we use our Monte Carlo procedure to draw 2000 sets of incidence matrices B and parameters θ from the posterior distribution of Eq. (5) (Fig. 1b). These samples vary in their structure: some edges, like the one connecting the plant N. vanhoutteanum and the pollinator A. mellifera, are present in nearly all samples, while others, like the one between M. sechellarum and A. mellifera, appear only a small fraction of the time. Some others never occur at all. Averaging over these sampled networks we can estimate the probability, Eq. (7), that each connection exists in the network of preferred interactions between plant and animal species—see Fig. 1c. Some connections have high probability, close to 1, meaning that we have a high degree of confidence that they exist. Others have probability close to 0, meaning we have a high degree of confidence that they do not exist. And some have intermediate probabilities, meaning we are uncertain about them (such as the M. sechellarum–A. mellifera connection, which has probability around 0.45). In the latter case the method is telling us that the data are not sufficient to reach a firm conclusion about these connections. Indeed, if we compare with the original data matrix M in Fig. 1a, we find that most of the uncertain connections are ones for which we have very few observations, relative to the total number of observations for these species—say Mij = 1 or 2 for species with dozens of total observations overall.As we have mentioned, we also need to check whether the model is a good fit to the data by performing a posterior–predictive test. Figure 2 shows the results of this test. The main plot in the figure compares the values of the 40 largest elements of the original data matrix M with the corresponding elements of the generated matrix (widetilde{{boldsymbol{M}}}). In each case, the original value is well within one standard deviation of the average value generated by the test, confirming the accuracy of the model. The inset of the figure shows the residue matrix ({boldsymbol{M}}-widetilde{{boldsymbol{M}}}), which reveals no systematic bias unaccounted for by the model. The discrepancy X2 of Eq. (12) takes the value 26.94 in this case, well below the expected value of npna = 168, which indicates that the good fit is not a statistical fluke.Fig. 2: Results of a posterior–predictive test on the data matrix M for the example data set analyzed in Fig. 1.The main plot shows the error on the 40 largest entries of M, while the inset shows the residue matrix ({boldsymbol{M}}-langle tilde{{boldsymbol{M}}}rangle). Because the actual data M are well within one standard deviation of the posterior–predictive mean, the test confirms that the model is a good fit in this case. Error bars correspond to one standard deviation and are computed with n = 2000 samples from the posterior distribution.Full size imageIn addition to inferring the structure of the network itself, our method allows us to estimate many other quantities from the data. There are two primary methods by which we can do this. One is to look at the values of the fitted model parameters, which represent quantities such as the preference r and species abundances σ, τ. The other is to compute averages of quantities that depend on the network structure or the parameters (or both) from Eq. (9).As an example of the former approach, consider the parameter ρ, which represents the average probability of an edge, also known as the connectance of the network. Figure 3a shows the distribution of values of this quantity over our set of Monte Carlo samples, and neatly summarizes our overall certainty about the presence or absence of edges. If we were certain about all edges in the network, then ρ would take only a single value and the distribution would be narrowly peaked. The distribution we observe, however, is somewhat broadened, indicating significant uncertainty. The most likely value of ρ, the peak of the distribution, turns out to be quite close to the value one would arrive at if one were simply to assume that every pair of species that interacts even once is connected in the network. This does not mean, however, that one could make this assumption and get good results. As we show below, the network one would derive by doing so would be badly in error in other ways.Fig. 3: Analyses that can be performed using samples from the posterior distribution of Eq. (5).a Distribution of the connectance ρ. Connectance values for binary networks obtained by thresholding the data matrix at Mij  > 0 and Mij ≥ 5 are shown as vertical lines for reference. b Distribution of the preference parameter r. The mean value of r is 〈r〉 = 45.9 and its mode close to 40, but individual values as high as 100 are possible. c Distribution of the nestedness measure NODF. Values obtained by thresholding the data matrix at Mij  > 0 and Mij  > 1 are shown for reference. d Measured and estimated abundances for each of the plant species (R2 = 0.54).Full size imageFigure 3b shows the distribution of another of the model parameters, the parameter r, which measures the extent to which pollinators prefer the plants they normally pollinate over the ones they do not. For this particular data set the most likely value of r is around 40, meaning that pollinators visit their preferred plant species about 40 times more often than non-preferred ones, indicating all other things being equal, an impressive level of selectivity on the part of the pollinators.For the calculation of more complicated network properties we can perform an average over the value of any function f(B, θ), as long as there is an algorithm to compute it. As an example, Fig. 3c shows a calculation of the quantity known as “Nestedness based on Overlap and Decreasing Fill” (NODF), a measure of the nestedness property discussed in the introduction. This quantity measures the extent to which specialist species—those with relatively few interactions—tend to interact with a subset of the partners of generalist species42. While it is complicated to compute NODF analytically, due to the fact that one must order the species by degrees22, it is straightforward to calculate it within our framework: we simply calculate the value for each sampled network B and plot the resulting distribution. Interestingly, the most likely value of NODF is significantly different from the one we would calculate had we assumed, as discussed above, that a single interaction is sufficient to consider two species connected. On the contrary, we find that the system is almost certainly more nested than this simple analysis would conclude.In Fig. 3d, we compare the values of our estimated floral abundance parameters σ to the measured abundances reported by Kaiser-Bunbury et al.41. These parameters are not measures of abundance in the usual sense, because they combine actual abundance (quantity or density) with other characteristics such as ease of observation. We do find a correlation between the estimated and observed abundances, but it is relatively weak (R2 = 0.54), signaling significant disagreement, on which we elaborate in the discussion section.Undersampled dataAs we have pointed out, the connections in the network about which we are most uncertain tend to be ones that are undersampled, i.e., those for which we have only a small amount of data. In an ideal world we could address this problem by taking more data, but it is rare that we have the opportunity to do this. More commonly the data have already been gathered and our task is to produce the best results we can with those data. There are nonetheless some remedies open to us, such as aggregating data over different geographical areas or time windows. In Fig. 4 we compare the edge probabilities estimated from data recorded individually at the four “restored” sites in the Mahé study during October 2012 to the edge probabilities we obtain when we aggregate these observations into a single data matrix and only then estimate the network. (We use restored sites observed during the same month because they are likely to be ecologically similar, meaning the data are measuring approximately the same system.) Comparison of the two distributions shows—as we would hope—that there are fewer uncertain edges in the aggregated network than in its disaggregated parts, i.e., there are fewer edges with probabilities in the middle of the distribution and more with probabilities close to zero or one.Fig. 4: Illustration of the effect of data aggregation on edge uncertainty.a Histogram of the edge probabilities P(Bij = 1∣M) for the four restored sites in the Mahé study as observed in October 2012 and analyzed individually. b Equivalent histogram after aggregating the data over the sites and then estimating a single network from the resulting data matrix. The horizontal lines, both drawn at fifty observations—are added merely as a guide to the eye. Note how the upper histogram has more mass near the middle of the plot, while the lower one has most of its mass close to probability zero or one, indicating greater certainty in the positions of the edges in the aggregated data.Full size imageIn other cases neither aggregation nor gathering more data is possible, for instance when reanalyzing a data set already collected by others or already maximally aggregated. Such data sets record the results of observational studies that are already over, and may contain too few observations, but our approach still allows us to perform rigorous inference in these circumstances.For instance, Jordano et al.43 used dozens of existing plant-pollinator and plant-frugivore data sets to argue that the degree distributions of mutualistic networks have a long tail, but this conclusion is undermined by issues with undersampling. As an example, one of the data sets they studied, originally gathered by Inouye and Pyke44, records 1314 individual interactions over a period of 3 months in Kosciusko National Park, Australia, between 40 plants and 85 pollinator species, which works out to an average of 0.386 unique observations per species pair. Is this sampling effort sufficient to establish edges with certainty? As a point of reference, the data analyzed in Fig. 1 comprises 201 observations between 8 plants and 21 pollinators species for an average of 1.196 observations per pair of species, and the aggregated data of Fig. 4 contain 1.420 observations for every pair. Nonetheless, there is uncertainty about some of the connections in these reconstructed networks; this suggests that the network of Inouye and Pyke, with less than a third as much data per species pair, will contain significant uncertainty.Even so, our method allows us to make inferences about this network. In Fig. 5, we show estimates of the degree distributions of both plant and pollinator nodes in the network obtained from the posterior distribution P(B∣M), along with naive estimates calculated by thresholding the (undersampled) data as in the study by Jordano et al.43. As the figure shows, the results derived from the two approaches are very different. The thresholded degree distributions were classified as scale-free by Jordano et al., but this classification no longer holds once we account for the issues with the data; the inferred degree distributions are in this case well-modeled as Poisson distributions of means 5.53 and 2.60 for plants and pollinators respectively and the power-law form is a poor fit. On the other hand, the abundance parameters of the model, shown in Fig. 5, do appear to have a broad distribution, an interesting finding that calls for a rethinking of the relationship between abundances and degree distributions. It is generally thought that interactions will tend to be evenly distributed under an even distribution of abundance13 but here the opposite seems to be true.Fig. 5: Distributions of species-level parameters for a network of plants and pollinators in Kosciusko National Park, Australia, from the study by Inouye and Pyke44.a Thresholded degree distributions calculated by connecting species i and j with an edge if Mij  > 0. Inferred degree distributions are calculated using the method of this paper, averaging the fraction pk of nodes with a given degree k over n = 2000 Monte Carlo samples. b Inferred distributions of abundances σ and τ, calculated as a histogram over n = 2000 Monte Carlo samples of the abundance parameters of the fitted model. Error bars correspond to one standard deviation in all cases.Full size image More

1.Pecl, G. T. et al. Biodiversity redistribution under climate change: impacts on ecosystems and human well-being. Science 355, eaai9214 (2017).PubMed 
Article 
CAS 

Google Scholar 
2.Diffenbaugh, N. S. & Field, C. B. Changes in ecologically critical terrestrial climate conditions. Science 341, 486–492 (2013).ADS 
CAS 
PubMed 
Article 

Google Scholar 
3.Viana, D. S., Santamaría, L. & Figuerola, J. Migratory birds as global dispersal vectors. Trends Ecol. Evol. 31, 763–775 (2016).PubMed 
Article 

Google Scholar 
4.Bauer, S. & Hoye, B. J. Migratory animals couple biodiversity and ecosystem functioning worldwide. Science 344, 1242552 (2014).CAS 
PubMed 
Article 

Google Scholar 
5.Corlett, R. T. & Westcott, D. A. Will plant movements keep up with climate change? Trends Ecol. Evol. 28, 482–488 (2013).PubMed 
Article 

Google Scholar 
6.Lenoir, J. & Svenning, J. C. Climate-related range shifts – a global multidimensional synthesis and new research directions. Ecography 38, 15–28 (2015).Article 

Google Scholar 
7.Lenoir, J. et al. Species better track climate warming in the oceans than on land. Nat. Ecol. Evol. 4, 1044–1059 (2020).PubMed 
Article 

Google Scholar 
8.Loarie, S. R. et al. The velocity of climate change. Nature 462, 1052–1055 (2009).ADS 
CAS 
PubMed 
Article 

Google Scholar 
9.Chen, I.-C., Hill, J. K., Ohlemüller, R., Roy, D. B. & Thomas, C. D. Rapid range shifts of species associated with high levels of climate warming. Science 333, 1024–1026 (2011).ADS 
CAS 
Article 

Google Scholar 
10.González-Varo, J. P., López-Bao, J. V. & Guitián, J. Seed dispersers help plants to escape global warming. Oikos 126, 1600–1606 (2017).Article 

Google Scholar 
11.Urban, M. C. et al. Improving the forecast for biodiversity under climate change. Science 353, aad8466 (2016).PubMed 
Article 
CAS 

Google Scholar 
12.Thuiller, W. et al. Predicting global change impacts on plant species’ distributions: future challenges. Perspect. Plant Ecol. Evol. Syst. 9, 137–152 (2008).Article 

Google Scholar 
13.Nadeau, C. P. & Urban, M. C. Eco-evolution on the edge during climate change. Ecography 42, 1280–1297 (2019).
Google Scholar 
14.Bacles, C. F. E., Lowe, A. J. & Ennos, R. A. Effective seed dispersal across a fragmented landscape. Science 311, 628 (2006).PubMed 
Article 

Google Scholar 
15.Jordano, P., García, C., Godoy, J. A. & García-Castaño, J. L. Differential contribution of frugivores to complex seed dispersal patterns. Proc. Natl Acad. Sci. USA 104, 3278–3282 (2007).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
16.Breitbach, N., Böhning-Gaese, K., Laube, I. & Schleuning, M. Short seed-dispersal distances and low seedling recruitment in farmland populations of bird-dispersed cherry trees. J. Ecol. 100, 1349–1358 (2012).Article 

Google Scholar 
17.Cain, M. L., Damman, H. & Muir, A. Seed dispersal and the Holocene migration of woodland herbs. Ecol. Monogr. 68, 325–347 (1998).Article 

Google Scholar 
18.Nathan, R. et al. Spread of North American wind-dispersed trees in future environments. Ecol. Lett. 14, 211–219 (2011).PubMed 
Article 

Google Scholar 
19.Nathan, R. et al. Mechanisms of long-distance seed dispersal. Trends Ecol. Evol. 23, 638–647 (2008).PubMed 
Article 

Google Scholar 
20.Viana, D. S., Gangoso, L., Bouten, W. & Figuerola, J. Overseas seed dispersal by migratory birds. Proc. R. Soc. Lond. B 283, 20152406 (2016).
Google Scholar 
21.Viana, D. S., Santamaría, L., Michot, T. C. & Figuerola, J. Migratory strategies of waterbirds shape the continental-scale dispersal of aquatic organisms. Ecography 36, 430–438 (2013).Article 

Google Scholar 
22.Carlquist, S. The biota of long-distance dispersal. V. Plant dispersal to Pacific islands. Bull. Torrey Bot. Club 94, 129–162 (1967).Article 

Google Scholar 
23.Esteves, C. F., Costa, J. M., Vargas, P., Freitas, H. & Heleno, R. H. On the limited potential of Azorean fleshy fruits for oceanic dispersal. PLoS ONE 10, e0138882 (2015).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
24.Viana, D. S., Santamaría, L., Michot, T. C. & Figuerola, J. Allometric scaling of long-distance seed dispersal by migratory birds. Am. Nat. 181, 649–662 (2013).PubMed 
Article 
PubMed Central 

Google Scholar 
25.Martínez-López, V., García, C., Zapata, V., Robledano, F. & De la Rúa, P. Intercontinental long-distance seed dispersal across the Mediterranean basin explains population genetic structure of a bird-dispersed shrub. Mol. Ecol. 29, 1408–1420 (2020).PubMed 
Article 

Google Scholar 
26.Newton, I. The Migration Ecology of Birds (Elsevier, 2010).27.Sorensen, A. E. Interactions between birds and fruit in a temperate woodland. Oecologia 50, 242–249 (1981).ADS 
CAS 
PubMed 
Article 

Google Scholar 
28.González-Varo, J. P., Arroyo, J. M. & Jordano, P. The timing of frugivore-mediated seed dispersal effectiveness. Mol. Ecol. 28, 219–231 (2019).PubMed 
Article 

Google Scholar 
29.Jordano, P. in Seeds: The Ecology of Regeneration of Plant Communities (ed. Gallagher, R. S.) 18–61 (CABI, 2014).30.Bascompte, J. & Jordano, P. Mutualistic Networks (Princeton Univ. Press, 2013).31.Gallinat, A. S. et al. Patterns and predictors of fleshy fruit phenology at five international botanical gardens. Am. J. Bot. 105, 1824–1834 (2018).PubMed 
Article 

Google Scholar 
32.Cadotte, M. W. Experimental evidence that evolutionarily diverse assemblages result in higher productivity. Proc. Natl Acad. Sci. USA 110, 8996–9000 (2013).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
33.Mitter, C., Farrell, B. & Futuyma, D. J. Phylogenetic studies of insect–plant interactions: insights into the genesis of diversity. Trends Ecol. Evol. 6, 290–293 (1991).CAS 
PubMed 
Article 

Google Scholar 
34.Siemann, E., Tilman, D., Haarstad, J. & Ritchie, M. Experimental tests of the dependence of arthropod diversity on plant diversity. Am. Nat. 152, 738–750 (1998).CAS 
PubMed 
Article 

Google Scholar 
35.Sanderson, F. J., Donald, P. F., Pain, D. J., Burfield, I. J. & van Bommel, F. P. J. Long-term population declines in Afro-Palearctic migrant birds. Biol. Conserv. 131, 93–105 (2006).Article 

Google Scholar 
36.Beresford, A. E. et al. Phenology and climate change in Africa and the decline of Afro-Palearctic migratory bird populations. Remote Sens. Ecol. Conserv. 5, 55–69 (2019).Article 

Google Scholar 
37.Nilsson, C., Bäckman, J. & Alerstam, T. Seasonal modulation of flight speed among nocturnal passerine migrants: differences between short- and long-distance migrants. Behav. Ecol. Sociobiol. 68, 1799–1807 (2014).Article 

Google Scholar 
38.Gaston, K. J. Valuing common species. Science 327, 154–155 (2010).ADS 
CAS 
PubMed 
Article 

Google Scholar 
39.Horton, K. G. et al. Phenology of nocturnal avian migration has shifted at the continental scale. Nat. Clim. Change 10, 63–68 (2020).ADS 
Article 

Google Scholar 
40.Miller-Rushing, A. J., Lloyd-Evans, T. L., Primack, R. B. & Satzinger, P. Bird migration times, climate change, and changing population sizes. Glob. Change Biol. 14, 1959–1972 (2008).ADS 
Article 

Google Scholar 
41.Brochet, A.-L. et al. Preliminary assessment of the scope and scale of illegal killing and taking of birds in the Mediterranean. Bird Conserv. Int. 26, 1–28 (2016).Article 

Google Scholar 
42.Kays, R., Crofoot, M. C., Jetz, W. & Wikelski, M. Terrestrial animal tracking as an eye on life and planet. Science 348, aaa2478 (2015).PubMed 
Article 
CAS 

Google Scholar 
43.Stiles, E. W. Patterns of fruit presentation and seed dispersal in bird-disseminated woody plants in the eastern deciduous forest. Am. Nat. 116, 670–688 (1980).Article 

Google Scholar 
44.Noma, N. & Yumoto, T. Fruiting phenology of animal-dispersed plants in response to winter migration of frugivores in a warm temperate forest on Yakushima Island, Japan. Ecol. Res. 12, 119–129 (1997).Article 

Google Scholar 
45.Lovas-Kiss, Á. et al. Shorebirds as important vectors for plant dispersal in Europe. Ecography 42, 956–967 (2019).Article 

Google Scholar 
46.Coughlan, N. E., Kelly, T. C., Davenport, J. & Jansen, M. A. K. Up, up and away: bird-mediated ectozoochorous dispersal between aquatic environments. Freshw. Biol. 62, 631–648 (2017).Article 

Google Scholar 
47.Olson, D. M. et al. Terrestrial ecoregions of the world: a new map of life on Earth: a new global map of terrestrial ecoregions provides an innovative tool for conserving biodiversity. Bioscience 51, 933–938 (2001).Article 

Google Scholar 
48.Rivas-Martínez, S., Penas, A. & Díaz, T. Bioclimatic Map of Europe, Thermoclimatic Belts (Cartographic Service, Univ. León, 2004).49.Olesen, J. M. et al. Missing and forbidden links in mutualistic networks. Proc. R. Soc. Lond. B 278, 725–732 (2011).
Google Scholar 
50.Snow, B. & Snow, D. Birds and Berries (T. and A. D. Poyser, 1988).51.Stiebel, H. & Bairlein, F. Frugivory in central European birds I: diet selection and foraging. Vogelwarte 46, 1–23 (2008).
Google Scholar 
52.González-Varo, J. P., Arroyo, J. M. & Jordano, P. Who dispersed the seeds? The use of DNA barcoding in frugivory and seed dispersal studies. Methods Ecol. Evol. 5, 806–814 (2014).Article 

Google Scholar 
53.Simmons, B. I. et al. Moving from frugivory to seed dispersal: incorporating the functional outcomes of interactions in plant–frugivore networks. J. Anim. Ecol. 87, 995–1007 (2018).PubMed 
PubMed Central 
Article 

Google Scholar 
54.Plein, M. et al. Constant properties of plant–frugivore networks despite fluctuations in fruit and bird communities in space and time. Ecology 94, 1296–1306 (2013).PubMed 
Article 

Google Scholar 
55.Albrecht, J. et al. Variation in neighbourhood context shapes frugivore-mediated facilitation and competition among co-dispersed plant species. J. Ecol. 103, 526–536 (2015).Article 

Google Scholar 
56.García, D. Birds in ecological networks: insights from bird–plant mutualistic interactions. Ardeola 63, 151–180 (2016).Article 

Google Scholar 
57.Farwig, N., Schabo, D. G. & Albrecht, J. Trait-associated loss of frugivores in fragmented forest does not affect seed removal rates. J. Ecol. 105, 20–28 (2017).Article 

Google Scholar 
58.Torroba Balmori, P., Zaldívar García, P. & Hernández Lázaro, Á. Semillas de Frutos Carnosos del Norte Ibérico: Guía de Identificación (Ediciones Univ. Valladolid, 2013).59.Stiebel, H. Frugivorie bei Mitteleuropäischen Vögeln. PhD thesis, Univ. Oldenburg (2003).60.Jordano, P. Data from: Angiosperm fleshy fruits and seed dispersers: a comparative analysis of adaptation and constraints in plant-animal interactions. Dryad https://doi.org/10.5061/dryad.9tb73 (2013).61.González-Varo, J. P., Carvalho, C. S., Arroyo, J. M. & Jordano, P. Unravelling seed dispersal through fragmented landscapes: frugivore species operate unevenly as mobile links. Mol. Ecol. 26, 4309–4321 (2017).PubMed 
Article 

Google Scholar 
62.Ratnasingham, S. & Hebert, P. D. N. bold: the Barcode of Life data system (http://www.barcodinglife.org). Mol. Ecol. Notes 7, 355–364 (2007).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
63.CBOL Plant Working Group et al. A DNA barcode for land plants. Proc. Natl Acad. Sci. USA 106, 12794–12797 (2009).PubMed Central 
Article 
PubMed 

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

Google Scholar 
65.González-Varo, J. P., Díaz-García, S., Arroyo, J. M. & Jordano, P. Seed dispersal by dispersing juvenile animals: a source of functional connectivity in fragmented landscapes. Biol. Lett. 15, 20190264 (2019).PubMed 
PubMed Central 
Article 

Google Scholar 
66.Fuentes, M. Latitudinal and elevational variation in fruiting phenology among western European bird-dispersed plants. Ecography 15, 177–183 (1992).Article 

Google Scholar 
67.Herrera, C. M. A study of avian frugivores, bird-dispersed plants, and their interaction in Mediterranean scrublands. Ecol. Monogr. 54, 1–23 (1984).Article 

Google Scholar 
68.Hampe, A. & Bairlein, F. Modified dispersal-related traits in disjunct populations of bird-dispersed Frangula alnus (Rhamnaceae): a result of its Quaternary distribution shifts? Ecography 23, 603–613 (2000).Article 

Google Scholar 
69.Thomas, P. A. & Mukassabi, T. A. Biological flora of the British Isles: Ruscus aculeatus. J. Ecol. 102, 1083–1100 (2014).Article 

Google Scholar 
70.Jordano, P. Biología de la reproducción de tres especies del género Lonicera (Caprifoliaceae) en la Sierra de Cazorla. An. Jardin Botanico Madr. 1979 48, 31–52 (1990).
Google Scholar 
71.Debussche, M. & Isenmann, P. A Mediterranean bird disperser assemblage: composition and phenology in relation to fruit availability. Rev. Ecol. 47, 411–432 (1992).
Google Scholar 
72.Jordano, P. Diet, fruit choice and variation in body condition of frugivorous warblers in Mediterranean scrubland. Ardea 76, 193–209 (1988).
Google Scholar 
73.Barroso, Á., Amor, F., Cerdá, X. & Boulay, R. Dispersal of non-myrmecochorous plants by a “keystone disperser” ant in a Mediterranean habitat reveals asymmetric interdependence. Insectes Soc. 60, 75–86 (2013).Article 

Google Scholar 
74.González-Varo, J. P. Fragmentation, habitat composition and the dispersal/predation balance in interactions between the Mediterranean myrtle and avian frugivores. Ecography 33, 185–197 (2010).Article 

Google Scholar 
75.Sánchez-Salcedo, E. M., Martínez-Nicolás, J. J. & Hernández, F. Phenological growth stages of mulberry tree (Morus sp.) codification and description according to the BBCH scale. Ann. Appl. Biol. 171, 441–450 (2017).Article 

Google Scholar 
76.García-Castaño, J. L. Consecuencias Demográficas de la Dispersión de Semillas por Aves y Mamíferos Frugívoros en la Vegetación Mediterránea de Montaña. PhD thesis, Univ. Sevilla (2001).77.Gilbert, O. L. Symphoricarpos albus (L.) S. F. Blake (S. rivularis Suksd., S. racemosus Michaux). J. Ecol. 83, 159–166 (1995).Article 

Google Scholar 
78.Billerman, S. M. et al. (eds) Birds of the World (Cornell Laboratory of Ornithology, 2020).79.Tellería, J., Asensio, B. & Díaz, M. Aves Ibéricas: II. Paseriformes (J. M. Reyero Editor, 1999).80.Díaz, M., Asensio, B. & Tellería, J. L. Aves Ibéricas: I. No paseriformes (J. M. Reyero Editor, 1996).81.SEO/Birdlife. La Enciclopedia de las Aves de España (SEO/Birdlife-Fundación BBVA, 2019).82.Spina, F. & Volponi, S. Atlante della Migrazione degli Uccelli in Italia. 2. Passeriformi (Ministero dell’Ambiente e della Tutela del Territorio e del Mare, Istituto Superiore per la Protezione e la Ricerca Ambientale (ISPRA), Tipografia SCR-Roma, 2008).83.Spina, F. & Volponi, S. Atlante della Migrazione degli Uccelli in Italia. 1. Non-Passeriformi (Ministero dell’Ambiente e della Tutela del Territorio e del Mare, Istituto Superiore per la Protezione e la Ricerca Ambientale (ISPRA), Tipografia CSR-Roma, 2008).84.Wernham, C. et al. The Migration Atlas: Movements of the Birds of Britain and Ireland (T. & A. D. Poyser, 2002).85.Cramp, S. The Complete Birds of the Western Paleartic (CD-ROM) (Oxford Univ. Press, 1998).86.Bairlein, F. et al. Atlas des Vogelzugs – Ringfunde deutscher Brut- und Gastvögel (Aula, 2014).87.Tomiałojć, L. & Stawarczyk, T. Awifauna Polski: Rozmieszczenie, Liczebność i Zmiany (PTPP pro. Natura, 2003).88.Busse, P., Gromadzki, M. & Szulc, B. Obserwacje przelotu jesiennego ptaków w roku 1960 w Górkach Wschodnich koło Gdańska (Observations on bird migration at Górki Wschodnie near Gdańsk Autumn 1960). Acta Ornithologica 7, 305–336 (1963).
Google Scholar 
89.Bobrek, R. et al. Międzysezonowa powtarzalność dynamiki jesiennej migracji wróblowych Passeriformes nad Jeziorem Rakutowskim. Ornis Polonica 57, 39–57 (2016).
Google Scholar 
90.Keller, M. et al. Ptaki Środkowej Wisły (M-ŚTO, 2017).91.Bocheński, M. et al. Awifauna przelotna i zimująca środkowego odcinka doliny Odry. Ptaki Śląska 16, 123–161 (2006).
Google Scholar 
92.BTO. BirdTrack. http://www.birdtrack.net (accessed October 2018).93.Brooks, M. E. et al. glmmTMB balances speed and flexibility among packages for zero-inflated generalized linear mixed modeling. R J. 9, 378–400 (2017).Article 

Google Scholar 
94.Fox, J. & Weisberg, S. An R Companion to Applied Regression 2nd edn (SAGE, 2011).95.Douma, J. C. & Weedon, J. T. Analysing continuous proportions in ecology and evolution: a practical introduction to beta and Dirichlet regression. Methods Ecol. Evol. 10, 1412–1430 (2019).Article 

Google Scholar 
96.Smith, S. A. & Brown, J. W. Constructing a broadly inclusive seed plant phylogeny. Am. J. Bot. 105, 302–314 (2018).PubMed 
Article 

Google Scholar 
97.Magallón, S., Gómez-Acevedo, S., Sánchez-Reyes, L. L. & Hernández-Hernández, T. A metacalibrated time-tree documents the early rise of flowering plant phylogenetic diversity. New Phytol. 207, 437–453 (2015).PubMed 
Article 
PubMed Central 

Google Scholar 
98.Paradis, E. & Schliep, K. ape 5.0: an environment for modern phylogenetics and evolutionary analyses in R. Bioinformatics 35, 526–528 (2019).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
99.Pagel, M. Inferring the historical patterns of biological evolution. Nature 401, 877–884 (1999).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
100.Freckleton, R. P., Harvey, P. H. & Pagel, M. Phylogenetic analysis and comparative data: a test and review of evidence. Am. Nat. 160, 712–726 (2002).CAS 
PubMed 
Article 

Google Scholar 
101.Molina-Venegas, R. & Rodríguez, M. Á. Revisiting phylogenetic signal; strong or negligible impacts of polytomies and branch length information? BMC Evol. Biol. 17, 53 (2017).PubMed 
PubMed Central 
Article 

Google Scholar 
102.Revell, L. J. phytools: an R package for phylogenetic comparative biology (and other things). Methods Ecol. Evol. 3, 217–223 (2012).Article 

Google Scholar 
103.Dormann, C. F., Fründ, J., Blüthgen, N. & Gruber, B. Indices, graphs and null models: analyzing bipartite ecological networks. Open Ecol. J. 2, 7–24 (2009).Article 

Google Scholar 
104.Bascompte, J., Jordano, P. & Olesen, J. M. Asymmetric coevolutionary networks facilitate biodiversity maintenance. Science 312, 431–433 (2006).ADS 
CAS 
PubMed 
Article 

Google Scholar 
105.Bates, D., Maechler, M. & Bolker, B. lme4: linear mixed-effects models using ‘Eigen’ and S4. R package version 1.1-19 https://CRAN.R-project.org/package=lme4 (2013). More

Study areaThe study was carried out in Bunnefjorden, Norway. The fjord froze over from January to April and we here analyze data from ice-free conditions in early winter (12 Nov to 2 Dec 2009). Bunnefjorden is a 150 m deep inner branch of the Oslofjord, and a 57 m deep sill at the entrance restricts water exchange with the outer part of the fjord. Klevjer and Kaartvedt24 provide a map of the study area. The fjord branch normally becomes hypoxic in the lower part of the water column. During the current study, oxygen contents were 2–3 ml l−1 between 15 and 60 m, while waters below 70–80 m were severely hypoxic and devoid of fishes9.Studies of overwintering sprat have been undertaken in Bunnefjorden during several winters, and the biology of sprat as well as the identity of the main acoustic targets in the fjord are well established9,18,24. In the winter of the current study, catches from 33 trawl samples were dominated by sprat; with ~ 40 times higher catches than the next most abundant species, herring (Clupea harengus)9.Study designSolberg and Kaartvedt9 and Solberg et al.18 provide details on methods, and we here only give a summary of the acoustic setup. In short, upward-looking Simrad EK 60 echosounders kept in pressure-proof casings were deployed at the bottom (150 m) and in buoys (80 and 30 m) for enhanced resolution in shallower part of the water column. Cables for electricity and transfer of data to a PC on shore enabled continuous operation of the systems. We here use the data from the shallowest echosounder (200 kHz) that provided superior resolution in near-surface water, though did not cover the full depth range of the population distribution. Echograms from the deeper located echosounders covering the whole (inhabited) water column and showing the full diel population behavior are given in Solberg and Kaartvedt9 and Solberg et al.18.Records of gas releaseReleased gas appeared as ascending lines in the echogram (Fig. 4). We quantified the release as explained by Solberg and Kaartvedt9. We only included ascending traces connected to the acoustic record of a fish, but without enumerating the release per individual fish. Since the same fish may release several bursts of bubbles within a short time interval, we here pooled any sequences of gas release within a 10-s period as one event. This procedure will also exclude cases with several different individuals releasing bubbles in the course of this short time interval, yet we chose this conservative approach not to generate an artificial high connection of gas releases between the fishes.Figure 4Echogram showing sprat releasing gas, with every oblique line representing one release event and lines with a different angle to the release events representing swimming sprat. Colors represent the volume-backscatter coefficient (Sv).Full size imageAnalyses of dataThe frequency of gas releases varied with time, both within a day and between the weeks. Such patterns compare to service systems like call centres and hospital emergency rooms25 that can be modelled as a Poisson process26,27. We therefore started our analysis with the statistical procedure suggested by Brown et al.28 in their influential analysis of the call dynamics in a banking call centre. The first step is to subdivide the day into time intervals, which are short enough to consider event rates as approximately constant. Here we chose to investigate alternative periods of respectively 1, 5, and 30 min, as well as 1, 2, 4, and 6 h. At the longest interval, the peaks in the gas release intensity are expected to be the result of a non-constant Poisson parameter, and therefore more likely to induce rejection of the null-hypothesis of a homogenous random process. In contrast, we expect to find higher concordance with a random process for the short intervals of 5 min. In assessing connectivity among gas bubble releases, we formulate a new model allowing for a formal test of non-randomness (summarized in Fig. 5). We name this approach the simulated connectivity test (S-CON test), which we implemented in R29, with the code being available in the Supplementary appendix.Figure 5Illustration of the steps related to the simulated connectivity test (S-CON test). Bubbles occurring within 10 s are pooled into single release events. We then determine the connectivity of each release event—aka the number of release events within the following 30 s time window. From these, we calculate the average connectivity for a specified period (1, 5, 30 min, 1-, 2-, 4- and 6-h intervals). In a final step, we compare the observed average connectivity to the critical value which is defined as the 95th percentile of 1000 simulations of random placements of the same number of releases. If the observed connectivity is larger than the critical value, we reject the null hypothesis of random gas releases for the specified time interval.Full size imageIf there is a common physiological reason or some form of communication among sprat, a burst of gas release is likely followed by subsequent releases. Thus, it is reasonable to assume that the total number of releases within a short time interval like 30 s would be effective in detecting dependencies between the releases. We therefore define the concept of connectivity as follows:Let the gas be released at times T1, T2, …, Tn and define the connectivity at each event as the numbers of records within the following 30 s. The average connectivity in any considered window of the investigated time-period (for example a window of 1 h) is defined as the average connectivity of all cases of connectivity within the considered window (see also Fig. 5).In order to test the null hypothesis of no dependency between gas bubble releases, we compare the measured average connectivity in the data set with the simulation of 1000 random placements of the total number of observations in a given time window. For example, if we consider a window of 30 min with 15 release events having an average connectivity of 2.1, we performed 1000 random placements of 15 points between 1 and 1800. In this way, we get 1000 simulated values of the average connectivity, from which we pick out the critical 95th percentile, following the common significance level of 0.05 in biology. If the observed average connectivity is larger than this critical value, we reject the null-hypothesis and conclude that the releases of gas bubbles are dependent random variables. Thus, if our example obtains a critical value of 1.7, the null-hypothesis of random arrival times of bubbles is rejected (because the observed value of 2.1 is larger than the critical value of 1.7).Since a dependency between the fish will induce a higher concentration of release events than produced by random releases, we expect the average connectivity to be quite sensitive to the alternative hypothesis of dependent arrival times. Also, note that the concept of connectivity has a combinatory nature, so we need only require that the considered window contains at least two releases of gas bubbles. In contrast, alternative approaches using Kolmogorov–Smirnov tests28 are based on the cumulative distribution function and therefore require at least five observed bubble releases.To test the dependency of the results on the chosen time interval, we also ran the analysis using connectivity intervals of 25 and 35 s, which revealed some variability to the estimates of non-random bubble release (Fig. 3) but did not influence the general pattern. We also tested whether the interval within which we consider subsequent bubbles to be part of one single release event influences our results. The more we consider sequential bubbles to be independent of each other, i.e. their own release event, the higher the proportion of non-random gas release and vice versa.Fish abundanceTo exclude the possibility that apparent connectivity would be a mere result of fluctuating fish abundance, we tested whether the number of released bubbles is a function of fish biomass. For this, we first calculated the total number of gas release events within 30-min periods. We then compared these values to the summed surface integrated acoustic scattering coefficient (SA) for the same periods and for the same depth interval (upper 30 m), assuming that the integrated scattering coefficient (SA) serves as a proxy for the total fish biomass9. We filtered the scattering data to remove noise from non-biological sources prior to use. Both variables were log-transformed prior to analysis. We then fitted a linear model of the two variables using generalized least squares. To account for temporal autocorrelation in the data, we also included a correlation structure of type 1 (corAR1). The analysis was done in R29 using the nlme package30.Ethics declarationsLive animals (fish) were not used in this study. More