Ecology

1.Baer, S. G. & Birgé, H. E. Soil ecosystem services: An overview. Manag. Soil Health Sustain. Agric. 1, 1–22 (2018).
Google Scholar 
2.Geisen, S., Wall, D. H. & van der Putten, W. H. Challenges and opportunities for soil biodiversity in the anthropocene. Curr. Biol. 29, R1036–R1044 (2019).CAS 
PubMed 

Google Scholar 
3.Guerra, C. A. et al. Tracking, targeting, and conserving soil biodiversity. Science 371, 239–241 (2021).ADS 
CAS 
PubMed 

Google Scholar 
4.Tsiafouli, M. A. et al. Intensive agriculture reduces soil biodiversity across Europe. Global Change Biol. 21, 973–985 (2015).ADS 

Google Scholar 
5.Bender, S. F., Wagg, C. & van der Heijden, M. G. A. An underground revolution: Biodiversity and soil ecological engineering for agricultural sustainability. Trends Ecol. Evol. 31, 440–452 (2016).PubMed 

Google Scholar 
6.Wagg, C., Bender, S. F., Widmer, F. & van der Heijden, M. G. A. Soil biodiversity and soil community composition determine ecosystem multifunctionality. PNAS 111, 5266–5270 (2014).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
7.Wall, D. H., Nielsen, U. N. & Six, J. Soil biodiversity and human health. Nature 528, 69–76 (2015).ADS 
CAS 
PubMed 

Google Scholar 
8.Smith, P. et al. Global change pressures on soils from land use and management. Global Change Biol. 22, 1008–1028 (2016).ADS 

Google Scholar 
9.Birkhofer, K., Smith, H. G. & Rundlöf, M. Environmental Impacts of Organic Farming. in eLS. 1–7 (John Wiley & Sons Ltd, 2016).10.Bengtsson, J., Ahnström, J. & Weibull, A.-C. The effects of organic agriculture on biodiversity and abundance: A meta-analysis: Organic agriculture, biodiversity and abundance. J. Appl. Ecol. 42, 261–269 (2005).
Google Scholar 
11.Abbott, L. K. & Manning, D. A. C. Soil health and related ecosystem services in organic agriculture. Sustain. Agric. Res. 4, 116 (2015).
Google Scholar 
12.de Graaff, M.-A., Hornslein, N., Throop, H. L., Kardol, P. & van Diepen, L. T. A. Effects of agricultural intensification on soil biodiversity and implications for ecosystem functioning: A meta-analysis. in Advances in Agronomy vol. 155 1–44 (Elsevier, 2019).13.Peters, M. K. et al. Climate–land-use interactions shape tropical mountain biodiversity and ecosystem functions. Nature 568, 88–92 (2019).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
14.Pokhrel, Y. et al. Global terrestrial water storage and drought severity under climate change. Nat. Clim. Change 11, 226–233 (2021).ADS 

Google Scholar 
15.Iglesias, A. & Garrote, L. Adaptation strategies for agricultural water management under climate change in Europe. Agric. Water Manage. 155, 113–124 (2015).
Google Scholar 
16.Pörtner, H. O. et al. IPBES-IPCC Co-sponsored Workshop Report Synopsis on Biodiversity and Climate Change. https://zenodo.org/record/4920414 (2021).17.Blankinship, J. C., Niklaus, P. A. & Hungate, B. A. A meta-analysis of responses of soil biota to global change. Oecologia 165, 553–565 (2011).ADS 
PubMed 

Google Scholar 
18.Holmstrup, M. et al. Long-term and realistic global change manipulations had low impact on diversity of soil biota in temperate heathland. Sci. Rep. 7, 41388 (2017).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
19.Fry, E. L. et al. Soil multifunctionality and drought resistance are determined by plant structural traits in restoring grassland. Ecology 99, 2260–2271 (2018).PubMed 

Google Scholar 
20.Zhou, Z., Wang, C. & Luo, Y. Meta-analysis of the impacts of global change factors on soil microbial diversity and functionality. Nat. Commun. 11, 3072 (2020).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
21.Schimel, J. P. Life in dry soils: Effects of drought on soil microbial communities and processes. Annu. Rev. Ecol. Evol. Syst. 49, 409–432 (2018).
Google Scholar 
22.Kundel, D. et al. Drought effects on nitrogen provisioning in different agricultural systems: Insights gained and lessons learned from a field experiment. Nitrogen 2, 1–17 (2021).
Google Scholar 
23.Abbasi, A. O. et al. Reviews and syntheses: Soil responses to manipulated precipitation changes: An assessment of meta-analyses. Biogeosciences 17, 3859–3873 (2020).ADS 
CAS 

Google Scholar 
24.Webber, H. et al. Diverging importance of drought stress for maize and winter wheat in Europe. Nat. Commun. 9, 4249 (2018).ADS 
PubMed 
PubMed Central 

Google Scholar 
25.Gomez-Zavaglia, A., Mejuto, J. C. & Simal-Gandara, J. Mitigation of emerging implications of climate change on food production systems. Food Res. Int. 134, 109256 (2020).CAS 
PubMed 
PubMed Central 

Google Scholar 
26.Yin, R. et al. Soil functional biodiversity and biological quality under threat: Intensive land use outweighs climate change. Soil Biol. Biochem. 147, 107847 (2020).CAS 
PubMed 
PubMed Central 

Google Scholar 
27.Rawls, W. J., Pachepsky, Y. A., Ritchie, J. C., Sobecki, T. M. & Bloodworth, H. Effect of soil organic carbon on soil water retention. Geoderma 116, 61–76 (2003).ADS 
CAS 

Google Scholar 
28.Lal, R. Soil health and carbon management. Food Energy Secur. 5, 212–222 (2016).
Google Scholar 
29.Iizumi, T. & Wagai, R. Leveraging drought risk reduction for sustainable food, soil and climate via soil organic carbon sequestration. Sci. Rep. 9, 19744 (2019).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
30.Fließbach, A., Oberholzer, H.-R., Gunst, L. & Mäder, P. Soil organic matter and biological soil quality indicators after 21 years of organic and conventional farming. Agric. Ecosyst. Environ. 118, 273–284 (2007).
Google Scholar 
31.Gattinger, A. et al. Enhanced top soil carbon stocks under organic farming. PNAS 109, 18226–18231 (2012).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
32.Schädler, M. et al. Investigating the consequences of climate change under different land-use regimes: A novel experimental infrastructure. Ecosphere 10, e02635 (2019).
Google Scholar 
33.Birkhofer, K. et al. Ecosystem services: Current challenges and opportunities for ecological research. Front. Ecol. Evol. 2, 87 (2015).
Google Scholar 
34.Birkhofer, K. et al. Relationships between multiple biodiversity components and ecosystem services along a landscape complexity gradient. Biol. Cons. 218, 247–253 (2018).
Google Scholar 
35.Chabert, A. & Sarthou, J.-P. Conservation agriculture as a promising trade-off between conventional and organic agriculture in bundling ecosystem services. Agric. Ecosyst. Environ. 292, 106815 (2020).CAS 

Google Scholar 
36.Felipe-Lucia, M. R. et al. Land-use intensity alters networks between biodiversity, ecosystem functions, and services. PNAS 117, 28140–28149 (2020).CAS 
PubMed 
PubMed Central 

Google Scholar 
37.Lori, M., Symnaczik, S., Mäder, P., De Deyn, G. & Gattinger, A. Organic farming enhances soil microbial abundance and activity: A meta-analysis and meta-regression. PLoS ONE 12, e0180442 (2017).PubMed 
PubMed Central 

Google Scholar 
38.Kundel, D. et al. Effects of simulated drought on biological soil quality, microbial diversity and yields under long-term conventional and organic agriculture. FEMS Microbiol. Ecol. 96, fiaa205 (2020).CAS 
PubMed 
PubMed Central 

Google Scholar 
39.Chen, Q.-L. et al. Rare microbial taxa as the major drivers of ecosystem multifunctionality in long-term fertilized soils. Soil Biol. Biochem. 141, 107686 (2020).CAS 

Google Scholar 
40.Garland, G. et al. Crop cover is more important than rotational diversity for soil multifunctionality and cereal yields in European cropping systems. Nat. Food 2, 28–37 (2021).
Google Scholar 
41.Tamburini, G. et al. Agricultural diversification promotes multiple ecosystem services without compromising yield. Sci. Adv. 6, eaba1715 (2020).ADS 
PubMed 
PubMed Central 

Google Scholar 
42.Vazquez, C., de Goede, R. G. M., Rutgers, M., de Koeijer, T. J. & Creamer, R. E. Assessing multifunctionality of agricultural soils: Reducing the biodiversity trade-off. Eur. J. Soil. Sci. 72, 1624–1639 (2020).
Google Scholar 
43.Zwetsloot, M. J. et al. Soil multifunctionality: Synergies and trade-offs across European climatic zones and land uses. Eur. J. Soil. Sci. 72, 1640–1654 (2020).
Google Scholar 
44.Delgado-Baquerizo, M. et al. Soil microbial communities drive the resistance of ecosystem multifunctionality to global change in drylands across the globe. Ecol. Lett. 20, 1295–1305 (2017).PubMed 

Google Scholar 
45.Bardgett, R. D. & Caruso, T. Soil microbial community responses to climate extremes: Resistance, resilience and transitions to alternative states. Phil. Trans. R. Soc. B 375, 20190112 (2020).CAS 
PubMed 
PubMed Central 

Google Scholar 
46.Meyer, S., Kundel, D., Birkhofer, K., Fliessbach, A. & Scheu, S. Soil microarthropods respond differently to simulated drought in organic and conventional farming systems. Ecol. Evol. 11, 10369–10380 (2021).PubMed 
PubMed Central 

Google Scholar 
47.De Smedt, P. et al. Linking macrodetritivore distribution to desiccation resistance in small forest fragments embedded in agricultural landscapes in Europe. Landscape Ecol. 33, 407–421 (2018).
Google Scholar 
48.Liu, W. P. A., Phillips, L. M., Terblanche, J. S., Janion-Scheepers, C. & Chown, S. L. An unusually diverse genus of Collembola in the Cape Floristic Region characterised by substantial desiccation tolerance. Oecologia 195, 873–885 (2021).ADS 
PubMed 

Google Scholar 
49.Birkhofer, K. et al. Long-term organic farming fosters below and aboveground biota: Implications for soil quality, biological control and productivity. Soil Biol. Biochem. 40, 2297–2308 (2008).CAS 

Google Scholar 
50.Mäder, P. Soil fertility and biodiversity in organic farming. Science 296, 1694–1697 (2002).ADS 
PubMed 

Google Scholar 
51.Birkhofer, K., Bezemer, T. M., Hedlund, K. & Setälä, H. Community composition of soil organisms under different wheat farming systems. in Microbial Ecology in Sustainable Agroecosystems 89–111 (CRC press Boca Raton, 2012).52.Birkhofer, K. et al. Soil fauna feeding activity in temperate grassland soils increases with legume and grass species richness. Soil Biol. Biochem. 43, 2200–2207 (2011).CAS 

Google Scholar 
53.Siebert, J. et al. Extensive grassland-use sustains high levels of soil biological activity, but does not alleviate detrimental climate change effects. Adv. Ecol. Res. 60, 25–58 (2019).
Google Scholar 
54.de Vries, F. T. et al. Land use alters the resistance and resilience of soil food webs to drought. Nat. Clim. Change 2, 276–280 (2012).ADS 

Google Scholar 
55.Torode, M. D. et al. Altered precipitation impacts on above-and below-ground grassland invertebrates: Summer drought leads to outbreaks in spring. Front. Plant Sci. 7, 1468 (2016).PubMed 
PubMed Central 

Google Scholar 
56.Jonas, J. L., Wilson, G. W. T., White, P. M. & Joern, A. Consumption of mycorrhizal and saprophytic fungi by Collembola in grassland soils. Soil Biol. Biochem. 39, 2594–2602 (2007).CAS 

Google Scholar 
57.Susanti, W. I., Pollierer, M. M., Widyastuti, R., Scheu, S. & Potapov, A. Conversion of rainforest to oil palm and rubber plantations alters energy channels in soil food webs. Ecol. Evol. 9, 9027–9039 (2019).PubMed 
PubMed Central 

Google Scholar 
58.Seres, A. et al. Collembola decrease the nitrogen uptake of maize through arbuscular mycorrhiza. ekol 28, 242–247 (2009).
Google Scholar 
59.Bender, S. F. & van der Heijden, M. G. A. Soil biota enhance agricultural sustainability by improving crop yield, nutrient uptake and reducing nitrogen leaching losses. J. Appl. Ecol. 52, 228–239 (2015).CAS 

Google Scholar 
60.Carson, J. K. et al. Low pore connectivity increases bacterial diversity in soil. Appl. Environ. Microbiol. 76, 3936–3942 (2010).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
61.Krause, H.-M., Fliessbach, A., Mayer, J. & Mäder, P. Implementation and management of the DOK long-term system comparison trial. in Long-Term Farming Systems Research 37–51, (Elsevier, 2020).62.Richner, W. et al. Grundlagen für die Düngung landwirtschaftlicher Kulturen in der Schweiz (GRUD 2017). Agrarforschung Schweiz 8, 47–66 (2017).
Google Scholar 
63.Kundel, D. et al. Design and manual to construct rainout-shelters for climate change experiments in agroecosystems. Front. Environ. Sci. 6, 14 (2018).
Google Scholar 
64.Garland, G. et al. A closer look at the functions behind ecosystem multifunctionality: A review. J. Ecol. 109, 600–613 (2021).
Google Scholar 
65.Anderson, M. J. Permutational Multivariate Analysis of Variance (PERMANOVA). in Wiley StatsRef: Statistics Reference Online 1–15.66.Fletcher, D. J. & Underwood, A. J. How to cope with negative estimates of components of variance in ecological field studies. J. Exp. Mar. Biol. Ecol. 273, 89–95 (2002).
Google Scholar 
67.Ho, J., Tumkaya, T., Aryal, S., Choi, H. & Claridge-Chang, A. Moving beyond P values: data analysis with estimation graphics. Nat. Methods 16, 565–566 (2019).CAS 
PubMed 
PubMed Central 

Google Scholar 
68.R Core Team. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing. Vienna, Austria. https://www.R-project.org.69.Revell, L. J. phytools: An R package for phylogenetic comparative biology (and other things). Methods Ecol. Evol. 3, 217–223 (2012).
Google Scholar  More

1.Swenson, N. G. et al. The biogeography and filtering of woody plant functional diversity in North and South America. Glob. Ecol. Biogeogr. 21, 798–808. https://doi.org/10.1111/j.1466-8238.2011.00727.x (2012).Article 

Google Scholar 
2.Wallace, A. R. Tropical Nature and Other Essays (Macmillan, 1878).
Google Scholar 
3.Connell, J. H. & Slatyer, R. O. Mechanisms of succession in natural communities and their role in community stability and organization. Am. Nat. 111(982), 1119–1144 (1977).Article 

Google Scholar 
4.Michalet, R. & Pugnaire, F. I. Facilitation in Communities: Underlying Mechanisms, Community and Ecosystem Implications (Wiley Online Library, 2016).
Google Scholar 
5.Connell, J. H. Diversity in tropical rain forests and coral reefs. Sci. Am. Nat. 199, 1302–1310 (1978).CAS 

Google Scholar 
6.Grime, J. P. Competitive exclusion in herbaceous vegetation. Nature 242(5396), 344–347 (1973).ADS 
Article 

Google Scholar 
7.Wilkinson, D. M. The disturbing history of intermediate disturbance. Oikos 84, 145–147 (1999).Article 

Google Scholar 
8.Al-Mufti, M. M., Sydes, C. L., Furness, S. B., Grime, J. P. & Band, S. R. A quantitative analysis of shoot phenology and dominance in herbaceous vegetation. J. Ecol. 65, 759–791 (1977).Article 

Google Scholar 
9.Huston, M. A. Disturbance, productivity, and species diversity: Empiricism vs. logic in ecological theory. Ecology 95(9), 2382–2396 (2014).Article 

Google Scholar 
10.Silvcrtown, J. Plant coexistence and the niche. Trends Ecol. Evol. 19, 605–611. https://doi.org/10.1016/j.tree.2004.09.003 (2004).Article 

Google Scholar 
11.Zobel, M. The relative of species pools in determining plant species richness: An alternative explanation of species coexistence?. Trends Ecol. Evol. 12(7), 266–269. https://doi.org/10.1016/S0169-5347(97)01096-3 (1997).CAS 
Article 
PubMed 

Google Scholar 
12.Cornwell, W. K., Schwilk, D. W. & Ackerly, D. D. A trait-based test for habitat filtering: Convex hull volume. Ecology 87(6), 1465–1471 (2006).Article 

