Ecology

1.Stal, L. J. et al. BASIC: Baltic Sea cyanobacteria. An investigation of the structure and dynamics of water blooms of cyanobacteria in the Baltic Sea—Responses to a changing environment. Cont. Shelf Res. 23, 1695–1714 (2003).Article 
ADS 

Google Scholar 
2.McGregor, G. B. et al. First report of a toxic Nodularia spumigena (nostocales/cyanobacteria) bloom in sub-tropical Australia. I. Phycological and public health investigations. Int. J. Env. Res. Public Health 9, 2396–2411 (2012).Article 

Google Scholar 
3.Popin, R. V. et al. Genomic and metabolomic analyses of natural products in Nodularia spumigena isolated from a shrimp culture pond. Toxins 12, 141 (2020).CAS 
PubMed Central 
Article 
PubMed 

Google Scholar 
4.Seaman, M., Ashton, P. & Williams, W. Inland salt waters of southern Africa. Hydrobiologia 210, 75–91 (1991).CAS 
Article 

Google Scholar 
5.Beutel, M. W., Horne, A. J., Roth, J. C. & Barratt, N. J. Saline Lakes 91–105 (Springer, 2001).Book 

Google Scholar 
6.Paerl, H. W. & Paul, V. J. Climate change: Links to global expansion of harmful cyanobacteria. Water Res. 46, 1349–1363 (2012).CAS 
PubMed 
Article 

Google Scholar 
7.Karjalainen, M. et al. Ecosystem consequences of cyanobacteria in the northern Baltic Sea. AMBIO J. Human Environ. 36, 195–202 (2007).CAS 
Article 

Google Scholar 
8.Sotton, B., Domaizon, I., Anneville, O., Cattanéo, F. & Guillard, J. Nodularin and cylindrospermopsin: A review of their effects on fish. Rev. Fish Biol. Fish. 25, 1–19 (2015).Article 

Google Scholar 
9.Mazur-Marzec, H., Bertos-Fortis, M., Toruńska-Sitarz, A., Fidor, A. & Legrand, C. Chemical and genetic diversity of Nodularia spumigena from the Baltic Sea. Mar. Drugs 14, 209. https://doi.org/10.3390/md14110209 (2016).CAS 
Article 
PubMed Central 
PubMed 

Google Scholar 
10.Voss, B. et al. Insights into the physiology and ecology of the brackish-water-adapted Cyanobacterium Nodularia spumigena CCY9414 based on a genome-transcriptome analysis. PLoS ONE 8, e60224–e60224. https://doi.org/10.1371/journal.pone.0060224 (2013).CAS 
Article 
PubMed 
PubMed Central 
ADS 

Google Scholar 
11.Le Manach, S. et al. Global metabolomic characterizations of Microcystis spp. highlights clonal diversity in natural bloom-forming populations and expands metabolite structural diversity. Front. Microbiol. https://doi.org/10.3389/fmicb.2019.00791 (2019).Article 
PubMed 
PubMed Central 

Google Scholar 
12.Welker, M. & von Döhren, H. Cyanobacterial peptides—Nature’s own combinatorial biosynthesis. FEMS Microbiol. Rev. 30, 530–563 (2006).CAS 
PubMed 
Article 

Google Scholar 
13.Kehr, J. C., Picchi, D. G. & Dittmann, E. Natural product biosyntheses in cyanobacteria: A treasure trove of unique enzymes. Beilstein J. Org. Chem. 7, 1622–1635. https://doi.org/10.3762/bjoc.7.191 (2011).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
14.Christiansen, G., Philmus, B., Hemscheidt, T. & Kurmayer, R. Genetic variation of adenylation domains of the anabaenopeptin synthesis operon and evolution of substrate promiscuity. J. Bacteriol. 193, 3822–3831 (2011).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
15.Ishida, K. et al. Biosynthesis and structure of aeruginoside 126A and 126B, cyanobacterial peptide glycosides bearing a 2-carboxy-6-hydroxyoctahydroindole moiety. Chem. Biol. 14, 565–576 (2007).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
16.Fewer, D.P. et al. The non-ribosomal assembly and frequent occurrence of the protease inhibitors spumigins in the bloom-forming cyanobacterium Nodularia spumigena. Mol. Microbiol. 73, 924–937. https://doi.org/10.1111/j.1365-2958.2009.06816.x (2009).CAS 
Article 
PubMed 

Google Scholar 
17.Portmann, C. et al. Isolation of aerucyclamides C and D and structure revision of microcyclamide 7806A: Heterocyclic ribosomal peptides from Microcystis aeruginosa PCC 7806 and their antiparasite evaluation. J. Nat. Prod. 71, 1891–1896 (2008).CAS 
PubMed 
Article 

Google Scholar 
18.Ersmark, K., Del Valle, J. R. & Hanessian, S. Chemistry and biology of the aeruginosin family of serine protease inhibitors. Angew. Chem. Int. Ed. 47, 1202–1223 (2008).CAS 
Article 

Google Scholar 
19.Liu, L. et al. Pseudoaeruginosins, nonribosomal peptides in Nodularia spumigena. ACS Chem. Biol. 10, 725–733 (2015).CAS 
PubMed 
Article 

Google Scholar 
20.Itou, Y., Suzuki, S., Ishida, K. & Murakami, M. Anabaenopeptins G and H, potent carboxypeptidase A inhibitors from the cyanobacterium Oscillatoria agardhii (NIES-595). Bioorg. Med. Chem. Lett. 9, 1243–1246 (1999).CAS 
PubMed 
Article 

Google Scholar 
21.Bister, B. et al. Cyanopeptolin 963A, a chymotrypsin inhibitor of Microcystis PCC 7806. J. Nat. Prod. 67, 1755–1757 (2004).CAS 
PubMed 
Article 

Google Scholar 
22.Neilan, B. A., Pearson, L. A., Muenchhoff, J., Moffitt, M. C. & Dittmann, E. Environmental conditions that influence toxin biosynthesis in cyanobacteria. Environ. Microbiol. 15, 1239–1253. https://doi.org/10.1111/j.1462-2920.2012.02729.x (2013).CAS 
Article 
PubMed 

Google Scholar 
23.Halstvedt, C. B., Rohrlack, T., Ptacnik, R. & Edvardsen, B. On the effect of abiotic environmental factors on production of bioactive oligopeptides in field populations of Planktothrix spp. (Cyanobacteria). J. Plankton Res. 30, 607–617 (2008).CAS 
Article 

Google Scholar 
24.Mazur-Marzec, H. et al. Diversity of peptides produced by Nodularia spumigena from various geographical regions. Mar. Drugs 11, 1–19. https://doi.org/10.3390/md11010001 (2012).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
25.Repka, S., Koivula, M., Harjunpa, V., Rouhiainen, L. & Sivonen, K. Effects of phosphate and light on growth of and bioactive peptide production by the Cyanobacterium anabaena strain 90 and its anabaenopeptilide mutant. Appl. Environ. Microbiol. 70, 4551–4560. https://doi.org/10.1128/aem.70.8.4551-4560.2004 (2004).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
26.Lehtimäki, J., Moisander, P., Sivonen, K. & Kononen, K. Growth, nitrogen fixation, and nodularin production by two Baltic Sea cyanobacteria. Appl. Environ. Microbiol. 63, 1647–1656 (1997).PubMed 
PubMed Central 
Article 

Google Scholar 
27.OECD. Test No. 201: Freshwater alga and cyanobacteria, growth inhibition test. OECD Guidelines for the Testing of Chemicals, Section 2. https://doi.org/10.1787/9789264069923-en (OECD
Publishing, Paris, 2011).28.Vaas, L. A. I., Sikorski, J., Michael, V., Göker, M. & Klenk, H.-P. Visualization and curve-parameter estimation strategies for efficient exploration of phenotype microarray kinetics. PLoS ONE 7, e34846. https://doi.org/10.1371/journal.pone.0034846 (2012).CAS 
Article 
PubMed 
PubMed Central 
ADS 

Google Scholar 
29.Higo, S., Yamatogi, T., Ishida, N., Hirae, S. & Koike, K. Application of a pulse-amplitude-modulation (PAM) fluorometer reveals its usefulness and robustness in the prediction of Karenia mikimotoi blooms: A case study in Sasebo Bay, Nagasaki, Japan. Harmful Algae 61, 63–70 (2017).Article 

Google Scholar 
30.Qi, H., Wang, J. & Wang, Z. A comparative study of maximal quantum yield of photosystem II to determine nitrogen and phosphorus limitation on two marine algae. J. Sea Res. 80, 1–11 (2013).Article 
ADS 

Google Scholar 
31.Briand, E., Bormans, M., Gugger, M., Dorrestein, P. C. & Gerwick, W. H. Changes in secondary metabolic profiles of Microcystis aeruginosa strains in response to intraspecific interactions. Environ. Microbiol. 18, 384–400. https://doi.org/10.1111/1462-2920.12904 (2016).CAS 
Article 
PubMed 

Google Scholar 
32.Koek, M. M., Muilwijk, B., van der Werf, M. J. & Hankemeier, T. Microbial metabolomics with gas chromatography/mass spectrometry. Anal. Chem. 78, 1272–1281 (2006).CAS 
PubMed 
Article 

Google Scholar 
33.R Development Core Team. R: A language and environment for statistical computing. (R Foundation for Statistical Computing, Vienna, Austria, 2020).34.Medlock, G. L. et al. Inferring metabolic mechanisms of interaction within a defined gut microbiota. Cell Syst. 7, 245-257.e247. https://doi.org/10.1016/j.cels.2018.08.003 (2018).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
35.Paul, C., Mausz, M. A. & Pohnert, G. A co-culturing/metabolomics approach to investigate chemically mediated interactions of planktonic organisms reveals influence of bacteria on diatom metabolism. Metabolomics 9, 349–359. https://doi.org/10.1007/s11306-012-0453-1 (2013).CAS 
Article 

Google Scholar 
36.Schatz, D. et al. Ecological implications of the emergence of non-toxic subcultures from toxic Microcystis strains. Environ. Microbiol. 7, 798–805 (2005).CAS 
PubMed 
Article 

Google Scholar 
37.Jensen, A., Rystad, B. & Skoglund, L. The use of dialysis culture in phytoplankton studies. J. Exp. Mar. Biol. Ecol. 8, 241–248 (1972).Article 

Google Scholar 
38.Kobayashi, K., Takata, Y. & Kodama, M. Direct contact between Pseudo-nitzschiaámultiseries and bacteria is necessary for the diatom to produce a high level of domoic acid. Fish. Sci. 75, 771–776 (2009).CAS 
Article 

Google Scholar 
39.McVeigh, I., & Brown, W. In vitro growth of chlamydomonas chlamydogama bold and haematococcus pluvialis flotow em. Wille in mixed cultures.
Bulletin of the Torrey Botanical Club, 81(3), 218–233. https://doi.org/10.2307/2481813 (1954).CAS 
Article 

Google Scholar 
40.Sieg, R. D., Poulson-Ellestad, K. L. & Kubanek, J. Chemical ecology of the marine plankton. Nat. Prod. Rep. 28, 388–399 (2011).CAS 
PubMed 
Article 

Google Scholar 
41.Yamasaki, A. An overview of CO2 mitigation options for global warming—Emphasizing CO2 sequestration options. J. Chem. Eng. Japan 36, 361–375 (2003).CAS 
Article 

Google Scholar 
42.Hajdu, S., Hoglander, H. & Larsson, U. Phytoplankton vertical distributions and composition in Baltic Sea cyanobacterial blooms. Harmful Algae 6, 189–205 (2007).Article 

Google Scholar 
43.Berman-Frank, I. & Dubinsky, Z. Balanced growth in aquatic plants: Myth or reality? Phytoplankton use the imbalance between carbon assimilation and biomass production to their strategic advantage. Bioscience 49, 29–37 (1999).Article 

Google Scholar 
44.Kruskopf, M. & Flynn, K. J. Chlorophyll content and fluorescence responses cannot be used to gauge reliably phytoplankton biomass, nutrient status or growth rate. New Phytol. 169, 525–536. https://doi.org/10.1111/j.1469-8137.2005.01601.x (2006).CAS 
Article 
PubMed 

Google Scholar 
45.Raven, J. A. & Beardall, J. Chlorophyll fluorescence and ecophysiology: Seeing red?. New Phytol. 169, 449–451. https://doi.org/10.1111/j.1469-8137.2006.01637.x (2006).CAS 
Article 
PubMed 

Google Scholar 
46.Li, Q. et al. A large-scale comparative metagenomic study reveals the functional interactions in six bloom-forming microcystis-epibiont communities. Front. Microbiol. https://doi.org/10.3389/fmicb.2018.00746 (2018).Article 
PubMed 
PubMed Central 

Google Scholar 
47.Harke, M. J. et al. A review of the global ecology, genomics, and biogeography of the toxic cyanobacterium Microcystis spp. Harmful Algae 54, 4–20 (2016).PubMed 
Article 

Google Scholar 
48.Caldwell, D. Associations between photosynthetic and heterotrophic prokaryotes in plankton. in Abstracts of the third International Symposium on Photosynthetic Prokaryotes (ed Nichols, J. M) (University of Liverpool, UK, 1979).49.Park, H. D. et al. Degradation of the cyanobacterial hepatotoxin microcystin by a new bacterium isolated from a hypertrophic lake. Environ. Toxicol. Int. J. 16, 337–343 (2001).CAS 
Article 
ADS 

Google Scholar 
50.Berg, C. et al. Dissection of microbial community functions during a cyanobacterial bloom in the Baltic Sea via metatranscriptomics. Front. Mar. Sci. 5, 55 (2018).Article 

Google Scholar 
51.Humbert, J.-F. et al. A tribute to disorder in the genome of the bloom-forming freshwater cyanobacterium Microcystis aeruginosa. PLoS ONE 8, e70747. https://doi.org/10.1371/journal.pone.0070747 (2013).CAS 
Article 
PubMed 
PubMed Central 
ADS 

