Philippa Kaur

1.Pakhomov, E. A., Froneman, P. W. & Perissinotto, R. Salp/krill interactions in the Southern Ocean: Spatial segregation and implications for the carbon flux. Deep Sea Res. II 49, 1881–1907 (2002).CAS 
Article 

Google Scholar 
2.Steinberg, D. K. et al. Long-term (1993–2013) changes in macrozooplankton off the Western Antarctic Peninsula. Deep Sea Res. I 101, 54–70 (2015).Article 

Google Scholar 
3.Whitehouse, M. J. et al. Role of krill versus bottom-up factors in controlling phytoplankton biomass in the northern Antarctic waters of South Georgia. Mar. Ecol. Prog. Ser. 393, 69–82 (2009).CAS 
Article 

Google Scholar 
4.Tarling, G. A. & Fielding, S. in Biology and Ecology of Antarctic krill (ed Siegel, V.) 279–319 (Springer International Publishing, 2016).5.Henschke, N., Everett, J. D., Richardson, A. J. & Suthers, I. M. Rethinking the role of salps in the ocean. Trends Ecol. Evol. 31, 720–733 (2016).PubMed 
Article 

Google Scholar 
6.Cavan, E. L. et al. The importance of Antarctic krill in biogeochemical cycles. Nat. Commun. 10, 4742 (2019).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
7.Belcher, A. et al. Krill faecal pellets drive hidden pulses of particulate organic carbon in the marginal ice zone. Nat. Commun. 10, 889 (2019).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
8.Phillips, B., Kremer, P. & Madin, L. P. Defecation by Salpa thompsoni and its contribution to vertical flux in the Southern Ocean. Mar. Biol. 156, 455–467 (2009).Article 

Google Scholar 
9.Siegel, V. & Watkins, J. L. in Biology and Ecology of Antarctic krill (ed Siegel, V.) 21–100 (Springer International Publishing, 2016).10.Atkinson, A. et al. Krill (Euphausia superba) distribution contracts southward during rapid regional warming. Nat. Clim. Change 9, 142–147 (2019).Article 

Google Scholar 
11.Vaughan, D. G. et al. Recent rapid regional climate warming on the Antarctic Peninsula. Clim. Change 60, 243–274 (2003).Article 

Google Scholar 
12.Ducklow, H. W. et al. Marine pelagic ecosystems: the West Antarctic Peninsula. Philos. Trans. R. Soc. B. 362, 67–94 (2007).Article 

Google Scholar 
13.Montes-Hugo, M. et al. Recent changes in phytoplankton communities associated with rapid regional climate change along the Western Antarctic Peninsula. Science 323, 1470–1473 (2009).CAS 
PubMed 
Article 

Google Scholar 
14.Clarke, A. et al. Climate change and the marine ecosystem of the western Antarctic Peninsula. Philos. T. R. Soc., B 362, 149–166 (2007).Article 

Google Scholar 
15.Atkinson, A., Siegel, V., Pakhomov, E. & Rothery, P. Long-term decline in krill stock and increase in salps within the Southern Ocean. Nature 432, 100–103 (2004).CAS 
PubMed 
Article 

Google Scholar 
16.Loeb, V. et al. Effects of sea-ice extent and krill or salp dominance on the Antarctic food web. Nature 387, 897–900 (1997).CAS 
Article 

Google Scholar 
17.Flores, H. et al. Impact of climate change on Antarctic krill. Mar. Ecol. Prog. Ser. 458, 1–19 (2012).Article 

Google Scholar 
18.Cox, M. J. et al. No evidence for a decline in the density of Antarctic krill Euphausia superba Dana, 1850, in the Southwest Atlantic sector between 1976 and 2016. J. Crust. Biol. 38, 656–661 (2018).19.Foxton, P. The Distribution and Life-history of Salpa thompsoni Foxton with observations on a Related Species, Salpa gerlachei Foxton (The University Press, 1966).20.Bernard, K. S., Steinberg, D. K. & Schofield, O. M. E. Summertime grazing impact of the dominant macrozooplankton off the Western Antarctic Peninsula. Deep Sea Res. I 62, 111–122 (2012).Article 

Google Scholar 
21.Condon, R. H. et al. Jellyfish blooms result in a major microbial respiratory sink of carbon in marine systems. Proc. Natl Acad. Sci. 108, 10225–10230 (2011).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
22.Meyer, M. A. & Elsayed, S. Z. Grazing of Euphausia superba Dana on natural phytoplankton populations. Polar Biol. 1, 193–197 (1983).Article 

Google Scholar 
23.Haberman, K. L., Ross, R. M. & Quetin, L. B. Diet of the Antarctic krill (Euphausia superba Dana): II. Selective grazing in mixed phytoplankton assemblages. J. Exp. Mar. Biol. Ecol. 283, 97–113 (2003).Article 

Google Scholar 
24.Schmidt, K. & Atkinson, A. in Biology and Ecology of Antarctic krill (ed Siegel, V.) 175–224 (Springer International Publishing, 2016).25.Andersen, V. in The Biology of Pelagic Tunicates (ed Bone, Q.) 125–137 (Oxford University Press, 1998).26.Mitra, A. et al. Bridging the gap between marine biogeochemical and fisheries sciences; configuring the zooplankton link. Prog. Oceanogr. 129, 176–199 (2014).Article 

Google Scholar 
27.Sailley, S. F., Polimene, L., Mitra, A., Atkinson, A. & Allen, J. I. Impact of zooplankton food selectivity on plankton dynamics and nutrient cycling. J. Plankton Res. 37, 519–529 (2015).CAS 
Article 

Google Scholar 
28.Hamner, W. M., Hamner, P. P., Strand, S. W. & Gilmer, R. W. Behavior of Antarctic krill, Euphausia superba: Chemoreception, feeding, schooling, and molting. Science 220, 433–435 (1983).CAS 
PubMed 
Article 

Google Scholar 
29.DeMott, W. R. in Behavioural Mechanisms of Food Selection (ed Hughes, R. N.) 569–594 (Springer, 1990).30.Le Fèvre, J., Legendre, L. & Rivkin, R. B. Fluxes of biogenic carbon in the Southern Ocean: Roles of large microphagous zooplankton. J. Mar. Syst. 17, 325–345 (1998).Article 

Google Scholar 
31.Moline, M. A., Claustre, H., Frazer, T. K., Schofield, O. & Vernet, M. Alteration of the food web along the Antarctic Peninsula in response to a regional warming trend. Glob. Change Biol. 10, 1973–1980 (2004).Article 

Google Scholar 
32.Frischer, M. E. et al. Selective feeding and linkages to the microbial food web by the doliolid Dolioletta gegenbauri. Limnol. Oceanogr. 66, 1993–2010 (2021).Article 

Google Scholar 
33.Dadon-Pilosof, A., Lombard, F., Genin, A., Sutherland, K. R. & Yahel, G. Prey taxonomy rather than size determines salp diets. Limnol. Oceanogr. 64, 1996–2010 (2019).Article 

Google Scholar 
34.Metfies, K., Nicolaus, A., von Harbou, L., Bathmann, U. & Peeken, I. Molecular analyses of gut contents: elucidating the feeding of co-occurring salps in the Lazarev Sea from a different perspective. Antarct. Sci. 26, 5545–5553 (2014).Article 

Google Scholar 
35.Cleary, A. C., Durbin, E. G. & Casas, M. C. Feeding by Antarctic krill Euphausia superba in the West Antarctic Peninsula: differences between fjords and open waters. Mar. Ecol. Prog. Ser. 595, 39–54 (2018).CAS 
Article 

Google Scholar 
36.Pompanon, F. et al. Who is eating what: Diet assessment using next generation sequencing. Mol. Ecol. 21, 1931–1950 (2012).CAS 
PubMed 
Article 

Google Scholar 
37.Passmore, A. J. et al. DNA as a dietary biomarker in Antarctic krill, Euphausia superba. Mar. Biotechnol. 8, 686–696 (2006).CAS 
Article 

Google Scholar 
38.von Harbou, L. et al. Salps in the Lazarev Sea, Southern Ocean: I. Feeding dynamics. Mar. Biol. 158, 2009–2026 (2011).Article 

Google Scholar 
39.Vernet, M. et al. Primary production throughout austral fall, during a time of decreasing daylength in the western Antarctic Peninsula. Mar. Ecol. Prog. Ser. 452, 45–61 (2012).CAS 
Article 

Google Scholar 
40.Moreau, S. et al. Variability of the microbial community in the western Antarctic Peninsula from late fall to spring during a low ice cover year. Polar Biol. 33, 1599–1614 (2010).Article 

Google Scholar 
41.Selz, V. et al. Distribution of Phaeocystis antarctica-dominated sea ice algal communities and their potential to seed phytoplankton across the western Antarctic Peninsula in spring. Mar. Ecol. Prog. Ser. 586, 91–112 (2018).CAS 
Article 

Google Scholar 
42.Nichols, D. S., Nichols, P. D. & Sullivan, C. W. Fatty acid, sterol and hydrocarbon composition of Antarctic sea ice diatom communities during the spring bloom in McMurdo Sound. Antarct. Sci. 5, 271–278 (1993).Article 

Google Scholar 
43.Fahl, K. & Kattner, G. Lipid Content and fatty acid composition of algal communities in sea-ice and water from the Weddell Sea (Antarctica). Polar Biol. 13, 405–409 (1993).Article 

Google Scholar 
44.Boyd, C. M., Heyraud, M. & Boyd, C. N. Feeding of the Antarctic krill Euphausia superba. J. Crust. Biol. 4, 123–141 (1984).Article 

Google Scholar 
45.Bone, Q., Carré, C. & Chang, P. Tunicate feeding filters. J. Mar. Biol. Assoc. U. K. 83, 907–919 (2003).Article 

Google Scholar 
46.Nelson, M. M., Phleger, C. F., Mooney, B. D. & Nichols, P. D. Lipids of gelatinous Antarctic zooplankton: Cnidaria and Ctenophora. Lipids 35, 551–559 (2000).CAS 
PubMed 
Article 

Google Scholar 
47.Huntley, M. E., Sykes, P. F. & Marin, V. Biometry and trophodynamics of Salpa thompsoni Foxton (Tunicata: Thaliacea) near the Antarctic Peninsula in austral summer, 1983–1984. Polar Biol. 10, 59–70 (1989).Article 

Google Scholar 
48.Hopkins, T. L. Food web of an Antarctic midwater ecosystem. Mar. Biol. 89, 197–212 (1985).Article 

Google Scholar 
49.Paffenhöfer, G. A. & Köster, M. Digestion of diatoms by planktonic copepods and doliolids. Mar. Ecol. Prog. Ser. 297, 303–310 (2005).Article 

Google Scholar 
50.von Harbou, L. Trophodynamics of Salps in the Atlantic Southern Ocean. PhD thesis, University of Bremen (2009).51.Hargraves, P. E. The ebridian flagellates Ebria and Hermesinum. Plankton Biol. Ecol. 49, 9–16 (2002).
Google Scholar 
52.Cavan, E. L. et al. Attenuation of particulate organic carbon flux in the Scotia Sea, Southern Ocean, is controlled by zooplankton fecal pellets. Geophys. Res. Lett. 42, 821–830 (2015).CAS 
Article 

Google Scholar 
53.Smith, K. L. Jr. et al. Large salp bloom export from the upper ocean and benthic community response in the abyssal northeast Pacific: day to week resolution. Limnol. Oceanogr. 59, 745–757 (2014).CAS 
Article 

Google Scholar 
54.Cadée, G. C., González, H. & Schnack-Schiel, S. B. Krill diet affects faecal string settling. Polar Biol. 12, 75–80 (1992).
Google Scholar 
55.Ploug, H., Iversen, M. H., Koski, M. & Buitenhuis, E. T. Production, oxygen respiration rates, and sinking velocity of copepod fecal pellets: Direct measurements of ballasting by opal and calcite. Limnol. Oceanogr. 53, 469–476 (2008).CAS 
Article 

Google Scholar 
56.Atkinson, A., Schmidt, K., Fielding, S., Kawaguchi, S. & Geissler, P. A. Variable food absorption by Antarctic krill: Relationships between diet, egestion rate and the composition and sinking rates of their fecal pellets. Deep Sea Res. II 59-60, 147–158 (2012).CAS 
Article 

Google Scholar 
57.Schmidt, K., Atkinson, A., Pond, D. W. & Ireland, L. C. Feeding and overwintering of Antarctic krill across its major habitats: The role of sea ice cover, water depth, and phytoplankton abundance. Limnol. Oceanogr. 59, 17–36 (2014).Article 

Google Scholar 
58.Cripps, G. C., Watkins, J. L., Hill, H. J. & Atkinson, A. Fatty acid content of Antarctic krill Euphausia superba at South Georgia related to regional populations and variations in diet. Mar. Ecol. Prog. Ser. 181, 177–188 (1999).CAS 
Article 

Google Scholar 
59.Schmidt, K., Atkinson, A., Petzke, K.-J., Voss, M. & Pond, D. W. Protozoans as a food source for Antarctic krill, Euphausia superba: Complementary insights from stomach content, fatty acids, and stable isotopes. Limnol. Oceanogr. 51, 2409–2427 (2006).CAS 
Article 

Google Scholar 
60.Hagen, W., Van Vleet, E. S. & Kattner, G. Seasonal lipid storage as overwintering strategy of Antarctic krill. Mar. Ecol. Prog. Ser. 134, 85–89 (1996).CAS 
Article 

Google Scholar 
61.Kawaguchi, S. & Takahashi, Y. Antarctic krill (Euphausia superba Dana) eat salps. Polar Biol. 16, 479–481 (1996).
Google Scholar 
62.Clarke, L. J., Bestley, S., Bissett, A. & Deagle, B. E. A globally distributed Syndiniales parasite dominates the Southern Ocean micro-eukaryote community near the sea-ice edge. ISME J. 13, 734–737 (2019).CAS 
PubMed 
Article 

Google Scholar 
63.Coats, D. W. & Park, M. G. Parasitism of photosynthetic dinoflagellates by three strains of Amoebophrya (Dinophyta): Parasite survival, infectivity, generation time, and host specificity. J. Phycol. 38, 520–528 (2002).Article 

Google Scholar 
64.Sutherland, K. R., Madin, L. P. & Stocker, R. Filtration of submicrometer particles by pelagic tunicates. Proc. Natl Acad. Sci. 107, 15129–15134 (2010).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
65.Gómez-Gutiérrez, J. & Morales-Avila, J. R. in Biology and Ecology of Antarctic krill (ed Siegel, V.) 351–387 (Springer International Publishing, 2006).66.Cleary, A. C., Casas, M. C., Durbin, E. G. & Gómez-Gutiérrez, J. Parasites in Antarctic krill guts inferred from DNA sequences. Antarct. Sci. 31, 16–22 (2019).Article 

Google Scholar 
67.Zamora-Terol, S., Novotny, A. & Winder, M. Molecular evidence of host-parasite interactions between zooplankton and Syndiniales. Aquat. Ecol. 55, 125–134 (2021).CAS 
Article 