Google Scholar 
13.Kraft, N. J., Valencia, R. & Ackerly, D. D. Functional traits and niche-based tree community assembly in an Amazonian forest. Science 322(5901), 580–582. https://doi.org/10.1126/science.1160662 (2008).ADS 
CAS 
Article 
PubMed 

Google Scholar 
14.Swenson, G. & Enquist, B. J. Opposing assembly mechanisms in a neotropical dry forest: Implications for phylogenetic and functional community ecology. Ecology 90(8), 2161–2170 (2009).Article 

Google Scholar 
15.Cavender-Bares, J., Kozak, K. H., Fine, P. V. & Kembel, S. W. The merging of community ecology and phylogenetic biology. Ecol. Lett. 12(7), 693–715 (2009).Article 

Google Scholar 
16.Cadotte, M. W. & Tucker, C. M. Should environmental filtering be abandoned?. Trends Ecol. Evol. 32(6), 429–437. https://doi.org/10.1016/j.tree.2017.03.004 (2017).Article 
PubMed 

Google Scholar 
17.Cariton, J. T. & Geller, J. B. Ecological roulette: The global transport of nonindigenous marine organisms. Science 261(5117), 78–82. https://doi.org/10.1126/science.261.5117.78 (1993).ADS 
CAS 
Article 
PubMed 

Google Scholar 
18.Connolly, L. M. T. S. R. Understanding diversity–stability relationships: Towards a unified model of portfolio effects. Ecol. Lett. 16, 140–150. https://doi.org/10.1111/ele.12019 (2013).Article 
PubMed 

Google Scholar 
19.Mccann, K. S. The diversity–stability debate. Nature 405, 228–233 (2000).CAS 
Article 

Google Scholar 
20.Baliddawa, C. W. Plant species diversity and crop pest control. An analytical review. Int. J. Trop. Insect Sci. 6(04), 479–487. https://doi.org/10.1017/s1742758400004306 (1985).Article 

Google Scholar 
21.Clara Nicholls, M. A. Plant biodiversity enhances bees and other insect pollinators in agroecosystems. A review. Agron. Sustain. Dev 33, 257–274. https://doi.org/10.1007/s13593-012-0092-y (2015).Article 

Google Scholar 
22.Haddad, N. M., Crutsinger, G. M., Gross, K., Haarstad, J. & Tilman, D. Plant diversity and the stability of foodwebs. Ecol. Lett. 14, 42–46. https://doi.org/10.1111/j.1461-0248.2010.01548.x (2011).Article 
PubMed 

Google Scholar 
23.Guyot, V., Castagneyrol, B., Vialatte, A. & Deconchat, M. Tree diversity reduces pest damage in mature forests across Europe. Biol. Lett. 12, 1–5 (2016).Article 

Google Scholar 
24.Fu, H. et al. Local and regional drivers of turnover and nestedness components of species and functional beta diversity in lake macrophyte communities in China. Sci. Total Environ. 687, 206–217. https://doi.org/10.1016/j.scitotenv.2019.06.092 (2019).ADS 
CAS 
Article 
PubMed 

Google Scholar 
25.MacDougall, A. S. et al. The Neolithic Plant Invasion Hypothesis: The role of preadaptation and disturbance in grassland invasion. New Phytol. 220(1), 94–103. https://doi.org/10.1111/nph.15285 (2018).Article 
PubMed 

Google Scholar 
26.Mouquet, N., Munguia, P., Kneitel, J. M. & Miller, T. E. Community assembly time and the relationship between local and regional species richness. Oikos 103(3), 618–626. https://doi.org/10.1034/j.1600-0706.2003.12772.x (2003).Article 

Google Scholar 
27.Sferra, C. O., Hart, J. L. & Howeth, J. G. Habitat age influences metacommunity assembly and species richness in successional pond ecosystems. Ecosphere 8(6), e01871 (2017).Article 

Google Scholar 
28.Macarthur, R. H. & Wilson, E. O. An equilibrium theory of insular zoogeography. Int. J. Org. Evol. 17(4), 373–387 (1963).Article 

Google Scholar 
29.Simberloff, D. S. Equilibrium theory of island biogeography and ecology. Annu. Rev. Ecol. Syst. 5(1), 161–182. https://doi.org/10.1146/annurev.es.05.110174.001113 (1974).Article 

Google Scholar 
30.Vijay, V., & Armsworth, P. R. Pervasive cropland in protected areas highlight trade-offs between conservation and food security. PNAS 118(4), e2010121118 (2021).31.Batisse, M. Action plan for biosphere reserves. Environmental conservation 12(1), 17–27 (1985).Article 

Google Scholar 
32.MAB. Criteria for Designation and Evaluation of Unesco Biosphere Reserves in Germany. German National Committee for the UNESCO Programme (2007).33.UNESCO. Biosphere Reserves. The Seville Strategy and the Statutory Framework of the World Network. Paris, France (1995).34.Hadush, M., Holden, S. T. & Tilahun, M. Does population pressure induce farm intensification? Empirical evidence from Tigrai Region, Ethiopia. Agric. Econ. 50(3), 259–277. https://doi.org/10.1111/agec.12482 (2019).Article 

Google Scholar 
35.Tilman, D. et al. Forecasting agriculturally driven global environmental change. Science 292(5515), 281–284. https://doi.org/10.1126/science.1057544 (2001).ADS 
CAS 
Article 
PubMed 

Google Scholar 
36.Holden, E., Linnerud, K., & Banister, D. Sustainable development: Our common future revisited. Global Environmental Change, 26, 130–139 (2014)Article 

Google Scholar 
37.M’Woueni, D., Gaoue, O. G., Balagueman, R. O., Biaou, H. S. & Natta, A. K. Road mediated spatio-temporal tree decline in traditional agroforests in an African biosphere reserve. Glob. Ecol. Conserv. 20, e00796. https://doi.org/10.1016/j.gecco.2019.e00796 (2019).Article 

Google Scholar 
38.Yaméogo, G., Yélémou, B., & Traoré, D. Pratique et perception paysannes dans la création de parc agroforestier dans le terroir de Vipalogo (Burkina Faso). Base (2005).39.Vodouhè, F. G., Adegbidi, A., Coulibaly, O. & Sinsin, B. Parkia biglobosa (Jacq.) R. Br. Ex Benth. Harvesting as a tool for conservation and source of income for local people in Pendjari Biosphere Reserve. Acta Botanica Gallica 158(4), 595–608 (2011).Article 

Google Scholar 
40.Vodouhè, F. G., Coulibaly, O., Biaou, G. & Sinsin, B. Traditional agroforestry systems and biodiversity conservation in Benin (West Africa). Agrofor. Syst. 82(1), 1–13 (2011).Article 

Google Scholar 
41.Bee, J. N., Kunstler, G. & Coomes, D. A. Resistance and resilience of New Zealand tree species to browsing. J. Ecol. 95(5), 1014–1026. https://doi.org/10.1111/j.1365-2745.2007.01261.x (2007).Article 

Google Scholar 
42.Hoffmann, W. A. The effects of fire and cover on seedling establishment in a neotropical savanna. J. Ecol. https://doi.org/10.2307/2261200 (1996).Article 

Google Scholar 
43.Gnangle, P. C. et al. Perceptions locales du changement climatique et mesures d’adaptation dans la gestion des parcs à karité au Nord-Bénin. Int. J. Biol. Chem. Sci. 6(1), 136–149. https://doi.org/10.4314/ijbcs.v6i1.13 (2012).Article 

Google Scholar 
44.Ouoba, H. Y., Bastide, B., Coulibaly-Lingani, P., Kabore, S. A. & Boussim, J. I. Connaissances et perceptions des producteurs sur la gestion des parcs à Vitellaria paradoxa CF Gaertn. (Karité) au Burkina Faso. Int. J. Biol. Chem. Sci. 12(6), 2766–2783 (2018).Article 

Google Scholar 
45.Alencar, N. L., de Sousa Araújo, T. A., de Amorim, E. L. C. & de Albuquerque, U. P. The inclusion and selection of medicinal plants in traditional pharmacopoeias—Evidence in support of the diversification hypothesis. Econ. Bot. 64(1), 68–79. https://doi.org/10.1007/s12231-009-9104-5 (2010).Article 

Google Scholar 
46.Gaoue, O. G. et al.. Theories and major hypotheses in ethnobotany. Economic Botany, 71(3), 269–287 (2017).Article 

Google Scholar 
47.Helm, J. et al. Recovery of Mediterranean steppe vegetation after cultivation: Legacy effects on plant composition, soil properties and functional traits. Appl. Veg. Sci. 22(1), 71–84. https://doi.org/10.1111/avsc.12415 (2019).Article 

Google Scholar 
48.Nash, K. L., Graham, N. A., Jennings, S., Wilson, S. K. & Bellwood, D. R. Herbivore cross-scale redundancy supports response diversity and promotes coral reef resilience. J. Appl. Ecol. 53(3), 646–655. https://doi.org/10.1111/1365-2664.12430 (2016).Article 

Google Scholar 
49.Lambin, E. F. et al. The causes of land-use and land-cover change: Moving beyond the myths. Glob. Environ. Change 11(4), 261–269. https://doi.org/10.1016/S0959-3780(01)00007-3 (2001).Article 

Google Scholar 
50.Camou-Guerrero, A., Reyes-García, V., Martínez-Ramos, M. & Casas, A. Knowledge and use value of plant species in a rarámuri community : A gender perspective for conservation. Hum. Ecol. 36(2), 259–272. https://doi.org/10.1007/s10745-007-9152-3 (2008).Article 

Google Scholar 
51.de Wet, H., Nkwanyana, M. N. & Vuuren, V. S. F. Medicinal plants used for the treatment of diarrhoea in northern Maputaland, KwaZulu-Natal Province, South Africa. J. Ethnopharmacol. J. 130, 284–289. https://doi.org/10.1016/j.jep.2010.05.004 (2010).Article 

Google Scholar 
52.Toledo, V. M., Ortiz-Espejel, B., Cortéz, L., Moguel, P., A., & Ordoñez, M. D. J. The multiple use of tropical forests by indigenous peoples in Mexico: A case of adaptive management. Conserv. Ecol. 7(3), 9. https://www.ecologyandsociety.org/vol7/iss3/art9/ (2003).53.Azihou, F. A. Elephants’ (Loxodonta africana) Impacts on Vegetation Structure and Availability of Plant Species that Other Animals Feed on in the Biosphere Reserve of Pendjari (University of Abomey-Calavi, 2008).
Google Scholar 
54.Faure, P. V. B. Some factors affecting regional differentiation of the soils in the Republic of Benin (West Africa). CATENA 32, 281–306 (1998).Article 

Google Scholar 
55.Tiomoko, D. Gestion de la Réserve de Biosphère de la Pendjari : Modes de gestion et proposition d’un modèle conceptuel de durabilité (Universite d’Abomey-Calavi, 2014).
Google Scholar 
56.ASECNA. Données climatiques, station de Natitingou. Bénin (2010).57.PNP. Plan d’Amenagement Participatif et de Gestion du Parc National de la Pendjari, Bénin 2004–2013 (2009).58.Assédé, E. P., Adomou, A. C. & Sinsin, B. Magnoliophyta, biosphere reserve of Pendjari, Atacora province, Benin. Check List 8(4), 642–661 (2012).Article 

Google Scholar 
59.Houinato, M. & Sinsin, B. La pression agro-pastorale sur la zone riveraine de la Réserve de la Biosphère de la Pendjari. Tropicultura 18(3), 112–117 (2000).
Google Scholar 
60.Gaoue, O. G. Determinant factors for the integrated management of Pendjari hunting reserve northern Benin (Université d’Abomey Calavi, 2000).61.Adomou, C. A. Vegetation patterns and environmental gradients in Benin: Implications for biogeography and conservation. PhD. Dissertation, Wageningen University. PhD. Dissertation, Wageningen University (2005).62.Inoussa, M. et al. Contrasting population structures of two keystone woodland species of W National Park, Niger. S. Afr. J. Bot. 112, 95–101 (2017).Article 

Google Scholar 
63.Shannon, C. E. A mathematical theory of communication. Bell Syst. Tech. J. 27(3), 379–423. https://doi.org/10.1002/j.1538-7305.1948.tb01338.x (1948).MathSciNet 
Article 
MATH 

Google Scholar 
64.R Core Team.R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. https://www.R-project.org/. (2015).65.Chao, A. & Lee, S.-M. Estimating the number of classes via sample coverage. J. Am. Stat. Assoc. 87(417), 210–217 (1992).MathSciNet 
Article 

Google Scholar 
66.Bernard, H. R. Research Methods in Anthropology: Qualitative and Quantitative Approaches (Rowman & Littlefield, 2017).
Google Scholar 
67.Ceuterick, M., Vandebroek, I., Torry, B. & Pieroni, A. Cross-cultural adaptation in urban ethnobotany: The Colombian folk pharmacopoeia in London. J. Ethnopharmacol. 120(3), 342–359 (2008).Article 

Google Scholar 
68.Cohen, N. & Arieli, T. Field research in conflict environments: Methodological challenges and snowball sampling. J. Peace Res. 48(4), 423–435 (2011).Article 

Google Scholar 
69.Tongco, M. D. C. Purposive sampling as a tool for informant selection. Ethnobot. Res. Appl. 5, 147–158 (2007).Article 

Google Scholar 
70.Phillips, O. & Gentry, A. H. The useful plants of Tambopata, Peru: I. Statistical hypotheses tests with a new quantitative technique. Econ. Bot. 47(1), 15–32 (1993).Article 

Google Scholar 
71.Whitney, C. EthnobotanyR: Calculate Quantitative Ethnobotany Indices. Package Version 0.1.8
https://CRAN.R-project.org/package=ethnobotanyR (2021).72.Oksanen, M. et al. Astrocyte alterations in neurodegenerative pathologies and their modeling in human induced pluripotent stem cell platforms. Cell. Mol. Life Sci. 76(14), 2739–2760. https://doi.org/10.1007/s00018-019-03111-7 (2019).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
73.Canty, A. J. Resampling methods in R: The boot package. R News 2(3): 2–7 (2002).74.Akaike, H. Maximum likelihood identification of Gaussian autoregressive moving average models. Biometrika 60(2), 255–265 (1973).MathSciNet 
Article 

Google Scholar 
75.Burnham, K. P. & Anderson, D. R. Multimodel inference: Understanding AIC and BIC in model selection. Sociol. Methods Res. 33(2), 261–304 (2004).MathSciNet 
Article 

Google Scholar  More