Google Scholar 
52.Toporowska, M., Mazur-Marzec, H. & Pawlik-Skowrońska, B. The effects of cyanobacterial bloom extracts on the biomass, Chl-a, MC and other oligopeptides contents in a natural Planktothrix agardhii population. Int. J. Env. Res. Public Health 17, 2881 (2020).CAS 
Article 

Google Scholar 
53.Grabowska, M., Kobos, J., Toruńska-Sitarz, A. & Mazur-Marzec, H. Non-ribosomal peptides produced by Planktothrix agardhii from Siemianówka Dam Reservoir SDR (northeast Poland). Arch. Microbiol. 196, 697–707 (2014).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
54.Penn, K., Wang, J., Fernando, S. C. & Thompson, J. R. Secondary metabolite gene expression and interplay of bacterial functions in a tropical freshwater cyanobacterial bloom. ISME J. 8, 1866–1878 (2014).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
55.Neilan, B. A. et al. Nonribosomal peptide synthesis and toxigenicity of cyanobacteria. J. Bacteriol. 181, 4089–4097 (1999).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
56.Long, B. M., Jones, G. J. & Orr, P. T. Cellular microcystin content in N-limited Microcystis aeruginosa can be predicted from growth rate. Appl. Environ. Microbiol. 67, 278–283 (2001).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
57.Qu, J. et al. Determination of the role of microcystis aeruginosa in toxin generation based on phosphoproteomic profiles. Toxins 10, 304 (2018).PubMed Central 
Article 
CAS 
PubMed 

Google Scholar 
58.Raven, J. A. Cyanotoxins: A poison that frees phosphate. Curr. Biol. 20, R850–R852 (2010).CAS 
PubMed 
Article 

Google Scholar 
59.Utkilen, H. & Gjølme, N. Iron-stimulated toxin production in Microcystis aeruginosa. Appl. Environ. Microbiol. 61, 797–800 (1995).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
60.Gan, N. et al. The role of microcystins in maintaining colonies of bloom-forming Microcystis spp. Environ. Microbiol. 14, 730–742 (2012).CAS 
PubMed 
Article 

Google Scholar 
61.Pomati, F., Rossetti, C., Manarolla, G., Burns, B. P. & Neilan, B. A. Interactions between intracellular Na+ levels and saxitoxin production in Cylindrospermopsis raciborskii T3. Microbiology 150, 455–461 (2004).CAS 
PubMed 
Article 

Google Scholar 
62.Seigler, D. & Price, P. W. Secondary compounds in plants: Primary functions. Am. Nat. 110, 101–105 (1976).CAS 
Article 

Google Scholar 
63.Zilliges, Y. et al. The cyanobacterial hepatotoxin microcystin binds to proteins and increases the fitness of Microcystis under oxidative stress conditions. PLoS ONE 6, e17615 (2011).CAS 
PubMed 
PubMed Central 
Article 
ADS 

Google Scholar 
64.Meissner, S., Fastner, J. & Dittmann, E. Microcystin production revisited: Conjugate formation makes a major contribution. Environ. Microbiol. 15, 1810–1820. https://doi.org/10.1111/1462-2920.12072 (2013).CAS 
Article 
PubMed 

Google Scholar 
65.Orr, P. T., Willis, A. & Burford, M. A. Application of first order rate kinetics to explain changes in bloom toxicity—The importance of understanding cell toxin quotas. J. Oceanol. Limnol. 36, 1063–1074. https://doi.org/10.1007/s00343-019-7188-z (2018).CAS 
Article 
ADS 

Google Scholar 
66.Rantala, A. et al. Phylogenetic evidence for the early evolution of microcystin synthesis. Proc. Natl. Acad. Sci. USA 101, 568–573 (2004).CAS 
PubMed 
Article 
ADS 

Google Scholar 
67.Orr, P. T. & Jones, G. J. Relationship between microcystin production and cell division rates in nitrogen-limited Microcystis aeruginosa cultures. Limnol. Oceanogr. 43, 1604–1614 (1998).CAS 
Article 
ADS 

Google Scholar 
68.Burford, M. A. et al. Understanding the winning strategies used by the bloom-forming cyanobacterium Cylindrospermopsis raciborskii. Harmful Algae 54, 44–53 (2016).PubMed 
Article 

Google Scholar 
69.Pierangelini, M. et al. Constitutive cylindrospermopsin pool size in Cylindrospermopsis raciborskii under different light and CO2 partial pressure conditions. Appl. Environ. Microbiol. 81, 3069–3076. https://doi.org/10.1128/aem.03556-14 (2015).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
70.Falkowski, P. G., Sukenik, A. & Herzig, R. Nitrogen limitation in Isochrysis galbana (Haptophyceae). II. Relative abundance of chloroplast proteins. J. Phycol. 25, 471–478 (1989).CAS 
Article 

Google Scholar 
71.Turpin, D. H. Effects of inorganic N availability on algal photosynthesis and carbon metabolism. J. Phycol. 27, 14–20 (1991).CAS 
Article 

Google Scholar 
72.Moffitt, M. C. & Neilan, B. A. Characterization of the nodularin synthetase gene cluster and proposed theory of the evolution of cyanobacterial hepatotoxins. Appl. Environ. Microbiol. 70, 6353–6362 (2004).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
73.Fewer, D. P. et al. New structural variants of aeruginosin produced by the toxic bloom forming cyanobacterium Nodularia spumigena. PLoS ONE 8, e73618 (2013).CAS 
PubMed 
PubMed Central 
Article 
ADS 

Google Scholar 
74.Fujii, K. et al. Comparative study of toxic and non-toxic cyanobacterial products: Novel peptides from toxic Nodularia spumigena AV1. Tetrahedron Lett. 38, 5525–5528 (1997).CAS 
Article 

Google Scholar 
75.Ishida, K. et al. Plasticity and evolution of aeruginosin biosynthesis in cyanobacteria. Appl. Environ. Microbiol. 75, 2017–2026 (2009).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
76.Suikkanen, S., Fistarol, G. O. & Granéli, E. Allelopathic effects of the Baltic cyanobacteria Nodularia spumdigena, Aphanizomenon flosaquae and Anabaena lemmermannii on algal monocultures. J. Exp. Mar. Biol. Ecol. 308, 85–101 (2004).Article 

Google Scholar 
77.Suikkanen, S., Engström-Öst, J., Jokela, J., Sivonen, K. & Viitasalo, M. Allelopathy of Baltic Sea cyanobacteria: No evidence for the role of nodularin. J. Plankton Res. 28, 543–550. https://doi.org/10.1093/plankt/fbi139 (2006).CAS 
Article 

Google Scholar 
78.Żak, A. & Kosakowska, A. The influence of extracellular compounds produced by selected Baltic cyanobacteria, diatoms and dinoflagellates on growth of green algae Chlorella vulgaris. Estuar. Coast. Shelf Sci. 167, 113–118 (2015).Article 
ADS 

Google Scholar 
79.Śliwińska-Wilczewska, S., Felpeto, A. B., Możdżeń, K., Vasconcelos, V. & Latała, A. Physiological effects on coexisting microalgae of the allelochemicals produced by the bloom-forming cyanobacteria Synechococcus sp. and Nodularia spumigena. Toxins 11, 712 (2019).PubMed Central 
Article 
CAS 
PubMed 

Google Scholar 
80.Gross, E. M. Allelopathy of aquatic autotrophs. Crit. Rev. Plant Sci. 22, 313–339 (2003).Article 

Google Scholar 
81.Legrand, C., Rengefors, K., Fistarol, G. O. & Graneli, E. Allelopathy in phytoplankton-biochemical, ecological and evolutionary aspects. Phycologia 42, 406–419 (2003).Article 

Google Scholar 
82.Leao, P. N., Vasconcelos, M. T. & Vasconcelos, V. M. Allelopathy in freshwater cyanobacteria. Crit. Rev. Microbiol. 35, 271–282. https://doi.org/10.3109/10408410902823705 (2009).CAS 
Article 
PubMed 

Google Scholar 
83.MacKintosh, C., Beattie, K. A., Klumpp, S., Cohen, P. & Codd, G. A. Cyanobacterial microcystin-LR is a potent and specific inhibitor of protein phosphatases 1 and 2A from both mammals and higher plants. FEBS Lett. 264, 187–192 (1990).CAS 
PubMed 
Article 

Google Scholar 
84.Pflugmacher, S. Possible allelopathic effects of cyanotoxins, with reference to microcystin-LR, in aquatic ecosystems. Environ. Toxicol. Int. J. 17, 407–413 (2002).CAS 
Article 
ADS 

Google Scholar 
85.Tilahun S. Exclusive partitioning of intra- and extra-cellular cyanotoxins: limitation of the conventional procedure. Environ. Sci. Pollut. Res. Int. 27(14), 17427–17428. https://doi.org/10.1007/s11356-020-08256-8 (2020).Article 
PubMed 

Google Scholar 
86.Park, H. D. et al. Temporal variabilities of the concentrations of intra-and extracellular microcystin and toxic Microcystis species in a hypertrophic lake, Lake Suwa, Japan (1991–1994). Environ. Toxicol. Water Qual. Int. J. 13, 61–72 (1998).CAS 
Article 
ADS 

Google Scholar 
87.Tsuji, K. et al. Stability of microcystins from cyanobacteria: Effect of light on decomposition and isomerization. Environ. Sci. Technol. 28, 173–177 (1994).CAS 
PubMed 
Article 
ADS 
PubMed Central 

Google Scholar 
88.Schatz, D. et al. Towards clarification of the biological role of microcystins, a family of cyanobacterial toxins. Environ. Microbiol. 9, 965–970 (2007).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
89.Makower, A. K. et al. Transcriptomics-aided dissection of the intracellular and extracellular roles of microcystin in Microcystis aeruginosa PCC 7806. Appl. Environ. Microbiol. 81, 544–554 (2015).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
90.Kaplan, A. et al. The languages spoken in the water body (or the biological role of cyanobacterial toxins). Front. Microbiol. 3, 138 (2012).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
91.Svercel, M. Negative allelopathy among cyanobacteria. in Cyanobacteria: Ecology, Toxicology and Management. (ed Ferrao-Filho, A. S.) 27–46 (Nova Science Publishers, New York, NY, USA, 2013).
Google Scholar 
92.Wiegand, C. & Pflugmacher, S. Ecotoxicological effects of selected cyanobacterial secondary metabolites a short review. Toxicol. Appl. Pharmacol. 203, 201–218. https://doi.org/10.1016/j.taap.2004.11.002 (2005).CAS 
Article 
PubMed 

Google Scholar 
93.Agrawal, M. & Agrawal, M. K. Cyanobacteria–herbivore interaction in freshwater ecosystem. J. Microbiol. Biotechnol. Res. 1, 52–66 (2011).
Google Scholar 
94.Sadler, T. & von Elert, E. Dietary exposure of Daphnia to microcystins: No in vivo relevance of biotransformation. Aquat. Toxicol. 150, 73–82. https://doi.org/10.1016/j.aquatox.2014.02.017 (2014).CAS 
Article 
PubMed 

Google Scholar 
95.Rohrlack, T., Christiansen, G. & Kurmayer, R. Putative antiparasite defensive system involving ribosomal and nonribosomal oligopeptides in cyanobacteria of the genus Planktothrix. Appl. Environ. Microbiol. 79, 2642–2647 (2013).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
96.Sivonen, K. et al. Occurrence of the hepatotoxic cyanobacterium Nodularia spumigena in the Baltic Sea and structure of the toxin. Appl. Environ. Microbiol. 55, 1990–1995 (1989).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
97.Burbage, C. D. & Binder, B. J. Relationship between cell cycle and light-limited growth rate in oceanic Prochlorococcus (MIT9312) and Synechococcus (WH8103) (Cyanobacteria). J. Phycol. 43, 266–274. https://doi.org/10.1111/j.1529-8817.2007.00315.x (2007).Article 

Google Scholar 
98.Lei, L., Dai, J., Lin, Q., Peng, L. Competitive dominance of Microcystis aeruginosa against Raphidiopsis raciborskii is strain-and temperature dependent. Knowl. Manag. Aquat. Ecosyst. 421, 36. https://doi.org/10.1051/kmae/2020023 (2020).Article 

Google Scholar  More

1.Schulthess, C. P. & Huang, C. P. Adsorption of heavy metals by silicon and aluminum oxide surfaces on clay minerals. Soil Sci. Soc. Am. J. 54, 679–688 (1990).ADS 
Article 

Google Scholar 
2.City of Joplin Health Department. Report to Jasper County EPA Superfund citizen’s task force. City of Joplin Health Department, Joplin, MO (1995).3.Beyer, W. N., Pattee, O. H., Sileo, L., Hoffman, D. J. & Mulhern, B. M. Metal contamination in wildlife living near two zinc smelters. Environ. Pollut. Ser. A38, 63–86 (1985).Article 

Google Scholar 
4.Khan, D. H. & Frankland, B. Effects of cadmium and lead on radish plants with particular reference to movement of metals through soil profile and plant. Plant Soil 70, 335–345 (1983).CAS 
Article 

Google Scholar 
5.Ruby, M. V., Davis, A., Kempton, J. H., Drexler, J. & Bergstrom, P. D. Lead bioavailability: Dissolution kinetics under simulated gastric conditions. Environ. Sci. Technol. 26, 1242–1248 (1992).ADS 
CAS 
Article 

Google Scholar 
6.Ruby, M. V., Davis, A. & Nicholson, A. In-situ formation of lead phosphates in soils as a method to immobilize lead. Environ. Sci. Tchnol. 28, 646–654 (1994).ADS 
CAS 
Article 

Google Scholar 
7.Nriagu, J. O. Lead orthophosphates-IV: Formation and stability in the environment. Geochim. Cosmochim. Acta 38, 887–898 (1974).ADS 
CAS 
Article 

Google Scholar 
8.Ma, Q. Y., Logan, T. J. & Traina, S. J. Lead immobilization from aqueous solutions and contaminated soils using phosphate rocks. Environ. Sci. Technol. 27, 1118–1126 (1995).ADS 
Article 

Google Scholar 
9.Xu, Y. & Schwartz, F. W. Lead immobilization by hydroxyapatite in aqueous solutions. J. Contam. Hydrol. 15, 187–206 (1994).ADS 
CAS 
Article 