Google Scholar 
68.Kawaguchi, S., Ichii, T. & Naganobu, M. Do krill and salps compete? Contrary evidence from the krill fisheries. CCAMLR Sci. 5, 205–216 (1998).
Google Scholar 
69.Fadeev, E. et al. Microbial communities in the east and west Fram Strait during sea ice melting season. Front. Mar. Sci. 5, 429 (2018).Article 

Google Scholar 
70.Vestheim, H. & Jarman, S. N. Blocking primers to enhance PCR amplification of rare sequences in mixed samples – a case study on prey DNA in Antarctic krill stomachs. Front. Zool. 5, 12 (2008).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
71.Ludwig, W. et al. ARB: a software environment for sequence data. Nucleic Acids Res. 32, 1363–1371 (2004).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
72.R Foundation for Statistical Computing. R: A Language and Environment for Statistical Computing (R Foundation for Statistical Computing, 2019).73.Callahan, B. J. et al. DADA2: High-resolution sample inference from Illumina amplicon data. Nat. Methods 13, 581–583 (2016).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
74.Callahan, B. DADA2 Pipeline Tutorial (1.16), available online: https://benjjneb.github.io/dada2/tutorial.html. Accessed: 3 Feb 2020.75.Martin, M. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet. J. 17, 3 (2011).Article 

Google Scholar 
76.Guillou, L. et al. The protist ribosomal reference database (PR2): A catalog of unicellular eukaryote small sub-unit rRNA sequences with curated taxonomy. Nucleic Acids Res. 41, D597–D604 (2013).CAS 
PubMed 
Article 

Google Scholar 
77.Gong, W. & Marchetti, A. Estimation of 18S gene copy number in marine eukaryotic plankton using a next-generation sequencing approach. Front. Mar. Sci. 6, 219 (2019).Article 

Google Scholar 
78.Metfies, K. et al. Uncovering the intricacies of microbial community dynamics at Helgoland Roads at the end of a spring bloom using automated sampling and 18S meta-barcoding. PLoS ONE 15, e0233921 (2020).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
79.Catlett, D. et al. Evaluation of accuracy and precision in an amplicon sequencing workflow for marine protist communities. Limnol. Oceanogr. Methods 18, 20–40 (2019).Article 

Google Scholar 
80.Kattner, G. & Fricke, H. S. G. Simple gas-liquid-chromatographic method for the simultaneous determination of fatty-acids and alcohols in wax esters of marine organisms. J. Chromatogr. 361, 263–268 (1986).CAS 
Article 

Google Scholar 
81.Gloor, G. B., Macklaim, J. M., Pawlowsky-Glahn, V. & Egozcue, J. J. Microbiome datasets are compositional: And this is not optional. Front. Microbiol. 8, 2224 (2017).PubMed 
PubMed Central 
Article 

Google Scholar 
82.Palarea-Albaladejo, J. & Martín-Fernández, J. A. zCompositions—R package for multivariate imputation of left-censored data under a compositional approach. Chemometrics Intell. Lab. Syst. 143, 85–96 (2015).CAS 
Article 

Google Scholar 
83.Fernandes, A. D., Macklaim, J. M., Linn, T. G., Reid, G. & Gloor, G. B. ANOVA-Like differential expression (ALDEx) analysis for mixed population RNA-Seq. PLoS ONE 8, e67019 (2013).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
84.Quinn, T. P., Richardson, M. F., Lovell, D. & Crowley, T. M. propr: an R-package for identifying proportionally abundant features using compositional data analysis. Sci. Rep. 7, 16252 (2017).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
85.Bian, G. et al. The gut microbiota of healthy aged chinese is similar to that of the healthy young. mSphere 2, e00327–00317 (2017).PubMed 
PubMed Central 
Article 

Google Scholar 
86.Gloor, G. B. & Reid, G. Compositional analysis: A valid approach to analyze microbiome high-throughput sequencing data. Can. J. Microbiol. 62, 692–703 (2016).CAS 
PubMed 
Article 

Google Scholar 
87.Borstein, S. R. dietr: An R package for calculating fractional trophic levels from quantitative and qualitative diet data. Hydrobiologia 847, 4285–4294 (2020).Article 

Google Scholar 
88.Lechowicz, M. J. The sampling characteristics of electivity indices. Oecologia 52, 22–30 (1982).PubMed 
Article 

Google Scholar 
89.Dalsgaard, J., St John, M., Kattner, G., Muller-Navarra, D. & Hagen, W. Fatty acid trophic markers in the pelagic marine environment. Adv. Mar. Biol. 46, 225–340 (2003).PubMed 
Article 

Google Scholar 
90.Graeve, M., Kattner, G. & Hagen, W. Diet-induced changes in the fatty acid composition of Arctic herbivorous copepods: Experimental evidence of trophic markers. J. Exp. Mar. Biol. Ecol. 182, 97–110 (1994).CAS 
Article 

Google Scholar 
91.Kharlamenko, V. I., Zhukova, N. V., Khotimchenko, S. V., Svetashev, V. I. & Kamenev, G. M. Fatty acids as markers of food sources in a shallow-water hydrothermal ecosystem (Kraternaya Bight, Yankich Island, Kurile Islands). Mar. Ecol. Prog. Ser. 120, 231–241 (1995).CAS 
Article 

Google Scholar 
92.Greenacre, M. Compositional Data Analysis in Practice (CRC Press, Taylor & Francis Group, 2018).93.Suh, H.-L. & Nemoto, T. Comparative morphology of filtering structure of five species of Euphausia (Euphausiacea, Crustacea) from the Antarctic Ocean. Proc. NIPR Symp. Polar Biol. 1, 72–83 (1987).
Google Scholar 
94.Alldredge, A. L. & Madin, L. P. Pelagic tunicates: unique herbivores in the marine plankton. Bioscience 32, 655–663 (1982).Article 

Google Scholar 
95.Kelly, P. S. The Ecological Role of Salpa Thompsoni in the Kerguelen Plateau Region of the Southern Ocean: A First Comprehensive Evaluation. PhD thesis, University of Tasmania (2019).96.Ericson, J. A. et al. Seasonal and interannual variations in the fatty acid composition of adult Euphausia superba Dana, 1850 (Euphausiacea) samples derived from the Scotia Sea krill fishery. J. Crust. Biol. 38, 662–672 (2018).
Google Scholar 
97.Martin, D. L., Ross, R. M., Quetin, L. B. & Murray, A. E. Molecular approach (PCR-DGGE) to diet analysis in young Antarctic krill Euphausia superba. Mar. Ecol. Prog. Ser. 319, 155–165 (2006).CAS 
Article 

Google Scholar 
98.Matsuoka, K. et al. Quantarctica, an integrated mapping environment for Antarctica, the Southern Ocean, and sub-Antarctic islands. Environ. Model. Softw. 140, 105015 (2021).Article 

Google Scholar  More