The binding of [2Fe–2S] and [3Fe–4S] in clMagRThree conserved cysteines (C60, C124, and C126) of clMagR in a CXnCGC sequence motif (n is 63–65 in most cases) play critical roles in iron–sulfur cluster binding18 (Fig. 1a), which has been further validated by alanines substitution mutant clMagR3M (C60A, C124A, and C126A mutation of clMagRWT). Strep-tagged clMagRWT and clMagR3M were freshly prepared (labeled as “as-isolated”) and purified to homogeneity under aerobic conditions. The clMagRWT protein showed brown color and clMagR3M appeared colorless in the solution, indicating the presence or absence of iron–sulfur cluster, respectively. Consistently, the Ultraviolet–visible (UV–Vis) spectrum of as-isolated clMagRWT showed absorption from 300-to-600-nm region, and with an absorption peak at 325 and 415 nm, and a shoulder at 470 nm, whereas these absorption peaks were abolished in clMagR3M (Fig. 1b). Circular dichroism (CD) spectroscopy was applied to characterize the types of iron–sulfur cluster and their protein environments during cluster maturation42,43,44. As shown in Fig. 1c, clMagRWT shows distinct positive peaks at 371 nm and 426 nm and three negative peaks at 324 nm, 396 nm, and 463 nm, respectively, suggesting the presence of [2Fe–2S] cluster45. However, it is worth pointing out that [4Fe–4S] or [3Fe–4S] clusters usually exhibit negligible CD intensity compared to [2Fe–2S] as shown previously in NifIscA45,46, thus CD spectroscopy cannot exclude the existence of [4Fe–4S] or [3Fe–4S]. Electron paramagnetic resonance (EPR) spectroscopy was then used to analyze different states of as-isolated clMagRWT. The oxidized clMagRWT was S = 1/2 species, characterized by a rhombic EPR signal with g values at g1 = 2.016, g2 = 2.002, and g3 = 1.997 (Fig. 1d) which disappeared at 45 K, suggesting the presence of [3Fe–4S]1+ cluster47,48. After reduced with sodium dithionite (Fig. 1e), EPR signal from [2Fe–2S] cluster can be observed until the temperature increased to 60K49,50,51. Thus, two distinct iron–sulfur clusters were assigned by EPR spectroscopy of clMagRWT. Figure 1Characterization of iron–sulfur clusters in as-isolated clMagR. (a) Sequence alignment of MagR in eight representative species. Predicted secondary structures are shown in the upper lines, with two alpha-helices (orange cylinders) and seven beta-strands (green arrows). Conserved residues with iron–sulfur cluster binding properties are shown in the red background (100% conserved), indicated by stars. Other conserved residues are shown in the gray background and bold fonts. Species’ common name, Latin name and sequence ID in NCBI are listed as follows: Pigeon (Columba livia), XP_005508102.1*; Zebra finch(Taeniopygia guttata), XP_002194930.1*; Fly(Drosophila melanogaster), NP_573062.1*; Monarch butterfly(Danaus plexippus), AVZ24723.1*; Salmon(Salmo salar), XP_013999046.1*; Octopus(Octopus bimaculoides), XP_014786756.1*; Little brown bat(Myotis lucifugus), XP_006102189.1*; Human(Homo sapiens), NP_112202.2*. (b) UV–Vis absorption spectrum of as-isolated pigeon MagR (clMagRWT, black) and C60AC124AC126A substitution mutant (clMagR3M, red), indicating three cysteines contribute to the iron–sulfur cluster binding. SDS-PAGEs of protein preparation are shown as inserts, theoretical mass of the clMagR monomer and clMagR3M monomer were 16.41 kDa, 16.31 kDa, respectively. (c) CD spectrum of as-isolated clMagRWT(black) and clMagR3M(red). (d, e) X-band EPR spectrum of as-isolated clMagRWT at oxidized (d) and reduced status (e). The samples were frozen in TBS buffer and the spectrums were recorded at various temperatures (10 K, 25 K, 45 K, 60 K). (f) Low-temperature resonance Raman spectra of as-isolated clMagRWT. Spectra were recorded at 17 K using 488 nm laser excitation.Full size imageConsidering some iron–sulfur clusters in proteins are diamagnetic and therefore EPR silent, low-temperature Resonance Raman (RR) spectroscopy was then utilized as a probe to characterize those clusters52. With 488 nm excitation, the RR spectra of clMagRWT in the iron–sulfur stretching region (240–450 cm−1) show the presence of [3Fe–4S]1+ cluster (represented by two bridging modes at 286 and 347 cm−1, and one terminal modes at 364 cm−1) and [2Fe–2S]2+ cluster (represented by three iron–sulfur bridging mode at 293, 308 and 330 cm−1 and two terminal modes at 407 and 422 cm−1, as shown in Fig. 1f)52,53,54,55,56. Taking together, we conclude that as-isolated clMagRWT contains both cystine-ligated [2Fe–2S] cluster and [3Fe–4S] cluster.The assembly and conversion of [2Fe–2S] and [3Fe–4S] in clMagRIron–sulfur cluster assembly of IscA, an clMagR homology protein in bacteria, is mediated by cysteine desulfurase IscS2. To elucidate how iron–sulfur cluster assembles in clMagR, time-course experiment was performed, and UV–Vis absorption and CD spectrum were used to monitor the IscS-catalyzed iron–sulfur cluster assembly in clMagR (Fig. 2). No signal of the iron–sulfur cluster was recorded when the reaction begins (0 min), and then the characteristic visible absorption peak and CD spectrum of clMagRWT appeared after 5 min, indicating [2Fe–2S] cluster assembled. As the reaction proceeds, the UV–Vis absorption intensity increased and after 180 min the signal was dominated by a broad shoulder centered at 415 nm (Fig. 2a). Concomitantly, the CD spectrum of the [2Fe–2S] center decreased and then almost disappeared after 180 min, indicating that [2Fe–2S] had been converted to [3Fe–4S] clusters and the reconstitution finished (Fig. 2b).Figure 2Iron–sulfur cluster assembly on clMagR. (a, b) IscS-mediated iron–sulfur cluster assembly on clMagR monitored as a function of time by UV–Vis absorption (a) and CD spectroscopy (b). The spectra shown were taken with samples of pretreated clMagR to remove iron–sulfur clusters before reconstitution (apo-clMagR, 0 min, light green), incubated with IscS after 5 min (green), and after 180 min (dark green). (c, d) chemical reconstitution-mediated iron–sulfur cluster assembly on clMagR monitored as a function of time by UV–Vis absorption (c) and CD spectroscopies (d). The spectra shown were taken with samples of pretreated clMagR to remove iron–sulfur clusters before reconstitution (apo-clMagR, light green) and chemically reconstituted clMagR (chem re clMagR, purple). (e) X-band EPR spectrum of chemically reconstituted clMagRWT. The spectrum was recorded at 10 K. (f) Low-temperature resonance Raman spectra of chemically reconstituted clMagR. Protein and reagent concentrations are described in the Methods. Spectra were recorded at 17 K using 488 nm laser excitation.Full size imageIron–sulfur cluster assembly can be achieved by chemical reconstitution as well, since iron–sulfur apo-proteins are able to spontaneously form iron–sulfur clusters in vitro when supplied with iron and sulfide under reducing conditions1,43,57. With this approach, started with apo-clMagRWT, we successfully reconstituted [3Fe–4S] cluster in clMagR protein, confirmed by UV–Vis absorption and CD spectrum result (Fig. 2c,d). To further validate if [3Fe–4S] is the sole type of iron–sulfur cluster in clMagR after chemical reconstitution, EPR and low-temperature Resonance Raman spectroscopy were applied (Fig. 2e,f). The chemically reconstituted clMagRWT was S = 1/2 species, characterized by a rhombic EPR signal with g values at g1 = 2.017, g2 = 2.002, and g3 = 1.994 (Fig. 2e). The signal is assigned to a S = 1/2 [3Fe–4S]1+ cluster. The Low-temperature Resonance Raman spectrum showed an intense band at 346 cm−1 and additional bands at 406 and 420 cm−1, which demonstrated that chemically reconstituted clMagRWT only contains [3Fe–4S]1+ cluster (Fig. 2f).We further investigated if clMagR could serve as an iron–sulfur carrier protein to accept [2Fe–2S] cluster from scaffold protein such as IscU58. Briefly, 400 µM holo-IscU was mixed with 400 µM strep-tagged apo-clMagRWT and incubated for 180 min under reduced condition, then, after desalting and strep-tactin affinity column separation, UV–Vis absorption and CD spectroscopy were applied the iron–sulfur cluster transfer process (Fig. 3a). The intensity of UV–Vis spectrum decreased in IscU (Fig. 3b) but significantly increased in clMagR after reaction (Fig. 3d), indicating [2Fe–2S] cluster was transferred from IscU to clMagR59. Consistently, CD spectrum of IscU and clMagR also confirmed that [2Fe–2S] transfer occurred between IscU and clMagR (Fig. 3c,e). The resulting spectrum is very similar to that of the [2Fe–2S] intermediate assembled on IscS mediated reconstituted apo-clMagR (Fig. 2b).Figure 3clMagR serve as carrier protein to accept [2Fe–2S] cluster from IscU in vitro. (a) A cartoon schematically illustrates the experimental procedures of in vitro iron–sulfur cluster transfer from IscU to clMagR. (b, c) The UV–Vis absorption (b) and CD spectra (c) of IscU. IscU protein samples were taken before mixing with apo-clMagR (holo-IscU, black lines) and after incubated with apo-clMagR for 180 min (pink lines). (d, e) The UV–Vis absorption (d) and CD spectra (e) of clMagR. clMagR samples were taken before mixing with holo-IscU (apo-clMagR, light green lines) and after incubated with holo-IscU for 180 min (holo-clMagR, brown lines).Full size imageCys-60 is essential for clMagR to bind [3Fe–4S] cluster, not [2Fe–2S] clusterThree conserved cysteines (C60, C124, and C126) of clMagR play critical roles in iron–sulfur cluster binding, and the substitute mutation of these three residues abolished iron–sulfur binding (Fig. 1b,c)18. To elucidate if three cysteines bind [2Fe–2S] and [3Fe–4S] differently, single Cys-to-Ala substitutions (C60A, C124A, and C126A) were made and their iron–sulfur binding properties were characterized.Freshly purified as-isolated clMagRC60A showed light brown color, and [2Fe–2S] cluster binding was verified by UV–Vis absorption and CD spectrum (Fig. 4a,b). A typical protein-bound [2Fe–2S] cluster absorption peak at 325 nm and a shoulder at 415 nm are visible in UV–Vis absorption (Fig. 4a, light orange line). Consistently, the CD spectrum of as-isolated clMagRC60A mutant had a negative peak at 397 nm and a positive peak at 451 nm (Fig. 4b, light orange line), confirmed the [2Fe–2S] cluster binding, similar to clMagRWT. However, in contrast to clMagRWT, chemical reconstitution failed to convert [2Fe–2S] cluster to [3Fe–4S] cluster in clMagRC60A. As shown in Fig. 4a,b (orange line), chemically reconstituted clMagRC60A showed similar and characteristic [2Fe–2S] UV–Vis absorption peaks and CD spectrum, but not [3Fe–4S] (Fig. 4a,b, orange lines), suggesting that C60A mutation abolished [3Fe–4S] cluster binding ability in clMagR.Figure 4Three conserved cysteines play different roles in iron–sulfur binding in clMagR. (a, b) Chemical reconstitution-mediated iron–sulfur cluster assembly on apo-clMagRC60A monitored by UV–Vis absorption (a) and CD spectroscopies (b). The samples of spectra shown are as-isolated clMagRC60A (light orange) and chemically reconstituted clMagRC60A (chem re clMagRC60A, orange). (c, d) chemical reconstitution-mediated iron–sulfur cluster assembly on clMagRC124A monitored by UV–Vis absorption (c) and CD spectroscopies (d). The samples of spectra shown are as-isolated clMagRC124A (light purple) and chemically reconstituted clMagRC124A (chem re clMagRC124A, purple). (e, f) chemical reconstitution-mediated iron–sulfur cluster assembly on pigeon clMagRC126A monitored by UV–Vis absorption (e) and CD spectroscopies (f). The samples of spectra shown are as-isolated clMagRC126A (light blue) and chemically reconstituted clMagRC126A (chem re clMagRC126A, blue). SDS-PAGE results were shown in the right of corresponding UV–Vis spectra as inserts (a, c, e). The theoretical mass of the clMagRC60A monomer, clMagRC124A monomer and clMagRC126A monomer were 16.38 kDa. (g, h) The UV–Vis absorption (c) and CD spectra (d) of clMagRC60A obtained by mixing apo-clMagRC60A and holo-IscU which was recorded before the addition of apo-clMagRC60A (dotted orange lines) and after incubation with apo-clMagRC60A for 180 min (orange lines). Protein and reagent concentrations are described in the Experimental procedures.Full size imageIn contrast, purified as-isolated clMagRC124A and clMagRC126A were colorless, and the binding of iron–sulfur clusters was barely detectable by UV–Vis and CD spectrum (Fig. 4c–f, light purple, and light blue lines, respectively). However, chemical reconstitution successfully reconstituted [3Fe–4S] cluster binding in both clMagRC124A and clMagRC126A (Fig. 4c–f, purple and blue lines, respectively). After chemical reconstitution, the UV–Vis absorption of both clMagRC124A and clMagRC126A mutants showed the signal of iron–sulfur cluster binding (Fig. 4c,e). Parallel CD spectrum studies confirmed both chemically reconstituted clMagRC124A and clMagRC126A have [3Fe–4S] cluster binding (Fig. 4d,f), similar to chemically reconstituted clMagRWT. The results demonstrated that Cys-124 and Cys-126 in clMagR play important roles in [2Fe–2S] cluster binding, thus, mutating these two residues lead to clMagR favors [3Fe–4S] binding.Considering clMagR can act as a carrier protein to accept iron–sulfur cluster from IscU (Fig. 3), it is worth testing if three cysteines play a different role in this process as well. Holo-IscU was mixed with apo-clMagR single cysteine mutants in a reduced state for 180 min. The apo status of all three mutants (labeled as apo-clMagRC60A, apo-clMagRC124A, and apo-clMagRC126A) had no iron–sulfur cluster binding before mixing with holo-IscU, as shown by negligible UV absorption and CD intensities (Fig. 4g,h and Supplementary Fig. 1a–d, dotted lines). After incubation with holo-IscU and separation of IscU and clMagR mutants, clMagRC60A showed distinct changes in UV–Vis absorption and CD spectrum (Fig. 4g,h). The UV–Vis absorption increased and showed better-resolved peaks at 322 nm, 410 nm, 504 nm (Fig. 4g, orange line), and parallel CD spectra had distinct positive peaks (319 nm, 355 nm, 445 nm, and 534 nm) and four negative peaks (333 nm, 392 nm, 477 nm, and 579 nm, Fig. 4h), indicating [2Fe–2S] cluster was transferred from IscU to clMagRC60A. Interestingly, clMagRC124A and clMagRC126A could also accept [2Fe–2S] cluster transferred from holo-IscU, though the binding efficiency is much lower than clMagRWT and clMagRC60A, as verified by UV–Vis and CD spectrum (Supplementary Fig. 1a–d). It seems that clMagRC60A accept [2Fe–2S] cluster from scaffold protein IscU more effectively compared with clMagRC124A and clMagRC126A. And after incubation with clMagR mutants, UV–Vis absorption of IscU significantly decreased, confirmed that iron–sulfur cluster transfer occurred in between holo-IscU and three clMagR mutants (Supplementary Fig. 1e).Again, our data demonstrated that three conserved cystines of clMagR played different roles on the iron–sulfur cluster binding, and especially Cys-60 is essential for clMagR to bind [3Fe–4S] cluster, not [2Fe–2S] cluster. Therefore, it is possible to obtain a [2Fe–2S] cluster binding only clMagR by mutating Cys-60. Thus, we labeled clMagR protein samples based on their iron–sulfur cluster in later experiments. For example, we labeled the chemically reconstituted clMagRWT as [3Fe–4S]-clMagRWT, and clMagRC60A that accepted [2Fe–2S] cluster from holo-IscU as [2Fe–2S]-clMagRC60A, to investigate the magnetic property of clMagR when it binds different iron–sulfur clusters.[3Fe–4S]-clMagR shows different magnetic properties from [2Fe–2S]-clMagRMagR has been reported as a putative magnetoreceptor and exhibits intrinsic magnetic moment experimentally and theoretically when forms complex with cryptochrome (Cry)18,20,21. To elucidate if different iron–sulfur clusters binding in clMagR have different magnetic features and respond to external magnetic fields differently, we obtained [3Fe–4S] and [2Fe–2S] bound only clMagR protein by chemical reconstitution of clMagRWT (as [3Fe–4S]-clMagRWT) and holo-IscU incubated and re-purified clMagRC60A (as [2Fe–2S]-clMagRC60A), respectively, and measured the magnetic moment of these proteins with Superconducting Quantum Interference Device (SQUID) magnetometry. SQUID is a highly sensitive magnetometry to measure extremely subtle magnetic fields and to study the magnetic properties of a range of samples, including extremely low magnetic moment biological samples. Therefore, it has been regularly used as a first test to identify the specific kind of magnetism of a given specimen, such as ferromagnetic, antiferromagnetic, paramagnetic or diamagnetic, by measuring at different temperatures and external magnetic field strength. For example, B-DNA was identified as paramagnetic under low temperature by SQUID60.Purified clMagR3M was utilized as a control since it had no iron–sulfur cluster binding due to lack of cysteine residues (Fig. 1b,c). The magnetic measurement was done at different temperatures (5 K and 300 K) and MH curves (magnetization (M) curves measured versus applied fields (H)) were generated for three proteins to reflect the protein magnetic anisotropy. The MH curves of clMagR3M clearly exhibited diamagnetic property at both 5 K and 300 K, suggesting that magnetism of clMagR is dependent on the iron–sulfur cluster (Fig. 5a,b, red lines). In contrast, [3Fe–4S]-clMagRWT showed superparamagnetic behavior at 5 K which has saturation magnetization (MS) at 2 T about 0.22771 emu/g protein (Fig. 5a, purple line), [2Fe–2S]-clMagRC60A is paramagnetic at 5 K (Fig. 5a, orange line). Interestingly, at higher temperature such as 300 K, [2Fe–2S]-clMagRC60A is diamagnetic while [3Fe–4S]-clMagRWT is paramagnetic (Fig. 5b, orange line and purple line). The different magnetism, as well as the different saturation magnetization of clMagR with different iron–sulfur binding, are clearly important features of this putative magnetoreceptor, and worth further investigation and validation in vivo in the future.Figure 5[3Fe–4S]-clMagRWT shows different magnetic properties from [2Fe–2S]-clMagRC60A. (a) Field-dependent magnetization curves (MH) at 5 K for [2Fe–2S]-clMagRC60A (orange), [3Fe–4S]-clMagRWT (chem re clMagRWT, purple), and clMagR3M (red). The magnetic susceptibility of [2Fe–2S]-clMagRC60A is 2.27749E−6 and the magnetic susceptibility of clMagR3M is − 4.0438E−7. (b) Field-dependent magnetization curves (MH) at 300 K for [2Fe–2S]-clMagRC60A (orange), [3Fe–4S]-clMagRWT (chem re clMagRWT, purple), and clMagR3M (red). And the magnetic susceptibility is − 1.83638E−7, 5.93483E−8, − 3.26432E−7, respectively.Full size image More