Google Scholar 
10.Zhang, P. C., Ryan, J. A. & Yang, J. In vitro soil Pb solubility in the presence of hydroxyapatite. Environ. Sci. Technol. 32, 2763–2768 (1998).ADS 
CAS 
Article 

Google Scholar 
11.Osborne, L. R., Baker, L. L. & Strawn, D. G. Lead immobilization and phosphorus availability in phosphate-amended, mine-contaminated soils. J. Environ. Qual. 44(1), 183–190 (2015).Article 

Google Scholar 
12.Yang, J. J., Mosby, D. E., Casteel, S. W. & Blanchar, R. W. Lead immobilization using phosphoric acid in smelter-contaminated urban soil. Environ. Sci. Technol. 35, 3553–3559 (2001).ADS 
CAS 
Article 

Google Scholar 
13.Tang, X., Yang, J., Goyne, K. W. & Deng, B. Long-term risk reduction of lead contaminated urban soil by phosphate treatment. Env. Eng. Sci. 26(12), 1747–1754 (2009).CAS 
Article 

Google Scholar 
14.Tang, X. & Yang, J. Long-term stability of risk assessment of lead mill waste treated by soluble phosphate. Sci. Total Environ. 438, 299–303 (2012).ADS 
CAS 
Article 

Google Scholar 
15.Brown, S. L. & Chaney, R. L. A rapid in-vitro procedure to a mimicked in-situ remediation study of metal-contaminated soils characterize the effectiveness of a variety of in-situ lead remediation with emphasis on Cd and Pb. J. Environ. Qual. 23, 58–63 (1997).
Google Scholar 
16.Li, Y. M., Chaney, R. L., Siebielec, G. & Kerschner, B. Response of four turf grass cultivars to limestone and biosolids-compost amendment of a zinc and cadmium contaminated soil at Palmerton, Pennsylvania. J. Environ. Qual. 29, 1440–1447 (2000).CAS 
Article 

Google Scholar 
17.Basta, N. T., Gradwhol, R., Snethen, K. L. & Schroder, J. L. Chemical immobilization of lead, zinc and cadmium in smelter-contaminated soils using biosolids and rock phosphate. J. Environ. Qual. 30, 1222–1230 (2001).CAS 
Article 

Google Scholar 
18.Brown, S., Chaney, R. L., Hallfrisch, J. G. & Xue, Q. Effects of biosolids processing on lead bioavailability in an urban soil. J. Environ. Qual. 32, 100–108 (2003).CAS 
Article 

Google Scholar 
19.Mosby, D.E., Miller, S., Bishop, C., Mehuys, J. Former and abandoned lead and zinc mines demonstration project. Final Report. Missouri Department of Nature Resources, Jefferson City, MO (2002).20.Singh, S. P., Ma, L. Q., Tack, F. G. & Verloo, M. G. Trace metal leachability of land-disposed dredged sediments. J. Environ. Qual. 29, 1124–1142 (2000).CAS 
Article 

Google Scholar 
21.Yang, J., Tang, X. & Wang, Z. Y. Water quality and ecotoxicity as influenced by phosphate and biosolid treatments in lead-contaminated soil and mine waste. J. Environ. Monit. Rest. 3, 21–33 (2007).
Google Scholar 
22.Ma, L. Q. & Rao, G. N. Chemical fractionation of cadmium, copper, nickel, and zinc in contaminated soils. J. Environ. Qual. 26, 259–264 (1997).CAS 
Article 

Google Scholar 
23.Bhattacharyya, P., Chakrabarti, K., Chakraborty, A., Tripathy, S. & Powell, M. A. Fractionation and bioavailability of Pb in municipal solid waste compost and Pb uptake by rice straw and grain under submerged condition in amended soil. Geosci. J. 12, 41–45 (2008).ADS 
CAS 
Article 

Google Scholar 
24.Brady, N. & Weil, R. The Nature and Property of Soils (Prentice Hall, 2002).
Google Scholar 
25.Fu, H. et al. Cadmium and lead speciation as affected by soil amendments in calcareous soil. Environ. Eng. Sci. https://doi.org/10.1089/ees.2017.0307 (2017).Article 

Google Scholar 
26.Xu, J. C., Huang, L. M., Chen, C., Wang, J. & Long, X. X. Effective lead immobilization by phosphate rock solubilization mediated by phosphate rock amendment and phosphate solubilizing bacteria. Chemosphere 237, 124540 (2019).ADS 
CAS 
Article 

Google Scholar 
27.Kungolos, A. et al. Toxic properties of metals and organotin compounds and their interactions on daphnia magna and vibrio fischeri. Water Ai, Soil Pollut. 4, 101–110 (2004).CAS 
Article 

Google Scholar 
28.Eighmy, T. T. et al. Heavy metal stabilization in municipal solid waste combustion dry scrubber residue using soluble phosphate. Environ. Sci. Technol. 31, 3330–3338 (1997).ADS 
CAS 
Article 

Google Scholar 
29.Crannell, B. S. et al. Heavy metal stabilization in municipal solid waste combustion bottom ash using soluble phosphate. Waste Manag. 20, 135–148 (2000).CAS 
Article 

Google Scholar 
30.Bubb, J. M. & Lester, J. N. Impact of heavy metals on low land rivers and implications for the man and the environment. Sci. Total Environ. 100, 207–258 (1991).ADS 
CAS 
Article 

Google Scholar 
31.Zhang, M. K. et al. Solubility of phosphorus and heavy metals in potting media amended with yard waste-biosolids compost. J. Environ. Qual. 33, 373–379 (2004).CAS 
Article 

Google Scholar 
32.Xian, X. Effect of chemical forms of cadmium, zinc, and lead in polluted soils on their uptake by cabbage plants. Plant Soil 113, 257–264 (1989).CAS 
Article 

Google Scholar 
33.Zhang, M. K. et al. Phosphorus and heavy metal attachment and release in sandy soil aggregate fractions. Soil Sci. Soc. Am. J. 67, 1158–1167 (2003).ADS 
CAS 
Article 

Google Scholar 
34.Pierzynski, G. M. & Schwab, A. P. Bioavailability of zinc, cadmium, and lead in a metal contaminated alluvial soil. J. Environ. Qual. 22, 247–254 (1993).CAS 
Article 

Google Scholar 
35.Samaras, V., Tsadilas, C.D. Distribution and availability of six heavy metals in a soil treated with sewage sludge. In Proceedings of International Conference Biogeochemistry of Trace Elements, Berkeley, CA (1997).36.Scheckel, K. G. & Ryan, J. A. Spectroscopic speciation and quantification of lead in phosphate-amended soils. J. Environ. Qual. 33, 1288–1295 (2004).CAS 
Article 

Google Scholar 
37.Lang, F. & Kaupenjohann, M. Effect of dissolved organic matter on the precipitation and mobility of the lead compound chloropyromorphite in solution. Eur. J. Soil Sci. 54, 139–147 (2003).CAS 
Article 

Google Scholar 
38.Shi, Q. et al. Lead immobilization by phosphate in the presence of iron oxides: Adsorption versus precipitation. Water Res. 179, 115853 (2020).CAS 
Article 

Google Scholar 
39.Andrunik, M., Wolowiec, M., Wojnarski, D., Zelek-Pogudz, S. & Bajda, T. Transformation of Pb, Cd, and Zn minerals using phosphates. Minerals. 10, 342 (2020).CAS 
Article 

Google Scholar 
40.Guo, J. H. et al. Stablizing lead bullets in shooting range soil by phosphate-based surface coating. AIMS Environ Sci. 3(3), 474–487 (2016).CAS 
Article 

Google Scholar  More

This section summarizes how we track a population’s biological resource demand and domestic availability. We also explain which income metrics we chose. A more complete discussion of the resource metric method is included in the Supplementary Methods.Measuring the biological resource balanceThe sustainable development literature has consistently recognized the importance of biological resource security. For example, the foundational Brundtland report expressed it as the need to live “within the planet’s ecological means” or “in harmony with the changing productive potential of the ecosystem”43.These principles call for comparing biological resource regeneration with a population’s demand on nature. Since people’s demands compete for nature’s products and services, one way of measuring this relationship between regeneration and human demand is by tracking how much mutually exclusive, biologically productive area is necessary to provide the resource flows that people demand. Humans demand biologically productive areas in several quantifiable ways: production of food, fibre and timber; physical infrastructure such as roads and buildings; and absorption of waste, particularly the carbon dioxide from fossil fuel combustion. The total demand for biologically productive surfaces can be compared with the productive areas available that provide regeneration. Since the productivity of areas varies, they need to be measured not in terms of their physical extension, but in terms of biological regeneration they represent. For example, one can use a biologically productive hectare with world-average productivity as the common measurement unit that then allows expression of both demand and availability of productive areas in units that become comparable across space and time.Ecological footprint accounting is a well-documented concept to measure the total supply and demand of biological regeneration. In ecological footprint accounting, the ecosystem capacity to regenerate biomass is called biocapacity. It is measured in standardized ‘global hectares’, which represent the productivity of a world-average biologically productive hectare. The human demand for biocapacity is called the population’s ‘ecological footprint’, and it is the sum of all the mutually exclusive demands on these bioproductive areas. Ecological footprints are also expressed in global hectares.The principles of ecological footprint accounting, and the derived methods for national and sub-national assessments, are documented extensively within scientific literature6,7,8,9,38,44,45,46. The national accounting methodology has also been reviewed and documented by numerous national government agencies47.The essence of the approach is that regeneration is used as the lens to analyse both availability and demand because biological assets are materially the most limiting factor of the human economy1,2. In addition, biocapacity and ecological footprint can be tracked and compared with each other on the basis of two principles:

1.

By scaling every area proportionally to its biological productivity, each biologically productive area becomes commensurable with any other one. This is the essence of the global hectare.

2.

By including only areas that exclude other uses, that is, by making sure that every area is counted only once, the areas can meaningfully be added up, both for all the competing demands on productive surfaces (the ecological footprint) and for the surfaces that contain the planet’s regenerative capacity (the biocapacity).

The country-level accounts, called the National Footprint and Biocapacity Accounts, show that humanity’s demand exceeds Earth’s biocapacity, and the gap has been increasing since the 1970s8,38,40. This is consistent with research on planetary boundaries or ecosystem health1,2,10,11.Countries’ resource demand can be analysed from a consumption or a production perspective. The consumption perspective, which is the one used in this study, adjusts for trade and indicates the total resource consumption demand of a population. The production perspective identifies how much demand activities within a country directly put on ecosystems. This could be interpreted as the demand associated with generating the country’s GDP.Countries that demand more than their domestic ecosystems regenerate run a biocapacity deficit. It is made possible by three mechanisms: (1) overuse of domestic ecosystems, or local overshoot; (2) net import of biocapacity; and (3) use of the global commons, as in the case of emitting CO2 from fossil fuel into the atmosphere or fishing international waters38.Global results indicate that as of 2017, Earth had about 12.1 billion biologically productive hectares, according to Food and Agriculture Organization land-use statistics48. This includes productive ocean areas. By definition, this equals 12.1 billion global hectares, as each global hectare represents the productive average of all these 12.1 billion hectares. By contrast, human demand in 2017 added up to 20.9 billion global hectares, 73% higher than the regeneration of all the planet’s ecosystems combined (in per-person numbers, an average footprint of 2.8 global hectares contrasted to 1.6 global hectares of biocapacity available per person worldwide). This 73% overshoot may have dropped to 56% in 2020 due to lockdowns during COVID-194. In 2017, ecological footprint country averages varied from 0.5 global hectares per person (Eritrea) to 14.7 global hectares per person (Qatar). Biocapacity averages among countries stretch from 0.1 global hectares per person (Singapore) to 84 global hectares per person (Suriname)40.The accounts include only human demands (including domesticated animals) and not those of the millions of other living species, which together make possible the continuous functioning of the global ecosystem. To maintain biodiversity, which is critical for the integrity of the global ecosystem, humanity’s footprint would need to be less than the planet’s total biocapacity. E.O. Wilson, for example, proposed to only use half the planet’s capacity to secure 85% of its current biodiversity49. Using this objective as reference would imply that humanity’s current biological metabolism would be three times too large. It also makes clear that zero biocapacity deficits are a necessary but not sufficient condition for planetary resource stability. Still, for simplicity, we use the zero biocapacity deficit line as the demarcation line.Currently, the single-largest competing demand on the biosphere is the need for carbon sequestration capacity to neutralize emissions from fossil fuel burning. In 2020, this demand made up 57% of humanity’s ecological footprint. To comply with the Paris Agreement’s stated goal (Article 2 of ref. 30), this portion of the footprint would need to fall rapidly to zero. This reduction may come at the cost of increasing other parts of the ecological footprint. For example, more forest or agricultural products may be used to substitute for fossil fuels. If the Paris Agreement is fully implemented, there will be legal pressure to eliminate the carbon-related part of the deficit. If it is not implemented, the reduction pressures will emerge more slowly, which will increase the likelihood that the biocapacity will become increasingly damaged by climate change. Taking either path forces a country to eliminate its carbon footprint one way or the other. Fossil fuel dependence is therefore turning into an ever-growing risk. The pressure of increased land use has historically been the leading factor in the extinction of biodiversity, but unless nations can effectively control climate change, it will soon predominate as the major factor responsible for the massive extinction event that we humans have already started as a result of our unsustainable consumption—the sixth such event in the history of our planet.Measuring ability to purchase resources from abroadAnnual value production of an economy is measured by its GDP. It can be calculated as the value add of all its produced goods and services, as the sum of all the incomes or as the sum of all expenditures. Therefore, GDP can be used as a measure for a country’s income50,51.The analysis here focuses on the relative purchasing power of countries’ economic actors on global markets. Therefore, we use nominal US$ (or for time series, constant US$) instead of purchasing-power-adjusted US$, which reflect purchasing power on local markets. As economic actors compete for global resources in the same global market, each dollar has approximately the same weight, independent of the dollar’s purchasing power in the actor’s domestic market (called purchasing power parity). While this simplifies the fact that many commodities do not have a single homogeneous global market, the price range for resources in international markets is much narrower than that between domestic markets.For the sake of this analysis, we use average country income. Although incomes within countries vary vastly, we assume that nominal per-capita GDP is a reasonable approximation for national purchasing power in international markets. As a medium of exchange, money gives its owner the option to trade it in for physical assets, including biological resources; hence, more money means access to more resources.Not all international resource transfers are traded on global markets. Purchases could be under the protection of government-to-government arrangements or long-term contracts. The more of the international resource exchanges that occur in global markets, the tighter the competition on the global market for the remaining resources. Such increased competition makes the implications of the analysis presented here even more dramatic.In the context of global ecological overshoot, biocapacity scarcity will increase; therefore, the competition for purchasing additional resources will become even fiercer. In this case, using world-average income as an approximation for the dividing line between those who can net-purchase from abroad and those who cannot is too lenient. This demarcation indicates only that statistically those above the line can net-purchase from abroad. It does not indicate, however, whether they can purchase enough from abroad to cover their biocapacity deficit. This means that even more national economies than those identified by the 72% in this paper are excluded from being able to purchase sufficient resources from abroad. More