We reviewed and mapped antimicrobial resistance in aquatic food animals in Asia during a period of substantial industry growth. Our findings indicate that between 2000 and 2018, antimicrobial resistance in bacteria from cultured aquatic food animals was stable (33%) while the resistance from wild-caught aquatic food animals decreased sharply (52% to 22%). These trends represent currently available evidence from point prevalence surveys, which serve as a surrogate in the absence of systematic surveillance and should be interpreted cautiously. Structured, systematic surveillance will be imperative to document trends in multi-drug resistance at the sub-national level in the future.Our results are consistent with an analysis of antimicrobial resistance in aquaculture-derived bacteria from forty countries, nearly half of which in Asia, which identified a global mean multi-antibiotic resistance index of .25, and a higher index ( >.35) in low-income and middle-income countries in Asia27. Although antimicrobial use in surveys from cultured animals was most frequently unspecified, in the limited surveys that recorded whether on-farm antimicrobials were either used or not used (n = 63; 11%), use was associated with higher multi-drug resistance than the absence of use (p  More

1.Food and Agriculture Organization of the United Nations. The State of World Fisheries and Aquaculture http://www.fao.org/3/i2727e/i2727e00.htm (2012).2.Isbell, F. et al. Biodiversity increases the resistance of ecosystem productivity to climate extremes. Nature 526, 574–577 (2015).ADS 
CAS 
PubMed 
Article 

Google Scholar 
3.IPBES. Summary for policymakers of the global assessment report on biodiversity and ecosystem services of the Intergovernmental Science-Policy Platform on Biodiversity and Ecosystem Services https://doi.org/10.5281/zenodo.3553579 (2019).4.Tittensor, D. P. et al. A mid-term analysis of progress toward international biodiversity targets. Science 346, 241–244 (2014).ADS 
CAS 
PubMed 
Article 

Google Scholar 
5.Dornelas, M. et al. Assemblage time series reveal biodiversity change but not systematic loss. Science 344, 296–299 (2014).ADS 
CAS 
PubMed 
Article 

Google Scholar 
6.Gonzalez, A. et al. Estimating local biodiversity change: a critique of papers claiming no net loss of local diversity. Ecology 97, 1949–1960 (2016).PubMed 
Article 

Google Scholar 
7.Elahi, R. et al. Recent trends in local-scale marine biodiversity reflect community structure and human impacts. Curr. Biol. 25, 1938–1943 (2015).CAS 
PubMed 
Article 

Google Scholar 
8.Blowes, S. A. et al. The geography of biodiversity change in marine and terrestrial assemblages. Science 366, 339–345 (2019).ADS 
CAS 
PubMed 
Article 

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

Google Scholar 
10.Loreau, M. Biodiversity and ecosystem functioning: current knowledge and future challenges. Science 294, 804–808 (2001).ADS 
CAS 
PubMed 
Article 

Google Scholar 
11.Covich, A. P. et al. The role of biodiversity in the functioning of freshwater and marine benthic ecosystems. Bioscience 54, 767–775 (2004).Article 

Google Scholar 
12.Hooper, D. U. et al. Effects of biodiversity on ecosystem functioning: a consensus of current knowledge. Ecol. Monogr. 75, 3–35 (2005).Article 

Google Scholar 
13.Widdicombe, C. E., Eloire, D., Harbour, D., Harris, R. P. & Somerfield, P. J. Long-term phytoplankton community dynamics in the Western English Channel. J. Plankton Res. 32, 643–655 (2010).Article 

Google Scholar 
14.Eloire, D. et al. Temporal variability and community composition of zooplankton at station L4 in the Western Channel: 20 years of sampling. J. Plankton Res. 32, 657–679 (2010).Article 

Google Scholar 
15.Hillebrand, H. et al. In Handbook on Marine Environment Protection (eds Salomon, M. & Markus, T.) 21 (Springer, 2018).16.Bindoff, N. L. et al. In IPCC Special Report on the Ocean and Cryosphere in a Changing Climate (eds H.-O. Pörtner, D. C. Roberts, V. Masson-Delmotte, P. Zhai, M. Tignor, E. Poloczanska, K. Mintenbeck, A. Alegría, M. Nicolai, A. Okem, J. Petzold, B. Rama, N. M. W) Cambridge University Press (2019).17.Barton, A. D., Irwin, A. J., Finkel, Z. V. & Stock, C. A. Anthropogenic climate change drives shift and shuffle in North Atlantic phytoplankton communities. Proc. Natl Acad. Sci. USA 113, 2964–2969 (2016).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
18.Pecuchet, L. et al. Spatio‐temporal dynamics of multi‐trophic communities reveal ecosystem‐wide functional reorganization. Ecography 43, 197–208 (2020).Article 

Google Scholar 
19.Poloczanska, E. S. et al. Responses of marine organisms to climate change across oceans. Front. Mar. Sci. 3, 62 (2016).20.Pennekamp, F. et al. Biodiversity increases and decreases ecosystem stability. Nature 563, 109–112 (2018).ADS 
CAS 
PubMed 
Article 

Google Scholar 
21.Duffy, J. E., Godwin, C. M. & Cardinale, B. J. Biodiversity effects in the wild are common and as strong as key drivers of productivity. Nature 549, 261–264 (2017).ADS 
CAS 
PubMed 
Article 

Google Scholar 
22.Worm, B. et al. Impacts of biodiversity loss on ocean ecosystem services. Science 314, 787–790 (2006).ADS 
CAS 
PubMed 
Article 

Google Scholar 
23.Blois, J. L., Zarnetske, P. L., Fitzpatrick, M. C. & Finnegan, S. Climate change and the past, present, and future of biotic interactions. Science 341, 499–504 (2013).ADS 
CAS 
PubMed 
Article 

Google Scholar 
24.Dossena, M. et al. Warming alters community size structure and ecosystem functioning. Proc. R. Soc. B Biol. Sci. 279, 3011–3019 (2012).Article 

Google Scholar 
25.Brander, K. & Kiørboe, T. Decreasing phytoplankton size adversely affects ocean food chains. Glob. Chang. Biol. https://doi.org/10.1111/gcb.15216 (2020).26.Mouw, C. B., Barnett, A., McKinley, G. A., Gloege, L. & Pilcher, D. Phytoplankton size impact on export flux in the global ocean. Glob. Biogeochem. Cycles 30, 1542–1562 (2016).ADS 
CAS 
Article 

Google Scholar 
27.Riahi, K. et al. RCP 8.5—a scenario of comparatively high greenhouse gas emissions. Clim. Change 109, 33–57 (2011).ADS 
CAS 
Article 

Google Scholar 
28.Magnan, A. K. et al. Implications of the Paris agreement for the ocean. Nat. Clim. Chang. 6, 732–735 (2016).ADS 
Article 

Google Scholar 
29.Kuhn, A. M. et al. Temporal and spatial scales of correlation in marine phytoplankton communities. J. Geophys. Res. Ocean. 124, 9417–9438 (2019).ADS 
Article 

Google Scholar 
30.Sonnewald, M., Dutkiewicz, S., Hill, C. & Forget, G. Elucidating ecological complexity: unsupervised learning determines global marine eco-provinces. Sci. Adv. 6, eaay4740 (2020).ADS 
PubMed 
PubMed Central 
Article 

Google Scholar 
31.Dutkiewicz, S., Boyd, P. W. & Riebesell, U. Exploring biogeochemical and ecological redundancy in phytoplankton communities in the global ocean. Glob. Chang. Biol. 27, 1196–1213 (2021).ADS 
PubMed 
PubMed Central 
Article 

Google Scholar 
32.Flombaum, P., Wang, W.-L., Primeau, F. W. & Martiny, A. C. Global picophytoplankton niche partitioning predicts overall positive response to ocean warming. Nat. Geosci. 13, 116–120 (2020).ADS 
CAS 
Article 

Google Scholar 
33.Righetti, D., Vogt, M., Gruber, N., Psomas, A. & Zimmermann, N. E. Global pattern of phytoplankton diversity driven by temperature and environmental variability. Sci. Adv. 5, eaau6253 (2019).ADS 
PubMed 
PubMed Central 
Article 

Google Scholar 
34.Ibarbalz, F. M. et al. Global trends in marine plankton diversity across kingdoms of life. Cell 179, 1084–1097.e21 (2019).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
35.Bopp, L. et al. Multiple stressors of ocean ecosystems in the 21st century: projections with CMIP5 models. Biogeosciences 10, 6225–6245 (2013).ADS 
Article 

Google Scholar 
36.Kwiatkowski, L. et al. Twenty-first century ocean warming, acidification, deoxygenation, and upper-ocean nutrient and primary production decline from CMIP6 model projections. Biogeosciences 17, 3439–3470 (2020).ADS 
CAS 
Article 

Google Scholar 
37.Cabré, A., Marinov, I. & Leung, S. Consistent global responses of marine ecosystems to future climate change across the IPCC AR5 earth system models. Clim. Dyn. 45, 1253–1280 (2015).Article 

Google Scholar 
38.Bopp, L., Aumont, O., Cadule, P., Alvain, S. & Gehlen, M. Response of diatoms distribution to global warming and potential implications: a global model study. Geophys. Res. Lett. 32, n/a−n/a (2005).Article 
CAS 

Google Scholar 
39.Dutkiewicz, S. et al. Dimensions of marine phytoplankton diversity. Biogeosciences 17, 609–634 (2020).ADS 
Article 

Google Scholar 
40.Dutkiewicz, S., Scott, J. R. & Follows, M. J. Winners and losers: ecological and biogeochemical changes in a warming ocean. Glob. Biogeochem. Cycles 27, 463–477 (2013).ADS 
CAS 
Article 

Google Scholar 
41.Marinov, I., Doney, S. C. & Lima, I. D. Response of ocean phytoplankton community structure to climate change over the 21st century: partitioning the effects of nutrients, temperature, and light. Biogeosciences 7, 3941–3959 (2010).ADS 
Article 

Google Scholar 
42.Dutkiewicz, S., Ward, B. A., Scott, J. R. & Follows, M. J. Understanding predicted shifts in diazotroph biogeography using resource competition theory. Biogeosciences 11, 5445–5461 (2014).ADS 
Article 

Google Scholar 
43.Dutkiewicz, S. et al. Impact of ocean acidification on the structure of future phytoplankton communities. Nat. Clim. Chang. 5, 1002–1006 (2015).ADS 
CAS 
Article 

Google Scholar 
44.Kooijman, S. A. L. M. & Troost, T. A. Quantitative steps in the evolution of metabolic organisation as specified by the dynamic energy budget theory. Biol. Rev. 82, 113–142 (2007).CAS 
PubMed 
Article 

Google Scholar 
45.Lévy, M., Jahn, O., Dutkiewicz, S., Follows, M. J. & D’Ovidio, F. The dynamical landscape of marine phytoplankton diversity. J. R. Soc. Interface 12, 20150481 (2015).PubMed 
PubMed Central 
Article 

Google Scholar 
46.Beaugrand, G., Edwards, M., Raybaud, V., Goberville, E. & Kirby, R. R. Future vulnerability of marine biodiversity compared with contemporary and past changes. Nat. Clim. Chang. 5, 695–701 (2015).ADS 
Article 

Google Scholar 
47.Thomas, M. K., Kremer, C. T., Klausmeier, C. A. & Litchman, E. A global pattern of thermal adaptation in marine phytoplankton. Science 338, 1085–1088 (2012).ADS 
CAS 
PubMed 
Article 

Google Scholar 
48.Hillebrand, H. et al. Biodiversity change is uncoupled from species richness trends: consequences for conservation and monitoring. J. Appl. Ecol. 55, 169–184 (2018).Article 

Google Scholar 
49.Litchman, E. & Klausmeier, C. A. Trait-based community ecology of phytoplankton. Annu. Rev. Ecol. Evol. Syst. 39, 615–639 (2008).Article 

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

Google Scholar 
51.Stock, C. A. et al. Reconciling fisheries catch and ocean productivity. Proc. Natl Acad. Sci. USA 114, E1441–E1449 (2017).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
52.Armengol, L., Calbet, A., Franchy, G., Rodríguez-Santos, A. & Hernández-León, S. Planktonic food web structure and trophic transfer efficiency along a productivity gradient in the tropical and subtropical Atlantic Ocean. Sci. Rep. 9, 2044 (2019).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
53.Cram, J. A. et al. The role of particle size, ballast, temperature, and oxygen in the sinking flux to the deep sea. Glob. Biogeochem. Cycles 32, 858–876 (2018).ADS 
CAS 
Article 

Google Scholar 
54.Kéfi, S., Dakos, V., Scheffer, M., Van Nes, E. H. & Rietkerk, M. Early warning signals also precede non-catastrophic transitions. Oikos 122, 641–648 (2013).Article 

Google Scholar 
55.Doncaster, C. P. et al. Early warning of critical transitions in biodiversity from compositional disorder. Ecology 97, 3079–3090 (2016).PubMed 
Article 

Google Scholar 
56.Gunderson, L. H. Ecological resilience—in theory and application. Annu. Rev. Ecol. Syst. 31, 425–439 (2000).Article 

Google Scholar 
57.Benedetti, F. et al. The seasonal and inter-annual fluctuations of plankton abundance and community structure in a North Atlantic Marine Protected Area. Front. Mar. Sci. 6, 214 (2019).58.Pannard, A., Bormans, M. & Lagadeuc, Y. Short-term variability in physical forcing in temperate reservoirs: effects on phytoplankton dynamics and sedimentary fluxes. Freshw. Biol. 52, 12–27 (2007).CAS 
Article 

Google Scholar 
59.Vidal, T., Calado, A. J., Moita, M. T. & Cunha, M. R. Phytoplankton dynamics in relation to seasonal variability and upwelling and relaxation patterns at the mouth of Ria de Aveiro (West Iberian Margin) over a four-year period. PLoS One 12, e0177237 (2017).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
60.Cermeño, P., de Vargas, C., Abrantes, F. & Falkowski, P. G. Phytoplankton biogeography and community stability in the ocean. PLoS One 5, e10037 (2010).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
61.Allen, S. et al. Interannual stability of phytoplankton community composition in the North-East Atlantic. Mar. Ecol. Prog. Ser. 655, 43–57 (2020).ADS 
Article 

Google Scholar 
62.Barton, A. D., Lozier, M. S. & Williams, R. G. Physical controls of variability in North Atlantic phytoplankton communities. Limnol. Oceanogr. 60, 181–197 (2015).ADS 
Article 

Google Scholar 
63.Collins, S., Rost, B. & Rynearson, T. A. Evolutionary potential of marine phytoplankton under ocean acidification. Evol. Appl. 7, 140–155 (2014).CAS 
PubMed 
Article 

Google Scholar 
64.Lohbeck, K. T., Riebesell, U. & Reusch, T. B. H. Adaptive evolution of a key phytoplankton species to ocean acidification. Nat. Geosci. 5, 346–351 (2012).ADS 
CAS 
Article 

Google Scholar 
65.Irwin, A. J., Finkel, Z. V., Müller-Karger, F. E. & Troccoli Ghinaglia, L. Phytoplankton adapt to changing ocean environments. Proc. Natl Acad. Sci. USA 112, 5762–5766 (2015).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
66.Cael, B. B. et al. Marine ecosystem changepoints spread under ocean warming in an Earth System Model. Geophys. Res. Lett.67.Cael, B. B., Dutkiewicz, S. & Henson, S. A. Abrupt shifts in 21st-century plankton communities. Sci. Adv.68.Parmesan, C. & Yohe, G. A globally coherent fingerprint of climate change impacts across natural systems. Nature 421, 37–42 (2003).ADS 
CAS 
Article 

Google Scholar 
69.Burrows, M. T. et al. The pace of shifting climate in marine and terrestrial ecosystems. Science 334, 652–655 (2011).ADS 
CAS 
PubMed 
Article 

Google Scholar 
70.Chivers, W. J., Walne, A. W. & Hays, G. C. Mismatch between marine plankton range movements and the velocity of climate change. Nat. Commun. 8, 14434 (2017).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
71.Poloczanska, E. S. et al. Global imprint of climate change on marine life. Nat. Clim. Chang. 3, 919–925 (2013).ADS 
Article 

Google Scholar 
72.Jonkers, L., Hillebrand, H. & Kucera, M. Global change drives modern plankton communities away from the pre-industrial state. Nature 570, 372–375 (2019).ADS 
CAS 
PubMed 
Article 

Google Scholar 
73.Pond, D. W., Tarling, G. A. & Mayor, D. J. Hydrostatic pressure and temperature effects on the membranes of a seasonally migrating marine copepod. PLoS One 9, e111043 (2014).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
74.Mayor, D. J., Sommer, U., Cook, K. B. & Viant, M. R. The metabolic response of marine copepods to environmental warming and ocean acidification in the absence of food. Sci. Rep. 5, 13690 (2015).ADS 
PubMed 
PubMed Central 
Article 

Google Scholar 
75.Richardson, D. M. & Pyšek, P. Elton, C.S. 1958: The ecology of invasions by animals and plants. London: Methuen. Prog. Phys. Geogr. Earth Environ. 31, 659–666 (2007).Article 

Google Scholar 
76.May, R. M. Qualitative stability in model ecosystems. Ecology 54, 638–641 (1973).Article 

Google Scholar 
77.Lotze, H. K. et al. Global ensemble projections reveal trophic amplification of ocean biomass declines with climate change. Proc. Natl Acad. Sci. USA 116, 12907–12912 (2019).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
78.Barange, M. et al. Impacts of climate change on marine ecosystem production in societies dependent on fisheries. Nat. Clim. Chang. 4, 211–216 (2014).ADS 
Article 

Google Scholar 
79.Marañón, E. et al. Unimodal size scaling of phytoplankton growth and the size dependence of nutrient uptake and use. Ecol. Lett. 16, 371–379 (2013).PubMed 
Article 

Google Scholar 
80.Sokolov, A. P. et al. The MIT Integrated Global System Model (IGSM) Version 2: Model Description and Baseline Evaluation Joint Program Report Series, pp. 40 https://globalchange.mit.edu/publication/14579 (2005).81.Dutkiewicz, S. et al. Ocean colour signature of climate change. Nat. Commun. 10, 578 (2019).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
82.Monier, E., Scott, J. R., Sokolov, A. P., Forest, C. E. & Schlosser, C. A. An integrated assessment modeling framework for uncertainty studies in global and regional climate change: the MIT IGSM-CAM (version 1.0). Geosci. Model Dev. 6, 2063–2085 (2013).ADS 
Article 

Google Scholar 
83.Moss, R. H. et al. The next generation of scenarios for climate change research and assessment. Nature 463, 747–756 (2010).ADS 
CAS 
PubMed 
Article 

Google Scholar 
84.Buitenhuis, E. T. et al. MAREDAT: towards a world atlas of MARine Ecosystem DATa. Earth Syst. Sci. Data 5, 227–239 (2013).ADS 
Article 

Google Scholar 
85.Ward, B. A. Temperature-correlated changes in phytoplankton community structure are restricted to polar waters. PLoS One 10, e0135581 (2015).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
86.Dutkiewicz, S. GUD IGSM depth integrated biomass https://doi.org/10.7910/DVN/LWHQNS (2021).87.Dutkiewicz, S. & Jahn, O. GUD IGSM numerical code and inputs https://doi.org/10.7910/DVN/UA8VNU (2021). More

The model is organized as follows. We simulated water striders as an array of the discrete “objects” which interact one with another according to more or less natural rules. The interaction includes strong short-range repulsion between the animals. Mathematically, short range repulsion means that one animal cannot penetrate inside a private area of other animal. Normally such repulsion appears at relatively short distances corresponding to a radius (R^{repuls}) of their private territory. Besides, there is mutual attraction at longer distances (R^{attract} > R^{repuls}), which biologically corresponds to a tendency to aggregation4. The tendency to aggregation often gives some competitive advantage due to possible collective reactions on the external challenges (for example, on the attacks of predators).The simplest way to simulate an interaction with regulated characteristic distance is to use Gaussian effective potential with some characteristic radius (R_{0}), since the Gaussian potential guaranty the stability in dynamics simulations with relatively large (Delta t)4. Repulsion force in this case looks as follows:$$f_{j}^{rep} = A_{0} left( {vec{r} – vec{r}_{j} } right)expleft[ { – left( {frac{{vec{r} – vec{r}_{j} }}{{R_{0} }}} right)^{2} } right],$$
(2)
where factor A0  R^{repuls}) attraction normally causes a minimum of the interaction energy at some intermediate distance (R^{min }). In infinite empty space, being used in the equations of motion, such a combination of the forces leads to an ordering of the objects with the equilibrium distance corresponding to the position (R_{min }) of the energy minimum. But, if the area ({ [0,Lx],[0,Ly]}) is limited, the equations of motion must be supplied by appropriate boundary conditions which do not allow the animals to leave this area.The simplest way to introduce the boundary conditions is to apply mathematically “soft” but extremely high and narrow walls around the system, which repulse the animals back to the internal space with exponentially growing force$$f_{j}^{bound} = B_{j}^{bound} exp left[ { – left| {frac{{overrightarrow {r}_{j} – overrightarrow {r}_{bound} }}{{R_{{}}^{bound} }}} right|} right]$$
(5)
acting in the direction opposite to that in which the animal occasionally crosses any of the boundaries. The words “extremely high” mean that the amplitude of this force should be supplied by the pre-factor (B_{j}^{bound}) which is much bigger than the amplitude of the repulsion force between the animals (regulated by the amplitude (B_{{_{jk} }}^{repuls})) to be able to overpower their mutual repulsion (B_{j}^{bound} > > B_{j}^{repuls}) pushing them out the boundaries. Therefore also the characteristic distance of the repulsion from the wall should be much shorter than the typical distance between the animals in the population (R^{bond} < < R^{min }). In this sense the boundary should be as narrow as possible.One can expect that mean density of the population is not extremely high and does not force the animals to stay exactly on the minimal distance (R^{min }) between them. In this case total combination of the attracting and repulsing interactions combined with the rejecting boundary conditions normally leads to the specific patterns where relatively dense groups of the animals spread on the distances close to the equilibrium ones are accompanied by almost empty voids between them5. We have checked this hypothesis numerically many times. Below, such patterns will be seen in all the particular visualizations of the numerical results.It should be noted that mathematically, if number of the animals increases inside of the same limited and already populated area, the individuals tend to fill all the voids with the equilibrium (sometimes, even higher than equilibrium) density. It happens in real population as well in that cases when the number of the animals grows quickly for a self-regulation of the population and they become simply forced to occupy every empty space inside the area.However, normally when the population grows too quickly the growth should be regulated by a number of factors. In particular, it will be restricted by a competition for the food and space accompanied by death of some of the participants of the process. From mathematical point of view it means that the equations of motion, which will be written below, must be accompanied by natural generation of the new individuals as well as by their (reasonable) disappearance from the system.It is almost impossible to write such a generation in analytical form, but in numerical model it can be formally presented as a set of natural rules. First of all the length of the array is supposed to be a variable integer. Let us suppose that at every step of calculation (Delta t) the length of the array can potentially increase to the new one (N to N + 1). The potential act of the generation becomes real if the numerically generated random number (varsigma) uniformly distributed in interval [0;1] is less than some threshold (varsigma_{thres} < 1). In the numerical procedure the rate of the generation is certainly regulated by (varsigma_{{{text{thres}}}}) and can be extremely low for example, if (varsigma_{thres} < < 1).New member of the array is generated with random coordinates within the prescribed area ({ [0,Lx],[0,Ly]}) with the mass ((m_{N + 1})) randomly distributed around a starting (basically small) one (m_{s}). Qualitatively it means that at every time moment the array can be supplemented with some probability by a new member which is physically placed in some arbitrary place inside the area ({ [0,Lx],[0,Ly]}). At the same time, we have to introduce a process of disappearance of some of the array members. Natural criterion for this will be established below from the checking of the accumulation or loss of the mass (e.g. reserve fat mass) (m_{j}) by every individual. It is supposed that, if an animal gets critical size (which can be both: maximal (m_{max }) or minimal (m_{min })) it disappears from the system and the length of the array decreases (N to N - 1).Both of such events must be adjusted numerically to some rate natural for a particular biological system. It is quite expected from the very beginning, that there should be a kind of balance between the creation and disappearance rates. If new animas statistically appear too often, the overpopulation will take place. In opposite limit the array will quickly shrink (N to 0) and it will cause a distinction of the population.Obviously both these rates should be naturally regulated by the available food resources. For the particular problem under consideration the resources are generated by a random deposition of the potentially available food onto the water surface. It means that we introduce a new array for the food with coordinates (overrightarrow {r}_{n}^{food}). The length of the array (N^{food}) is also variable and index (n) running in the interval (n = 0,..,N^{food}). In principle, it is possible, and very often happened in our simulations, that at some particular moment available food inside the area can completely disappear. In this case the length of the food array reduces to zero (N^{food} = 0).The food is generated by the random deposition of the food portions distributed inside the area of the system ({ [0,Lx],[0,Ly]}). Generally, it is organized in the same manner as the deposition of the new individuals. If the random number (zeta) uniformly distributed in the interval [0;1] is less than some threshold (varsigma_{thres} < 1) the food portion is deposed at an arbitrary place of the area, in principle at any given time step (Delta t).It is obvious that if the threshold is much smaller than unit: (zeta_{thres}^{food} < < 1) the food portions are produced very rarely. Certainly, the rate of the “food production” in the frames of the model must be properly regulated to make its behaviour natural. Of course, it relates to the size of the portions, but also to the intervals between the food depositions. In any case, these intervals are expected to be much longer than discrete time steps (Delta t) used to solve numerically the equations of motion. At the same time, it must be much shorter than other (biologically reasonable) time-scales of the problem. To control the stability and reasonability of the simulations we have varied this rate in very wide intervals. It was found that at all the reasonable rates food balance in the system tends to the stationary scenario, which corresponds to the expected scenarios in nature.When a portion of the food falls onto the surface the animals which occasionally appear relatively close to it are attracted to this food portion and “eat” it with some characteristic rate. In the particular simulation it was simulated by the additional term of the attraction force:$$f_{{_{jn} }}^{food} (overrightarrow {r}_{j} ,overrightarrow {r}_{n} ) = B_{jn}^{food} (overrightarrow {r}_{j} - overrightarrow {r}_{n} )exp left[ { - left( {frac{{overrightarrow {r}_{j} - overrightarrow {r}_{n} }}{{R_{{_{j} }}^{food} }}} right)^{2} } right]$$ (6) As it is seen from the interactions in the model the animals compete for the food, repulsing one another and reacting faster on the force (f_{{_{jk} }}^{food} (overrightarrow {r}_{j} ,overrightarrow {r}_{k} )). It’s why we apply here non-uniform coefficient (B_{{_{jk} }}^{food}), which is different for different index (j) and in fact depends on some power of mass (B_{jk}^{food} sim m_{j}^{alpha }) with some exponent (alpha). Scaling estimation gives the value of the exponent (alpha) = 2/3. Another form of the competition is related to their simultaneous consumption the food with different rate depending on the size of the different individuals. We suppose that the consumption is proportional to already accumulated mass of the individual: (partial m_{j} /partial t = mu m_{j}). Effective dumping (eta_{j}) on the water surface also depends on the mass of animal (eta_{j} sim m_{j}^{beta }), where estimated exponent (beta) = 1/3.Let us remind that the consumption of the food portion with index (n) is possible only when the animal is close enough to this portion. In other words, one more threshold has to be incorporated: (left| {overrightarrow {r}_{j} - overrightarrow {r}_{n} } right| < R_{thres}^{food}). As we see, the model consists of the dynamic equations of motion and a number of the procedures, which work in parallel and essentially affect both: the dynamic behavior and the results.The equations of motion can be formally written accumulating all the above mentioned forces of the problem:$$m_{j} partial^{2} r_{j} /partial t^{2} + eta_{j} partial r_{j} /partial t = sumlimits_{k} {left[ {f_{{_{jk} }}^{repuls} (overrightarrow {r}_{j} ,overrightarrow {r}_{k} ) + f_{{_{jk} }}^{attract} (overrightarrow {r}_{j} ,overrightarrow {r}_{k} )} right]} + sumlimits_{n} {f_{{_{jn} }}^{food} (overrightarrow {r}_{j} ,overrightarrow {r}_{n} )} + sumlimits_{Boundaries} {f_{j}^{bound} }$$ (7) with their solution combined with the procedures described above. All the initial velocities are zero. This combination makes the solution nontrivial because it is performed for the arrays with changeable lengths and with the varied sizes and forces of the participants. Nevertheless, expected scenarios can be generally predicted and verified later by the large set of numerical experiments at different combinations of the parameters.The longer the time an individual spends near the portion of food and the larger the size it has to the current moment the faster it eats up and faster accumulates additional mass. In general, the larger animals consume more food. However, one should remember that large animals may also have disadvantages. For example, due to inertia of motion faster animals can occasionally jump over a good position near the food. We have observed many such local events during the simulations. Such events partially can be observed in the movies presented below. It is also possible that small animals can appear near the randomly deposed food earlier than the big ones and start consuming it.In some cases such possibilities even emerges due to not just probabilistic but quite regular reason. For example, it can happen due to a further described effect of the “wave of fear”. By influence of such a wave (being scared by somebody moving along the shore), the animals almost synchronically start to run away from one of the boundaries. Big, strong and fast animals escape far away, while the small ones still remain near to the shore. As result, the food deposed into this region will be partially (or even completely) consumed by the weak individuals before the strong ones will return. This effect will be reported further on.Basically, dynamic balance of the population is determined by the relation between a set of the time constants of the problem: how often the food appears, how quick the animals consume it and how efficiently they grow. However, the population structure is also important for the total balance. The more an animal consumes the bigger and stronger it becomes, than faster moves to the new portion of the food and more eats. The individuals are born small and more or less equal. They grow initially depending on their “luck” only. However, population quickly forms some non-uniform distribution of the sizes.In terms of distributions one can say that with the time the individual members of the array shift along the distribution to the larger or smaller size depending on the personal balance of accumulation or losing the mass. At every instant moment the distribution actually demonstrates two opposite fluxes: to the limit of the small and large sizes. Strongest animals evolve to the predefined critical large size (m_{max }), at which they leave the population. Weakest ones tend to the critically small weight mmin and also leave the population. The balance is maintained by the fact that at every time moment the population is incorporated by new young members who either grow or die.In this respect, it is interesting to observe the dynamics of the animal’s sizes distribution. Obviously, the initial population consisting of young animals has histogram localized around small sizes. Later the distribution becomes wider and its edge moves to the larger sizes. The asymptotic histogram shape is determined by the two opposite processes described above. This shape will be compared further with a real distribution found from the field observations.Typical behavior of the model evolution for a population inside a limited area is presented in illustrative movie (video_S2). For a convenience, the population is formally divided into 3 subgroups, which are plotted by the circles having different sizes (small, medium, and big) and colors (blue, green, and red, respectively). The randomly deposed portions of food are shown by large black circles.We start from relatively young population, which contains 100 individuals and allow them to move, grow and disappear according to the rules of game. This behavior demonstrates correlative motions, as well as variation in the population composition, which look self-consistent and rather natural. The interactions between the animals with other animals and with food cause quite predictable motions of the individuals and the changes of their sizes (and colors, respectively) when they leave one of the three groups and join to another. The process stabilizes with the time and seems to become stationary. Below we will study this process by the quantitative time-dependencies and statistical histograms.The bigger an animal the faster it moves in average. Therefore, it is convenient to sort the animals according to the masses and velocities. Such representation is reproduced in the third movie “video_S3.mp4”. The separation between subgroups, their evolution from the initial to a stationary state is clearly seen from the movie in dynamics. Moreover, one can observe even how group of the small animals divides by itself into two subgroups with quite well pronounced gap between them. This separation is caused by the mentioned above two fluxes of the sub-populations, which either grow from the initially small sizes to the medium ones or “in unlucky case” decrease and disappear.Actually, the rate of evolution with fast changes of the masses after very few acts of the interaction with deposed food represented in both movies is strongly overestimated in contrast to the reality. Therefore, as a next step we reduced the rates of accumulation and loss of the mass and proportionally prolonged the simulation time. This process was also recorded in two movies “video_S4.mp4” and “video_S5.mp4”. To reduce the lengths of the videos the time intervals between the frames were made 10 times longer. Because of this reduction the movies look almost stroboscopic. Nevertheless, the movies illustrate quite well long-time dynamics of the system at realistic rate of food deposition and consumption.Information about the system accumulated during long runs including typical pattern formed by the moving animals at some intermediate stage and other plots is presented in Figs. 1, 2, 3 and 4. The spatial pattern in Fig. 1 is taken at some intermediate stage of evolution. Blue, green and red circles correspond to the small medium and large sizes of the individuals. Large black circles mark the food portions, which are recently deposed and not eaten yet.Figure 1Typical spatial pattern formed by the population at intermediate stage of evolution. Blue, green and red circles correspond to the small medium and large sizes of the individuals. One can see the groups packed by the individuals with the distances close to the equilibrium Rmin(R^{min }) and voids. Black circles show food portions remaining to the current moment.Full size imageFigure 2Histograms of the distances between nearest neighbors. Blue line shows some instant distribution. Black line presents histogram accumulated during long run. Maximum of the histogram corresponds to the distance Rmin close to the equilibrium. Long tail of the black curve corresponds to the presence of some individuals inside the voids.Full size imageFigure 3Mass depending values: (a) instant configuration of the velocities of different subgroups marked by the same symbols as in Fig. 1; (b) instant and long run averaged histograms of the masses for complete population; (c) instant and accumulated distributions of the velocities monotonously increasing with masses (volumes) of the individuals.Full size imageFigure 4The comparison between the numerically (curves with points) and experimentally (bars) obtained data. The subplot (a) shows the histogram of the distances between nearest neighbors and the subplot (b) presents the histograms of the volumes of the animals (which are supposed to be approximately proportional to their masses).Full size imageOne can well distinguish in Fig. 1 relatively dense domains packed with the distances between the individuals close to the equilibrium radius (R^{min }) determined by the relation between repulsive (f_{{_{jk} }}^{repuls} (overrightarrow {r}_{j} ,overrightarrow {r}_{k} )) and attraction (f_{{_{jk} }}^{attract} (overrightarrow {r}_{j} ,overrightarrow {r}_{k} )) forces. These domains are separated by voids. To characterize such patterns quantitatively the histograms of the distances between nearest neighbors6 were accumulated during a long simulation run. These histograms are shown in Fig. 2. Thin blue line here reproduces some arbitrary instant distribution. Black line with the circles presents the histogram accumulated during a long run. Maximum of the histogram is close to the equilibrium distance (R^{min }) determined by the balance of interactions. Long tail of the black curve, which extends to the few times longer distances than (R^{min }), exists, since some individuals are located inside voids far-away from all the neighbors.Typical size distributions to which the population evolves with the time are reproduced in Fig. 3. It shows three important size depending values. Figure 3a demonstrates an instant configuration of the velocities. Three different subgroups are marked by the same symbols as in Fig. 1. Figure 3b shows instant histogram of the animal’s masses in complete population and the histogram averaged over long run. They are plotted by the blue and black lines respectively. How mean velocities of the animals depend on their mass is shown in Fig. 3c. As above, the instant and accumulated distributions are plotted by the blue and black lines respectively.It can be noticed from the subplot (c), Fig. 3, that the averaged velocities correlate with the mass of the individuals. In principle, it could be expected according to the model rules, namely due to the nonlinear dependence of strength (f_{{_{jk} }}^{repuls} (overrightarrow {r}_{j} ,overrightarrow {r}_{k} )) as well as effective damping (eta_{j}) for the individual animals on their size (mass (m_{j})). However, taken into account multiple interactions in the system the result was not obvious in advance. The velocity increase is especially pronounced when the animals become large and are close to the final stage of its growth.Further, we would like to compare the numerical and experimentally obtained data. This comparison is presented in Fig. 4 by the curves and bar-plots respectively. From the field videos it is difficult to determine the masses (m_{j}) of the individuals, but we can estimate their volumes (V_{j}) which are approximately proportional to their masses (V_{j} sim m_{j}). Therefore, to compare numerical results with the measurements from videos we have plotted the distributions along the volume coordinate. The subplot (a) in Fig. 4 shows the histograms of the distances between nearest neighbors, and the subplot (b) presents the histograms of the volumes of the animals. The histograms obtained from simulations match very well to the histograms calculated using experimental data. Two sample Kolmogorov–Smirnov test does not reject the hypothesis, that the distributions obtained in experiment and in numerical simulations are from the same continues distribution for both the animal volumes (p = 0.59, D55,500 = 0.105) and for the distances between nearest neighbors (p = 0.84, D55,500 = 0.083).It is important to note that due to stability of the system the distributions independent on their initial shape converge to the same quasi-stationary histograms. To check this convergence, we started from two extreme populations. The first one consisted almost from the large individuals only and another one contained only small individuals. Time depending volume histograms were accumulated into the volume distribution over time gray-scale maps. In Fig. 5 the results for the two cases are shown in the subplots (a) and (b) respectively. Lighter gray color corresponds to the higher density. Both distributions attract to the practically the same final one with the time. The shapes of distributions, which have started from the two extreme initial distributions, almost coincide after long runs, while the distributions flanks shift in different directions as marked by the red and blue arrows in Fig. 6 and visualized in the supplementary movies “video_S6.mp4” and “video_S7.mp4” respectively. As in the static figures the instant and time averaged histograms are shown by the blue and black lines respectively. For the supplementary video “video_S6.mp4” the distribution is initially localized mainly near the right side of the interval and after some quick transient period “jumps” to the distribution close to the final histogram. Comparison with second video “video_S7.mp4” shows how both averaged distributions are attracted after a long simulation run to the similar distribution.Figure 5Gray-scale maps of time-volume histograms accumulated for two limit cases: starting from large and small animals, are shown in the subplots (a,b) respectively. Lighter gray color corresponds to the higher density. Dashed line marks common position of distribution maximums to which they converge during extremely long runs.Full size imageFigure 6Final volume distributions after long runs. It is seen that the maximums coincide already. Remaining trends of the histogram alterations for the distribution started from the large (open circles) and small (closed circles) individuals are marked by the red and blue arrows respectively.Full size imageDuring the system evolution some of the animals leave the population and some new appear. One can accumulate the numbers of died and survived animals for each given time moment. Accumulated numbers of died and survived animals are plotted in Fig. 7a by black and red curves respectively. The balance between these numbers depends on the relation of all the model parameters. Here the same set of the parameters was used as for simulations presented in Figs. 1, 2 and 3.Figure 7Time depending integral populations: (a) accumulated numbers of died (N1) and survived (N2) animals (black and red curves); (b) variations of the numbers of animals in subpopulations of small (Nc1), medium (Nc2) and large (Nc3) animals shown by blue, green and red curves respectively.Full size imageThe curves presenting the integral number of animals are relatively smooth, Fig. 7a, yet the instant numbers of animals in small, medium and large subpopulations shown by blue, green and red curves respectively in Fig. 7b vary much stronger, because at every particular time moment, their numbers are relatively small and the fluctuations are strong. From the figure, it is also obvious that during initial transient time interval the number of small individuals do not exceeds the number of two other populations. It is natural and was expected, because all newborn animals are small. Therefore, subpopulation with small size prevails at stationary stage.As we already announced nontrivial transformation of the spatial pattern and modification of all the resulting distributions can appear due to the “waves of fear”. Such a wave can appear when the animals almost synchronously try to escape from the shoreline if they are being scared by somebody moving along. Mathematically such a wave can be provoked by an additional repulsion force acting from one side of the simulation area to the distance much longer than already incorporated short range reflection from the boundary (Eq. 6). This force influences all the members of the array but appears for relatively limited time. It can be added to the model interaction occurring either randomly or periodically. Second variant is more regular and preferable and allows accumulate statistics faster.Analytically this force has the same form as (f_{j}^{bound}) with the same boundary position (overrightarrow {r}_{bound}):$$f_{j}^{wave} = B_{j}^{wave} exp left[ { - left| {frac{{overrightarrow {r}_{j} - overrightarrow {r}_{bound} }}{{R_{{}}^{wave} }}} right|} right],$$ (8) but it has larger amplitude (B_{j}^{wave} > B_{j}^{bound}) and much longer distance of the exponential decay (R_{{}}^{wave} > > R_{{}}^{bound}). So, it influences not only the individuals occasionally attempting to cross the boundary but almost all others too, even if they are currently relatively far from the boundary.The video in supplementary material “video_S8.mp4” illustrates typical behavior of the system at presence of the “waves of fear” occasionally observed in original sequences (video_S1.avi). One can see how bigger, stronger and faster animals escape quicker and further from the dangerous boundary, while some smallest ones do not react so quickly and remain almost near to it. Food deposition is not correlated with the “waves of fear”. As result, the food deposed near the dangerous boundary is consumed mostly by small animals, since that boundary region is practically depopulated by bigger animals. While stronger individuals return to the shore, such food may be already consumed by the weaker individuals.Static image of such “shifted” to one side pattern for a moment when food and weak members of the array are presented in one side of the area while absolute majority of the population is collected in a center or closer to another side of the area is reproduced in Fig. 8.Figure 8Shifted pattern at the presence of the “wave of fear”. Some food was just deposed inside the depopulated area near to the dangerous boundary. Some small individuals are in the vicinity to the food, while absolute majority of the population is localized in a center or near opposite boundary. The symbols have the same meaning as in Fig. 1.Full size imageBeing regularly applied the “waves of fear” in principle leads to a redistribution of the food and animals in the space. It is quite expected that the density of the animals will spend some time in the areas distanced from the dangerous wall. As result, the food will longer remain available in the places close to the dangerous wall. We have accumulated these densities during long runs at different parameters and obtained such shifted distributions. Typical distributions of the food and animals for the set of parameters which were used in simulations presented in the previous figures are shown in Fig. 9.Figure 9Density distributions of the food and animals at presence of regularly applied “waves of fear”. The densities of the food and animals are shown by the solid black and dashed magenta curves respectively. The curves for the densities of large, middle and small individuals are plotted by the same colors as above. The density shift with increasing size is marked by the arrow.Full size imageMean density of the food is shown by the solid black curve. Total density of the animals is presented by the dashed magenta curve. These densities are anti-correlated. Besides, there is a fine structure of the partial densities of the animals. The curves for the densities of large, middle and small individuals are plotted by the same colors as above. Strongest correlation between the size of the animals and their spatial density is observed in interval between two dash-dotted lines. The density shift with increasing size is marked by the arrow. In close proximity to the dangerous wall, where mean density of the food is abnormally high, the density shift dependence on the individuals size change the direction. Number of the small individuals reduces in average, because they quickly move to the category of middle or even large individuals. More