Our results show that Chl a alone is not an adequate proxy for prey growth rates in dilution grazing experiments when mixoplankton are present5,10. Chlorophyll is, in any case, a poor proxy for phototrophic plankton biomass31 because of inter-species variations, and also for the photoacclimation abilities of some species (for which very significant changes can occur within a few hours). The problem extends to the involvement of mixoplanktonic prey and grazers. Nevertheless, even very recent studies continue to rely on this parameter for quantifications of grazing despite acknowledging the dominance, both in biomass and abundance, of mixoplanktonic predators in their system30. Moreover, the detailed analysis of the species-specific dynamics revealed that different prey species are consumed at very different rates. In our experiments, and contrary to expectations (see32,33, and Fig. S1 in the Supplementary Information), C. weissflogii was only actively ingested in the ciliate experiment and, according to the results from the control bottles (Table 2), not by M. rubrum (see Fig. 4 and Fig. S1a).Certainly, it is not the first time that a negative selection against diatoms has been seen; for example, Burkill et al.34 noticed that diatoms were less grazed by protist grazers than other phytoplankton species, as assessed by a dilution technique paired with High-Performance Liquid Chromatography for pigment analysis. Using the same method, Suzuki et al.35 reported that diatoms became the dominant phytoplankton group, which suggests that other groups were preferentially fed upon. Calbet et al.36, in the Arctic, also found only occasional grazing over the local diatoms. In our study, diatoms were not only not consumed, but the presence of dinoflagellates appeared to contribute to their growth (Fig. 4), this relationship being partly dependent on the concentration of the predator (see Fig. 2c, d). This result could be a direct consequence of assimilation and use of compounds (e.g.,37,38) released by microplankton such as ammonium (e.g.,39,40) and urea (e.g.,41), which were not supplied in the growth medium, but which would have supported prey growth. Alternatively, this unexpected outcome may have been a consequence of the selective ingestion of R. salina by the two predators, relieving the competition for nutrients and light and resulting in a higher growth rate of the diatom in the presence of the predators. We cannot rule out the fact that diatoms sink faster than flagellates which, as the bottles were not mixed during most of the incubation period (although gently mixed at every sampling point), may have also involuntarily decreased ingestion rates on C. weissflogii. Still, one C. weissflogii cell contains, on average, ca. 2.5 times more Chl a than one R. salina cell (initial value excluded, see Table 3). Taken together with the preference for R. salina it is not surprising that the proportion of total Chl a represented by the diatoms increased over time, in particular in the L/D treatment (Figs. 6a, c and 7a, c).Table 3 Chl a content (pg Chl a cell−1) of the target species at each sampling point as calculated from the control bottles.Full size tableAnother factor clearly highlighted by our experiments, is that protozooplankton themselves contribute a significant portion of the total chlorophyll of the system (due to ingested Chl a), in particular at the beginning of the incubation (see Figs. 6 and 7); this being invariably ignored in a traditional dilution experiment. The high Chl a detected inside the protozooplanktonic grazers at the beginning of the incubations could suggest that the system was initially not in equilibrium, and that this was the result of superfluous feeding (e.g.,42). This would, nevertheless, be surprising since we required ca. 1 h to collect the initial samples (t = 0 h) after joining all the organisms together (see the section “Dilution grazing experiments” in the “Methods” section); previous studies, like the one on G. dominans and Oxyrrhis marina by Calbet et al.42, showed that the hunger response and consequent vacuole replenishment occurred in ca. 100 min for very high prey concentrations and it is expected to decrease at lower prey concentrations as the ones used in our study. Therefore, even if one assumes that the first 4 h of incubation are a result of superfluous feeding, after 24 h, the “estimated”, “observed”, and “from dilution slope” grazing estimates are not significantly different to those displayed in Fig. 5 (P  > 0.05 in all instances) and, therefore, we can assume that the hunger response was likely irrelevant (e.g.,43) and did not mask our results. In any case, as stated before, an actual field grazing dilution experiment also suffers from similar problems, because grazers and prey are suddenly diluted and not pre-adapted to distinct food concentrations. Nevertheless, this is not novel information, since Chl a and its degradation products have been found inside several protozooplankton species from different phylogenetic groups immediately after feeding44 and even after some days without food45. An increase in intracellular Chl a concentrations immediately after feeding has also been found in mixoplankton46,47, on which this increase is derived both from ingested prey as well as from new synthesis of their own Chl a. Additionally, several experiments with Live Fluorescently Labelled Algae (LFLA) show that predators (irrespective of their trophic mode) seem to maximise the concentration of intracellular prey shortly after the initiation of the incubation (e.g.,48; Ferreira et al., submitted). Indeed, some authors have even been able to measure photosynthesis in protozooplankton, like the ciliates Mesodinium pulex49 and Strombidinopsis sp.50.The fact that Chl a is a poor indicator of phytoplankton biomass and the inherent consequences discussed so far can be solved by the quantification of the prey community abundance (e.g.,51) by microscopy or by the use of signature pigments for each major phytoplankton group. The latter method, however, is not as thorough as the former, since rare are the cases where one pigment is exclusively associated with a single group of organisms (see52 and references therein). In any case, any pigment-based proxy is subject to the same problems, as identified by Kruskopf & Flynn31. Irrespective of the quantification method, it has been made evident that the different algae are consumed at different rates (e.g., pigments10,34,35; microscopy5,36).Prey selection in protistan grazers is a common feature (e.g.,23,26,27,28). Given the diversity of grazers in natural communities and the array of preferred prey that each particular species possesses, it is logical to think that dilution experiments will capture the net community response properly. Likewise, grazers interact with each other through toxins, competition, and intraguild predation among other factors. An example of intraguild predation could be the observed on K. armiger by G. dominans (see Figs. 2f and 4 and Table 1), which caused an average loss of ca. 18.72 pg of K. armiger carbon per G. dominans per hour in the D treatment. Interestingly, in the same treatment, a slight negative effect of K. armiger on its predator G. dominans can also be deduced (i.e., positive g, Table 1), resulting in an average loss of ca. 0.33 pg G. dominans carbon per K. armiger per hour. This could be a consequence of algal toxins, since K. armiger is a known producer of karmitoxin22, whose presence may have negative effects even on metazoan grazers21. Regarding ciliates, none of the species used is a known producer of toxic compounds, which suggests that the average loss of ca. 1.25 pg M. rubrum carbon per hour in the D treatment was due to S. arenicola predation. Altogether, it seems clear from our data that intraguild predation cannot be ignored when analysing dilution experiments (Fig. 4). Furthermore, our results clearly show that single functional responses cannot be used to extrapolate community grazing impacts, as evidenced by the differences in estimated and measured ingestion rates based on the disappearance of prey in combined grazers experiments (Fig. 5). Nevertheless, this is a relatively common procedure (e.g.,53 and references therein). Often in modelling approaches, individual predator’s functional responses have been used to extrapolate prey selectivity and community grazing responses27; in reality complex prey selectivity functions are required to satisfactorily describe prey selectivity and inter-prey allelopathic interactions54.It is, however, also evident that the measured ingestion rates in combined grazers experiments were not the same as those calculated from the slope of the dilution grazing experiment. This raises the question of why was that the case. It is well known that phytoplankton cultures, when extremely diluted, show a lag phase of different duration55 which has been attributed to the net leakage of metabolites56. Assuming that the duration of the lag phase will be dependent on the level of dilution, it seems reasonable to deduce that after ca. 24 h the instantaneous growth rates (µ) in the most diluted treatments will be lower than that of the undiluted treatments. This has consequences, not only for the estimated prey growth rates but also for the whole assessment of the grazing rate, due to the flattening of the regression line (i.e., the decrease in the computed growth rate). This artefact may be more evident in cultures acclimated to very particular conditions (as the laboratory cultures used in this study) than in nature.Another important finding of our research is the importance of light on the correct expression of the feeding activity by both mixoplankton and protozooplankton. We noticed that irrespective of the light conditions, all species exhibited a diurnal feeding rhythm (R. salina panels in Figs. 2 and 3), which is in accordance with earlier observations on protists (e.g.,29,57,58). The presence of light typically increased the ingestion rates. Additionally, the ingestion rates differed during the night period between L/D and D treatments, which implies that receiving light during the day is also vital in modulating the night behaviour of protoozoo- and mixoplankton. In particular, mixoplankton grazing is usually affected by light conditions, typically increasing (e.g.,32,59), but also sometimes decreasing(e.g.,60) in the presence of light. Different irradiance levels can also affect the magnitude of ingestion rates both in protozoo- and mixoplankton (see61 and references therein).For those reasons, we hoped for a rather consistent pattern among our protists that would help us discriminate mixoplankton in dilution grazing experiments. As a matter of fact, based on the results from Arias et al.29, we expected that in the dinoflagellate experiment, the D treatment would have inhibited only the grazing of K. armiger, enabling a simple discrimination between trophic modes. The reality did not meet the expectations since the day and night-time carbon-specific ingestion rates (as assessed using the control bottles, Table 2) of K. armiger were respectively higher and equal than those of G. dominans. Conversely, in the ciliate experiment, protozooplankton were the major grazers in our incubations regardless of the day period and light conditions. This response was not as straightforward as one would expect it to be because M. rubrum has been recently suggested to be a species complex containing at least 7 different species (62 and references therein), which hinders any possible conjecture on their grazing impact. Indeed, the uneven responses found between and within trophic modes precluded such optimistic hypothetical procedure.The D treatment in the present paper illustrated the importance of mimicking natural light conditions, a factor also addressed in the original description of the technique by Landry and Hassett1. It is crucial for the whole interpretation of the dilution technique that incubations should be conducted in similar light (and temperature) conditions as the natural ones to allow for the continued growth of the phototrophic prey. However, here we want to stress another aspect of the incubations: should they start during the day or the night? Considering our (and previous) results on diel feeding rhythms, and on the contribution of each species to the total Chl a pool, it is clear that different results will be obtained if the incubations are started during the day or the night. Besides, whether day or night, organisms are also likely to be in a very different physiological state (either growing or decreasing). Therefore, we recommend that dilution experiments conducted in the field should always be started at the same period of the day to enable comparisons (see also Anderson et al.14 for similar conclusions on bacterivory exerted by small flagellates). Ideally, incubations would be started at different times of the day to capture the intricacies of the community dynamics on a diel cycle. Nevertheless, should the segmented analysis be impossible, we argue that the right time to begin the incubations would be during the night, as this is the time where ingestion rates by protozooplankton are typically lower (e.g.,29,57,58, this study) and would, consequently, reduce their quota of Chl a in the system.Lastly, we want to stress that we are aware that our study does not represent natural biodiversity because our experiments were conducted in the laboratory with a few species. Nevertheless, we attempted to use common species of wide distribution for each major group of protists to provide a better institutionalisation of our conclusions. Further to the choice of predator and prey is their concentrations and proportions. Being a laboratory experiment designed to understand fundamental mechanisms within a dilution grazing experiment, we departed from near saturating food conditions from where we started the dilution series. In nature, the concentrations that we used may be high but are not unrealistic, and actually lower than in many bloom scenarios. We included diatoms at high concentrations, even knowing that they are not the preferred prey of most grazers34, because diatoms are very abundant in many natural ecosystems and to stress the point of food selection within the experiment. For sure, using different proportions of prey would have rendered different results. However, as previously mentioned, our aim was not to seek flaws in the dilution technique, but to understand the role of mixoplankton in these experiments and the complex trophic interactions that may occur within. Ultimately, with our choice of prey and their concentrations, we have proven that when there is no selection for a massively abundant prey, the use of Chl a as a proxy for community abundances may underestimate actual grazing rates.Some other aspects of our experiments may also be criticised because they do not fully match a standard dilution experiment. For instance, we manipulated light, adding complexity to the study. However, this manipulation enabled the deepening into the drivers of the mixoplanktonic and protozooplanktonic grazing responses. Another characteristic, perhaps awkward, of our study is that we allowed the grazers to deplete their prey before starting the experiment. One may argue this procedure does not mimic the natural previous trophic history a grazer may have in nature. Yet, in nature, when facing a dilution experiment, it is impossible to ascertain whether the organisms are encountering novel prey or not. Indeed, they (prey and predator) could have just migrated into such conditions, or be subject to famine, or just moved from a food patch. In any case, it is true that a consistent “hunger response” would have affected our initial grazing values, biasing grazing rate estimates. To overcome this artefact, we let the grazers feed for about one hour before starting the actual dilution assay (see the “Methods” section). From that point on, any dilution is, in fact, an abrupt alteration of the food scenario, which is likely more important than the previous trophic history of the grazer.In summary, with these laboratory experiments, we have presented evidence calling for a revision of the use of chlorophyll in dilution grazing experiments5,10, and we have highlighted the need to observe the organismal composition of both initial and final communities to better understand the dynamics during the dilution grazing experiments51. This approach will not incorporate mixoplanktonic activity into the dilution technique per se however if combined with LFLA (see5,17), a semi-quantitative approach to disentangle the contribution of mixoplankton to community grazing could be achieved (although not perfect). An alternative (and perhaps more elegant) solution could be the integration of the experimental technique with in silico modelling. The modelling approaches of the dilution technique have already been used, for example, to disentangle niche competition63 and to explore nonlinear grazer responses20. We believe that our experimental design and knowledge of the previously indicated data could be of use for the configuration of a dilution grazing model, which could then be validated in the field (and, optimistically, coupled to the ubiquitous application of the dilution technique across the globe). We cannot guarantee that having a properly constructed model that mimics the dilution technique will be the solution to the mixoplankton paradigm. However, it may provide a step towards that goal as it could finally shed much-needed light on the mixo- and heterotrophic contributions to the grazing pressure of a given system. To quote from the commentary of Flynn et al.6, it could provide the answer to the question of whether mixoplankton are de facto “another of the Emperor’s New Suit of Clothes” or, “on the other hand (…) collectively worthy of more detailed inclusion in models”. More