In the case of the first studied facility—Facility (I), thermo-visual examinations did not demonstrate any signs of thermal activity in the majority of the area. Above all, the new layered sections of the facility are free from thermal phenomena. However, the central section of the dump is characterized by rather intense thermal phenomena. The cone is also thermally active, yet in this section the thermal activity is not very intensive (see Fig. 1). In the section with the strongest thermal activity, i.e. the southern part of the top, the recorded surface temperature exceeded 600 °C (see Fig. 2). Outside the thermally active zone, the surface temperature did not exceed 25–35 °C.Figure 1Thermo-visual examination of the extractive waste dump (Facility I); October, 2017.Full size imageFigure 2Subsidence (the so-called crater) at the top of Facility I; October, 2017.Full size imageThe tests conducted on the premises of Facility II demonstrated the lack of thermal activity in the majority of the area. Only in the central section, on the top of the older part of the dump, the measurements showed the occurrence of thermal anomalies (see Fig. 3). These were minor areas located along the edge of the scarp on the northern and southern sides where the measured surface temperature was above 80 °C. It is worth mentioning that in 2010, the recorded temperatures for this section were at the level of several hundred degrees Celsius. Outside the sections, the temperature did not exceed 25–40 °C. The relatively high temperature in the area without thermal activity may be explained by the fact that in the morning hours it was exposed to the sun.Figure 3Thermo-visual examination of the extractive waste dump (Facility II); June, 2017.Full size imageDuring the course of the thermo-visual examination of the extractive waste dump in Facility III, no thermal activity was observed in the majority of the area. Only in the central section of the Facility, in the area which is still in operation, the measurements demonstrated the occurrence of thermal anomalies (see Fig. 4). These were minor areas where the local measured surface temperature was above 200 °C. Outside these areas, the temperature did not exceed 20–30 °C.Figure 4Thermo-visual examination of the extractive waste dump (Facility III); October, 2017.Full size imageThe simplest method of analyzing the differences and similarities between the studied objects is their visualization in the space of measured parameters. When two or three parameters are measured, such visualization is not problematic. However, in our study there are nine measured parameters. The graphic presentation of a 9-dimensional space is not possible. This is the rationale of applying the HCA which enables to analyze the similarities between the studied objects (three different dumping facilities in different periods of time) as well as the similarities between the parameters in object space. Nevertheless, the HCA does not allow to simultaneously analyze the relationships between the objects and the measured parameters. This problem was solved by the use of a color map of the experimental data, which enabled an in-depth interpretation of the data structure. In addition, the application of the color map facilitates highlighting the differences and similarities among the clusters showed in the dendrograms, and, in consequence, it helps to distinguish the facilities which are characterized by the highest or the lowest values of the measured parameters. Figure 5 demonstrates the dendrogram for 21 objects representing the studied dumping facilities in different periods of time in the space of 9 measured parameters, the dendrogram for the measured parameters in the object space as well as a color map presenting the values of the measured parameters for particular objects.Figure 5Dendrograms for (a) 21 objects representing the studied dumping facilities (see Table 1) in the space of 9 measured parameters; (b) parameters in the object space; (c) a color map presenting the values of measured parameters for particular dumping facilities.Full size imageBased on the dendrogram presented in Fig. 5a, a clear distinction of the examined objects representing the discussed dumping facilities into two clusters—A and B can be observed. Cluster A includes all samples representing Facility II in the whole period of the monitoring as well as samples representing Facility I taken in Quarter 3, 2017 and in Quarters 1–4, 2018 (objects nos. 1 and 3–6). All samples representing Facility III as well as two samples representing Facility I taken in Quarter 4, 2017 and in Quarter 4, 2018, respectively (objects nos. 2 and 6) are collected in Cluster B. In addition, within each of the clusters certain sub-clusters may be distinguished. Within Cluster A, the following sub-groups can be observed:

sub-group A1 collecting all samples taken from Facility II during the whole period of the monitoring (objects nos. 8–14);

sub-group A2 including samples taken from Facility I in Quarter III, 2017 and in Quarters 1–4, 2018 (objects nos. 1 and 3–6).

In turn, Cluster B contains the following three sub-groups:

sub-group B1 collecting two samples taken from Facility I in Quarter 4, 2017 and 2018, respectively (objects nos. 2 and 6) as well as samples representing Facility III in Quarter 3, 2017, Quarters 1–3, 2018 and Quarter 1, 2019 (objects nos. 15, 17–19 and 21);

sub-group B2 encompassing the remaining two samples taken from Facility III in Quarter 4, 2017 and in Quarter 4, 2018 (objects nos. 16 and 20).

The dendrogram obtained by means of Ward’s linkage method for the measured parameters in the space of 21 objects (see Fig. 5b) enables to distinguish the three principal clusters of parameters listed below:

class A collecting parameters nos. 2, 3 and 4 (describing the concentrations of acenaphtene, fluorene and phenanthrene);

class B containing parameters nos. 1 and 9 (describing the concentrations of naphthalene and chrysene);

class C including the remaining parameters nos. 5, 6, 7 and 8 (describing the concentrations of anthracene, pyrene, fluoranthene and B(a)anthracene).

The PAHs emissions from a burning mine waste dump must be carefully monitored due to their potential toxicity and genotoxicity38,39. PAHs have two different roots; one is incomplete combustion of organic matter, whereas the other one is their production in the geological formation when organic sediments were chemically transformed into fossil fuels. In our study, the first path, namely spontaneous coal waste combustion is observed40,41. It is also worth mentioning that PAHs pose significant human health hazards. The exposure to PAHs may result in skin, lung or stomach cancers in the human organism28.As mentioned before, a considerable drawback of the Hierarchical Clustering Analysis lies in that it does not allow for simultaneous interpretation of the dendrograms describing objects which represent particular samples taken from the extractive waste facility in different periods in the space of measured parameters and parameters in the object space. The lack of the possibility of interpreting the above relationships significantly limits the knowledge of the studied phenomena because the purpose of the analysis is to determine not only the differences among particular samples but also the underlying cause of such differences. Therefore, the dendrogram presenting the examined samples taken from the three selected facilities in different periods of time (see Fig. 5a) was juxtaposed with a color map of the experimental data (see Fig. 5c) demonstrating the values of the measured parameters arranged according to the order of the objects and parameters organization. The juxtaposition enables to determine the reason why the examined samples were distributed in such a way. In addition, the interpretation of the dendrogram for the objects in parameters space complemented with the color map of experimental data allows distinguishing samples which are characterized by the highest values of the measured parameters.Analyzing the dendrogram presented in Fig. 5a together with the color map of the experimental data, it can be observed that all samples within Cluster A were characterized by relatively lower values of the measured parameters. Moreover, sub-group A1 identified within Cluster A was characterized by relatively lowest concentrations of acenaphtene, fluorene, phenanthrene and pyrene (parameters nos. 2, 3, 4, and 7) among all of the examined samples taken from the three selected facilities in the whole period of the monitoring. It confirms the observations which were made previously with the use of a thermo-visual camera which demonstrated that there were no thermal phenomena in the majority of the waste dump areas. The observed sparse emissions of gases result from certain thermal anomalies occurring in the central section of Facility II as well as some small areas located along the northern and southern side of the scarp; however, in any of the places the temperature did not exceed 80 °C (see Fig. 3). The thermo-visual examinations of the remaining sections of the facility showed that the surface temperature did not exceed the level of 25–40 °C and it is a result of sun exposure rather than the occurrence of any thermal phenomena.Sub-group A2 which includes all samples taken from Facility I in Quarter 3, 2017 and Quarters 1–4, 2018 (objects nos. 1 and 3–6) differs from A1 objects in terms of relatively higher concentrations of acenaphtene, fluorene, phenanthrene and pyrene (parameters nos. 2, 3, 4, and 7). Additionally, within sub-group A2, the uniqueness of sample 1 taken in Quarter 3, 2018 (object no. 5) can be observed; it was characterized by the lowest concentration of chryzene (parameter no. 9) of all the examined samples.Despite the fact that the results of the thermo-visual analysis of Facility III in general did not show any major signs of thermal activity, the Hierarchical Clustering Analysis by means of which the samples were examined in terms of the emission of the 9 parameters indicated that some thermal phenomena take place at this dump.Furthermore, the analysis of the color map of the experimental data for the objects grouped in Cluster B makes it abundantly clear that thermal activity takes place in the majority of the monitored areas located at Facility III. What is more, for two samples taken from Facility I in Quarter 4, 2017 and Quarter 4, 2018 (objects nos. 2 and 6), thermal activity was observed. For the two objects, the concentrations of acenaphtene and fluorene are the highest among all of the examined samples (parameters nos. 2 and 3); similarly, the concentrations of phenanthrene, pyrene and B(a)anthracene (parameters nos. 4, 7, and 8) are also high. It confirms the observations made by means of a thermo-visual camera which demonstrated that there were no thermal phenomena if one considers the facility in its entirety. However, there are certain areas in the central section of the facility where thermal phenomena occur.In the case of the samples taken from Facility III in Quarter 3, 2017 as well as in Quarter 2 and Quarter 3, 2018 (objects nos. 15, 18, and 19) within sub-group B1, the values of all measured parameters are only slightly increased, which can indicate the occurrence of thermal phenomena. Yet, it must be kept in mind that the concentrations are decidedly lower than for the remaining samples taken from Facility III. It may be also an indication that self-heating of the coal extractive waste takes place. The other two samples taken from Facility III in Quarter 1, 2018 and Quarter 1, 2019 (objects nos. 17 and 21) are characterized by relatively highest concentrations of pyrene and B(a)anthracene (parameters nos. 7 and 8) as well as high concentrations of anthracene and fluorene (parameters nos. 5 and 6), which can result from an intense thermal activity.A similar observation can be made for the two samples taken from Facility III in Quarter 4, 2017 and Quarter 4, 2018 (objects nos. 16 and 20) which are characterized by the highest concentrations of naphthalene, anthracene, fluoranthene and chrysene of all the examined samples (parameters nos. 1, 5, 6 and 9) as well as high concentrations of phenanthrene (parameter no. 4). Although the results of the thermo-visual analysis of the whole area of Facility III did not demonstrate signs of thermal activity in most of its sections, there exist certain spots with intense thermal processes, which is confirmed by relatively high values of all the measured parameters. More

1.Baker, B. J. et al. Diversity, ecology and evolution of Archaea. Nat. Microbiol. 5, 887–900 (2020).PubMed 
Article 
CAS 

Google Scholar 
2.Baker, B. J., Appler, K. E. & Gong, X. New microbial biodiversity in marine sediments. Ann. Rev. Mar. Sci. 13, 161–175 (2020).PubMed 
Article 

Google Scholar 
3.Hug, L. A. et al. A new view of the tree of life. Nat. Microbiol. 1, 16048 (2016).CAS 
PubMed 
Article 

Google Scholar 
4.Kozubal, M. A. et al. Geoarchaeota: a new candidate phylum in the Archaea from high-temperature acidic iron mats in Yellowstone National Park. ISME J. 7, 622–634 (2013).CAS 
PubMed 
Article 

Google Scholar 
5.Jay, Z. J. et al. Marsarchaeota are an aerobic archaeal lineage abundant in geothermal iron oxide microbial mats. Nat. Microbiol. 3, 732–740 (2018).CAS 
PubMed 
Article 

Google Scholar 
6.Hua, Z. S. et al. Genomic inference of the metabolism and evolution of the archaeal phylum Aigarchaeota. Nat. Commun. 9, 2832 (2018).ADS 
PubMed 
PubMed Central 
Article 
CAS 

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

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

Google Scholar 
9.Orsi, W. D. et al. Metabolic activity analyses demonstrate that Lokiarchaeon exhibits homoacetogenesis in sulfidic marine sediments. Nat. Microbiol. 5, 248–255 (2020).CAS 
PubMed 
Article 

Google Scholar 
10.Eme, L., Spang, A., Lombard, J., Stairs, C. W. & Ettema, T. J. G. Archaea and the origin of eukaryotes. Nat. Rev. Microbiol. 15, 711–723 (2017).CAS 
PubMed 
Article 

Google Scholar 
11.Zaremba-Niedzwiedzka, K. et al. Asgard archaea illuminate the origin of eukaryotic cellular complexity. Nature 541, 353–358 (2017).ADS 
CAS 
PubMed 
Article 

Google Scholar 
12.Williams, T. A. et al. Integrative modeling of gene and genome evolution roots the archaeal tree of life. Proc. Natl Acad. Sci. U.S.A. 114, E4602 –E4611 (2017).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
13.Spang, A., Caceres, E. F. & Ettema, T. J. G. Genomic exploration of the diversity, ecology, and evolution of the archaeal domain of life. Science 357, eaaf3883 (2017).14.Trembath-Reichert, E. et al. Methyl-compound use and slow growth characterize microbial life in 2-km-deep subseafloor coal and shale beds. Proc. Natl Acad. Sci. USA 114, E9206–E9215 (2017).CAS 
PubMed 
Article 