1.Stork, N. E., McBroom, J., Gely, C. & Hamilton, A. J. New approaches narrow global species estimates for beetles, insects, and terrestrial arthropods. PNAS 112, 7519–7752 (2015).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
2.Sømme, L. & Block, W. Adaptations to alpine and polar environments in insects and other terrestrial arthropods. In Insects at Low Temperature (eds Lee, R. E. & Denlinger, D. L.) (Springer, 1991).
Google Scholar 
3.Springer, M. Marine insects. In: Marine Biodiversity of Costa Rica, Central America. (ed. Wehrtmann, I. S. & Cortés, J.) Monographiae Biologicae, vol 86. (Springer, 2009).4.Cebeci, H. H. et al. The wood boring insects (Coleoptera: Cerambycidae and Buprestidae) recorded as the new pests for acer Undulatum Pojark from the Babadag mountain (SW Turkey). Fresenius Environ. Bull. 27, 9325–9328 (2018).CAS 

Google Scholar 
5.Becker, G. Untersuchungen über die Ernährungsphysiologie der Hausbockkäferlarven. Z. Vgl. Physiol. 29, 315–388 (1942).CAS 
Article 

Google Scholar 
6.Nardi, J. B., Mackie, R. I. & Dawson, J. O. Could microbial symbionts of arthropod guts contribute significantly to nitrogen fixation in terrestrial ecosystems?. J. Insect Physiol. 48, 751–763 (2002).CAS 
PubMed 
Article 