Our results showed that covariates at different geographical scales (national and local) may have important effects on the risk of snakebite, both for humans and animals. The results indicate that the risk of snakebite in the Terai varies at national scale between clusters and at local scale between households. The evaluation of the final models without spatial random components and the worsening of the models’ goodness of fit as a result highlighted how snakebite risk and its determining factors are indeed spatially structured.A strong association between high snakebite incidence and mortality, and poverty was established from the analysis of 138 countries affected by the disease32. In this study, we identified the PPI, an indicator for poverty, as a highly influential risk-increasing factor for humans. This not only confirms the critical role of poverty as a driver for this Neglected Tropical Disease, but also offers the possibility to use a standardized index at individual household scale for similar studies. Chaves et al.33 used the Poverty Gap, which is a simpler index expressing how far a person is from the average national poverty line, but to our knowledge, no study has used PPI for snakebite in any way. Applying PPI as a snakebite risk predictor also addresses previous expert calls for an Ecohealth approach to consider the relationship between the structural characteristics of houses, poverty, and snakebite34.Three of the survey covariates had significant effects on the odds of snakebite. Food storage and straw storage increased them, while sleeping on the floor reduce them. The effect of the first two covariates is likely to be related to prey availability, represented by rodents, which are attracted by food and shelter sources. Both food and straw are very often stored near dwellings, which in the end multiply the number of possible encounters between humans, domestic animals, and the hunting snakes20. The expected snakebite risk reduction effect by sleeping on the floor is more complex though. Previously, a higher snakebite incidence was reported among rural Hindus in Maharashtra, India, due to their custom of sleeping on the floor35,36, while in Nepal, Chappuis et al. did not find any protective effect or significant difference in snakebite cases between sleeping on a cot or on the floor37. This result, nevertheless, might be influenced again by regional customs that make sleeping on the floor more common in eastern Terai (71.1% of all affirmative answers to this question), and second, by the commonly acknowledged prevalence of kraits (Bungarus spp.) in western Terai, which are the species most commonly linked to bites to people sleeping on the floor while hunting at night inside houses22,38. This geographic separation, between the human behaviour and the distribution of the species considered to cause most bites linked to it, could explain the observed shift in the odds towards a reduction effect. This effect should be further explored in localized studies designed to capture behavioural differences in humans and snakes.For both the general human risk model and its equivalent prediction model, the covariate Distance to water had a significant risk-increasing effect. For each additional km in distance from permanent water sources, the odds of snakebite increased by 1.38 and 1.51 times, respectively. From a human perspective and in this socio-economic framework, it would be important to consider not only the distance to water, but also the path taken to get the water (or any other resource). If this path would lead a person through grasslands and open fields, this could imply an increased risk of snakebite. From an ecological perspective, there are two important aspects to consider in relation to water sources. One is, as in this study, the distance from large, constant water sources, which usually represent stable environments subject to less hydric stress. The second (not considered here) are the human-made water sources, such as ponds, reservoirs, and paddy fields that change often, are usually closer to human dwellings, and are known to attract some medically important venomous snakes (MIVS)5. Studies on snake migration and home range use have concluded that depending on species and ecological conditions, snakes can move between a few tens of meters per day and more than 10 km between seasons, while searching for water and prey resources38,39,40,41. In sub-tropical regions like the Terai, snakes living closer to continuous sources of water and vegetation should have easier access to a wider variety of prey. On the contrary, those living in agricultural areas might need to scout farther in the search for resources, encountering human-made waterbodies and prey, such as rodents42 and amphibians, abundant in this region10. Further studies considering all sources of water, and species ecology, biology and richness would be necessary to completely understand the effect of this and similar eco-physiological covariates.Another important factor was the NDVI, which is a commonly used value to express photosynthetic activity, leaf production and in summary the ‘greenness’ of the environment43. As is the case for other covariates, its interpretation depends on the study circumstances. In Iran, it was considered an indicator of prey availability for snakes and linked to snake habitat suitability14. Elevated NDVI values have been associated with higher number of hospitalizations in Nigeria and northern Ghana, in particular during the periods of high agricultural activity, which is also related to higher snake-human contact and higher snakebite incidence43. In our study, its ‘protective’ effect can indeed be the consequence of better access to prey associated with healthier ecosystems, explained in the Terai by the higher NDVI values of the multiple dense forests distributed along the region. In addition, the averaged NDVI values for agricultural areas should be lower than those for perennial forests, because they include the highs and lows of production and harvest.Environmental drivers like temperature and precipitation are common factors in geospatial analyses of snakebite13,14,17,44. They are found in many cases to be the main factors modulating the incidence or risk of snakebite, while varying in importance according to study conditions. For example, in Iran, precipitation seasonality was the most prevalent climatic covariate determining the habitat suitability leading to snakebite risk14, while in Mozambique, temperature seasonality was the predominant covariate13. Despite the Terai’s sub-tropical climate, the range of the average minimum temperature of the coldest month (BIO6) was 1.8–10.9 °C. For our snakebite risk analysis in animals, an increase of 10 °C of BIO6 between any two points represented an increase in the odds of snakebite of 23.41 times. For snakes, this range could be the difference between total lethargy and partial activity45, which could lead to increased numbers of snakebites. In addition, and according to the production and holding practices of domestic animals in the Terai, this temperature range can also represent the difference between animals (mainly ruminants) being kept in sheds when at the lower range limits, or being let out of them at the upper limits, which would again increase the chances of encounters with snakes.Similarly, for the animal model, pig density and sheep density, significantly influenced the variation in the risk of snakebite for animals in the Terai. This could be due to the conditions in which the animals and their feed are kept, favouring environments that are beneficial for either snakes or their prey. At more local scales, rather than the distribution, the presence of other animal species could instead be the factor associated with higher snakebite rates12. However, since the available data on domestic animal density was produced more than 10 years ago, and the animal population has grown substantially in the last years in Nepal, this outcome should be interpreted with caution.For the animal risk, the possession of an animal shed also significantly increased the odds of snakebite. Similar to straw storage, animal sheds and similar constructions offer some shelter and at the same time attract small (prey) animals, both of which are likely to attract snakes, increasing snakebite risk for the animals using the shed. If in addition, the sheds function as poultry coops, the snake hunting behaviour might be instead targeted towards chicks and chickens12. Mitigation measures such as raising the coop’s floor or securing openings with fine metal mesh have been suggested to reduce this risk12.The human modification of terrestrial systems was the only non-significant covariate in the animal risk model. However, as its strong, risk-reducing effect still seems to explain a lot of the response variation, it was retained. Its change in one unit, i.e., going from a pristine to fully modified environment, decreased the odds of snakebite by 0.13 (equivalent to 7.69 times), which agrees with previous national survey results from Sri Lanka21.For our human risk prediction model, four covariates were either significant or helped to explain the changes in the response. Distance to water and NDVI were clearly significant, and precipitation of the driest quarter (BIO17) and the mean annual temperature (BIO1) helped to explain some of the response variation with convincing, unambiguous effects. For BIO17, an increase of 100 mm of rain during the driest months of the year represented an odds-reduction effect equivalent to 8.33 times. This agrees with the results of distance to water, suggesting that the additional availability of resources during water shortage periods, i.e., almost four times more rain (BIO17 range: 18–71 mm), could locally improve ecological conditions for snakes also leading to less scouting and fewer human encounters. Previous studies have analysed the multilevel ecological effects of droughts, e.g., reducing snake prey and leading snakes to engage in riskier behaviours46,47. For BIO1, the protective effect was weaker. An increase of 10 °C represented a reduction of the odds of snakebite equivalent to 3.57 times. Average temperatures for specific locations are difficult to interpret, since they might depend on mild highs and lows, strong highs and lows, or relative combinations of both. Thus, despite having a relatively important effect on the response, this effect still might be the consequence of confounding and unknown interactions.Several other evaluated covariates, for both humans and animals, showed a negligible effect on describing the response, were not significant while having very large uncertainties, or both. Consequently, they were discarded as predicting factors. For the list of baseline covariates evaluated, see supplementary Table S1. For a complete list of available survey covariates, see Alcoba et al.27.Some of our discarded covariates have been important in other studies, for example, to quantify snakebite risk based on reclassification methods of covariates such as habitat suitability, species presence, or envenoming severity13,14,17,44,48. These methods are especially relevant when one species (or very few) is the cause of most snakebite cases, and has differentiated optimal and sub-optimal habitats. In Nepal, and particularly in the Terai, there are at least two, and sometimes more than 10 MIVS with overlapping distributions49. Thus, it could be said that practically the whole region offers suitable habitat for multiple MIVS. In addition, the impossibility of reliably identifying the species having bitten the surveyed victims hindered the use of single species in the analysis. In our analysis, species richness was removed, as it showed almost no effect on the response. A recent meta-analysis reported an equivalent result at global scale, finding no significant difference between the number of venomous snake species in tropical and temperate locations, while the number of snakebites is clearly higher in tropical areas50. These results suggested that high incidence of snakebite is unrelated to species richness, but instead related to other factors like the number of people working in agricultural environments21,32,50. Another important driver of snakebite incidence has been population density50. In our study, however, any possible effect from population density on the risk was diminished by the random selection of households at specific numbers during study design. This was later confirmed by the minimal effect of population density as covariate in the human risk analysis.This study presents a few limitations. For instance, despite the capacity of the INLA method to borrow strength from neighbouring observations and areas, the selection of adequate covariates with enough explanatory power still depends greatly on the number of snakebite cases, which even for a national scale study like this remains small. Also, some of the covariates with the strongest explanatory power came from our household survey, which prevented their use for generalized spatial prediction models. Concerning the animal risk analysis, due to the small number of snakebite cases we opted to aggregate all animal species and consider a grouped response. Thus, for a spatial analysis of animal risk, it was not worth it to consider each species, since that would dilute further an already sparse dataset in individual models and selection processes. Moreover, the data gathered for animals was dependent on the random selection of (human) households and unrelated to the current distribution of animal populations. This, in addition to the possible number of dry bites that go unnoticed, might be responsible for the low number of animal victims recorded (even combined across all species), making a more detailed analysis unfeasible.Despite the large number of covariates examined during our analysis, very few were useful to predict snakebite risk along the Terai. It is possible that confounders or other difficult-to-measure covariates could better explain the complex relationship between the ecology and biology of MIVS, socio-economic factors, human behavioural traits, and the circumstances around domestic animal keeping. This needs to be further explored, following a recent call for an overarching One Health and Ecohealth approach to better understand the drivers for snakebite risk, incidence, and mortality under specific situations34.In conclusion, snakebite is a multi-factorial disease and there is no possible universal approach to model its risk. Each model should be individually designed for each set of socio-economical, geographic, ecological, cultural, and environmental circumstances19. To better understand and address the snakebite problem, it is necessary to approach it, whenever possible, with local data collected at a national scale, so that the conclusions drawn can fuel appropriate national public health policies and actions. As long as people work, live, and keep their domestic animals in close contact with natural environments with MIVS, the risk of snakebite will be present. However, better understanding of the factors influencing that risk at the most granular scale possible, and the estimation of the populations at risk, can help to better target prevention and mitigation measures. For humans, this evidence can channel efforts towards improved access to treatment through the optimized stockpiling of antivenom, and the improvement, relocation or new construction of treating facilities, but more importantly, towards community education and sensitization in preventive campaigns51. Part of that preventive and educative efforts can be done at household level, by promoting and facilitating the use of protective equipment such as rubber boots, or the guidance on how to improve and adapt their immediate surroundings to make them ecologically less attractive for snakes and their prey. For domestic animals, this information could help better target awareness-raising activities for animal owners and implement mitigation strategies. For animals at higher risk, tailored interventions such as the improvement of housing conditions, regular cleaning of sheds and surrounding areas (e.g., from food waste and surrounding vegetation), and using light when animals are walked out of the enclosure at night could be deployed specifically as snakebite prevention measures52. It is also important to highlight that many of the factors analysed in this study affect most directly the snakes themselves, not only as snakebite agents, but also as a diverse group of species, differently affected by ecological, climatic and environmental factors in a multiplicity of settings shared with humans and domestic animals. It is therefore necessary to further investigate how those factors influence the behavioural and ecological traits of MIVS in order to truly understand this disease from a One Health viewpoint. At stake is the reduction of snakebite envenoming incidence rates in humans and animals, and of its possible long-term sequelae on human populations. More

1.Muok AR, Briegel A. Intermicrobial hitchhiking: how nonmotile microbes leverage communal motility. Trends Microbiol. 2021;29:542–50.CAS 
PubMed 

Google Scholar 
2.Kohlmeier S, Smits THM, Ford RM, Keel C, Harms H, Wick LY. Taking the fungal highway: mobilization of pollutant-degrading bacteria by fungi. Environ Sci Technol. 2005;39:4640–6.CAS 
PubMed 

Google Scholar 
3.Simon A, Bindschedler S, Job D, Wick LY, Filippidou S, Kooli WM, et al. Exploiting the fungal highway: development of a novel tool for the in situ isolation of bacteria migrating along fungal mycelium. FEMS Microbiol Ecol. 2015;91:fiv116.PubMed 

Google Scholar 
4.Deveau A, Bonito G, Uehling J, Paoletti M, Becker M, Bindschedler S, et al. Bacterial–fungal interactions: ecology, mechanisms and challenges. FEMS Microbiol Rev. 2018;42:335–52.CAS 
PubMed 

Google Scholar 
5.Harms H, Schlosser D, Wick LY. Untapped potential: exploiting fungi in bioremediation of hazardous chemicals. Nat Rev Microbiol. 2011;9:177.CAS 
PubMed 

Google Scholar 
6.Otten W, Hall D, Harris K, Ritz K, Young IM, Gilligan CA. Soil physics, fungal epidemiology and the spread of Rhizoctonia solani. N. Phytol. 2001;151:459–68.
Google Scholar 
7.Sun B, Chen X, Zhang X, Liang A, Whalen JK, McLaughlin NB. Greater fungal and bacterial biomass in soil large macropores under no-tillage than mouldboard ploughing. Eur J Soil Biol. 2020;97:103155.CAS 

Google Scholar 
8.Otto S, Bruni EP, Harms H, Wick LY. Catch me if you can: dispersal and foraging of Bdellovibrio bacteriovorus 109J along mycelia. ISME J. 2017;11:386–93.PubMed 

Google Scholar 
9.Kjeldgaard B, Listian SA, Ramaswamhi V, Richter A, Kiesewalter HT, Kovács ÁT. Fungal hyphae colonization by Bacillus subtilis relies on biofilm matrix components. Biofilm. 2019;1:100007.PubMed 
PubMed Central 

Google Scholar 
10.Narr A, Nawaz A, Wick LY, Harms H, Chatzinotas A. Soil viral communities vary temporally and along a land use transect as revealed by virus-like particle counting and a modified community fingerprinting approach (fRAPD). Front Microbiol. 2017;8:1975.PubMed 
PubMed Central 

Google Scholar 
11.Rosner A, Gutstein R. Adsorption of actinophage Pal 6 to developing mycelium of Streptomyces albus. Can J Microbiol. 1981;27:254–7.CAS 
PubMed 

Google Scholar 
12.Ghanem N, E. Stanley C, Harms H, Chatzinotas A,Y, Wick L. Mycelial effects on phage retention during transport in a microfluidic platform. Environ Sci Technol. 2019;53:11755–63.CAS 
PubMed 

Google Scholar 
13.Dennehy JJ. What ecologists can tell virologists. Annu Rev Microbiol. 2014;68:117–35.CAS 
PubMed 

Google Scholar 
14.Hurst CJ, Gerba CP, Cech I. Effects of environmental variables and soil characteristics on virus survival in soil. Appl Environ Microbiol. 1980;40:1067–79.CAS 
PubMed 
PubMed Central 

Google Scholar 
15.Yeager JG, Brien RT. Enterovirus inactivation in soil. Appl Environ Microbiol. 1979;38:694–701.CAS 
PubMed 
PubMed Central 

Google Scholar 
16.Schwartz DA, Lindell D. Genetic hurdles limit the arms race between Prochlorococcus and the T7-like podoviruses infecting them. ISME J. 2017;11:1836–51.PubMed 
PubMed Central 

Google Scholar 
17.Shan J, Ramachandran A, Thanki AM, Vukusic FBI, Barylski J, Clokie MRJ. Bacteriophages are more virulent to bacteria with human cells than they are in bacterial culture; insights from HT-29 cells. Sci Rep. 2018;8:5091.PubMed 
PubMed Central 

Google Scholar 
18.Chaudhry W, Lee E, Worthy A, Weiss Z, Grabowicz M, Vega NM, et al. Mucoidy, a general mechanism for maintaining lytic phage in populations of bacteria. FEMS Microbiology Ecology. 2020;96:fiaa162.19.Yu Z, Schwarz C, Zhu L, Chen L, Shen Y, Yu P. Hitchhiking behavior in bacteriophages facilitates phage infection and enhances carrier bacteria colonization. Environ Sci Technol. 2020;55:2462–72.PubMed 

Google Scholar 
20.Tarafder AK, von Kügelgen A, Mellul AJ, Schulze U, Aarts DGAL, Bharat TAM. Phage liquid crystalline droplets form occlusive sheaths that encapsulate and protect infectious rod-shaped bacteria. Proc Natl Acad Sci. 2020;117:4724–31.CAS 
PubMed 
PubMed Central 

Google Scholar 
21.Callaway RM, Ridenour WM. Novel weapons: invasive success and the evolution of increased competitive ability. Front Ecol Environ. 2004;2:436–43.
Google Scholar 
22.Granato ET, Meiller-Legrand TA, Foster KR. The evolution and ecology of bacterial warfare. Curr Biol. 2019;29:521–37.
Google Scholar 
23.Gama JA, Reis AM, Domingues I, Mendes-Soares H, Matos AM, Dionisio F. Temperate Bacterial viruses as double-edged swords in bacterial warfare. PLoS One. 2013;8:e59043.CAS 
PubMed 
PubMed Central 

Google Scholar 
24.Dragoš A, Andersen AJC, Lozano-Andrade CN, Kempen PJ, Kovács ÁT, Strube ML. Phages carry interbacterial weapons encoded by biosynthetic gene clusters. Curr Biol. 2021;31:3479–89.PubMed 

Google Scholar 
25.Pyšek P, Bacher S, Kühn I, Novoa A, Catford JA, Hulme PE, et al. Macroecological framework for invasive aliens (MAFIA): disentangling large-scale context dependence in biological invasions. NeoBiota. 2020;62:407–61.
Google Scholar 
26.Blackburn TM, Pyšek P, Bacher S, Carlton JT, Duncan RP, Jarošík V, et al. A proposed unified framework for biological invasions. Trends Ecol Evol. 2011;26:333–9.PubMed 