Google Scholar 
15.Zhuang, G. C., Peña-Montenegro, T. D., Montgomery, A., Hunter, K. S. & Joye, S. B. Microbial metabolism of methanol and methylamine in the Gulf of Mexico: insight into marine carbon and nitrogen cycling. Environ. Microbiol. 20, 4543–4554 (2018).CAS 
PubMed 
Article 

Google Scholar 
16.Chistoserdova, L. Modularity of methylotrophy, revisited. Environ. Microbiol. 13, 2603–2622 (2011).CAS 
PubMed 
Article 

Google Scholar 
17.Chistoserdova, L. & Kalyuzhnaya, M. G. Current trends in methylotrophy. Trends Microbiol. 26, 703–714 (2018).CAS 
PubMed 
Article 

Google Scholar 
18.Sun, J., Mausz, M. A., Chen, Y. & Giovannoni, S. J. Microbial trimethylamine metabolism in marine environments. Environ. Microbiol. 21, 513–520 (2018).PubMed 
Article 

Google Scholar 
19.Zhuang, G.-C., Montgomery, A. & Joye, S. B. Heterotrophic metabolism of C1 and C2 low molecular weight compounds in northern Gulf of Mexico sediments: controlling factors and implications for organic carbon degradation. Geochim. Cosmochim. Acta 247, 243–260 (2019).ADS 
CAS 
Article 

Google Scholar 
20.Vanwonterghem, I. et al. Methylotrophic methanogenesis discovered in the archaeal phylum Verstraetearchaeota. Nat. Microbiol. 1, 16170 (2016).CAS 
PubMed 
Article 

Google Scholar 
21.Lazar, C. S. et al. Genomic evidence for distinct carbon substrate preferences and ecological niches of Bathyarchaeota in estuarine sediments. Environ. Microbiol. 18, 1200–1211 (2016).CAS 
PubMed 
Article 

Google Scholar 
22.Zhuang, G. Methylotrophic methanogenesis and potential methylated substrates in marine sediment. (University of Bremen, 2014).23.Richards, M. A. et al. Exploring hydrogenotrophic methanogenesis: a genome scale metabolic reconstruction of Methanococcus maripaludis. J. Bacteriol. 198, 3379–3390 (2016).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
24.Sousa, D. Z. et al. The deep-subsurface sulfate reducer Desulfotomaculum kuznetsovii employs two methanol-degrading pathways. Nat. Commun. 9, 239 (2018).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
25.Dombrowski, N., Teske, A. P. & Baker, B. J. Extensive metabolic versatility and redundancy in microbially diverse, dynamic Guaymas Basin hydrothermal sediments. Nat. Commun. 9, 4999 (2018).26.Fricke, W. F. et al. The genome sequence of Methanosphaera stadtmanae reveals why this human intestinal archaeon is restricted to methanol and H2 for methane formation and ATP synthesis †. J. Bacteriol. 188, 642–658 (2006).27.McKay L., et al. Co-occurring genomic capacity for anaerobic methane and dissimilatory sulfur metabolisms discovered in the Korarchaeota. Nat. Microbiol.4, 614–622 (2019).28.Muñoz-Velasco, I. et al. Methanogenesis on early stages of life: ancient but not primordial. Orig. Life Evol. Biosph. 48, 407–420 (2019).29.Adam, P. S., Borrel, G. & Gribaldo, S. An archaeal origin of the Wood–Ljungdahl H4MPT branch and the emergence of bacterial methylotrophy. Nat. Microbiol. 4, 2155–2163 (2019).30.Swan, B., Reifel, K. & Valentine, D. Periodic sulfide irruptions impact microbial community structure and diversity in the water column of a hypersaline lake. Aquat. Microb. Ecol. 60, 97–108 (2010).Article 

Google Scholar 
31.Adam, P. S., Borrel, G. & Gribaldo, S. Evolutionary history of carbon monoxide dehydrogenase/acetyl-CoA synthase, one of the oldest enzymatic complexes. PNAS 115, E5837 (2018).Article 
CAS 

Google Scholar 
32.Orita, I. et al. The ribulose monophosphate pathway substitutes for the missing pentose phosphate pathway in the archaeon Thermococcus kodakaraensis. J. Bacteriol. 188, 4698–4704 (2006).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
33.Urschel, M. R., Kubo, M. D., Hoehler, T. M., Peters, J. W. & Boyd, E. S. Carbon source preference in chemosynthetic hot spring communities. Appl. Environ. Microbiol. 81, 3834–3847 (2015).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
34.Yokohama, H., Wagner, I. D. & Wiegel, J. Caldicoprobacter oshimai gen. nov., sp. nov., an anaerobic, xylanolytic, extremely thermophilic bacterium isolated from sheep faeces, and proposal of Caldicoprobacteraceae fam. nov. Int. J. Syst. Evol. Microbiol. 60, 67–71 (2010).Article 
CAS 

Google Scholar 
35.Zhang, X. et al. Petroclostridium xylanilyticum gen. Nov., sp. nov., a xylan-degrading bacterium isolated from an oilfield, and reclassification of clostridial cluster iii members into four novel genera in a new hungateiclostridiaceae fam. nov.Int. J. Syst. Evol. Microbiol. 68, 3197–3211 (2018).CAS 
PubMed 
Article 

Google Scholar 
36.Girbal, L., Croux, C., Vasconcelos, I. & Soucaille, P. Regulation of metabolic shifts in Clostridium acetobutylicum ATCC 824. FEMS Microbiol. Rev. 17, 287–297 (1995).CAS 
Article 

Google Scholar 
37.Qi, F. et al. Improvement of butanol production in Clostridium acetobutylicum through enhancement of NAD(P)H availability. J. Ind. Microbiol. Biotechnol. 45, 993–1002 (2018).CAS 
PubMed 
Article 

Google Scholar 
38.Branduardi, P., Longo, V., Berterame, N. M., Rossi, G. & Porro, D. A novel pathway to produce butanol and isobutanol in Saccharomyces cerevisiae. Biotechnol. Biofuels 6, 68 (2013).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
39.Johnsen, U. & Schönheit, P. Novel xylose dehydrogenase in the halophilic archaeon Haloarcula marismortui. J. Bacteriol. 186, 6198–6207 (2004).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
40.Ravachol, J. et al. Mechanisms involved in xyloglucan catabolism by the cellulosome-producing bacterium Ruminiclostridium cellulolyticum. Sci. Rep. 6, 22770 (2016).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
41.Macdonald, S. S., Blaukopf, M. & Withers, S. G. N-acetylglucosaminidases from CAZy family GH3 are really glycoside phosphorylases, thereby explaining their use of histidine as an acid/Base catalyst in place of glutamic acid. J. Biol. Chem. 290, 4887–4895 (2015).CAS 
PubMed 
Article 

Google Scholar 
42.Wang, Y. et al. Environmental conditions constrain the distribution and diversity of Archaeal merA in Yellowstone National Park, Wyoming, U.S.A. Microb. Ecol. 62, 739–752 (2011).CAS 
PubMed 
Article 

Google Scholar 
43.Nunes, C. I. P. et al. ArsC3 from Desulfovibrio alaskensis G20, a cation and sulfate-independent highly efficient arsenate reductase. J. Biol. Inorg. Chem. 19, 1277–1285 (2014).CAS 
PubMed 
Article 

Google Scholar 
44.Silver, S. & Phung, L. T. Genes and enzymes involved in bacterial oxidation and reduction of inorganic arsenic. Appl. Environ. Microbiol. 71, 599–608 (2005).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
45.Colman, D. R., Lindsay, M. R. & Boyd, E. S. Mixing of meteoric and geothermal fluids supports hyperdiverse chemosynthetic hydrothermal communities. Nat. Commun. 10, 681 (2019).46.Zhou, Z. et al. Genome- and community-level interaction insights into carbon utilization and element cycling functions of Hydrothermarchaeota in hydrothermal sediment. mSystems 5, e00795-19 (2020).PubMed 
PubMed Central 
Article 

Google Scholar 
47.Rabus, R., Venceslau, S. S., Lars, W., Wall, J. D. & Pereira, I. A. C. A post-genomic view of the ecophysiology, catabolism and biotechnological relevance of sulphate-reducing prokaryotes. Adv. Micro. Physiol. 66, 55–321 (2015).CAS 
Article 

Google Scholar 
48.Tóth, A., Takács, M., Groma, G., Rákhely, G. & Kovács, K. L. A novel NADPH-dependent oxidoreductase with a unique domain structure in the hyperthermophilic Archaeon, Thermococcus litoralis. FEMS Microbiol. Lett. 282, 8–14 (2008).PubMed 
Article 
CAS 

Google Scholar 
49.Ma, K., Weiss, R. & Adams, M. W. W. Characterization of hydrogenase II from the hyperthermophilic archaeon Pyrococcus furiosus and assessment of its role in sulfur reduction. J. Bacteriol. 182, 1864–1871 (2000).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
50.Jenney, F. E. & Adams, M. W. W. Hydrogenases of the model hyperthermophiles. Ann. N. Y. Acad. Sci. 1125, 252–266 (2008).ADS 
CAS 
PubMed 
Article 

Google Scholar 
51.Van Haaster, D. J., Silva, P. J., Hagedoorn, P. L., Jongejan, J. A. & Hagen, W. R. Reinvestigation of the steady-state kinetics and physiological function of the soluble NiFe-hydrogenase I of Pyrococcus furiosus. J. Bacteriol. 190, 1584–1587 (2008).PubMed 
Article 
CAS 

Google Scholar 
52.Greening, C. et al. Genomic and metagenomic surveys of hydrogenase distribution indicate H2 is a widely utilised energy source for microbial growth and survival. ISME J. 10, 761–777 (2016).CAS 
PubMed 
Article 

Google Scholar 
53.Stetter, K. O. Hyperthermophiles in the history of life. Philos. Trans. R. Soc. B Biol. Sci. 361, 1837–1842 (2006).CAS 
Article 

Google Scholar 
54.Collins, T., Gerday, C. & Feller, G. Xylanases, xylanase families and extremophilic xylanases. FEMS Microbiol. Rev. 29, 3–23 (2005).CAS 
PubMed 
Article 

Google Scholar 
55.Schädel, C., Richter, A., Blöchl, A. & Hoch, G. Hemicellulose concentration and composition in plant cell walls under extreme carbon source-sink imbalances. Physiol. Plant. 139, 241–255 (2010).PubMed 

Google Scholar 
56.Chen, S. et al. The Great Oxidation Event expanded the genetic repertoire of arsenic metabolism and cycling. 117, 10414–10421 (2020).57.Rogers, K. L. & Schulte, M. D. Organic sulfur metabolisms in hydrothermal environments. Geobiology 10, 320–332 (2012).CAS 
PubMed 
Article 

Google Scholar 
58.Baker, B. J. et al. Genomic inference of the metabolism of cosmopolitan subsurface Archaea, Hadesarchaea. Nat. Microbiol. 1, 16002 (2016).CAS 
PubMed 
Article 

Google Scholar 
59.Hatzenpichler, R., Krukenberg, V., Spietz, R. L. & Jay, Z. J. Next-generation physiology approaches to study microbiome function at single cell level. Nat. Rev. Microbiol. 18, 241–256 (2020).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
60.Eren, A. M. et al. Anvi’o: an advanced analysis and visualization platform for ‘omics data. PeerJ 3, e1319 (2015).PubMed 
PubMed Central 
Article 

Google Scholar 
61.Kang, D. D., Froula, J., Egan, R. & Wang, Z. MetaBAT, an efficient tool for accurately reconstructing single genomes from complex microbial communities. PeerJ 3, e1165 (2015).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
62.Dick, G. J. et al. Community-wide analysis of microbial genome sequence signatures. Genome Biol. 10, R85 (2009).63.Darling, A. E. et al. PhyloSift: phylogenetic analysis of genomes and metagenomes. PeerJ 2, e243 (2014).PubMed 
PubMed Central 
Article 

Google Scholar 
64.Katoh, K. & Standley, D. M. MAFFT multiple sequence alignment software version 7: improvements in performance and usability. Mol. Biol. Evol. 30, 772–780 (2013).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
65.Nguyen, L. T., Schmidt, H. A., Von Haeseler, A. & Minh, B. Q. IQ-TREE: a fast and effective stochastic algorithm for estimating maximum-likelihood phylogenies. Mol. Biol. Evol. 32, 268–274 (2015).CAS 
PubMed 
Article 

Google Scholar 
66.Hyatt, D. et al. Prodigal: prokaryotic gene recognition and translation initiation site identification. BMC Bioinform. 11, 119 (2010).Article 
CAS 

Google Scholar 
67.Aramaki, T. et al. KofamKOALA: KEGG ortholog assignment based on profile HMM and adaptive score threshold. Bioinformatics 36, 2251–2252 (2020).CAS 
PubMed 
Article 

Google Scholar 
68.Jones, P. et al. InterProScan 5: genome-scale protein function classification. Bioinformatics 30, 1236–1240 (2014).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
69.Søndergaard, D., Pedersen, C. N. S. & Greening, C. HydDB: a web tool for hydrogenase classification and analysis. Sci. Rep. 6, 34212 (2016).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
70.Zhang, H. et al. DbCAN2: a meta server for automated carbohydrate-active enzyme annotation. Nucleic Acids Res. 46, W95–W101 (2018).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
71.De Anda, V. et al. MEBS, a software platform to evaluate large (meta)genomic collections according to their metabolic machinery: unraveling the sulfur cycle. Gigascience 6, 1–17 (2017).ADS 
PubMed 
PubMed Central 
Article 

Google Scholar 
72.Zhichao, Z. et al METABOLIC: High-throughput. profiling of microbial genomes for functional traits, biogeochemistry, and community-scale metabolic networks. Preprint at bioRxiv 761643 (2019).73.Rawlings, N. D., Morton, F. R., Kok, C. Y., Kong, J. & Barrett, A. J. MEROPS: the peptidase database. Nucleic Acids Res. 38, D227–D233 (2010).CAS 
PubMed 
Article 