Google Scholar 
7.Wermelinger, B., Duelli, P. & Obrist, M. K. Dynamics of saproxylic beetles (Coleoptera) in windthrow areas in alpine spruce forests. For. Snow Landsc. 77, 133–148 (2002).
Google Scholar 
8.Evans, H. F., Moraal, L. G. & Pajares, J. A. Biology, ecology and economic importance of Buprestidae and Cerambycidae. In: Bark and Wood Boring Insects in Living Trees in Europe, a Synthesis (ed. Lieutier, F., Day, K. R., Battisti, A., Grégoire, J. C. & Evans, H. F.) 447–474 (Springer, 2004).9.FAO (Food and Agriculture Organization of the United Nations) Global review of forest pests and diseases. FAO Forestry Paper 156 (2009).10.Offenberg, J., Nielsen, J. S. & Damgaard, V. Wood Ant (Formica polyctena) services and disservices in a danish apple plantation. Sociobiology 66, 247–256 (2019).Article 

Google Scholar 
11.Cline, A., Ivie, M. A., Bellamy, C. L. & Scher, J. A Resource for Wood Boring Beetles of the World: Wood Boring Beetle Families, Lucid v. 3.4. USDA/APHIS/PPQ Center for Plant Health Science and Technology, Montana State University, and California Department of Food and Agriculture. http://www.lucidcentral.org/keys/v3/WBB (2009).12.Rudinski, J. A. Ecology of Scolytidae. Annu. Rev. Entomol. 7, 327–348 (1962).Article 

Google Scholar 
13.Huber, D. P. W., Aukema, B. H., Hodgkinson, R. S. & Lindgren, B. S. Successful colonization, reproduction, and new generation emergence in live interior hybrid spruce, Picea engelmannii x glauca, by mountain pine beetle, Dendroctonus ponderosae. Agric. For. Entomol. 11, 83–89 (2009).Article 

Google Scholar 
14.O’neill, K. M., Fultz, J. E. & Ivie, M. A. Distribution of adult Cerambycidae and Buprestidae (Coleoptera) in a subalpine forest under shelterwood management. Coleopt. Bull. 62(1), 27–36 (2008).Article 

Google Scholar 
15.Leather, S. R., Day, K. R. & Salisbury, A. N. The biology and ecology of the large pine weevil, Hylobius abietis (Coleoptera: Curculionidae): A problem of dispersal?. Bull. Entomol. Res. 89, 3–16 (1999).Article 

Google Scholar 
16.Nordlander, G., Bylund, H., Örlander, G. & Wallertz, K. Pine weevil population density and damage to coniferous seedlings in a regeneration area with and without shelterwood. Scand. J. For. Res. 18, 438–448 (2003).Article 

Google Scholar 
17.Nordlander, G., Hellqvist, C., Johansson, K. & Nordenhem, H. Regeneration of European boreal forests: Effectiveness of measures against seedling mortality caused by the pine weevil Hylobius abietis. For. Ecol. Manag. 262(12), 2354–2363 (2011).Article 

Google Scholar 
18.Bale, J. S. et al. Herbivory in global climate change research: Direct effects of rising temperature on insect herbivores. Glob. Change Biol. 8, 1–16 (2002).ADS 
Article 

Google Scholar 
19.Kahuthia-Gathu, R., Kirubi, D. T. & Gitonga, D. Composition and abundance ofwood-boring beetles of Acacia xanthophloea and their associated natural enemies in Thika, Kenya. J. Asia Pac. Biodivers. 11, 248–254 (2018).Article 

Google Scholar 
20.Holmes, T. & Koch, F. Bark beetle epidemics, life satisfaction, and economic well-being. Forests 10, 696 (2019).Article 

Google Scholar 
21.Adlung, K. G. A critical evaluation of the European research on use of red wood ants (Formica rufa group) for the protection of forests against harmful insects. Z. Angew. Entomol. 57, 167–189 (1966).Article 

Google Scholar 
22.Goldazarena, A., Romón, P. & López, S. Bark beetles control in forests of Northern Spain. Integrated Pest Management and Pest Control—Current and Future Tactics, Dr. Sonia Soloneski (Ed), ISBN: 978-953-51-0050-8, In Tech. http://www.intechopen.com/books/integrated-pest-management-and-pest-control-current-and-future-tactics/bark-beetles-control-in-forests-of-northern-spain (2012).23.IUCN 2020. The IUCN Red List of Threatened Species. Version 2020-1. https://www.iucnredlist.org (2020).24.Abdullah, H., Darvishzadeh, R., Skidmore, A. K., Groen, T. A. & Heurich, M. European spruce bark beetle (Ips typographus, L.) green attack affects foliar reflectance and biochemical properties. Int. J. Appl. Earth Obs. 64, 199–209 (2018).Article 

Google Scholar 
25.Bentz, B. J. & Jӧnsson, A. M. Modeling bark beetle responses to climate change. In Bark Beetles: Biology and Ecology of Native and Invasive Species (eds Vega, F. & Hofstetter, R.) 533–553 (Elsevier Academic Press, 2015).Chapter 

Google Scholar 
26.Linder, M. et al. Climate change impacts, adaptive capacity, and vulnerability of European forest ecosystems. For Ecol. Manag. 259, 698–709 (2010).Article 

Google Scholar 
27.Marini, L. et al. Climate drivers of bark beetle outbreak dynamics in Norway spruce forests. Ecography 40, 1426–1435 (2017).Article 

Google Scholar 
28.Schowalter, T. D. Insect ecology: an ecosystem approach, 3rd ed. (ed. Schowalter, T. D.) (Academic Press, 2011).29.Frizzi, F., Masoni, A., Quilghini, G., Ciampelli, P. & Santini, G. Chronicle of an impact foretold: The fate and effect of the introduced Formica paralugubris ant. Biol. Invasions 20, 3575–3589 (2018).Article 

Google Scholar 
30.Leong, M. et al. The habitats humans provide: Factors affecting the diversity and composition of arthropods in houses. Sci. Rep. 7, 15347 (2017).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
31.Schweiger, O. et al. Quantifying the impact of environmental factors on arthropod communities in agricultural landscapes across organisational levels and spatial scales. J. Appl. Ecol. 42, 1129–1139 (2005).Article 

Google Scholar 
32.Trigos-Peral, G. et al. Three categories of urban green areas and the effect of their different management on the communities of ants, spiders and harvestmen. Urban Ecosyst. 23, 803–818 (2020).Article 

Google Scholar 
33.Bellard, C., Bertelsmeier, C., Leadley, P., Thuiller, W. & Courchamp, F. Impacts of climate change on the future of biodiversity. Ecol. Lett. 15, 365–377 (2012).PubMed 
PubMed Central 
Article 

Google Scholar 
34.Czechowski, W., Radchenko, A., Czechowska, W. & Vepsäläinen, K. The ants of Poland with reference to the myrmecofauna of Europe. Fauna Poloniae 4 (Natura optima dux Foundation, 2012).35.Seidl, R. et al. Forest disturbances under climate change. Nat. Clim. Chang. 7, 395–402 (2017).ADS 
PubMed 
PubMed Central 
Article 

Google Scholar 
36.Sommerfeld, A. et al. Patterns and drivers of recent disturbances across the temperate forest biome. Nat. Commun. 9, 4355 (2018).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
37.Sorvari, J. & Hakkarainen, H. Wood ants are wood ants: Deforestation causes population declines in the polydomous wood ant Formica aquilonia. Ecol. Entomol. 32, 707–711 (2007).Article 

Google Scholar 
38.Roura-Pascual, N. et al. Geographic potential of Argentine ants (Linepithema humile Mayr) in the face of global climate change. Proc. R. Soc. Lond. B. 271(1557), 2527–2535 (2004).Article 

Google Scholar 
39.Kroll, J. C. & Fleet, R. R. Impact of woodpecker predation on over-wintering within-tree populations of the southern pine beetle (Dendroctonus frontalis). In The Role of Insectivorous Birds in Forest Ecosystems (eds Dickson, J. G. et al.) 269–281 (Academic Press, 1979).Chapter 

Google Scholar 
40.Khanday, A. L., Sureshan, P. M., Buhroo, A. A., Ranjith, A. P. & Tselikh, E. Pteromalid wasps (Hymenoptera: Chalcidoidea) associated with bark beetles, with the description of a new species from Kashmir, India. J. Asia Pac. Biodivers. 12, 262–272 (2019).Article 

Google Scholar 
41.Reeve, J. D. Predation and bark beetle dynamics. Oecologia 112, 48–54 (1997).ADS 
PubMed 
Article 

Google Scholar 
42.Curtsdotter, A. et al. Ecosystem function in predator–prey food webs confronting dynamic models with empirical data. J. Anim. Ecol. 88, 196–210 (2018).PubMed 
Article 

Google Scholar 
43.Choate, B. & Drummond, F. A. Ants as biological control agents in agricultural cropping systems. Terr. Arthropod Rev. 4, 157–180 (2011).Article 

Google Scholar 
44.Maňák, V., Björklund, N., Lenoir, L. & Nordlander, G. The effect of red wood ant abundance on feeding damage by the pine weevil Hylobius abietis. Agric. For. Entomol. 17, 57–63 (2015).Article 

Google Scholar 
45.Robinson, E. J. H., Stockan, J. A. & Glenn, R. I. Wood ants and their interaction with other organisms. In Wood Ant Ecology and Conservation (eds Stockan, A. & Robinson, E. J. H.) 177–206 (Cambridge Universsity Press, 2016).
Google Scholar 
46.Sorvari, J. Foraging distances and potentiality in forest pest insect control: An example with two candidate ants (Hymenoptera: Formicidae). Myrmecol. News 12, 211–215 (2009).
Google Scholar 
47.Gößwald, K. Die Waldameise im Ökosystem Wald, ihr Nutzen und ihre Hege. Aula Verlag, Wiesbaden, 510 pp (1990). (shortened editon is 978-3-89104-755-2 (ISBN)).48.Juhász, O. et al. Large- and small-scale environmental factors drive distributions of ant mound size across a latitudinal gradient. Insects 11, 350 (2020).PubMed Central 
Article 
PubMed 