Google Scholar 
27.Richardson DM, Pyšek P. Plant invasions: merging the concepts of species invasiveness and community invasibility. Prog Phys Geogr Earth Environ. 2006;30:409–31.
Google Scholar 
28.Williamson M. Explaining and predicting the success of invading species at different stages of invasion. Biol Invasions. 2006;8:1561–8.
Google Scholar 
29.Demerec M, Adelberg EA, Clark AJ, Hartman PE. A proposal for a uniform nomenclature in bacterial genetics. Genetics 1966;54:61–76.CAS 
PubMed 
PubMed Central 

Google Scholar 
30.Dechesne A, Wang G, Gülez G, Or D, Smets BF. Hydration-controlled bacterial motility and dispersal on surfaces. Proc Natl Acad Sci. 2010;107:14369–72.CAS 
PubMed 
PubMed Central 

Google Scholar 
31.Maurhofer M, Keel C, Schnider U, Voisard C, Haas D, Defao G. Influence of enhanced antibiotic production in Pseudomanas fluorescens strain CHA0 on its disease suppressive capacity. Phytopathol. 1992;82:190–5.CAS 

Google Scholar 
32.Schamfuß S, Neu TR, van der Meer JR, Tecon R, Harms H, Wick LY. Impact of mycelia on the accessibility of fluorene to PAH-degrading bacteria. Environ Sci Technol. 2013;47:6908–15.PubMed 

Google Scholar 
33.Bichet MC, Chin WH, Richards W, Lin Y-W, Avellaneda-Franco L, Hernandez CA, et al. Bacteriophage uptake by mammalian cell layers represents a potential sink that may impact phage therapy. iScience. 2021;24:102287.CAS 
PubMed 
PubMed Central 

Google Scholar 
34.Lu F, Wu S-H, Hung Y, Mou C-Y. Size effect on cell uptake in well-suspended, uniform mesoporous silica nanoparticles. Small. 2009;5:1408–13.CAS 
PubMed 

Google Scholar 
35.Fortier L-C, Moineau S Phage production and maintenance of stocks, including expected stock Lifetimes. In: Clokie MRJ, Kropinski AM, editors. Bacteriophages: Methods and Protocols, Volume 1: Isolation, Characterization, and Interactions. Totowa: Humana Press; 2009. p. 203–19.36.Mazzocco A, Waddell TE, Lingohr E, Johnson RP Enumeration of Bacteriophages Using the Small Drop Plaque Assay System In: Clokie MRJ, Kropinski AM, editors. Bacteriophages: Methods and Protocols, Volume 1: Isolation, Characterization, and Interactions Totowa: Humana Press; 2009. p. 81–85.37.Kropinski AM, Mazzocco A, Waddell TE, Lingohr E, Johnson RP Enumeration of Bacteriophages by Double Agar Overlay Plaque Assay In: Clokie MRJ, Kropinski AM, editors. Bacteriophages: Methods and Protocols, Volume 1: Isolation, Characterization, and Interactions. Totowa: Humana Press; 2009. p. 69–76.38.Thanki AM, Taylor-Joyce G, Dowah A, Yakubu Nale J, Malik D, Rebecca Jane Clokie M. Unravelling the Links between Phage Adsorption and Successful Infection in Clostridium difficile. Viruses. 2018;10:441.39.Nair RR, Fiegna F, Velicer GJ. Indirect evolution of social fitness inequalities and facultative social exploitation. Proc R Soc B Biol Sci. 2018;285:20180054.
Google Scholar 
40.Postma M, Goedhart J. PlotsOfData—A web app for visualizing data together with their summaries. PLOS Biol. 2019;17:e3000202.CAS 
PubMed 
PubMed Central 

Google Scholar 
41.Wood M. Statistical inference using bootstrap confidence intervals. Significance. 2004;1:180–2.
Google Scholar 
42.Cumming G, Finch S. Inference by eye: confidence interval and how to read pictures of data. Am Psychol. 2005;60:170–80.PubMed 

Google Scholar 
43.Frada MJ, Schatz D, Farstey V, Ossolinski JE, Sabanay H, Ben-Dor S, et al. Zooplankton may serve as transmission vectors for viruses infecting algal blooms in the ocean. Curr Biol. 2014;24:2592–7.CAS 
PubMed 

Google Scholar 
44.Frada MJ, Vardi A. Algal viruses hitchhiking on zooplankton across phytoplankton blooms. Commun Integr Biol. 2015;8:e1029210.PubMed 
PubMed Central 

Google Scholar 
45.Totsche KU, Kögel-Knabner I. Mobile organic sorbent affected contaminant transport in soil: numerical case studies for enhanced and reduced mobility. Vadose Zo J. 2004;3:352–67.CAS 

Google Scholar 
46.Reche I, D’Orta G, Mladenov N, Winget DM, Suttle CA. Deposition rates of viruses and bacteria above the atmospheric boundary layer. ISME J. 2018;12:1154–62.CAS 
PubMed 
PubMed Central 

Google Scholar 
47.Lehmann K, Lehmann R, Totsche KU. Event-driven dynamics of the total mobile inventory in undisturbed soil account for significant fluxes of particulate organic carbon. Sci Total Environ. 2021;756:143774.CAS 
PubMed 

Google Scholar 
48.Storms ZJ, Sauvageau D. Modeling tailed bacteriophage adsorption: insight into mechanisms. Virology. 2015;485:355–62.CAS 
PubMed 

Google Scholar 
49.Shan Y, Harms H, Wick LY. Electric field effects on bacterial deposition and transport in porous media. Environ Sci Technol. 2018;52:14294–301.CAS 
PubMed 

Google Scholar 
50.Junier P, Cailleau G, Palmieri I, Vallotton C, Trautschold OC, Junier T, et al. Democratization of fungal highway columns as a tool to investigate bacteria associated with soil fungi. FEMS Microbiol Ecol. 2021;97:fiab003.CAS 
PubMed 
PubMed Central 

Google Scholar 
51.Furuno S, Remer R, Chatzinotas A, Harms H, Wick LY. Use of mycelia as paths for the isolation of contaminant-degrading bacteria from soil. Micro Biotechnol. 2012;5:142–8.CAS 

Google Scholar 
52.Jiang F, Zhang L, Zhou J, George TS, Feng G. Arbuscular mycorrhizal fungi enhance mineralisation of organic phosphorus by carrying bacteria along their extraradical hyphae. N. Phytol. 2021;230:304–15.CAS 

Google Scholar 
53.Jansa J, Hodge A. Swimming, gliding, or hyphal riding? On microbial migration along the arbuscular mycorrhizal hyphal highway and functional consequences thereof. N. Phytol. 2021;230:14–16.
Google Scholar 
54.Zhang Y, Kastman EK, Guasto JS, Wolfe BE. Fungal networks shape dynamics of bacterial dispersal and community assembly in cheese rind microbiomes. Nat Commun. 2018;9:336.PubMed 
PubMed Central 

Google Scholar 
55.Ping D, Wang T, Fraebel DT, Maslov S, Sneppen K, Kuehn S. Hitchhiking, collapse, and contingency in phage infections of migrating bacterial populations. ISME J 2020;14:2007–18.PubMed 
PubMed Central 

Google Scholar 
56.Testa S, Berger S, Piccardi P, Oechslin F, Resch G, Mitri S. Spatial structure affects phage efficacy in infecting dual-strain biofilms of Pseudomonas aeruginosa. Commun Biol. 2019;2:405.PubMed 
PubMed Central 

Google Scholar 
57.May T, Tsuruta K, Okabe S. Exposure of conjugative plasmid carrying Escherichia coli biofilms to male-specific bacteriophages. ISME J. 2011;5:771–5.CAS 
PubMed 

Google Scholar 
58.Abedon ST. Phage “delay” towards enhancing bacterial escape from biofilms: a more comprehensive way of viewing resistance to bacteriophages. AIMS Microbiol. 2017;3:186.PubMed 
PubMed Central 

Google Scholar 
59.Adams MH Bacteriophages (Interscience Publishers, Inc., New York – London, 1959)60.Schrader HS, Schrader JO, Walker JJ, Bruggeman NB, Vanderloop JM, Shaffer JJ, et al. Effects of host starvation on bacteriophage dynamics. Bact Oligotrophic Environ Starvation-Survival Lifestyle. 1997; 368-85.61.Schrader HS, Schrader JO, Walker JJ, Wolf TA, Nickerson KW, Kokjohn TA. Bacteriophage infection and multiplication occur in Pseudomonas aeruginosa starved for 5 years. Can J Microbiol. 1997;43:1157–63.CAS 
PubMed 

Google Scholar 
62.Łoś M, Golec P, Łoś JM, Węglewska-Jurkiewicz A, Czyż A, Węgrzyn A, et al. Effective inhibition of lytic development of bacteriophages λ, P1 and T4 by starvation of their host, Escherichia coli. BMC Biotechnol. 2007;7:13.PubMed 
PubMed Central 

Google Scholar 
63.Bryan D, El-Shibiny A, Hobbs Z, Porter J, Kutter EM. Bacteriophage T4 infection of stationary phase E. coli: life after log from a phage perspective. Front Microbiol. 2016;7:1391.PubMed 
PubMed Central 

Google Scholar 
64.Yin J. A quantifiable phenotype of viral propagation. Biochem Biophys Res Commun. 1991;174:1009–14.CAS 
PubMed 

Google Scholar 
65.Chatterjee A, Willett JLE, Dunny GM, Duerkop BA. Phage infection and sub-lethal antibiotic exposure mediate Enterococcus faecalis type VII secretion system dependent inhibition of bystander bacteria. PLOS Genet. 2021;17:e1009204.CAS 
PubMed 
PubMed Central 

Google Scholar 
66.Berthold T, Centler F, Hübschmann T, Remer R, Thullner M, Harms H, et al. Mycelia as a focal point for horizontal gene transfer among soil bacteria. Sci Rep. 2016;6:36390.CAS 
PubMed 
PubMed Central 

Google Scholar 
67.Lee KL, Hubbard LC, Hern S, Yildiz I, Gratzl M, Steinmetz NF. Shape matters: the diffusion rates of TMV rods and CPMV icosahedrons in a spheroid model of extracellular matrix are distinct. Biomater Sci. 2013;1. https://doi.org/10.1039/C3BM00191A.68.Hudson P, Greenman J. Competition mediated by parasites: biological and theoretical progress. Trends Ecol Evol. 1998;13:387–90.CAS 
PubMed 

Google Scholar 
69.Sax DF, Stachowicz JJ, Brown JH, Bruno JF, Dawson MN, Gaines SD, et al. Ecological and evolutionary insights from species invasions. Trends Ecol Evol. 2007;22:465–71.PubMed 

Google Scholar 
70.Wagner PL, Waldor MK. Bacteriophage control of bacterial virulence. Infect Immun. 2020;70:3985–93.
Google Scholar 
71.Chantrey J, Dale TD, Read JM, White S, Whitfield F, Jones D, et al. European red squirrel population dynamics driven by squirrelpox at a gray squirrel invasion interface. Ecol Evol. 2014;4:3788–99.PubMed 
PubMed Central 

Google Scholar 
72.Essl F, Bacher S, Genovesi P, Hulme PE, Jeschke JM, Katsanevakis S, et al. Which taxa are alien? Criteria, applications, and uncertainties. Bioscience 2018;68:496–509.
Google Scholar 
73.Seebens H, Gastner MT, Blasius B. The risk of marine bioinvasion caused by global shipping. Ecol Lett. 2013;16:782–90.CAS 
PubMed 

Google Scholar 
74.Seebens H, Essl F, Blasius B. The intermediate distance hypothesis of biological invasions. Ecol Lett. 2017;20:158–65.PubMed 

Google Scholar 
75.Hulme PE, Bacher S, Kenis M, Klotz S, Kühn I, Minchin D, et al. Grasping at the routes of biological invasions: a framework for integrating pathways into policy. J Appl Ecol. 2008;45:403–14.
Google Scholar 
76.Liebhold AM, Brockerhoff EG, Garrett LJ, Parke JL, Britton KO. Live plant imports: the major pathway for forest insect and pathogen invasions of the US. Front Ecol Environ. 2012;10:135–43.
Google Scholar  More

1.Wall, D. H., Bardgett, R. D. & Kelly, E. Biodiversity in the dark. Nat. Geosci. 3(5), 297–298 (2010).ADS 
CAS 

Google Scholar 
2.Eisenhauer, N., Bonn, A. & Guerra, C. A. Recognizing the quiet extinction of invertebrates. Nat. Commun. 10(1), 1–3 (2019).
Google Scholar 
3.Koch, A. et al. Soil security: Solving the global soil crisis. Global Pol. 4(4), 434–441 (2013).
Google Scholar 
4.Wall, D. H., Nielsen, U. N. & Six, J. Soil biodiversity and human health. Nature 528(7580), 69–76 (2015).ADS 
CAS 
PubMed 

Google Scholar 
5.Guerra, C. A. et al. Blind spots in global soil biodiversity and ecosystem function research. Nat. Commun. 11(1), 1–13 (2020).
Google Scholar 
6.Bardgett, R. D. & van der Putten, W. H. Belowground biodiversity and ecosystem functioning. Nature 515, 505–511 (2014).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
7.Zou, K., Thébault, E., Lacroix, G. & Barot, S. Interactions between the green and brown food web determine ecosystem functioning. Funct. Ecol. 30(8), 1454–1465 (2016).
Google Scholar 
8.Lavelle, P. et al. Soil invertebrates and ecosystem services. Eur. J. Soil Biol. 42, S3–S15 (2006).
Google Scholar 
9.de Vries, F. T. et al. Soil food web properties explain ecosystem services across European land use systems. Proc. Natl. Acad. Sci. 110(35), 14296–14301 (2013).ADS 
PubMed 
PubMed Central 

Google Scholar 
10.Adhikari, K. & Hartemink, A. E. Linking soils to ecosystem services—A global review. Geoderma 262, 101–111 (2016).ADS 
CAS 

Google Scholar 
11.Cameron, E. K. et al. Global mismatches in aboveground and belowground biodiversity. Conserv. Biol. 33(5), 1187–1192 (2019).PubMed 

Google Scholar 
12.Phillips, H. R., Heintz-Buschart, A. & Eisenhauer, N. Putting soil invertebrate diversity on the map. Mol. Ecol. 29(4), 655–657 (2020).PubMed 
PubMed Central 

Google Scholar 
13.El Mujtar, V., Muñoz, N., Mc Cormick, B. P., Pulleman, M. & Tittonell, P. Role and management of soil biodiversity for food security and nutrition; where do we stand?. Glob. Food Sec. 20, 132–144 (2019).
Google Scholar 
14.Schuldt, A. et al. Biodiversity across trophic levels drives multifunctionality in highly diverse forests. Nat. Commun. 9(1), 2989 (2018).ADS 
PubMed 
PubMed Central 

Google Scholar 
15.Eisenhauer, N. et al. Priorities for research in soil ecology. Pedobiologia 63, 1–7 (2017).PubMed 
PubMed Central 

Google Scholar 
16.Brose, U. & Scheu, S. Into darkness: Unravelling the structure of soil food webs. Oikos 123(10), 1153–1156 (2014).
Google Scholar 
17.Phillips, H. R. et al. Red list of a black box. Nat. Ecol. Evol. 1(4), 1–1 (2017).
Google Scholar 
18.Hairston, N. G., Smith, F. E. & Slobodkin, L. B. Community structure, population control, and competition. Am. Nat. 94(879), 421–425 (1960).
Google Scholar 
19.Vidal, M. C. & Murphy, S. M. Bottom-up vs top-down effects on terrestrial insect herbivores: A meta-analysis. Ecol. Lett. 21(1), 138–150 (2018).PubMed 

Google Scholar 
20.Wagg, C., Bender, S. F., Widmer, F. & van der Heijden, M. G. Soil biodiversity and soil community composition determine ecosystem multifunctionality. Proc. Natl. Acad. Sci. 111(14), 5266–5270 (2014).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
21.Soliveres, S. et al. Biodiversity at multiple trophic levels is needed for ecosystem multifunctionality. Nature 536(7617), 456–459 (2016).ADS 
CAS 
PubMed 

Google Scholar 
22.Holling, C. S. Resilience and stability of ecological systems. Annu. Rev. Ecol. Syst. 4(1), 1–23 (1973).
Google Scholar 
23.Allesina, S. & Tang, S. Stability criteria for complex ecosystems. Nature 483(7388), 205–208 (2012).ADS 
CAS 
PubMed 