Google Scholar 
74.Taboada, B., Estrada, K., Ciria, R. & Merino, E. Operon-mapper: a web server for precise operon identification in bacterial and archaeal genomes. Bioinformatics 34, 4118–4120 (2018).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
75.Lombard, V., Golaconda Ramulu, H., Drula, E., Coutinho, P. M. & Henrissat, B. The carbohydrate-active enzymes database (CAZy) in 2013. Nucleic Acids Res. 42, D490–D495 (2014).CAS 
PubMed 
Article 

Google Scholar 
76.Yu, N. Y. et al. PSORTb 3.0: improved protein subcellular localization prediction with refined localization subcategories and predictive capabilities for all prokaryotes. Bioinformatics 26, 1608–1615 (2010).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
77.Boyd, J. A., Woodcroft, B. J. & Tyson, G. W. GraftM: a tool for scalable, phylogenetically informed classification of genes within metagenomes. Nucleic Acids Res. 46, e59 (2018).78.Hua, Z. S. et al. Insights into the ecological roles and evolution of methyl-coenzyme M reductase-containing hot spring Archaea. Nat. Commun. 10, 4574 (2019).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
79.Chaumeil, P., Mussig, A. J., Parks, D. H. & Hugenholtz, P. Genome analysis GTDB-Tk: a toolkit to classify genomes with the Genome Taxonomy Database. Bioinformatics 36, 1925–1927 (2020).80.Kahle, D. & Wickham, H. ggmap: spatial visualization with ggplot2. R. J. 5, 144–161 (2013).Article 

Google Scholar  More

In this study we successfully applied a DNA metabarcoding approach to identify consumed food items of plant and animal origin in human stomach content samples, even when digestion was advanced and macroscopic inspection no longer possible. A wide panel of common and less common edible food items were found, including meat, fish, legumes, cereals, nuts, fruits and spices. So far, gastric content analyses in a forensic context are typically based on microscopic and macroscopic identification of food items (reviewed e.g. in1). However, this approach is characterised by low taxonomic resolution, low sensitivity, and proves ineffective when meal leftovers are rendered unidentifiable due to chewing and digestive processes. In the field of molecular ecology, studies on animals have shown that morphological identification of prey items in the stomach underestimates prey diversity, which is particularly true when digestion is advanced (e.g.25). The only study to date applying DNA metabarcoding to infer human diet was based on faecal samples and did not assess any animal components of diet, although including a controlled feeding trial of an animal-based diet12. The comparison of the obtained plant DNA sequences to self-reporting indicated that, while some items were not reported but detected by DNA metabarcoding, all but one self-reported items were detected (the only exception being coffee), thus highlighting the sensitivity of the method. The present study, based on a random sampling of 48 human stomach contents collected during routine autopsies, includes a higher number of vegetal items and shows for the first time the successful detection of dietary items of animal origin. We found no correlation between the diversity of species detected and the time since death or digestion degree, which advocates for the utility of this methodology. The Vert01 primer set, highly specific to vertebrates, enables to distinguish between commonly eaten animal taxa and is clearly advantageous over morphological identification. In line with regional eating habits and previously published diet surveys26, we found within the 48 samples mainly pig, cattle/dairy and OTUs assigned to the plant families Poaceae, Rosaceae and Asteraceae (likely cereals, fruits, lettuces; Fig. 1). We did not detect coffee (Coffea spp.) in any of the stomach content samples, in line with12, which might be due to a degrading effect of roasting procedures on DNA, the absence of this popular beverage in all of the stomach samples being unlikely. Similarly, although common in Swiss eating habits, we also did not detect potato, which is usually eaten boiled or baked. Note that additional edible plant species, not listed in Fig. 2 since not constituting at least 10% of RRA but with 100% match with the database, were also detected (e.g. buckwheat, citrus fruits, flax, mangoes, sesame; Supplementary Table S2). Because we could obviously not compare our results to self-reported diets, we applied very stringent filtering parameters to avoid the occurrence of false positives (see Bioinformatic data treatment). It is beyond the approach of this study to distinguish between the animal source and a final processed food item (e.g. dairy or egg products) based on the obtained DNA sequences. However, this could be achieved by complementing the primer set with a bacterial marker (to e.g. identify the presence of a particular cheese27) or using proteomics (see below).Overall, the Vert01 metabarcode is able to discriminate well among commonly eaten genera. However, owing to its limited taxonomic resolution (72.4% at the species level, based on in silico testing11), species-level distinction is not always possible (e.g. between perch and pikeperch) or between potentially-eaten wild species and their conspecific domestic counterparts (e.g. wild boar and pig). In Fig. 2, we present the taxonomical assignation done using ObiTools together with a common name, selected after manually inspecting each sequence using BLAST and only considering 100% matches with edible species. In some cases, the common name refers to a group of species because the barcode was not specific enough to distinguish between genera or species. This is more relevant concerning plants, as the Sper01 metabarcode length ranges from 10 to 220 bp, implying that some items with shorter metabarcode and/or closely related phylogenetically could not be distinguished to genus or species level due to limited resolutive power. This is related to the nature of this universal plant marker, which has been designed to target a region of the trnL intron of chloroplast DNA which lacks taxonomic resolution within several plant families (only 21.5% resolution at the species level9,11) but has wide taxonomic coverage. This trade-off meant for our study that we could genetically not distinguish between some close species which are clearly different morphologically (e.g. stone fruits, cucurbits). To overcome this issue and increase the taxonomic resolution of the results, it is possible to envisage multiplexing within the same PCR of additional primers specifically targeting groups of species that cannot be identified at the species level by the P6 loop of the trnL intron. Such a strategy has already been implemented to distinguish between Carpinus betulus and Corylus avellana in bison diet28. Furthermore, it must be outlined that by using these primer sets only, diet assessment is not comprehensive as it does not target all possibly present food products. Even so-called universal primers may result in preferential amplification of some taxa over others and non-amplification of target taxa29,30. For this pilot study, we chose to use two universal PCR primer pairs with wide taxonomic coverage but limited specific resolution, in order to detect a broad range of items. To gain resolution for specific vertebrate or plant taxonomic groups (e.g. fish, birds, cereals) or target taxa not covered by these primers and which could be of forensic interest (e.g. marine crustaceans and molluscs, algae, fungi), it is possible to complement Vert01 and Sper01 with additional, taxonomically-restricted PCR metabarcoding primers described in the literature (e.g.31; examples reviewed in11). Taxonomic assignation of an unknown DNA sequence strongly depends on the exhaustiveness and quality of a reference database, either public as e.g. GenBank or custom-made/local (reviewed in32). In case of a priori knowledge of the overall consumed diet in samples, local databases may be restrained to the expected DNA sequences, which subsequently improves taxonomic assignment. For this study we in silico compiled databases containing all possible sequences amplified by our markers, but restricted these to vertebrates and spermatophytes (i.e. seed plants), respectively.The duration of stomach emptying has been estimated by the percentage of a meal present in a stomach3, but this process is influenced by several variables including the type and volume of consumed food, lifestyle and health, and can therefore last from few hours to days2. While one could argue that plant items usually remain longer in the stomach, our findings do not allow to draw robust conclusions about correlations of certain food items and digestion times. In order to establish hypotheses useful for time-frame estimations, additional experiments are necessary. In a controversial case of death, MS-based proteomics provided additional information through the analysis of food-derived proteins and peptides in the gastric content sampled at autopsy, indicating a last breakfast of milk and bread. While this method is certainly promising, it might reveal difficult if digestion is in an advanced stage, and has a less comprehensive scope than a DNA metabarcoding assay33. Furthermore, the effect of food processing techniques on DNA quality must be taken into account since cooking denatures e.g. proteins which in turn renders DNA amplification preferential to immunological approaches1. Different cooking treatments (variable duration of boiling, frying, baking) of tomato seeds showed that DNA extraction yielded in good quality DNA only for fresh seeds34, while digestion did not destroy DNA21. Hence, there might be an implicit bias of DNA metabarcoding to preferentially detect non-processed food (i.e. raw versus cooked). Another issue of environmental DNA-based methods is that it is not possible to distinguish between different states of food products based on DNA sequences. As mentioned before, we could not discriminate between e.g. grapes/wine, fruits/juices, beef meat/dairy products or chicken meat/eggs, since the DNA sequence of a derived product is identical to the DNA sequence of its source. While it is less common to encounter such biases for plants, mainly in cereal-derived products, it has to be taken into account when extrapolating diet patterns from DNA metabarcoding results.Stomach content sampling is invasive, but advantageous or even required with certain animal species and in particular circumstances, including definitely the human forensic context. An advantage of stomach content over faecal samples is that food is in an early stage of digestion before passing through the pyloric sphincter into the intestines, thus the effects of inhibition by bacteria or enzymes and degradation of DNA are less significant11,18. While some food particles such as seeds sometimes remain identifiable, even morphologically, after passing through the digestive system21, others do not and the same applies to DNA which is degraded by the digestive processes taking place in the intestinal tract. In a controlled feeding experiment on insects, the detectability of food DNA in different types of dietary samples showed that regurgitates and entire animals (including stomach content) outperformed faeces regarding detectability of prey DNA13. While food journals in dietary surveys may contain errors or deliberate omissions12, they are a comprehensive and easily accessible method of human diet assessment. However, in case of deceased persons that option is no longer available.Stomach content analyses provided crucial information for criminal investigations about cases of sudden and unexplained death on numerous occasions in recent years, enabling investigators to interpret perimortem events in detail (case examples reviewed in2). The results of this pilot study show that human stomach content analyses by DNA metabarcoding can be used as a complementary tool to traditional forensic macro- and microscopic approaches, with clear advantages such as an almost unlimited flexibility in terms of nature and range of taxa targeted, as well as high sensitivity and taxonomic resolution. Consequently, information that might otherwise remain undetected can be revealed, highlighting timings and circumstances surrounding the last hours of a person and his/her food intake. In a broader perspective, taking into account the potential improvements and refinements described above, and the growing amount of research literature available for wildlife species (i.e. environmental DNA-based studies), our results open up promising and novel prospects in the broader framework of human biomedical investigations of dietary patterns, based on partially or fully digested food found in the gastrointestinal tract or in faecal samples. More