Google Scholar 
49.Juhász, O. et al. Consequences of climate change-induced habitat conversions on red wood ants in a Central European mountain: a case study. Animals 10, 1677 (2020).PubMed Central 
Article 
PubMed 

Google Scholar 
50.Risch, A., Ellis, S. & Wiswell, H. Where and why? Wood ant population ecology. In Wood Ant Ecology and Conservation (eds Stockan, A. & Robinson, E. J. H.) 81–105 (Cambridge Universsity Press, 2016).
Google Scholar 
51.R Core Team. R: A Language and Environment for Statistical Computing. https://www.R-project.org/ (2020).52.Oksanen, J., Blanchet, F. G., Friendly, M., Kindt, R., Legendre, P., McGlinn, D., Minchin, P. R., O’Hara, R. B., Simpson, G. L., Solymos, P., Stevens, M. H. H., Szoecs, E., Wagner, H. Vegan: Community Ecology Package (2019).53.Roberts, D.W. (2010) labdsv: Ordination and Multivariate Analysis for Ecology. http://cran.r-project.org/package=labdsv54.Bates, D., Maechler, M., Bolker, B. & Walke, S. Fitting linear mixed-effects models using lme4. J. Stat. Softw. 67(1), 1–48 (2015).Article 

Google Scholar 
55.Hartig, F. DHARMa: residual diagnostics for hierarchical (multi-level/mixed) regression models. R package version 0.1. 0. CRAN/GitHub (2016).56.Lindén, A. & Mäntyniemi, S. Using the negative binomial distribution to model overdispersion in ecological count data. Ecology 97, 1414–1421 (2011).Article 

Google Scholar 
57.Wickham, H. ggplot2: Elegant Graphics for Data Analysis (Springer-Verlag, 2016).MATH 
Book 

Google Scholar 
58.Schowalter, T. D. Ecology and management of bark beetles (Coleoptera: Curculionidae: Scolytinae) in Southern Pine Forests. J. Integ. Pest Mngmt. 3, 1–7 (2012).Article 

Google Scholar 
59.Hölldobler, B. & Wilson, E. O. The Ants (Harvard University Press, 1990).Book 

Google Scholar 
60.Siemann, E., Haarstad, J. & Tilman, D. Dynamics of plant and arthropod diversity during old field succession. Ecography 22, 406–414 (1999).Article 

Google Scholar 
61.Schlyte, F. & Anderbrandt, O. Competition and niche separation between two bark beetles: Existence and mechanisms. Oikos 68, 437–447 (1993).Article 

Google Scholar 
62.Raffa, K. F. et al. Cross-scale drivers of natural disturbances prone to anthropogenic amplification: The dynamics of bark beetle eruptions. BioSci 58, 501–517 (2008).Article 

Google Scholar 
63.Drees, B. M., Jackman, J. A. & Merchant, M. E. Wood-Boring Insects of Trees and Shrub. Texas Agricultural Extension Service. http://hdl.handle.net/1969.1/160397 (1994).64.Ray, C. et al. Patterns of woodboring beetle activity following fires and bark beetle outbreaks in montane forests of California, USA. Fire Ecol. 15, 21 (2019).Article 

Google Scholar 
65.Risch, S. J. & Carroll, C. R. The ecological role of ants in two Mexican agroecosystems. Oecologia 55, 114–119 (1982).ADS 
PubMed 
Article 

Google Scholar 
66.Maák, I. et al. Behaviours indicating cannibalistic necrophagy in ants are modulated by the perception of pathogen infection level. Sci. Rep. 10, 17906 (2020).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
67.Doležal, P. & Sehnal, F. A simple method for the detection of imaginal diapause in beetles. J. Appl. Entomol. 131(3), 221–223 (2007).Article 

Google Scholar 
68.Horstmann, K. Untersuchungen ueber den Nahrungserwerb der Waldameisen (Formica polyctena Foerster) im Eichenwald II Abhaengigkeit vom Jahresverlauf und vom Nahrungsangebot. Oecologia (Berlin) 8, 371–390 (1972).ADS 
Article 

Google Scholar 
69.Overbeck, M. & Schmidt, M. Modelling infestation risk of Norway spruce by Ips typographus (L.) in the Lower Saxon Harz Mountains (Germany). For. Ecol. Manag. 266, 115–125 (2012).Article 

Google Scholar 
70.Ayre, G. Response to movement by Formica polyctena Forst. Nature 199, 405–406 (1963).ADS 
Article 

Google Scholar 
71.Thom, D., Seidl, R., Steyrer, G., Krehan, H. & Formayer, H. Slow and fast drivers of the natural disturbance regime in Central European forest ecosystems. For. Ecol. Manag. 307, 293–302 (2013).Article 

Google Scholar 
72.Czechowski, W. & Vepsäläinen, K. Territory size of wood ants (Hymenoptera: Formicidae): A search for limits of existence of Formica polyctena Först., an inherently polygynic and polycalic species. Ann. Zool. 59, 179–187 (2009).Article 

Google Scholar  More