Google Scholar 
24.Crowther, T. W. et al. Biotic interactions mediate soil microbial feedbacks to climate change. Proc. Natl. Acad. Sci. 112(22), 7033–7038 (2015).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
25.Maran, A. M. & Pelini, S. L. Predator contributions to belowground responses to warming. Ecosphere 7(9), e01457 (2016).
Google Scholar 
26.Geisen, S., Wall, D. H. & van der Putten, W. H. Challenges and opportunities for soil biodiversity in the Anthropocene. Curr. Biol. 29(19), R1036–R1044 (2019).CAS 
PubMed 

Google Scholar 
27.Rooney, N., McCann, K., Gellner, G. & Moore, J. C. Structural asymmetry and the stability of diverse food webs. Nature 442(7100), 265–269 (2006).ADS 
CAS 
PubMed 

Google Scholar 
28.Murphy, S. M., Lewis, D. & Wimp, G. M. Predator population size structure alters consumption of prey from epigeic and grazing food webs. Oecologia 192(3), 791–799 (2020).ADS 
PubMed 

Google Scholar 
29.Scheu, S. Plants and generalist predators as links between the below-ground and above-ground system. Basic Appl. Ecol. 2, 3–13 (2001).
Google Scholar 
30.Wardle, D. A. et al. Ecological linkages between aboveground and belowground biota. Science 304(5677), 1629–1633 (2004).ADS 
CAS 
PubMed 

Google Scholar 
31.de Vries, F. T. & Wallenstein, M. D. Below-ground connections underlying above-ground food production: A framework for optimising ecological connections in the rhizosphere. J. Ecol. 105(4), 913–920 (2017).
Google Scholar 
32.Wu, T., Ayres, E., Bardgett, R. D., Wall, D. H. & Garey, J. R. Molecular study of worldwide distribution and diversity of soil animals. Proc. Natl. Acad. Sci. 108(43), 17720–17725 (2011).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
33.Symondson, W. O. C., Sunderland, K. D. & Greenstone, M. H. Can generalist predators be effective biocontrol agents?. Annu. Rev. Entomol. 47(1), 561–594 (2002).CAS 
PubMed 

Google Scholar 
34.Dainese, M. et al. A global synthesis reveals biodiversity-mediated benefits for crop production. Sci. Adv. 5(10), eaax0121 (2019).ADS 
PubMed 
PubMed Central 

Google Scholar 
35.Karp, D. S. et al. Crop pests and predators exhibit inconsistent responses to surrounding landscape composition. Proc. Natl. Acad. Sci. 115(33), E7863–E7870 (2018).CAS 
PubMed 
PubMed Central 

Google Scholar 
36.Johnson, S. N. et al. New frontiers in belowground ecology for plant protection from root-feeding insects. Appl. Soil. Ecol. 108, 96–107 (2016).
Google Scholar 
37.Veen, C. et al. Applying the aboveground-belowground interaction concept in agriculture: Spatio-temporal scales matter. Front. Ecol. Evol. 7, 300 (2019).
Google Scholar 
38.Birkhofer, K., Wise, D. H. & Scheu, S. Subsidy from the detrital food web, but not microhabitat complexity, affects the role of generalist predators in an aboveground herbivore food web. Oikos 117(4), 494–500 (2008).
Google Scholar 
39.Birkhofer, K. et al. Organic farming affects the biological control of hemipteran pests and yields in spring barley independent of landscape complexity. Landsc. Ecol. 31(3), 567–579 (2016).
Google Scholar 
40.van der Putten, W. H. et al. Empirical and theoretical challenges in aboveground–belowground ecology. Oecologia 161(1), 1–14 (2009).ADS 
PubMed 
PubMed Central 

Google Scholar 
41.Kleijn, D. et al. Ecological intensification: Bridging the gap between science and practice. Trends Ecol. Evol. 34(2), 154–166 (2019).PubMed 

Google Scholar 
42.Bender, S. F., Wagg, C. & van der Heijden, M. G. An underground revolution: Biodiversity and soil ecological engineering for agricultural sustainability. Trends Ecol. Evol. 31(6), 440–452 (2016).PubMed 

Google Scholar 
43.Gagic, V. et al. Combined effects of agrochemicals and ecosystem services on crop yield across Europe. Ecol. Lett. 20(11), 1427–1436 (2017).PubMed 

Google Scholar 
44.Briones, M. J. The serendipitous value of soil fauna in ecosystem functioning: The unexplained explained. Front. Environ. Sci. 6, 149 (2018).
Google Scholar 
45.Kaya, H. K. & Gaugler, R. Entomopathogenic nematodes. Annu. Rev. Entomol. 38(1), 181–206 (1993).
Google Scholar 
46.Ferris, H. & Tuomisto, H. Unearthing the role of biological diversity in soil health. Soil Biol. Biochem. 85, 101–109 (2015).CAS 

Google Scholar 
47.Tsiafouli, M. A. et al. Intensive agriculture reduces soil biodiversity across Europe. Glob. Change Biol. 21, 973–985 (2015).ADS 

Google Scholar 
48.Bender, S. F. & van der Heijden, M. G. Soil biota enhance agricultural sustainability by improving crop yield, nutrient uptake and reducing nitrogen leaching losses. J. Appl. Ecol. 52(1), 228–239 (2015).CAS 

Google Scholar 
49.De Vries, F. T. et al. Land use alters the resistance and resilience of soil food webs to drought. Nat. Clim. Change 2, 276–280 (2012).ADS 

Google Scholar 
50.Bastida, F. et al. Climatic vulnerabilities and ecological preferences of soil invertebrates across biomes. Mol. Ecol. 29(4), 752–761 (2020).PubMed 

Google Scholar 
51.Pereira, H. M., Navarro, L. M. & Martins, I. S. Global biodiversity change: The bad, the good, and the unknown. Annu. Rev. Environ. Resour. 37, 25–50 (2012).
Google Scholar 
52.Polis, G. A. Complex trophic interactions in deserts: An empirical critique of food-web theory. Am. Nat. 138(1), 123–155 (1991).
Google Scholar 
53.Polis, G. A. & Strong, D. R. Food web complexity and community dynamics. Am. Nat. 147(5), 813–846 (1996).
Google Scholar 
54.Lavelle, P. et al. Ecosystem engineers in a self-organized soil: A review of concepts and future research questions. Soil Sci. 181(3/4), 91–109 (2016).ADS 
CAS 

Google Scholar 
55.Nielsen, U. N. et al. The enigma of soil animal species diversity revisited: The role of small-scale heterogeneity. PLoS ONE 5(7), e11567 (2010).ADS 
PubMed 
PubMed Central 

Google Scholar 
56.Heinen, R., van der Sluijs, M., Biere, A., Harvey, J. A. & Bezemer, T. M. Plant community composition but not plant traits determine the outcome of soil legacy effects on plants and insects. J. Ecol. 106(3), 1217–1229 (2018).
Google Scholar 
57.Ramirez, K. S., Geisen, S., Morriën, E., Snoek, B. L. & van der Putten, W. H. Network analyses can advance above-belowground ecology. Trends Plant Sci. 23(9), 759–768 (2018).CAS 
PubMed 

Google Scholar 
58.Boyer, S., Snyder, W. E. & Wratten, S. D. Molecular and isotopic approaches to food webs in agroecosystems. Food Webs 9, 1–3 (2016).
Google Scholar 
59.Casey, J. M. et al. Reconstructing hyperdiverse food webs: Gut content metabarcoding as a tool to disentangle trophic interactions on coral reefs. Methods Ecol. Evol. 10(8), 1157–1170 (2019).
Google Scholar 
60.Choate, B. A. & Lundgren, J. G. Invertebrate communities in spring wheat and the identification of cereal aphid predators through molecular gut content analysis. Crop Prot. 77, 110–118 (2015).
Google Scholar 
61.Furlong, M. J. Knowing your enemies: Integrating molecular and ecological methods to assess the impact of arthropod predators on crop pests. Insect Sci. 22(1), 6–19 (2015).PubMed 

Google Scholar 
62.Eitzinger, B., Rall, B. C., Traugott, M. & Scheu, S. Testing the validity of functional response models using molecular gut content analysis for prey choice in soil predators. Oikos 127(7), 915–926 (2018).
Google Scholar 
63.Barberán, A., Bates, S. T., Casamayor, E. O. & Fierer, N. Using network analysis to explore co-occurrence patterns in soil microbial communities. ISME J. 6(2), 343–351 (2012).PubMed 

Google Scholar 
64.Morriën, E. Understanding soil food web dynamics, how close do we get?. Soil Biol. Biochem. 102, 10–13 (2016).
Google Scholar 
65.Digel, C., Curtsdotter, A., Riede, J., Klarner, B. & Brose, U. Unravelling the complex structure of forest soil food webs: Higher omnivory and more trophic levels. Oikos 123(10), 1157–1172 (2014).
Google Scholar 
66.Toscano, B. J., Hin, V. & Rudolf, V. H. Cannibalism and intraguild predation community dynamics: Coexistence, competitive exclusion, and the loss of alternative stable states. Am. Nat. 190(5), 617–630 (2017).PubMed 

Google Scholar 
67.Coleman, D. C. & Wall, D. H. Soil fauna: Occurrence, biodiversity, and roles in ecosystem function. Soil Microbiol. Ecol. Biochem. 4, 111–149 (2015).
Google Scholar 
68.Brussaard, L. Biodiversity and ecosystem functioning in soil. Ambio 26, 563–570 (1997).
Google Scholar 
69.Briar, S. S. et al. The distribution of nematodes and soil microbial communities across soil aggregate fractions and farm management systems. Soil Biol. Biochem. 43, 905–914 (2011).CAS 

Google Scholar 
70.Oelbermann, K. & Scheu, S. Trophic guilds of generalist feeders in soil animal communities as indicated by stable isotope analysis (15N/14N). Bull. Entomol. Res. 100(5), 511 (2010).CAS 
PubMed 

Google Scholar 
71.Cohen, J. E., Pimm, S. L., Yodzis, P. & Saldaña, J. Body sizes of animal predators and animal prey in food webs. J. Anim. Ecol. 62, 67–78 (1993).
Google Scholar 
72.Nielsen, U. N., Wall, D. H. & Six, J. Soil biodiversity and the environment. Annu. Rev. Environ. Resour. 40, 63–90 (2015).
Google Scholar 
73.Veresoglou, S. D., Halley, J. M. & Rillig, M. C. Extinction risk of soil biota. Nat. Commun. 6(1), 1–10 (2015).
Google Scholar 
74.Ruf, A. A maturity index for predatory soil mites (Mesostigmata: Gamasina) as an indicator of environmental impacts of pollution on forest soils. Appl. Soil. Ecol. 9(1–3), 447–452 (1998).
Google Scholar 
75.Zak, D. R., Holmes, W. E., White, D. C., Peacock, A. D. & Tilman, D. Plant diversity, soil microbial communities, and ecosystem function: Are there any links?. Ecology 84(8), 2042–2050 (2003).
Google Scholar 
76.Leach, J. E., Triplett, L. R., Argueso, C. T. & Trivedi, P. Communication in the phytobiome. Cell 169(4), 587–596 (2017).CAS 
PubMed 

Google Scholar 
77.Barnes, A. D. et al. Energy flux: The link between multitrophic biodiversity and ecosystem functioning. Trends Ecol. Evol. 33(3), 186–197 (2018).PubMed 
PubMed Central 

Google Scholar 
78.Heinen, R., Biere, A., Harvey, J. A. & Bezemer, T. M. Effects of soil organisms on aboveground plant-insect interactions in the field: Patterns, mechanisms and the role of methodology. Front. Ecol. Evol. 6, 106 (2018).
Google Scholar 
79.Rillig, M. C. et al. The role of multiple global change factors in driving soil functions and microbial biodiversity. Science 366(6467), 886–890 (2019).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
80.Wardle, D. A., Hyodo, F., Bardgett, R. D., Yeates, G. W. & Nilsson, M. C. Long-term aboveground and belowground consequences of red wood ant exclusion in boreal forest. Ecology 92(3), 645–656 (2011).PubMed 

Google Scholar 
81.Preisser, E. L. & Strong, D. R. Climate affects predator control of an herbivore outbreak. Am. Nat. 163(5), 754–762 (2004).PubMed 

Google Scholar 
82.Hamilton, J. et al. Elevated atmospheric CO2 alters the arthropod community in a forest understory. Acta Oecol. 43, 80–85 (2012).ADS 

Google Scholar 
83.Zaller, J. G. et al. Future rainfall variations reduce abundances of aboveground arthropods in model agroecosystems with different soil types. Front. Environ. Sci. 2, 44 (2014).
Google Scholar 
84.Koltz, A. M., Classen, A. T. & Wright, J. P. Warming reverses top-down effects of predators on belowground ecosystem function in Arctic tundra. Proc. Natl. Acad. Sci. 115(32), E7541–E7549 (2018).CAS 
PubMed 
PubMed Central 

Google Scholar 
85.Santonja, M. et al. Plant litter mixture partly mitigates the negative effects of extended drought on soil biota and litter decomposition in a Mediterranean oak forest. J. Ecol. 105(3), 801–815 (2017).
Google Scholar 
86.Garratt, M. P. et al. Enhancing soil organic matter as a route to the ecological intensification of European arable systems. Ecosystems 21(7), 1404–1415 (2018).CAS 

Google Scholar 
87.Smith-Ramesh, L. M. Predators in the plant–soil feedback loop: Aboveground plant-associated predators may alter the outcome of plant–soil interactions. Ecol. Lett. 21(5), 646–654 (2018).PubMed 

Google Scholar 
88.Gurr, G. M., Wratten, S. D., Landis, D. A. & You, M. Habitat management to suppress pest populations: Progress and prospects. Annu. Rev. Entomol. 62, 91–109 (2017).CAS 
PubMed 

Google Scholar 
89.Rypstra, A. L., Carter, P. E., Balfour, R. A. & Marshall, S. D. Architectural features of agricultural habitats and their impact on the spider inhabitants. J. Arachnol. 27, 371–377 (1999).
Google Scholar 
90.Von Berg, K., Thies, C., Tscharntke, T. & Scheu, S. Changes in herbivore control in arable fields by detrital subsidies depend on predator species and vary in space. Oecologia 163(4), 1033–1042 (2010).ADS 

Google Scholar 
91.Rowen, E., Tooker, J. F. & Blubaugh, C. K. Managing fertility with animal waste to promote arthropod pest suppression. Biol. Control 134, 130–140 (2019).
Google Scholar 
92.Perović, D. J. et al. Managing biological control services through multi-trophic trait interactions: Review and guidelines for implementation at local and landscape scales. Biol. Rev. 93(1), 306–321 (2018).PubMed 

Google Scholar 
93.Roger-Estrade, J., Anger, C., Bertrand, M. & Richard, G. Tillage and soil ecology: Partners for sustainable agriculture. Soil Tillage Res. 111(1), 33–40 (2010).
Google Scholar 
94.Dias, T., Dukes, A. & Antunes, P. M. Accounting for soil biotic effects on soil health and crop productivity in the design of crop rotations. J. Sci. Food Agric. 95(3), 447–454 (2015).CAS 
PubMed 

Google Scholar 
95.Tamburini, G., De Simone, S., Sigura, M., Boscutti, F. & Marini, L. Conservation tillage mitigates the negative effect of landscape simplification on biological control. J. Appl. Ecol. 53(1), 233–241 (2016).
Google Scholar 
96.Pretty, J. et al. Global assessment of agricultural system redesign for sustainable intensification. Nat. Sustain. 1(8), 441–446 (2018).
Google Scholar 
97.Swift, M. J., Heal, O. W., Anderson, J. M. & Anderson, J. M. Decomposition in Terrestrial Ecosystems Vol. 5 (University of California Press, 1979).
Google Scholar 
98.van Straalen, N. M., Butovsky, R. O., Pokarzhevskii, A. D., Zaitsev, A. S. & Verhoef, S. C. Metal concentrations in soil and invertebrates in the vicinity of a metallurgical factory near Tula (Russia). Pedobiologia 45(5), 451–466 (2001).
Google Scholar 
99.Birkhofer, K. et al. Methods to identify the prey of invertebrate predators in terrestrial field studies. Ecol. Evol. 7(6), 1942–1953 (2017).PubMed 
PubMed Central 

Google Scholar 
100.Potapov, A. M., Tiunov, A. V. & Scheu, S. Uncovering trophic positions and food resources of soil animals using bulk natural stable isotope composition. Biol. Rev. 94(1), 37–59 (2019).
Google Scholar  More