The δ-MAPS analysis is performed onto monthly mean SST anomalies from the Mediterranean Sea Physical Reanalysis (CMEMS MED-Physics25) over the period 1987–2017. The advantage of using a reanalysis resides in the availability of a velocity field consistent with the SSTs that allows us to confirm the coupling between network domains and ocean currents within the euphotic layer.Validation of the δ-MAPS frameworkThe proposed ecoregionalization is first applied to the 2007–2010 period, when the domains, representing ecoregions, can be compared to those identified by Ref.2 using Lagrangian methods. The details of this validation are reported below and relevant figures can be found in the Supplementary Information.The 2007–2010 ecoregions in Figure S1 are consistent with the ones derived in Ref.2 through computationally intensive simulations. The name of each domain corresponds with those used in Ref.2 to ease comparison. It is worthwhile remarking that this work and the one of Ref.2 not only use very different methods to define connectivity, but also different data sources. Our study uses velocity and SST output fields from CMEMS MED Physics reanalysis, while Ref.2 uses the configuration PSY2V3 of the operational system MERCATOR OCEAN with a resolution of 8 km in the horizontal downscaled to a connectivity grid of 50 × 50 km. The data assimilation and clustering algorithms are different and Ref.2 employs a cut-off in addition to the clustering grid downscaling. These differences unavoidably translate into slightly different shapes and patterns of the domains inferred. For example, the D + V area in panel (a) of Figure S1 is effectively two separate ecoregions in Ref.2, in which the Messina Strait is not resolved at the connectivity grid level. However, this separation appears inconsistent with the surface kinetic energy (K.E.) of panel (b) in Figure S1, computed from the horizontal currents, e.g. zonal (u) and meridional (v) velocity components, as K.E. = 1/2 |V|2 where |V|= (u2 + v2)0.5. Indeed, there is no clear separation between the regions north and south at Messina Strait in our dataset. Having detailed this example and acknowledged that some differences should be expected, the overall basin eco-regionalization using δ-MAPS is consistent with that in Ref.2. The spatial accuracy is enough to well separate the main ecological areas, despite small-scale differences (i.e. some km, due to resolution choices).By and large, the SST anomaly domains in Figure S1 are bounded by ocean currents, in agreement with Ref.2. This is due to the dominance of advective forcing by ocean currents on the SSTs at equatorial and mid latitudes, on monthly timescales and spatial scales of few hundreds kilometers24. This link, which is foundational to the proposed methodology, is further quantified as follows: First, we calculate the surface K.E. per unit mass averaged over the time slot of interest; second, we select the points in the validation period (2007–2010) that exceed the 50th percentile of surface K.E. computed for the entire basin over 1987–2017 (e.g. 0.004 m2/s2); third, we compute the domain-boundary matrix augmented by 1 grid point in each direction; finally, we count which fraction of the domain boundaries computed in the boundary matrix overlaps with the K.E. fronts (above the 50th percentile threshold). The fraction obtained is high and equal to 0.73, and remains elevated when increasing the threshold to the 60th percentile (0.66). This procedure was repeated for all the time slots with Δ = 7 years used next in this study, obtaining high and very stable values in each case (mean ± variance = 0.73 ± 0.01 for the 50th percentile threshold, and 0.65 ± 0.01 for the 60th percentile threshold).Additionally, the correlation between the surface K.E. and K.E. at 50 m, 100 m, and 150 m over the whole 1987–2017 period (Figure S2) remains positive and significant, with coefficients for the whole domain (field mean c.c.  ± variance) of 0.83 ± 0.04 at 50 m, 0.68 ± 0.05 at 100 m, and 0.54 ± 0.06 at 150 m, indicating that the link extends to the whole euphotic layer.Mediterranean Sea ecoregions: long-term changesThe space-averaged (e.g. averaged on the whole basin) SSTs over the 1987–2017 period are characterized by a linear warming trend of about 0.04 °C per year, stronger in the eastern portion of the basin (Figure S3 in Supplementary Information). Over the same period, the K.E. per unit mass is characterized by different trends over decadal or quasi-decadal periods (Fig. 2, shown for surface only but the trend extends similarly to 50 m and 100 m depths) and no clear east–west contrast. A positive trend is found in the first part of the curve (1987–2001, 2.3 × 10–4 m2/s2 per year, green line in figure), followed by a central decade without statistically significant changes (2001–2010, blue line), and a steep negative trend afterward (2010–2017, – 4.1 10–4 m2/s2 per year, red line). We refer to 1987–2001, 2001–2010 and 2010–2017, as the UP, MAX and DOWN periods. The dynamical changes associated with the strengthening and weakening of ocean currents are hypothesized to coincide with a reshaping of the sub-basin ecoregions and reciprocal connectivity. The ecoregionalization inference is therefore performed considering time slots of varying length, so that yrend = yrini + Δ with yrini = y0 + n, n = 0,1,…,N, where y0 is the initial year of the dataset (1987) and N is the total number of time slots, each of duration Δ years, between 6 and 8. Time slots overlapping by more than one year among different trends periods are excluded. The choice of Δ = 7 years represents the best trade-off for having enough time slots to quantify the evolution of ecoregions and a sufficiently large number of data points in each time slot for statistical inference. We will focus on this case, but results are verified also for the other Δ values (see Supplementary Information).Figure 2Mean surface kinetic energy timeseries. Monthly time series of deseasonalized surface kinetic energy per unit mass (m2/s2), averaged over the whole Mediterranean Sea between 1987 and 2017. The shaded areas indicate the 1987–1993 (during the UP period), 2004–2010 (during MAX) and 2011–2017 (during DOWN) time slots used in Fig. 3.Full size imageStrength maps for three representative time slots are presented in Fig. 3a,c,e while maps of domain strengths for all Δ = 7 time slots can be found in Figure S4. The mean surface kinetic energy averaged within each timeslot is next compared to the number of ecoregions in corresponding timeslots. The fragmentation level, or the total number of ecoregions, and the mean surface kinetic energy content are highly correlated (Figure S5b in Supplementary Information), with a Pearson’s coefficient of 0.79 for the whole Mediterranean Sea, and 0.8 (0.65) for the eastern (western) basin. The fact that time slots are not independent does not invalidate the analysis, and a large correlation (c.c = 0.73) is retained even when using four non-overlapping time slots. A higher fragmentation occurs whenever the upper ocean layer is more energetic, and this relationship is robust to changes of Δ (see Supplementary Information). The domain strength is next compared to the mean K.E. content. For each timeslot, the domain strength is spatially averaged over the eastern and western basin separately. The correlations between the averaged strengths and the corresponding time slot mean surface K.E. values, both varying as the time slots change, are then calculated for eastern and western basins separately. No linkage is found in the western basin, but a strong anticorrelation describes the relationship in the eastern Mediterranean (c.c. − 0.74). This anticorrelation remains high (− 0.73) also when the eastern basin strengths are related to the whole basin surface K.E. averaged over each timeslot.Figure 3Domains and connectivity networks for the domain containing the Suez Canal. The three 7-year timeslots selected as representative of the UP (a,b), MAX (c,d) and DOWN (e,f) periods. The color of the domains represents their strength (left column), and the red dot shows the location of the Suez Canal. Links in the connectivity nets (right column) are colored according to the correlation between (the domain containing) the Suez Canal and other domains as labeled. Only correlations stronger than 0.35 are plotted. (Domains maps visualization produced with Matlab R2018a, https://www.mathworks.com/).Full size imageWe hypothesize that the amount of K.E. associated with semi-permanent jets, currents or large mesoscale eddies, grouped here together and named KE fronts, can be used as an indicator of their role as connectivity modulators. We identify KE fronts applying a pattern recognition algorithm on the K.E. fields for each time slot. The resulting pictures are processed by an image segmentation technique, based on K-means clustering, to separate the K.E. in four clusters of increasing energy content. The maximum-intensity group is selected as indicators for KE fronts and the number of pixels contained in each cluster is counted and used to estimate the size or abundance of each one. The maximum-intensity cluster well represents the energy-containing structures as measured by the correlation between the mean surface K.E. content in each time slot and the pixels within the corresponding cluster (c.c.  > 0.99). The more pixels reside within each cluster, the larger the KE fronts-populated areas that this cluster approximates. This estimation is carried out for the whole basin, and separately in the eastern and western parts. The number of pixels is then correlated to the number of inferred ecoregions for the whole Mediterranean (c.c. = 0.81), and for eastern (c.c. = 0.81) and western (c.c. = 0.69) basins. Figure S6 in the Supplementary Information compares the clustering maps of a low energy time slot (1987–1993, in panel (a)) and a higher one (2004–2010 in panel (b)), for the whole Mediterranean Sea for the maximum cluster. The number of ecoregions is highly correlated with the KE fronts everywhere and especially in the eastern Mediterranean Sea. The higher level of fragmentation found in the MAX period is thus associated with more abundant and/or larger surface KE fronts, acting as eco-dynamical barriers.To further strengthen this assessment, we consider that energy fronts can act as modulators for SSTa-derived domains. The ecoregionalization over a certain time slot characterizes that time range in one single ecoregion-map but stems from data known at several time points (i.e. monthly SSTa in our case). The resulting domains account therefore for the inherent physical variability of the system over time. A higher (lower) ecoregions fragmentation may therefore by associated with dynamical fronts occurring at different times and not necessarily in the same place, over a certain time range. If this is plausible, we expect to count more (less) occurrences of higher energy in broad areas where the domains are more (less) fragmented. For each time slot, the number of occurrences of a front in each pixel is therefore counted. Specifically, having defined a front as a K.E. realization above the 50th percentile of the overall (1987–2017) time varying surface K.E., we count how many times a front appears in the considered time slot at each pixel. In Figure S7 pixels are colored according to the number of occurrences in each time slot. The result is consistent with the domain fragmentation evolution. The higher fragmentation occurring in timeslots from 2001 to 2010 in the eastern basin is associated with more frequent fronts. Similar considerations hold for the other sub-basins, clearly distinguishing low energy periods from higher ones.Mediterranean ecoregions connectivity networksChanges in functional networks or connectivity among ecoregions can be assessed by comparing a network from each energy period (UP: 1987–1993, MAX: 2004–2010 and DOWN: 2011–2017) (Fig. 3b,d,f for the eastern basin and Figure S8 in the Supplementary Information for the western basin).In 1987–1993 the western basin was characterized by a high mean positive correlation of 0.73, with a strong, non-directional connectivity among the Tyrrhenian and Ligurian-Algero Provençal domains. In 2004–2010 the connectivity was overall weaker, and in particular reduced among Tyrrhenian waters. The connectivity between the Balearic domain (Bal) and the Tyrrhenian ones was also reduced. In 2011–2017 the connectivity was mostly recovered, especially in Tyrrhenian waters. In this period, the Algero-Provençal domain separated from the Ligurian Sea (Lig), enforcing its connectivity with the Balearic and the Alboran ecoregions.In the eastern basin we focus our attention on the ecoregion immediately offshore the Suez Canal (Fig. 3), the major anthropogenic corridor for the introduction of non-indigenous marine species in the Mediterranean Sea, the so-called Lessepsian immigrants32. According to δ-MAPS, connectivity from the domain surrounding Suez was high in the first decade, decreased approaching MAX, remained small until about 2010–2011 with fewer statistically significant links, and increased again in the more recent time slot considered. During the UP and DOWN periods, the strongest connections were with the eastern Levantine (domain N), followed by that with the Aegean, Ionian and Tunisian Seas. During UP the connectivity extended to the Provençal and Algerian Seas, in the western basin, while in DOWN these links were absent and replaced by a connection with the Adriatic Sea.The 1987–1993 and 2011–2017 periods, while not too dissimilar in energy levels, differed indeed for the phase of the Ionian-Adriatic Bimodal Oscillating System or BiOS33,34. The BiOS is a mode of variability characterized by a decadal reversal of the Northern Ionian Gyre (NIG) from cyclonic to anticyclonic, and vice versa. In its anticyclonic spinning the NIG deviates the inflowing Modified Atlantic Water (MAW) from the Sicily Channel towards the northern Ionian, entering the Adriatic Sea and decreasing its salinity and temperature. This prevents a portion of the MAW from reaching the Levantine basin, and enhances the outflow of Levantine waters into the western basin, along a pathway that follows the African coastline. The anticyclonic NIG co-occurs with higher concentrations of Atlantic and Western Mediterranean organisms in the Adriatic Sea. When the NIG is cyclonic, on the other hand, Levantine waters enter the Adriatic Sea, whereas the MAW preferably flows toward the Levantine35 and Lessepsian migrations influence the Adriatic Sea at various latitudes, affecting also phytoplankton phenology33,36,37. The corresponding regions and connectivity networks in the two opposite NIG periods are detailed in Figure S9 in the Supplementary Information. More

1.Becker, G. S. & Lewis, H. G. On the interaction between the quantity and quality of children. Journal of Political Economy 81(2, pt2), s279–s288 (1973).Article 

Google Scholar 
2.Bulatao, R. A. & Lee, R. D. Determinants of Fertility in Developing Countries (Academic Press, 1983).
Google Scholar 
3.Caldwell, J. C. The mechanisms of demographic change in historical perspective. Popul. Stud. 35(1), 1–27 (1981).Article 

Google Scholar 
4.Carlsson, G. The decline of fertility: innovation or adjustment process. Popul. Stud. 20(2), 149–174 (1966).CAS 
Article 

Google Scholar 
5.Easterlin, R. A. & Crimmins, E. M. The Fertility Revolution (University of Chicago Press, 1985).
Google Scholar 
6.Winterhalder, B. & Leslie, P. Risk-sensitive fertility: The variance compensation hypothesis. Evol. Hum. Behav. 23(1), 59–82. https://doi.org/10.1016/S1090-5138(01)00089-7 (2002).Article 

Google Scholar 
7.Sear, R., Lawson, D. W., Kaplan, H. & Shenk, M. K. Understanding Variation in Human Fertility: What Can We Learn from Evolutionary Demography? (The Royal Society, 2016).
Google Scholar 
8.Wood, J. W. Dynamics of Human Reproduction (Aldine de Gruyter, 1994).
Google Scholar 
9.Hruschka, D. J. & Burger, O. How does variance in fertility change over the demographic transition?. Philos. Trans. R. Soc. B 371, 20150155. https://doi.org/10.1098/rstb.2015.0155 (2016).Article 

Google Scholar 
10.Henry, L. Some data on natural fertility. Eugen. Q. 8(2), 81–91 (1961).CAS 
PubMed 
Article 

Google Scholar 
11.Campbell, K. L. & Wood, J. W. Fertility in traditional societies. In Natural Human Fertility Social and Biological Determinants (eds Diggory, P. et al.) 39–61 (MacMillan Press, 1988).
Google Scholar 
12.Ellison, P. T. On Fertile Ground (Harvard University Press, 2001).
Google Scholar 
13.Colleran, H., Jasienska, G., Nenko, I., Galbarczyk, A. & Mace, R. 2015 Fertility decline and the changing dynamics of wealth, status and inequality. Proc. R. Soc. B: Biol. Sci. 282, 20150287 (1806).Article 

Google Scholar 
14.Colleran, H. The cultural evolution of fertility decline. Philos. Trans. R. Soc. B: Biol. Sci. 371(1692), 20150152 (2016).Article 

Google Scholar 
15.Knight, F. H. Risk, Uncertainty, and Profit (Houghton Mifflin, 1921).
Google Scholar 
16.Ellison, P. T. Energetics and reproductive effort. Am. J. Hum. Biol. 15, 342–351 (2003).PubMed 
Article 
PubMed Central 

Google Scholar 
17.Jasienska, G. & Ellison, P. Energetic factors and seasonal changes in ovarian function in women from rural Poland. Am. J. Hum. Biol. 16, 563–580 (2004).PubMed 
Article 
PubMed Central 

Google Scholar 
18.Jasienska, G. & Ellison, P. T. Physical work causes suppression of ovarian function in women. Proc. R. Soc. Lond. B 265, 1847–1851 (1998).CAS 
Article 

Google Scholar 
19.Kramer, K. L. & McMillan, G. P. The effect of labor saving technology on longitudinal fertility changes. Curr. Anthropol. 47(1), 165–172 (2006).Article 

Google Scholar 
20.Panter-Brick, C. Lactation, birth spacing and maternal workloads among two cases in rural Nepal. J. Biosoc. Sci. 23, 137–154 (1991).CAS 
PubMed 
Article 

Google Scholar 
21.Sear, R., Steele, F., McGregor, I. A. & Mace, R. The effects of kin on child mortality in Gambia. Demography 39(1), 43–63 (2002).PubMed 
Article 

Google Scholar 
22.Valeggia, C. R. & Ellison, P. T. Lactation, energetics, and postpartum fecundity. In Reproductive Ecology and Human Evolution (ed. Ellison, P. T.) 85–105 (Aldine de Gruyter, 2001).
Google Scholar 
23.Gibson, M. & Mace, R. An energy-saving development initiative increases birth rate and childhood malnutrition in rural Ethiopia. PloS Med. 3, 476–484 (2006).Article 

Google Scholar 
24.Kramer, K. L. & McMillan, G. P. Women’s labor, fertility, and the introduction of modern technology in a rural Maya village. J. Anthropol. Res. 55(4), 499–520 (1999).CAS 
PubMed 
Article 

Google Scholar 
25.Low, B., Simon, C. & Anderson, K. An evolutionary ecological perspective on demographic transitions: modeling multiple currencies. Am. J. Hum. Biol. 14(2), 149–167 (2002).PubMed 
Article 

Google Scholar 
26.Low, B. S., Simon, C. S. & Anderson, K. G. The biodemography of modern women: tradeoffs when resources become limiting. In The Biodemography of Human Reproduction and Fertility (ed. Rodgers, J. L.) 105–134 (Kuwer Academic Publishers, 2003).
Google Scholar 
27.Dyson, T. & Murphy, M. The onset of fertility transition. Popul. Dev. Rev. 11(3), 399–440 (1985).Article 

Google Scholar 
28.Early, J. & Headland, T. N. Population Dynamics of a Philippine Rain Forest People (University of Florida Press, 1998).
Google Scholar 
29.Hill, K. & Hurtado, A. M. Ache Life History (Aldine de Gruyter, 1996).
Google Scholar 
30.Kramer, K. L. & Greaves, R. D. Changing patterns of infant mortality and fertility among Pumé foragers and horticulturalists. Am. Anthropol. 109(4), 713–726 (2007).Article 

Google Scholar 
31.Goldstein, J. R. & Klüsener, S. Spatial analysis of the causes of fertility decline in Prussia. Popul. Dev. Rev. 40(3), 497–525 (2014).Article 

Google Scholar 
32.Montgomery, M. R. & Casterline, J. B. The diffusion of fertility control in Taiwan: Evidence from pooled cross-section time-series models. Popul. Stud. 47(3), 457–479 (1993).CAS 
Article 

Google Scholar 
33.Schmertmann, C. P., Assunção, R. M. & Potter, J. E. Knox meets Cox: Adapting epidemiological space-time statistics to demographic studies. Demography 47(3), 629–650 (2010).PubMed 
PubMed Central 
Article 

Google Scholar 
34.Galloway, P. R., Hammel, E. A. & Lee, R. D. Fertility decline in Prussia, 1875–1910: A pooled cross-section time series analysis. Popul. Stud. 48(1), 135–158 (1994).CAS 
Article 

Google Scholar 
35.Schmertmann C. P., Potter J. E. & Assunção R. M. 2011 An innovative methodology for space-time analysis with an application to the 1960–2000 Brazilian mortality transition. In Navigating Time and Space in Population Studies 19–36 (Dordrecht, 2011).36.Bongaarts, J., Cleland, J., Townsend, J. W., Bertrand, J. T. & Gupta, M. D. Family Planning Programs for the 21st Century (Population Council, 2012).
Google Scholar 
37.Casterline J. B. Diffusion processes and fertility transition: Introduction. In Diffusion processes and fertility transition: Selected perspectives (ed. Population N.R.C.C.o.) (National Academies Press, US, 2011).38.Cleland, J. The effects of improved survival on fertility: A reassessment. In Global Fertility Transitions Population and Development Review Supplement to Vol 27 (eds. Bulatao R.A., Casterline J.B.) (Population Council, 2001).39.Cleland, J. & Wilson, C. Demand theories of the fertility transition: an iconoclastic view. Popul. Stud. 41(1), 5–30 (1987).Article 

Google Scholar 
40.Montgomery, M. R. & Casterline, J. B. Social learning, social influence, and new models of fertility. Popul. Dev. Rev. 22, 151–175 (1996).Article 

Google Scholar 
41.Sear, R. Evolutionary contributions to the study of human fertility. Popul. Stud. 69(sup1), S39–S55 (2015).Article 

Google Scholar 
42.Knodel, J. & Van de Walle, E. Lessons from the past: Policy implications of historical fertility studies. Popul. Dev. Rev. 5(2), 217–245 (1979).Article 

Google Scholar 
43.Watkins, S. C. From local to national communities: The transformation of demographic regimes in Western Europe, 1870–1960. Popul. Dev. Rev. 16(2), 241–272 (1990).Article 

Google Scholar 
44.Bongaarts, J. & Watkins, S. C. Social interactions and contemporary fertility transitions. Popul. Dev. Rev. 22(4), 639–682 (1996).Article 

Google Scholar 
45.Boyd, R. & Richerson, P. Culture and the Evolutionary Process (Univ. Press, 1985).
Google Scholar 
46.Cavalli-Sforza, L. L. & Feldman, M. W. Cultural Transmission and Evolution: A Quantitative Approach (Princeton University Press, 1981).MATH 

Google Scholar 
47.Colleran, H., Jasienska, G., Nenko, I., Galbarczyk, A. & Mace, R. Community-level education accelerates the cultural evolution of fertility decline. Proc. R. Soc. B Biol. Sci. https://doi.org/10.1098/rspb.2013.2732 (2014).Article 

Google Scholar 
48.Conrad, C., Lechner, M. & Werner, W. East German fertility after unification: crisis or adaptation?. Popul. Dev. Rev. 22(2), 331–358. https://doi.org/10.2307/2137438 (1996).Article 

Google Scholar 
49.Easterlin, R. A. Towards a socio-economic theory of fertility: a survey of recent research on economic factors in American fertility. In Fertility and Family Planning: A World View (ed. Behrman, S. J.) 127–156 (University of Michigan Press, 1969).
Google Scholar 
50.Easterlin, R. A. An economic framework for fertility analysis. Stud. Fam. Plann. 6, 54–63 (1975).CAS 
PubMed 
Article 

Google Scholar 
51.Galloway, P. R., Lee, R. D. & Hammel, E. A. Infant mortality and the fertility transition: Macro evidence from Europe and new findings from Prussia. In From Death to Birth Mortality Decline and Reproductive Change (eds Montgomery, M. R. & Cohen, B.) 182–226 (National Academy Press, 1998).
Google Scholar 
52.Kaplan, H. A theory of fertility and parental investment in traditional and modern human societies. Yearb. Phys. Anthropol. 39, 91–135 (1996).Article 

Google Scholar 
53.Lee, R. D. & Bulatao, R. A. The demand for children: a critical essay. In Determinants of Fertility in Developing Countries (eds Bulatao, R. A. & Lee, R. D.) 233–287 (Academic Press, 1983).
Google Scholar 
54.Lesthaeghe R. & Wilson C. Modes of production secularization and the pace of the fertility decline in Western Europe 1870–1930 (1986).55.Turke, P. Evolution and demand for children. Popul. Dev. Rev. 15(1), 61–90 (1989).MathSciNet 
Article 

Google Scholar 
56.Colleran, H. Farming in transition: land and property inheritance in a rural Polish population. Soc. Biol. Hum. Aff. 78, 7–19 (2014).
Google Scholar 
57.González-Bailón, S. & Murphy, T. E. The effects of social interactions on fertility decline in nineteenth-century France: an agent-based simulation experiment. Popul. Stud. 67(2), 135–155 (2013).Article 

Google Scholar 
58.Shenk, M. K., Towner, M. C., Kress, H. C. & Alam, N. A model comparison approach shows stronger support for economic models of fertility decline. Proc. Natl. Acad. Sci. 110(20), 8045–8050. https://doi.org/10.1073/pnas.1217029110 (2013).ADS 
Article 
PubMed 

Google Scholar 
59.Alvergne, A., Gurmu, E., Gibson, M. A. & Mace, R. Social transmission and the spread of modern contraception in rural Ethiopia. PLoS ONE 6, e22515. https://doi.org/10.1371/journal.pone.0022515 (2011).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
60.Mace, R. & Colleran, H. Kin influence on the decision to start using modern contraception: a longitudinal study from rural Gambia. Am. J. Hum. Biol. 21, 472–477. https://doi.org/10.1002/ajhb.20940 (2009).Article 
PubMed 

Google Scholar 
61.Montgomery, M., Casterline, J. B. & Heiland, F. Social Networks and the Diffusion of Fertility Control (Population Council, 1998).Book 

Google Scholar 
62.Veile, A. & Kramer, K. L. Pregnancy, birth and babies: motherhood and modernization in a Yucatec village. In Maternal Health, Pregnancy-Related Morbidity and Death among Indigenous Women of Mexico & Central America (ed. Schwartz, D.) 205–224 (Springer, Berlin, 2018).
Google Scholar 
63.Snopkowski, K., Towner, M. C., Shenk, M. K. & Colleran, H. Pathways from education to fertility decline: a multi-site comparative study. Philos. Trans. R. Soc. B: Biol. Sci. 371(1692), 20150156 (2016).Article 

Google Scholar 
64.Schultz T. P. The Fertility Transition: Economic Explanations. Economic Growth Center Discussion Paper No. 833. Available at SSRN: https://ssrn.com/abstract=286291 (2001).65.Becker, S. O., Cinnirella, F. & Woessmann, L. Does women’s education affect fertility? Evidence from pre-demographic transition Prussia. Eur. Rev. Econ. Hist. 17(1), 24–44 (2013).Article 

Google Scholar 
66.Gandrud, C. simPH: an R package for illustrating estimates from cox proportional hazard models including for interactive and nonlinear effects. J. Stat. Softw. 65(3), 1–20 (2015).Article 

Google Scholar 
67.Seiber, E. E., Bertrand, J. T. & Sullivan, T. M. Changes in contraceptive method mix in developing countries. Int. Fam. Plan. Perspect. 33(3), 117–123 (2007).PubMed 
Article 

Google Scholar 
68.Leite, I. D. C., Gupta, N. & Rodrigues, R. D. Female sterilization in Latin America: cross-national perspectives. J. Biosoc. Sci. 36(6), 683 (2004).Article 

Google Scholar 
69.Bertrand, J. T., Sullivan, T. M., Knowles, E. A., Zeeshan, M. F. & Shelton, J. D. Contraceptive method skew and shifts in method mix in low-and middle-income countries. Int. Perspect. Sexual Reprod. Health 40(3), 144–153 (2014).Article 

Google Scholar 
70.Leslie, P. & Winterhalder, B. Demographic consequences of unpredictability in fertility outcomes. Am. J. Hum. Biol. 14(2), 168–183 (2002).PubMed 
Article 
PubMed Central 

Google Scholar 
71.Gibson, M. & Mace, R. Labor-saving technology and fertility increase in rural Africa. Curr. Anthropol. 43(4), 631–637 (2002).PubMed 
Article 
PubMed Central 

Google Scholar 
72.Alvergne, A., Lawson, D. W., Clarke, P. M. R., Gurmu, E. & Mace, R. Fertility, parental investment, and the early adoption of modern contraception in rural Ethiopia. Am. J. Hum. Biol. 25(1), 107–115 (2013).PubMed 
Article 
PubMed Central 

Google Scholar 
73.Mace, R., Allal, N., Sear, R. & Prentice, A. M. The uptake of modern contraception in a Gambian community: the diffusion of an innovation over 25 years. In Social Information Transmission and Human Biology (eds Wells, J. C. K. et al.) 191–206 (Taylor & Francis Group, 2006).
Google Scholar 
74.Lerner, I.M. Heredity, evolution and society. San Francisco: W.H. Freeman (1968)75.Lewontin, R. C. & Levins, R. Biology Under the Influence (Monthly Review Press, 2007).
Google Scholar 
76.Donaldson-Matasci, M. C., Lachmann, M. & Bergstrom, C. T. Phenotypic diversity as an adaptation to environmental uncertainty. Evol. Ecol. Res. 10(4), 493–515 (2008).
Google Scholar 
77.Meyers, L. A. & Bull, J. J. Fighting change with change: adaptive variation in an uncertain world. Trends Ecol. Evol. 17(12), 551–557 (2002).Article 

Google Scholar 
78.Sermonti, G. The butterfly and the lion. In Organisms, Genes and Evolution: Evolutionary Theory at the Crossroads; Proceedings of the 7th International Senckenberg Conference (eds Peters, S. T. & Weingarten, M.) 103 (Franz Steiner Verlag, 2000).
Google Scholar 
79.Boone, J. L. & Kessler, K. L. More status or more children? Social status, fertility reduction and long-term fitness. Evol. Hum. Behav. 20, 257–277 (1999).Article 

Google Scholar 
80.Nolin, D. A. & Ziker, J. P. Reproductive responses to economic uncertainty. Hum. Nat. 27(4), 351–371 (2016).PubMed 
Article 
PubMed Central 

Google Scholar 
81.Jensen, R. The (perceived) returns to education and the demand for schooling. Q. J. Econ. 125(2), 515–548 (2010).Article 

Google Scholar 
82.Borgerhoff, M. M. The demographic transition: are we any closer to an evolutionary explanation?. Trends Ecol. Evol. 13, 266–270 (1998).Article 

Google Scholar 
83.Skirbekk, V. Fertility trends by social status. Demogr. Res. 18, 145–180 (2008).Article 

Google Scholar 
84.Vining, D. R. J. Social verses reproductive success: the central theoretical problem of human sociobiology. Behav. Brain Sci. 9(167), 216 (1986).
Google Scholar 
85.Kramer, K. L. Maya Children: Helpers at the Farm (Harvard University Press, 2005).
Google Scholar 
86.Kramer, K. L. & Boone, J. L. Why intensive agriculturalists have higher fertility: a household labor budget approach to subsistence intensification and fertility rates. Curr. Anthropol. 43(3), 511–517 (2002).Article 

Google Scholar 
87.Lee, R. D. & Kramer, K. L. Children’s economic roles in the Maya family life cycle: Cain, Caldwell and Chayanov revisited. Popul. Dev. Rev. 28(3), 475–499 (2002).Article 

Google Scholar 
88.Kramer, K. L. Reconsidering the cost of childbearing: the timing of children’s helping behavior across the life cycle of Maya families. In SocioEconomic Aspects of Human Behavioral Ecology (ed. Alvard, M.) 335–353 (Elsevier, 2004).
Google Scholar 
89.Kramer, K. L., Veile, A. & Otárola-Castillo, E. Sibling competition, growth tradeoffs. Biological vs. statistical significance. PLoS ONE 11(3), e0150126. https://doi.org/10.1371/journal.pone.0150126 (2016).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
90.Veile, A. & Kramer, K. L. Shifting weanling’s optimum: breastfeeding ecology and infant health in Yucatan. In Anthropology and Breastfeeding (eds Tomori, C. et al.) Chapter 12 (Routledge Press, 2018).
Google Scholar 
91.Feltz, C. J. & Miller, G. E. An asymptotic test for the equality of coefficients of variation from k populations. Stat. Med. 15(6), 647–658 (1996).Article 

Google Scholar 
92.Marwick, B. & Krishnamoorthy. K. Cvequality: Tests for the Equality of Coefficients of Variation from Multiple Groups. R software package version 0.1.3 (2019). Retrieved from https://github.com/benmarwick/cvequality, on 05/01/2019.93.Cahoy, D. O. A bootstrap test for equality of variances. Comput. Stat. Data Anal. 54(10), 2306–2316 (2010).MathSciNet 
MATH 
Article 

Google Scholar 
94.Therneau, T. A Package for Survival Analysis in S. version 2.38 (2015), https://CRAN.R-project.org/package=survival. More