Field survey and samplingField data were collected from 130 study sites spanning a latitudinal gradient of 35.89−50.70° N and a longitudinal gradient of 76.62−122.41°E and covering five provinces across the temperate region in northern China (Xinjiang Autonomous Region, Qinghai Province, Gansu Province, Ningxia Autonomous Region, and Inner Mongolia Autonomous Region; Fig. 2a). Locations for the field study target natural drylands, delineated as areas with aridity level above 0.35 (ref. 30), and represent a large aridity gradient including dry-subhumid (N = 12), semiarid (N = 42), arid (N = 56), and hyperarid (N = 20) regions (Fig. 2a), which are highly vulnerable to expected increases in aridity with human activity and climate change33,71. The aridity level of each site was calculated as 1 – AI, where AI is the ratio of precipitation to potential evapotranspiration38. We obtained AI from the Global Aridity Index and Potential Evapotranspiration Climate database (https://cgiarcsi.community/). The selection of the field sites aimed to minimize the potential impacts of human activity and other disturbances on soil, vegetation, and geomorphological characteristics based on the following three criteria: (i) sites were at least 1 km away from major roads and >50 km from human habitations; (ii) sites were under pristine or unmanaged conditions without visible signs of domestic animal grazing, grass/wood collection, engineering restoration plantings, and infrastructure construction; and (iii) the soil was dry without experiencing rainfall events for at least 3 days prior to sampling. Collectively, our field survey involved a wide range of the abiotic and biotic features of dryland ecosystems across northern China. These sites encompass the 14 soil types, i.e., arenosols, calcisols, cambisols, chernozems, fluvisols, gleysols, greyzems, gypsisols, kastanozems, leptosols, luvisols, phaeozems, solonchaks, and solonetz, and the four main vegetation types44, i.e., typical grassland (dominated by Stipa spp., Leymus spp., Cleistogenes spp., and Agropyron spp.), desert grassland (dominated by Stipa spp., Cleistogenes spp., Suaeda spp., and Artemisia spp.), alpine grassland (dominated by Stipa spp., Leymus spp., Carex spp., and Festuca spp.), and desert (dominated by Reaumuria spp., Salsola spp., Calligonum spp., and Nitraria spp.). Elevation, mean annual temperature, and mean annual precipitation (1970–2000; https://www.worldclim.org/) of the sites varied from 204 to 3,570 m a.s.l. (mean, 1,294 m a.s.l.), from –4.3 to 12.8 °C (mean, 5.0 °C), and from 21 to 453 mm (mean, 195 mm), respectively (Supplementary Table 1).Field sampling was conducted between June and September from 2015 to 2017 (each site was visited once over this period) following well-established standardized protocols as described in refs. 13,34. In brief, three 30 m × 30 m quadrats were established at each site to represent the local vegetation and soil types that covered an area of no less than 10,000 m2. The cover of perennial vegetation was estimated and all perennial plant species were listed by walking steadily along four 1.5 m × 30 m parallel transects (spaced 8 m apart) located within each quadrat using the belt transect method72. Site-level estimate for perennial plant cover was obtained by averaging the values measured in the 12 transects established. After vegetation survey, we located five 1 m × 1 m (for typical grassland, desert grassland, and alpine grassland) or five 5 m × 5 m (for desert) plots within each quadrat (at each corner and the center of the quadrat) to measure site-level plant aboveground and root biomass (g m−2). In each 1 m × 1 m plot, all grasses and dwarf shrubs were harvested to ground level for measurement of aboveground biomass. Five soil cores (7 cm diameter; 0–40 cm depth) per 1-m2 plot were collected randomly, and the roots were removed using a 1-mm sieve and washed cleanly to measure root biomass. All shoot and root samples were dried to constant weight at 65 °C. In each 5 m × 5 m plot, we recorded the number of individuals per dominant shrub species and canopy cover and height of each individual, thereby estimating aboveground and root biomass according to the allometric models developed in previous studies that were conducted in the same regions as sampled here (see Supplementary Table 9 for details). Based on these measurements, we further estimated BNPP. However, BNPP is typically difficult to observe and measure, especially over large spatial scales and environmental gradients as in this study, because the root system is subject to simultaneous growth and turnover73,74. Across our survey areas, ~77–98% of the precipitation occurs between June and September (during the peak-growing season) corresponding to the period of the highest plant above- and belowground biomass34,35,41,75. Therefore, we argue that BNPP can be estimated approximately at each site by the following equation:$$frac{{{{{{rm{Aboveground}}}}}},{{{{{rm{biomass}}}}}}}{{{{{{rm{Root}}}}}},{{{{{rm{biomass}}}}}}}cong frac{{{{{{rm{Aboveground}}}}}},{{{{{rm{net}}}}}},{{{{{rm{primary}}}}}},{{{{{rm{productivity}}}}}},({{{{{rm{ANPP}}}}}})}{{{{{{rm{BNPP}}}}}}}$$
(1)
where both aboveground and root biomass are site-level measurements (g m−2). We used normalized difference vegetation index (NDVI) as a metric for ANPP as explained in recent studies in drylands14,33,70. NDVI data were obtained from the moderate resolution imaging spectroradiometer aboard NASA’s Terra satellites (https://neo.sci.gsfc.nasa.gov/). We used the average NDVI values during our sampling dates as a proxy for ANPP at the site level as described in ref. 14.Five soil cores (0–20 cm depth) per quadrat were then taken randomly under the canopies of the dominant perennial plant species and in bare areas devoid of perennial vegetation, respectively, and then were mixed as one sample for vegetation areas and the other sample for bare ground. When more than one dominant perennial plant species was observed, another three composite samples were collected under the canopies of co-dominant perennial plant species. All vegetation and soil surveys were carried out during the wet season (June to September) when biological activity and productivity are maximal; as such, we do not expect the different sampling times and years or seasonality to be a major factor influencing our conclusions. Collectively, 6–21 soil samples per site were collected, and in total 864 samples were taken and analyzed for each of the seven individual soil functions (see below) and multifunctionality. All soil functions evaluated in the field study were calculated at site level by using a weighted average of the mean values observed in vegetated areas and bare ground by their respective cover13,14,38. After field sampling, the visible pieces of plant material were removed carefully from the soil, which was sieved and divided into three portions. The first portion was air-dried and used for soil organic C, total N, total P, available P, and pH analyses. The second portion was immediately mixed with 2 M KCl and stored at 4 °C for soil ammonium and nitrate analyses. The third portion was immediately frozen at –80 °C for assessing soil microbial diversity.Microcosm experimentIn addition to the large-scale field study described above, we manipulated soil water availability in a microcosm experiment to evaluate the linkages between moisture content, soil microbial diversity, and multifunctionality. It is important to note that our intention is not to directly compare results between these two different approaches [i.e., in the field, measures of soil functions are related to nutrient pools, which we use to associate soil multifunctionality with both plant and soil microbial diversity, whereas in the microcosm experiment the measures of soil functions are related to process rates such as respiration rate and key enzyme activities (see below), which we use to associate soil multifunctionality with microbial diversity in the absence of plants]. Rather, by using an experimental microcosm approach, we aimed to complement the field study and thus further verify the potential increases in aridity to alter the relationship between soil microbial diversity and multifunctionality in the absence of plants. In parallel with the sampling protocols described above, we collected a greater mass of soil (c. 30 kg) under vegetation canopies from one site [i.e., Jingtai country (37.40°N, 104.26°E; Gansu Province, China)]. Soil type, mean annual temperature, mean annual precipitation, and aridity level (1970–2000; https://www.worldclim.org/) of the site is calcisols, 7.9 °C, 205 mm and 0.81, respectively. Following field sampling, the soil was stored immediately at 4 °C until subsequent processing in the laboratory.In brief, a total of 30 experimental microcosms composed of 10 moisture levels with three replicates were established under sterile conditions in a closed incubation chamber (Supplementary Fig. 1a). Each microcosm was filled with 1 kg of soil. These microcosms were incubated at 18.5 °C [the annual mean land surface temperature (1981–2010) for the sampling site; http://data.cma.cn/en], and moisture contents were adjusted and artificially maintained at the ten levels respectively equivalent to 3, 5, 8, 10, 20, 40, 60, 80, 100, and 120% field capacity (27.6%) during the duration of the experiment for 30 days. The corresponding moisture content (%) measured at the end of the experiment varied from 2.03 ± 0.034 to 33.57 ± 1.94, which matched well with differences in moisture conditions among a subset of field soil samples (N = 521; Supplementary Fig. 1b). After incubation, the soil was removed from each microcosm; a portion of the soil was immediately frozen at –80 °C for molecular analysis, and the other fraction was air-dried, sieved, and stored at –20 °C for assessing multiple soil functions as described below.DNA extraction, PCR amplification, and amplicon sequencingFor both the field and experimental studies, we assessed the diversity of soil archaea, bacteria, and fungi using Illumina-based sequencing. Genomic DNA was extracted from 0.5 g of each defrosted soil sample (N = 864 for the field study and N = 30 for the experimental study) using the PowerSoil® DNA Isolation Kit (MO BIO Laboratories, USA) according to the manufacturer’s instructions. For our field study, extracted DNA was pooled at site level, ultimately resulting in 130 composite DNA samples under canopies of vegetation and in bare ground, respectively. Pooling DNA samples may outperform the commonly used method that extracts genomic DNA from mixed soil samples, which could remove large amounts of information on the diversity of soil microorganisms14,22. Negative controls (deionized H2O in place of soil) underwent identical procedures during the extraction to ensure zero contamination in downstream analyses.The V3−V5 regions of the archaeal 16S rRNA gene were amplified using the primer pair Arch344F and Arch915R. Thermal conditions were composed of an initial denaturation of 3 min at 95 °C, ten cycles of touchdown PCR (95 °C for 30 s, annealing temperatures starting at 60 °C for 30 s then decreasing 0.5 °C per cycles, and 72 °C for 1 min), followed by 25 cycles at 95 °C for 30 s, 55 °C for 30 s, and 72 °C for 1 min, with a final extension at 72 °C for 10 min. The primer pair 338F and 806R was used for amplification of the V3−V4 regions of the bacterial 16S rRNA gene. Thermocycling conditions consisted of 3 min at 95 °C and then subjected to 30 amplification cycles of 30 s denaturation at 95 °C, 30 s annealing at 55 °C, followed by 72 °C for 45 s, and a final extension of 72 °C for 10 min. The fungal internal transcribed spacer (ITS) region 1 was amplified using the primer pair ITS1F and ITS2. The amplification conditions involved denaturation at 95 °C for 3 min, 35 cycles of 94 °C for 1 min, 51 °C for 1 min, and 72 °C for 1 min and a final extension at 72 °C for 10 min. Details of primers for each microbial taxa were given in Supplementary Table 10. These primers contained variable length error-correcting barcodes unique to each sample. All amplification reactions were performed in a total volume of 20 μl containing 4 μl of 5× FastPfu Buffer, 2 μl of 2.5 mM dNTPs, 0.8 μl of both the forward and reverse primers, 10 ng of template DNA, and 0.4 μl of FastPfu DNA Polymerase (TransGen Biotech., China). To mitigate individual PCR reaction biases each sample was amplified in triplicate and pooled together. All PCRs were done with the ABI GeneAmp® 9700 Thermal Cycler (Thermo Fisher Scientific, USA). PCR products were evaluated on 2.0% agarose gel with ethidium bromide staining to ensure correct amplicon length, and were gel-purified using the AxyPrep DNA Gel Extraction Kit (Axygen Biosciences, USA). Purified amplicons were combined at equimolar concentrations and paired-end sequenced (2 × 300 bp) on an Illumina MiSeq platform (Illumina, USA) at the Majorbio Bio-pharm Technology Co., Ltd. (Shanghai, China) according to standard protocols.Sequence processingInitial sequence processing was conducted with the QIIME pipeline76. Briefly, reads were quality-trimmed with a threshold of an average quality score higher than 20 over 10 bp moving-window sizes and a minimum length of 50 bp. Paired-end reads with at least 10 bp overlap and 2 indicate that the models are different; Supplementary Table 2]. We further assessed whether soil multifunctionality responded more rapidly to aridity than did any individual soil functions. To this end, we explored the presence of aridity thresholds for those relationships that were better fitted by nonlinear regressions (Fig. 2b–i) using the standard protocols developed in ref. 33. The presence of an aridity threshold means that once an aridity level is reached, a given variable either changes abruptly its value (i.e., discontinuous threshold) or its relationship with aridity (i.e., continuous threshold). Hence, a lower aridity threshold indicates that a given variable is more vulnerable to increasing aridity than are others33. We further fitted step (a linear regression that modifies only intercept at a given aridity level) and stegmented (showing changes both in intercept and slope at a given aridity level) regressions for the determination of discontinuous thresholds, and segmented (exhibiting changes only in slope at a given aridity level) regressions for continuous thresholds. Each of these models yields a change point (i.e., threshold) describing the aridity level that evidences the shift in a given nonlinear relationship evaluated. We also used AIC to choose the best threshold model and the corresponding threshold in each case (Supplementary Table 2).We then employed analysis of variance based on type-I sum of squares in a linear mixed-effects model (Eq. (2); Table 1) to test the relationships between the multiple biotic (BNPP, plant species richness, and the soil microbial diversity index) and abiotic (aridity, soil pH, and soil clay content) factors and soil multifunctionality:$${{{{{rm{Soil}}}}}},{{{{{rm{multifunctionality}}}}}} sim; {{{{{rm{Year}}}}}}+{{{{{rm{Plant}}}}}},{{{{{rm{species}}}}}},{{{{{rm{richness}}}}}}\ quad+,{{{{{rm{Soil}}}}}},{{{{{rm{microbial}}}}}},{{{{{rm{diversity}}}}}},{{{{{rm{index}}}}}}\ quad+{{{{{rm{Plant}}}}}},{{{{{rm{species}}}}}},{{{{{rm{richness}}}}}}times {{{{{rm{Soil}}}}}},{{{{{rm{microbial}}}}}},{{{{{rm{diversity}}}}}},{{{{{rm{index}}}}}}+{{{{{rm{Aridity}}}}}}\ quad+,{{{{{rm{Aridity}}}}}}times {{{{{rm{Plant}}}}}},{{{{{rm{species}}}}}},{{{{{rm{richness}}}}}}\ quad+,{{{{{rm{Aridity}}}}}}times {{{{{rm{Soil}}}}}},{{{{{rm{microbial}}}}}},{{{{{rm{diversity}}}}}},{{{{{rm{index}}}}}}+{{{{{rm{BNPP}}}}}}+{{{{{rm{Soil}}}}}},{{{{{rm{pH}}}}}}+{{{{{rm{Soil}}}}}},{{{{{rm{clay}}}}}},{{{{{rm{content}}}}}}\ quad+,{{{{{rm{Elevation}}}}}}+{{{{{rm{Latitude}}}}}}+{{{{{rm{Longitude}}}}}}+(1|{{{{{rm{Soil}}}}}},{{{{{rm{type}}}}}})+(1|{{{{{rm{Vegetation}}}}}},{{{{{rm{type}}}}}})$$
(2)
where × indicates an interaction term. We obtained information on soil clay content (%) from the SoilGrids system (https://soilgrids.org/), and eliminated variation due to different sampling years by first entering the term “Year” into the statistical model41. The elevation, latitude, and longitude of the study sites were included to account for the spatial structure of our dataset13,70. To account for the similarities of soil and vegetation types among study sites we included “Soil type” and “Vegetation type” as random terms.We further simplified the Eq. (2) to focus only on the relationships between aridity, biodiversity, and soil multifunctionality (Eq. (3); Supplementary Fig. 5). We did so because excluding additional biotic and abiotic factors did not change qualitatively the main results presented here (Table 1 and Supplementary Fig. 5), and therefore we used the simplest model to test our hypotheses more clearly. Our simplified model was:$${{{{{rm{Soil}}}}}},{{{{{rm{multifunctionality}}}}}} sim {{{{{rm{Year}}}}}}+{{{{{rm{Plant}}}}}},{{{{{rm{species}}}}}},{{{{{rm{richness}}}}}}\ quad+,{{{{{rm{Soil}}}}}},{{{{{rm{microbial}}}}}},{{{{{rm{diversity}}}}}},{{{{{rm{index}}}}}}\ quad+,{{{{{rm{Aridity}}}}}}+{{{{{rm{Aridity}}}}}}times {{{{{rm{Plant}}}}}},{{{{{rm{species}}}}}},{{{{{rm{richness}}}}}}\ quad+,{{{{{rm{Aridity}}}}}}times {{{{{rm{Soil}}}}}},{{{{{rm{microbial}}}}}},{{{{{rm{diversity}}}}}},{{{{{rm{index}}}}}}\ quad+,{{{{{rm{Aridity}}}}}}times {{{{{rm{Plant}}}}}},{{{{{rm{species}}}}}},{{{{{rm{richness}}}}}}times {{{{{rm{Soil}}}}}},{{{{{rm{microbial}}}}}},{{{{{rm{diversity}}}}}},{{{{{rm{index}}}}}}\ quad+,(1|{{{{{rm{Soil}}}}}},{{{{{rm{type}}}}}})+(1|{{{{{rm{Vegetation}}}}}},{{{{{rm{type}}}}}})$$
(3)
To evaluate how the biodiversity–multifunctionality relationships varied along aridity gradients, we conducted a moving-window analysis as detailed in ref. 69. Briefly, we performed the linear mixed-effects model described in Eq. (3) for a subset window of 60 study sites with the lowest aridity values (this number of sites provided sufficient statistical power for our model), and repeated the same calculations as many times as sites remained (i.e., 70). We then bootstrapped the standardized coefficients of each fixed term within each subset window, which was matched to the average value of aridity across the 60 sites. We fitted linear and nonlinear regressions to the bootstrapped coefficients of biodiversity and its interaction with aridity along aridity gradients (Fig. 3a, b and Supplementary Table 2), and identified the aridity thresholds for the changes in the coefficients of biodiversity (Fig. 3a and Supplementary Table 2) using the same procedure already described above. To provide further support for the aridity thresholds identified here, we also assessed the significance of the bootstrapped standardized coefficients of biodiversity and its interaction with aridity at 95% confidence intervals for each subset window (Fig. 3e). Before fitting threshold regressions, we evaluated whether the variables followed either a unimodal or bimodal distribution using the fitgmdist function in MATLAB (The MathWorks Inc., USA). Our results showed that all variables used for threshold detection presented unimodal distributions (Supplementary Table 11), suggesting that the three threshold regressions mentioned above (i.e., segmented, step, and stegmented) are appropriate in all cases33. We used the chngpt and gam packages in R (http://cran.r-project.org/) to fit segmented/step/stegmented and GAM regressions, respectively. To further check the validity of the thresholds identified, we bootstrapped linear regressions at both sides of each threshold for each variable. We then used the nonparametric Mann–Whitney U-test to compare the slope and the predicted value evaluated before and after each threshold. In all cases, we found significant differences in both of these two parameters (Fig. 3c, d and Supplementary Figs. 2, 3, 6).Given a clear shift in the relationships between plant or microbial diversity and soil multifunctionality occurring at a threshold around an aridity level of 0.8 (Fig. 3), we further used OLS regressions to clarify the relationships between each component of plant or microbial diversity and soil multifunctionality in less and more arid regions separately, as well as across all sites (Fig. 4). To do so, we split our study sites into two groups: sites with aridity 0.8 (more arid regions; N = 76). Moreover, we fitted the mixed-effects model described in Eq. (2) for less and more arid regions separately to ensure the robustness of these bivariate correlations when accounting for multiple biotic and abiotic factors simultaneously, with the exception of using all components of microbial diversity metrics (i.e., soil archaeal, bacterial, and fungal richness) instead of the soil microbial diversity index in the models (Supplementary Table 4). All linear mixed-effects models were performed using the R package lme4. We used a variance inflation factor (VIF) to evaluate the risk of multicollinearity, and selected variables with VIF  0.05). Finally, we also used SEMs to compare the hypothesized direct and indirect relationships between moisture content, microbial diversity, and soil multifunctionality at low and high moisture levels (see an a priori model in Supplementary Figs. 17b and 31c, d). Test of goodness-of-fit for SEMs were same as described above. All the SEM analyses were conducted using AMOS 21.0 (IBM SPSS Inc., USA). Data and code used to perform above analyses are available in figshare98.Reporting SummaryFurther information on research design is available in the Nature Research Reporting Summary linked to this article. More

Collection of the samplesTen (10) wood pellet samples were purchased from a different location in B&H, of known suppliers from the market (supermarkets, garden shops, and gas stations). The samples were accompanied by a declaration describing that nine of them were originated from B&H, and one of them was from Italy. Characteristics of collected wood pellet samples (type of wood, energetic value, declared moisture, declared and determined ash amount) are listed in Table 1. All of the samples were analyzed for moisture and ash content. Additionally, in ash samples of mentioned wood pellets, heavy metal concentration (Cd, Co, Cr, Cu, Fe, Mn, Ni, Pb, and Zn) was determined.Table 1 Characteristic of analyzed samples wood pellets.Full size tableAll pellet samples were originated from B&H, purchased from different cities, often used for house heating, instead of sample S3 which was from Italy.Ash determination of wood biomass samplesThe wood pellet samples were oven-dried at 105 °C for 24 h. The content of ash was determined by gravimetric method according to the procedure published by Pan and Eberhardt18 as follows: pellet samples, 1 g (± 0.1 mg) of each was weighed into a previously annealed ceramic pot (m1) and burned in a muffle furnace (Nabertherm) for one hour at 300 °C, following by increasing the temperature to 400 °C for one hour more and then burning the samples for next six hours at 550 °C. The procedure is repeated until a constant mass (m2) was reached. The ash content is determined by the Eq. (1):$${text{Ash content}}, % = frac{{{text{(m}}_{2} – {text{m}}_{{1}} {)}}}{{{text{m}}_{{{text{sample}}}} }} times {100 }{text{.}}$$
(1)
Preparation of samplesThe chemical determinations of the heavy metals in wood pellet ashes (Table 2) were made by wet digestion by soaking the samples in 25 mL of 65% HNO3 in polytetrafluoroethylene (PTFE) vessels. After evaporation of the nitrogen oxides, the vessels were closed and allowed to react for 14 h at 80 °C, following by cooling to room temperature. Then, the digest was filtered, transferred to a 25 mL volumetric flask, and filled up with redistilled water to the mark. All samples and blank were prepared in three replicates19,20,21.Table 2 Heavy metal concentrations (mg kg−1 d.w.) in the wood pellet ashes.Full size tableHeavy metal analysisMetal analyses in ash samples of mentioned wood pellets were performed using a flame atomic absorption spectrometry (Varian AA240FS) for Mn, Fe, Pb, and Zn and graphite furnace (Varian AA240Z) for Cd, Co, Cr, Cu, and Ni. A blank probe was prepared using the same digestion method to avoid the matrix effect. Standard metal solutions used for the calibration graphs were prepared by diluting 1000 mg L−1 stock single-element atomic absorption standard solutions of Cd, Co, Cr, Cu, Fe, Mn, Ni, Pb, and Zn (Certipur Grade, Merck, Germany). Linear calibration graphs with correlation coefficients  > 0.99 were obtained for all analyzed metals. The accuracy of the method was evaluated using the standard reference materials: Fine Fly Ash (CTA-FFA-1, Institute of Nuclear Chemistry and Technology Poland) and Fly Ash from pulverized coal (BCR-038, Institute of reference materials and measurements-IRMM, Belgium). The obtained results were in the range of the reference materials. The detection limit (LOD) and limit of quantification (LOQ) for the nine analyzed metals were calculated based on Xb + 3 SDb and Xb + 10 SDb, respectively, where Xb is the mean concentration of the blank sample (n = 8) and SDb is the standard deviation of the blank for eight readings22. The values of the LOD were: Cd (0.61 µg L−1), Co (0.49 µg L−1), Cr (0.67 µg L−1), Cu (20.10 µg L−1), Fe (83.85 µg L−1), Mn (6.42 µg L−1), Ni (1.12 µg L−1), Pb (23.77 µg L−1), Zn (58.68 µg L−1), and LOQ values were: Cd (1.25 µg L−1), Co (1.41 µg L−1), Cr (1.42 µg L−1), Cu (47.66 µg L−1), Fe (111.2 µg L−1), Mn (16.14 µg L−1), Ni (2.70 µg L−1), Pb (47.73 µg L−1) and Zn (71.05 µg L−1).Pollution evaluationThe metal pollution index (MPI) as the geometric mean of the concentration of all metals found in ashes of wood samples was calculated by the following Eq. (2)23:$${text{MPI}} = left( {{text{C}}_{1} cdot {text{C}}_{2} cdot cdots {text{C}}_{{text{k}}} } right)^{{1/{text{k}}}} ,$$
(2)
where C1 is the concentration value of the first metal, C2 is the concentration value of the second metal, Ck is the concentration value of the kth metal.Evaluation of the presence and the grade of anthropogenic activity were demonstrated through the calculation of the enrichment factor (EF), widely used in environmental issues24. To understand which elements were relatively enriched in the different wood pellet ash samples, the heavy metal enrichment factor was calculated relative to soil values according to Eq. (3)25.$${text{EF}} = frac{{left( {frac{{{text{C}}_{{text{k}}} }}{{{text{E}}_{{{text{ref}}}} }}} right)_{{{text{ashes}}}} }}{{left( {frac{{{text{C}}_{{text{k}}} }}{{{text{E}}_{{{text{ref}}}} }}} right)_{{{text{soil}}}} }},$$
(3)
where Ck is the concentration of the element in the sample or the soil, Eref the concentration of the reference element used for normalization. A reference element is an element commonly stable in the soil characterized by the absence of vertical mobility and/or degradation phenomena. As in many studies as a reference element were Fe, Al, Mn, Sc, or total organic carbon used26,27. Therefore Fe has been chosen as reference material in this study. Iron is one of the major constituents of soil, as well as the average chemical constituent of the upper continental crust26.Health risk assessmentThe general exposure equations used in this study were adapted according to the US Environmental Protection Agency guidance28,29,30. The daily exposure (D) to heavy metals via wood pellet ash was calculated for the three main routes of exposure: (i) direct ingestion of ash particles (Ding); (ii) inhalation of suspended particles via mouth and nose (Dinh); and (iii) dermal absorption to skin adhered ash particles (Ddermal). Equations (4) to (6) were used to calculate exposure via ingestion, inhalation, and dermal route, respectively22,31.$${text{D}}_{{{text{ing}}}} = {text{ C }} cdot frac{{{text{ IngR }} cdot {text{ EF }} cdot {text{ ED}}}}{{{text{BW }} cdot {text{ AT}}}}{ } cdot {text{CF}}1{, }$$
(4)
$${text{D}}_{{{text{inh}}}} = {text{ C }} cdot frac{{{text{ InhR}} cdot {text{ EF }} cdot {text{ ED}}}}{{{text{PEF }} cdot {text{ BW }} cdot {text{ AT}}}}{, }$$
(5)
$${text{D}}_{{{text{dermal}}}} = {text{ C }} cdot frac{{{text{ SA }} cdot {text{ SL }} cdot {text{ABS }} cdot {text{EF }} cdot {text{ ED}}}}{{{text{BW }} cdot {text{ AT}}}}{ } cdot {text{CF}}1{, }$$
(6)

where c (mg kg−1) is the heavy metals concentrations in ash samples; IngR (mg day−1) is the conservative estimates of dust ingestion rates, 50 for adults, 200 for children30,32; InhR (m3 h−1) is the inhalation rate, 2.15 for adults, 1.68 for children32; EF (h year−1) is the exposure frequency, 1225 for adults and children22; ED (years) is the exposure duration, 70 for adults, 6 for children22; BW (kg) is the body weight, 80 for adults, 18.60 for children32; AT (days) is the averaging time, 25,550 for adults, 2190 for children22; PEF is the particle emission factor (m3 kg−1), 6.80 × 108 for adults and children31; SA (cm3) is the exposed skin area, 6840 for adults, 2550 for children32; SL (mg cm−2) is the skin adherence factor, 0.22 for adults, 0.27 for children32; ABS is the dermal absorption factor, 0.001 for adults and children31; CF1 is the unit conversation factor, 10–6 for adults and children22.The potential non-carcinogenic risk for each metal was estimated using the Hazard coefficient (HQ), as suggested by US EPA33. The HQ under various routes of exposure such as ingestion (HQing), inhalation (HQinh), and dermal (HQdermal) was calculated as a ratio of daily exposure (D) to reference dose of each metal (RfD) according to Eq. (7)32.$${text{HQ}}_{{text{k}}} = frac{{{text{D}}_{{text{k}}} }}{{{text{RfD}}}},$$
(7)

where k is ingestion, inhalation, or dermal route. The total hazard index (HI) of heavy metal for all routes of exposure was calculated as a sum of HQing, HQinh, and HQdermal as given in Eq. (8)34.$${text{HI}} = {text{ HQ}}_{{text{ing }}} + {text{ HQ}}_{{text{inh }}} + {text{ HQ}}_{{text{dermal }}} .$$
(8)
The carcinogenic risk (Risk) for potential carcinogenic metals was calculated by multiplying the doses by the corresponding slope factor (SF), as given in Eq. (9)35. The carcinogenic oral, inhalation, and dermal SF, as well as dermal absorption toxicity values, were provided from the Integrated Risk Information System30. The reference doses for Pb were taken from the Guidelines for Drinking Water Quality published by the World Health Organization36.$${text{Risk}} = { }mathop sum limits_{{{text{k}} = 1}}^{{text{n}}} {text{D}}_{{text{k}}} cdot {text{ SF}}_{{text{k}}} ,$$
(9)
where SF is the cancer slope factor for individually metal and k route of exposure (ingestion, inhalation, or dermal path). The total cancer risk (Risktotal) of potential carcinogens was calculated as the sum of the individual risk values using the following Eq. (10).$${text{Risk}}_{{{text{total}}}} = {text{Risk}}_{{{text{ing}}}} + {text{Risk}}_{{{text{inh}}}} + {text{Risk}}_{{{text{dermal}}}} .$$
(10) More

Regular structures and isothermal theoremFor networks where each node represents a single individual, the isothermal theorem of evolutionary graph theory shows that the fixation probability is the same as the fixation probability of a well-mixed population if the temperature distribution is homogeneous across the whole population1. The temperature of a node defined as the sum over all the weights leads to that node. This theorem extends to structured meta-populations for any migration probability (lambda ): If the underlying structure of the meta-population that connects the patches is a regular network and the local population size is identical in each patch, the temperature of all individuals is identical, regardless of the value of the migration probability. Therefore, the fixation probability in a population with such a structure is the same as the fixation probability in a well-mixed population of the same total population size (N=sum _{j=1}^M N_j), given by ( phi _{mathrm{wm}}^N(r)).Small migration regimeIf the migration probability is small enough such that the time between two subsequent migration events (( sim frac{1}{lambda } )) is much longer than the absorption time within any patch, then at the time of each migration event we may suppose that the meta-population is in a homogeneous configuration22,28. In other words, the low migration regime is an approximation in which we neglect the probability that the meta-population is not in a homogeneous configuration at the time of migration events. We define a homogeneous configuration of the meta-population as a configuration in which in all patches either all individuals are mutants, or all are wild-types.Therefore, instead of having (2^N) states, where N is the population size, the system has only (2^M) states, where M is the number of patches. Thus, we can calculate the fixation probability exactly as in the case of a standard evolutionary graph model where each node represents a single individual but with a modified transition probabilities.In a network with homogeneous patches, in order to increase the number of homogeneous mutant-patches one individual mutant needs to migrate to one of its neighbouring homogeneous wild-type-patches and reaches fixation there. For example if node j is occupied by mutants and one of its neighbouring patches, node k, is occupied by wild-types, the probability that one mutant individual from patch j migrates to patch k and reaches fixation there is (frac{lambda }{mathrm{deg} (j)}phi _{mathrm{wm}}^{N_{k}}(r) ), where (mathrm{deg} (j) ) is the degree of node j to take into account that the mutant can move to different patches. This is analogous to the probability that one mutant in node j replaces one wild-type in node k ,(T^{jrightarrow k}), in the network of individuals.Similarly, if node j is occupied by wild-types and one of its neighbouring patches, node j, is occupied by mutants the probability that one wild-type individual from patch j migrates to patch k and reaches fixation there equals to (frac{lambda }{mathrm{deg} (j)}phi _{mathrm{wm}}^{N_{k}}(1/r) ) where (mathrm{deg} (j) ). Overall, we can move from network of individuals to the network of homogeneous patches by replacing the transition probabilities with the product of migration and fixation probabilities.Two-patch meta-populationThe simplest non-trivial case is the fixation probability in a two-patch meta-population with different local size for small migration probability (lambda ). If the migration probability (lambda ) is very small, a new mutant first needs to take over its own patch and only then the first migrant arrives in the second patch. To be more precise, the time between two migration events has to be much higher than the typical time that it takes for the migrant to take over the patch or go extinct again38. In this case, we can divide the dynamics into two phases: A first phase in which a mutant invades one patch and a second phase in which a homogeneous patch of mutants invades the whole meta-population. Assume a new mutation arises in patch 1. Only if this mutant reaches fixation in patch 1, it also has a chance to reach fixation in patch 2. When patch 1 consists of only mutants and patch 2 consists of only wild-types, there are two possibilities for the ultimate fate of the mutant:

(i)

Eventually, the offspring of one mutant selected from patch 1 for reproduction will migrate to patch 2 and reach fixation there. The wild-type goes extinct. This happens with probability ( frac{N_1 r}{N_1 r+N_2} phi _{mathrm{wm}}^{N_2}(r)).

(ii)

Eventually, the offspring of one wild-type selected from patch 2 for reproduction will migrate to patch 1 and the mutant goes extinct. This occurs with probability ( frac{N_2}{N_1r+N_2} phi _{mathrm{wm}}^{N_1}(tfrac{1}{r})).

Therefore, the probability that a single mutant arising in patch 1 reaches fixation in the entire population is $$begin{aligned} phi _{mathrm{wm}}^{N_1}(r) frac{frac{N_1 r}{N_1 r+N_2} phi _{mathrm{wm}}^{N_2}(r)}{frac{N_1 r}{N_1 r+N_2} phi _{mathrm{wm}}^{N_2}(r)+frac{N_2}{N_1r+N_2} phi _{mathrm{wm}}^{N_1}left( tfrac{1}{r}right) }=phi _{mathrm{wm}}^{N_1}(r) phi _{mathrm{wm}}^{N_2}(r) frac{1 }{ phi _{mathrm{wm}}^{N_2}(r) +frac{N_2}{N_1} frac{1}{r}phi _{mathrm{wm}}^{N_1} left( tfrac{1}{r}right) }. end{aligned}$$
(3a)
Similarly the probability that a mutant arising in patch 2 takes over the whole population equals$$begin{aligned} phi _{mathrm{wm}}^{N_2}(r) phi _{mathrm{wm}}^{N_1}(r) frac{1 }{phi _{mathrm{wm}}^{N_1}(r)+frac{N_1}{N_2} frac{1}{r} phi _{mathrm{wm}}^{N_2}left( tfrac{1}{r}right) }. end{aligned}$$
(3b)
If we assume that the mutant arises in a patch with a probability proportional to the patch size, the average fixation probability (phi _{bullet !!-!!bullet }) in a two patch population for small migration probability is the weighted sum of Eqs. (3a) and (3b),$$begin{aligned} phi _{bullet !!-!!bullet }&= phi _{mathrm{wm}}^{N_1}(r) phi _{mathrm{wm}}^{N_2}(r) nonumber \&quad times left( frac{frac{N_1}{N_1+N_2} }{ phi _{mathrm{wm}}^{N_2}(r) +frac{N_2}{N_1} frac{1}{r}phi _{mathrm{wm}}^{N_1}left( tfrac{1}{r}right) } +frac{frac{N_2}{N_1+N_2} }{ phi _{mathrm{wm}}^{N_1}(r) +frac{N_1}{N_2} frac{1}{r} phi _{mathrm{wm}}^{N_2}left( tfrac{1}{r}right) }right) . end{aligned}$$
(4)
In the case of neutrality, (r=1), we recover (phi _{bullet !!-!!bullet } = frac{1}{N_1+N_2})—the fixation probability in a population of the total size of the two patches. For identical patch sizes, ( N_1=N_2 ), Eq. (4) simplifies to$$begin{aligned} phi _{bullet !!-!!bullet } = left( phi _{mathrm{wm}}^{N_1}(r)right) ^2 frac{1}{phi _{mathrm{wm}}^{N_1}(r)+frac{1}{r} phi _{mathrm{wm}}^{N_1}left( tfrac{1}{r}right) } = phi _{mathrm{wm}}^{2 N_1}(r), end{aligned}$$
(5)
where the simplification to the fixation probability within a single population of size (2N_1) reflects the validity of the isothermal theorem.For (N_1 ne N_2), we approximate Eq. (4) for weak and strong selection. Let us first consider highly advantageous mutants, (r gg 1). In this case, we have (phi _{mathrm{wm}}^{N_1}(r) gg phi _{mathrm{wm}}^{N_1}(tfrac{1}{r})) and thus we can neglect the possibility that a wild-type takes over a mutant patch if patch sizes are sufficiently large. The probability (phi _{bullet !!-!!bullet } ) then becomes a weighted average reflecting patch sizes. For identical patch size (N_1=N_2 = N/2), it reduces to (phi _{bullet !!-!!bullet } approx phi _{mathrm{wm}}^{N_1}(r)=phi _{mathrm{wm}}^{N/2}(r)). In other words, taking over the first patch is sufficient to make fixation in the entire population certain. For patches of very different size, (N_1 gg N_2), we have (N approx N_1) and find (phi _{bullet !!-!! bullet } approx phi _{mathrm{wm}}^{N}(r), ) which implies that fixation is driven by the fixation process in the larger patch, regardless of where the mutant arises. Note that there is a difference between the case of identical patch size and very different patch size . The case of highly disadvantageous mutants, (r ll 1), can be handled in a very similar way.Next, we consider weak selection, (r approx 1). We can approximate the fixation probability as (phi _{mathrm{wm}}^{N}(r^{pm 1}) approx frac{1}{N} pm frac{N-1}{2N} (r-1)). With this, we find$$begin{aligned} phi _{bullet !!-!!bullet } approx frac{1}{N_1+N_2} +frac{1}{2} left( 1 – frac{1}{N_1+N_2} -frac{(N_1-N_2)^2}{(N_1^2+N_2^2)^2} N_1 N_2right) (r-1). end{aligned}$$
(6)
For identical patch size (N_1=N_2 = N/2), this reduces to$$begin{aligned} phi _{bullet !!-!!bullet } approx tfrac{1}{N} +tfrac{N-1}{2N} (r-1), end{aligned}$$
(7)
which is the known result for a single population of size (N=N_1+N_2). When patches have very different size, (N_1 gg N_2) such that (N approx N_1), we recover the same result. Thus, the difference between the fixation probability of a two-patch meta-population with identical patch size and the fixation probability of a two-patch meta-population with very different patch size that we found for highly advantageous mutants is no longer observed for weak selection.When migration probabilities become larger, our approximation is no longer valid and we need to rely on numerical approaches. Figure 2 illustrates the difference between the fixation probability of a two-patch structure meta-population and the equivalent well-mixed population of size (N_1+N_2 ) when migration is low using Eq. (4) and comparing with the numerical approach in Ref.39.While the fixation probability of the two-patch meta-population is very close to the fixation probability of the well-mixed population40, a close inspection reveals an interesting property: For low migration probabilities and (N_1 ne N_2), the two patch structure is a suppressor of selection in the original sense of Lieberman et al.1: For advantageous mutations, (r >1), it decreases the fixation probability, whereas for disadvantageous mutations, (r1) and negative for (r1 ) the minimum fixation probability occurs when the two patch sizes are identical, ( N_1=N_2=N/2 ). Similarly, for fitness values ( r More