Study site and populationThe study was conducted between July and October 2019 in the North of Peru. We focused on populations of coexisting Morpho species present in the regional park of the Cordillera Escalera (San Martin Department) near the city of Tarapoto. Both the capture-recapture and the dummy experiment were performed at the exact same location, on the bank of the Shilcayo river (06°27′14.364″S, 76°20′45.852″W).DNA extraction and RAD-SequencingThirty-one wild males caught on the study site were sequenced to perform population genomic analyses (M. achilles—n = 13, M. helenor—n = 10 and M. deidamia—n = 8). DNA was extracted from each sample from a slice of the thorax, using Qiagen kit DNeasy Blood & Tissue. DNA quantification (using the microfluorimetric method) and quality controls (using electrophoresis and spectrophotometric method) were performed prior to sequencing. RAD-library preparation and sequencing were performed at the MGX-Montpellier GenomiX platform (Montpellier, France). DNA was digested with the Pst1 enzyme and the library was prepared according to Baird and Etter’ protocol47 in a slightly modified version. Paired-end RAD-sequencing was performed on a 2 lanes flow cell of an Illumina HiSeq2500 in a rapid mode so that reads (125 bp) were expected to be of high quality with no missing base (N content). We obtained 299 million sequences, comprising R1 and R2 reads for each sequenced fragment. Adapters were removed from the reads.Read quality control, alignment and dataset generationRead quality was assessed with FastQC v0.11.9 (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). The per base sequence quality was high across all reads (no lower than 36 for R1 and 32 for R2) with an average quality score of 39 (40 being the maximum). Overall, FastQC highlighted the high quality of the sequencing data, allowing us to skip the step of read trimming.The data were demultiplexed, assigning each sequence to its sample ID and the reads were aligned using Stacks V2.5 (http://catchenlab.life.illinois.edu/stacks/)48,49. Parameters were set following the 80% polymorphic (r80) loci rule, which only considers loci shared by at least 80% of the samples50. The optimised parameters are ‘max distance between stacks’ (inside each sample) and ‘number of mismatches between stacks’ (between samples). Every other parameter was kept to default values. After aligning all reads, we selected 2740 biallelic loci shared by all samples, including 88,889 SNPs in total. Each locus had a length of 463.12 bp on average (range [343; 908]). These loci are assumed to be evenly distributed throughout the genome but cover only a limited portion of the genome (around 0.5%). Datasets were stored in a VCF file (containing all the SNPs found in the alignment) and a fasta file (containing the two alleles found at every locus for each sample). To run DILS-ABC inferences, Stacks fasta files were converted to another fasta format compatible with DILS (https://github.com/CoBiG2/RAD_Tools).Demographic inferencesEight categories of demographic models were compared, according to temporal patterns of introgression. This was done to answer two questions on gene flow in Morpho: (1) is there ongoing migration between M. helenor and M. achilles? (2) do M. helenor and/or M. achilles exchange alleles with M. deidamia? This was assessed by an ABC approach using a version of DILS adapted to samples of three populations/species32. Since Stacks does not report monomorphic RAD loci, the ABC analysis was conditioned in the same way, by excluding monomorphic loci from the simulations. Focusing on polymorphic loci may only limit our ability to estimate the absolute values of parameters (i.e. population sizes expressed in numbers of individuals, and ages of past events expressed in numbers of generations); nevertheless, this framework excluding monomorphic loci still allows reliable comparisons of models51 and estimations of relative parameter values, as performed to investigate the human history51.A generalist model was studied (Supplementary Fig. 12). This model describes an ancestral population subdivided in two populations: the ancestor of M. deidamia and the common ancestor of M. helenor/M. achilles. The latter population was further subdivided into the three species/populations currently sampled. Each split event is accompanied by a change in demographic size, the value of which is independent of the ancestral size. In addition, given clear genomic signatures for recent demographic changes with largely negative Tajima’s D, we implemented variations for the effective sizes of the three modern lineages at independent times. Finally, migration can occur between each pair of species/populations. Migration affecting the M. helenor/M. achilles pair can either be the result of secondary contact after a period of isolation (ongoing migration), or of ancestral migration (current isolation) as in50,52.As this model is over-parameterised, our general strategy is to investigate the above two questions by comparing variations of this generalist model. Thus, to test the gene flow between M. helenor and M. achilles, we compared two categories of models. (1) With random parameter values for all model parameters including the ongoing migration between M. helenor and M. achilles (gene flow resulting from a secondary contact between them); (2) as above, but with the migration between M. helenor and M. achilles set to zero after a randomly drawn number of generations following their split. An overlap between ‘current isolation’ and ‘ongoing migration’ models can occur when the transition time (from ancestral migration to current isolation forward in time for a ‘current isolation’ model; or from ancestral isolation to ongoing migration forward in time for an ‘ongoing migration’ model) tends towards the extreme values 0 or Tsplit hel-ach (Supplementary Fig. 12). To reduce this effect, the transition times were drawn in a Beta distribution with parameters (α = 5, β = 1) when migration has to be restricted to a past period, and in a Beta distribution with parameters (α = 1, β = 5) when migration is assumed to occur after a recent secondary contact.When two broad categories of models are statistically compared, each category is represented by simulations performed under the four sub-models allowing or not allowing genomic heterogeneities for effective sizes (Ne) and for migration rates (N.m). For instance, to test for gene flow between M. helenor and M. achilles, the model of ‘ongoing migration’ is actually represented by simulations with the four possible combinations of homogeneity/heterogeneity, all labelled as being ‘ongoing migration’.As for any inferential analysis, it is important to recognise that the best-supported model is based on a classification of models within a studied set. Intermediate models, with more subtle cycles of genetic isolation and secondary contact could produce a better fit to the data, but it would be surprising to detect a strong support for the model assuming a lack of recent gene flow, if the most recent secondary contact of such cyclicity induced elevated gene flow.For each model, 50,000 simulations using random combinations of parameters were performed. Parameters were drawn from uniform prior distributions. Population sizes were sampled from the uniform prior [0–1,000,000] (in diploid individuals); the older time of split was sampled from the uniform prior [0–8,000,000] (generations); ages of the subsequent demographic events were sampled in a uniform prior between 0 and the sampled time of split. Migration rates 4.N.m were sampled from the uniform prior [0–50]. Both migration rates and effective population sizes are allowed to vary throughout the genomes as a result of linked selection, following refs. 53,54,55.On each simulated dataset, we calculated a vector of means and standard deviations for different summary statistics: intraspecific statistics (π for M. helenor, π for M. achilles, π for M. deidamia, θW for M. helenor, θW for M. achilles, θW for M. deidamia, Tajima’s D for M. helenor, Tajima’s D for M. achilles, Tajima’s D for M. deidamia) and interspecific statistics (gross divergence, net divergence and FST for all three possible pairs; ABBA-BABA D). Our version of DILS includes part of the DaDi56 and Moments57 strategy involving the identification of the best model proposed demographic model from the molecular patterns of polymorphism and divergence (proportion of shared polymorphisms, fixed differences between species, exclusive polymorphisms, etc.), excluding monomorphic loci. Thus, only loci containing at least one SNP in an alignment of the three species studied are considered, including singletons. Importantly, each locus carrying at least one SNP in a tri-specific alignment is associated with a mutation rate assumed to be 3 · 10−9 mutations per generation and per base pair to convert demographic parameters into demographic units from coalescence units.We first conditioned the mutations occurring during coalescent simulations by using theta (=4 · N · µ · Li; where N is the effective population size, µ the mutation rate per nucleotide and per generation; Li the length of locus i). The number of simulated segregating sites for a given locus strongly depends on the coalescent history (i.e the total length of the simulated coalescent tree), occasionally generating monomorphic loci. To confirm that the inferences are not impacted by differences in the number of monomorphic loci in the simulated datasets, we then used an alternative simulation approach, by randomly placing in simulated coalescent trees a fixed number of mutations corresponding to the observed number of SNPs for each locus. Thus, a randomly simulated dataset consists of 2740 loci whose lengths (ranging from 339 to 894 nucleotides) and number of SNPs (ranging from 1 to 91) individually match the properties of the observed loci in the actual dataset. Since the results drawn from both approaches were similar, we report only the estimations provided by the simulations based on the actual number of SNPs. Comparisons between the two approaches can be found in supplementary (Supplementary Tables 8, 9).Statistical comparisons between simulated and observed statistics were performed using the R package abcrf version 1.8.158,59.Mark-recapture experimentTo estimate the timing of patrolling activity among Morpho species, we performed capture-mark-recapture between 9 a.m. and 2 p.m. (flight activity in Morpho is drastically reduced in the afternoons at this site) during 17 sunny days. Although on a few days, capture was cancelled because of bad weather annihilating butterfly activity, the 17 capture sessions were mostly consecutives, as they were performed in a 22 days period (Supplementary Table 1 and Supplementary Fig. 15). All butterflies were captured with hand-nets, identified at the species level, and numbered on their dorsal wing surface using a black marker. The exact time of each capture was annotated. Butterflies captured while inactive, such as those laying on a branch or on the ground were excluded from the analysis to focus exclusively on actively patrolling individuals. We measured patrolling time for a total of 295 occasions, including 78 recaptures (i.e. 217 individuals were captured at least once). All captured individuals were males. Individuals M. achilles were the most frequently captured (n = 121), followed by M. helenor (n = 95). Individuals M. deidamia were about half less captured (n = 48), and individual M. menelaus were the least captured (n = 34). Because striking differences in patrolling time were observed among Morpho species, we used time of the day as a predictor of species identity in order to distinguish between M. helenor and M. achilles in the below-described experiment because butterflies from these two species are morphologically too similar to be identified while flying (Supplementary Fig. 13). After the 17 nearly-consecutive days of capture, one day of capture was repeated every 2 weeks during 2 months in parallel to the dummy experiment (described below), to verify that temporal activity was stable over time (Supplementary Fig. 13).Estimating population size from mark-recapture dataBased on capture-recapture histories, we estimated individual abundance for each species using a loglinear model implemented in the R package Rcapture version 1.4.360 (Supplementary Fig. 15). Given the short duration the sampling period (22 days) relative to the longevity of adult Morpho butterflies (several months61), we used a closed-population model assuming no effect of births, deaths, immigration and emigration. Abundance was estimated in Morpho helenor and M. achilles only, as capture and recapture events were too few in the other species (M. deidamia and M. menelaus) to allow estimating population size (Supplementary Table 1).Experiment with dummy butterfliesWe investigated the response of patrolling males to sympatric conspecifics, congeners and of exotic conspecifics, using dummies placed on their flight path. Dummies were built with real wings dissected and washed with hexane to remove volatile compounds and cuticular hydrocarbons, ensuring to test only the visual aspect of the dummies. We mounted the wings on a solar-powered fluttering device (Butterfly Solar Héliobil R029br) that mimics a flying butterfly, thereby increasing the attractiveness of the dummy. The fluttering dummy was positioned on the riverbank, and placed at the centre of a 1 m3 space delimitated with four vertical stacks (Fig. 1a). The set-up was continuously monitored by a human observer and filmed using a camera (Gopro Hero5 Black set at 120 images per second) mounted on a tripod. Patrolling Morpho butterflies that deviated from their flight path to approach the dummy but did not enter the cubic space were categorised as approaching. Any Morpho butterfly entering the cubic space was considered as interacting with the dummy. Those passing without showing interest to the setup were categorised as passing. The category of behaviour and the exact time of the butterfly responses were annotated on site by the human observer. Patrolling individuals were mainly identified at the species level by the observer on the site: M. menelaus can be easily distinguished from M. deidamia, and these two species are also quite different from M. helenor and M. achilles. However, the sister species M. helenor and M. achilles cannot be discriminated during flight, and we thus rely on an indirect method, based on flight hours, to infer the species identity of wild visitors looking as a M. helenor/M. achilles (Supplementary Fig. 13). Note that removing data with the highest levels of uncertainty in species identity (i.e. when discarding visits performed in the period where M. helenor/M. achilles temporally overlap) does not quantitively affect our results (Supplementary Fig. 14 and Supplementary Tables 5, 6). Using the recorded video, we also measured the duration of the interactions (i.e. the time spent in the cubic space) occurring between patrolling male and the dummy. The ten dummies were each tested during 4 sunny days from 9 a.m. to 2 p.m. (i.e. during 5 h). This resulted in 40 days of experiment over which each dummy was left fluttering on the river bank for a combined duration of 20 h. Dummies were randomly attributed to each day of the experiment. Mark-recapture data suggested a very low rate of individuals passing through the site several times per day (mean percentage of recapture within the same day = 0.95%), thus limiting potential pseudoreplication within each dummy replicate. We recently showed that intraspecific variation in wing colour pattern within the locality is very low in these species25. Using a single dummy per sex and species, as done here, should thus have little impact on the observed behaviours.In order to control for variation in weather (affecting both the activity of patrolling butterflies and of the solar-powered device), we collected hourly data on the percentage of cloud cover for the period and location of our experiment (available at https://www.visualcrossing.com/). A percentage of cloud cover was then associated with all the behavioural observations, and used as a control variable in all statistical analyses.Three-dimensional kinematics of flight interaction with the dummiesTo test whether Morpho males showed different flight behaviours when interacting with the male and female dummy, we filmed the flight interactions using two orthogonally positioned video cameras (Gopro Hero5 Black, recording at 120 images per second) around the dummy setup (Fig. 1a). Stereoscopic video sequences obtained from the two cameras were synchronised with respect to a reference frame (here using a clapperboard). Prior to each filming session, the camera system was calibrated with the direct linear transformation (DLT) technique62 by digitising the positions of a wand moved around the dummy. Wand tracking was done using DLTdv863, and computation of the DLT coefficients was performed using easyWand64. After spatial and temporal calibration, we also used DLTdv8 to digitise the three-dimensional positions of both the visiting (real) butterfly and the dummy butterfly at each video frame by manually tracking the body centroid in each camera view. Butterfly positions throughout the flight trajectory were post-processed using a linear Kalman filter65, providing smoothed temporal dynamics of spatial position, velocity and acceleration of the body centroid. Based on these data, we investigated how spatial position, speed and acceleration of the visitor butterfly varied over the course of the interaction. We proceeded by dividing space into 10 cm spherical intervals around the dummy position ranging from 0 to 1.2 m distance (this step standardises interactions of different durations), and computed the proportion of time spent, the mean speed and acceleration of the interacting butterfly within each distance interval (Fig. 2). We analysed a total of 28 interactions performed by individual Morpho achilles male, including 14 with the dummy of its conspecific male and 14 with the dummy of its conspecific female. Analysed interactions lasted in average 1.44 ± 0.87 (mean ± sd) s.Statistical analysis of behavioural experimentsDifferences in patrolling time were assessed by testing the effect of species on time of capture using Kruskal–Wallis test. To test the effect of visitor identity and dummy characteristics on the number of approaches and interactions, we performed logistic regressions. Approach was treated as a binary variable, where 0 meant ‘passing without approaching’ and 1 meant ‘approaching the dummy setup’. For the interactions, we only considered individuals approaching the setup, such as 0 meant ‘approaching without entering the cubic space’ and 1 meant ‘entering the cubic space’. This allowed getting rid of the uncertainties on whether passing individuals had actually seen the setup or not. We first tested the effect of visiting species on approach and interaction while controlling for dummy’s characteristics to test for intrinsic differences in territoriality (or ‘curiosity’) among species. We then tested the effect of the dummy sex and identity on approach and interaction separately in Morpho helenor and M. achilles. The percentage of cloud cover was also included in the models to control for variation in dummy movements (generated by the solar-powered device), potentially affecting the butterfly response (Supplementary Tables 3 and 4). We further tested if variation in wing area and proportion of iridescent blue among dummies affected the frequency of approach and interaction, again using logistic regression analyses (Supplementary Fig. 7). Statistical significance of each variables was assessed using likelihood ratio tests comparing logistic regression models66. Finally, we tested the effect of dummy sex and identity on the duration of interaction using Kruskal–Wallis tests.Based on the flight kinematic data, we investigated whether flight behaviour during the interaction differed with male vs. female dummies. We ran a mixed-effects model testing the effect of (1) the sex of the dummy and of (2) the distance from dummy (fixed effects), on the proportion of time spent (fixed effects), using the flight ID as a random effect. The flight ID corresponds to the behaviour of a single wild males flying within the ‘interaction space’. Specifically, we tested for the statistical interaction between the sex of the dummy and distance from dummy on the proportion of time spent in the different distance intervals. We then similarly tested for difference in acceleration over the course of the flight interaction, by testing the effect of (1) the sex of the dummy and of (2) the distance from dummy (fixed effects), on the acceleration, with the flight ID as a random effect. We focused on the statistical interaction between the sex of the dummy and the distance from dummy on the mean acceleration in the different distance intervals.Reporting summaryFurther information on research design is available in the Nature Research Reporting Summary linked to this article. More