Ecology

Chemotaxonomic characteristicsThe predominant respiratory quinone for all novel strains was ubiquinone 8 (Q-8), consistent with other Ramlibacter species. C16:0 and summed feature 3 (consisting of C16:1 ω7c and/or C16:1 ω6c) were identified as the common major fatty acids ( > 10%) of the novel strains USB13T, AW1T, GTP1T, and HM2T. Other than the aforementioned fatty acids, strain USB13T had C10:0 3-OH additionally as its major fatty acid, whereas strains AW1T and HM2T shared C17:0 cyclo and summed feature 8 (consisting of C18:1 ω7c and/or C18: 1 ω6c) as its additional fatty acids. Detailed comparisons of the fatty acid profiles of the novel strains and their reference strains are summarized in Table S1.Strains USB13T, AW1T, GTP1T, and HM2T shared major polar lipids diphosphatidylglycerol (DPG), phosphatidylglycerol (PG), and phosphatidylethanolamine (PE), which was consistent with the major polar lipids of the reference strains. Additionally, the polar lipid profile of USB13T consisted of one unidentified phosphoaminolipid, two unidentified phosphoglycoaminolipids, and six unidentified polar lipids while the polar lipid profile of AW1T had one unidentified lipid, one unidentified phosphoglycolipid, and three unidentified glycolipids in addition. The polar lipid profile of strain GTP1T additionally consisted of two unidentified phosphoaminolipids, and that of strain HM2T additionally had one unidentified phosphoaminolipid, one unidentified phosphoglycolipid, one unidentified phosphoglycoaminolipid, and two unidentified phospholipids. Polar lipid profiles of the novel strains USB13T, AW1T, GTP1T, and HM2T are shown in Figure S1.Physiological, morphological characteristics, and screening of cellulose-degrading strainsWhen grown on R2A agar, strain USB13T produced reddish white and flat colonies while strain AW1T produced orange, convex colonies, strain GTP1T produced white, convex colonies, and strain HM2T produced cream-colored, flat, transparent colonies. Under TEM, monotrichous flagella were observed only in strain HM2T, and when tested for motility, strain USB13T and AW1T showed gliding motility, whereas strain GTP1T was non-motile. Strains USB13T and HM2T showed positive results for both catalase and oxidase activities; strain AW1T showed positive results for catalase and negative results for oxidase activity, and strain GTP1T showed negative results for catalase and positive results for oxidase activity. All strains were identified to be strictly aerobic, while showing negative results for urea, gelatin, starch, chitin, and DNA hydrolysis and positive results for hydrolysis of Tween 80. In addition, strain USB13T was the only strain to produce iron-chelating siderophores. When tested for NaCl tolerance, growth of strain USB13T was observed in NaCl concentrations of 0–7% (w/v), possibly due to the fact the strain was isolated from a marine environment. A detailed comparison of physiological and morphological characteristics between the novel species and its closely related Ramlibacter strains is presented in Table 1, while TEM images of the novel strains are shown in Figure S2. Results of the reference strains in Table 1 coincided with the data from the original literature1,3,4,5,7,8.Table 1 Characteristics differentiating strains USB13T, AW1T, GTP1T, and HM2T from closely related strains of the genus Ramlibacter.Full size tableStrains: 1, USB13T; 2, AW1T; 3, GTP1T; 4, HM2T; R. monticola KACC 19175T; 6, R. alkalitolerans KACC 19305T; 7, R. ginsenosidimutans KACC 17527T; 8, R. humi KCTC 52922T; 9, R. henchirensis KACC 11925T; 10, R. tataouinensis KACC 11924T; 11, R. rhizophilus KCTC 52083T. All strains are positive for esterase lipase (C8), while all strains are negative for chitin hydrolysis. All data were obtained from this study unless indicated otherwise. + , Positive; w + , weakly positive; -, negative.R2A agar plates supplemented with 1% (w/v) CMC were stained with Congo red dye after 7 days of incubation. Clear zones only formed around colonies of strain USB13T, indicating that strain USB13T solely possessed CMC-hydrolyzing activity among the four novel strains. When inoculated in basal salt medium, filter paper from the USB13T sample underwent degradation, whereas samples containing strains AW1T, GTP1T, and HM2T did not show any signs of degradation.Phylogenetic and genomic analysesEzBioCloud search results and BLASTn searches revealed that the novel strains belonged to the family Comamonadaceae and genus Ramlibacter. Using BLASTn, 16S rRNA gene sequence similarities were determined where strain USB13T was closest to strain GTP1T (98.5%), followed by strain HM2T (98.1%) and strain AW1T (97.1%). Strain AW1T shared the highest similarity with strain GTP1T (97.3%), followed by strain HM2T (97.1%), while strain GTP1T shared a similarity of 98.2% with strain HM2T. Phylogenetic analysis based on the MP method (Fig. 1) showed the clustering of the novel strains USB13T, AW1T, GTP1T, and HM2T with strains such as R. monticola G-3-2T, R. ginsenosidimutans BXN5-27T, R. alkalitolerans CJ661T, and R. rhizophilus YS3.2.7T. Similar topologies were observed in trees reconstructed by ML (Figure S3) and MP methods. The UBCG phylogenomic tree (Fig. 2), which was reconstructed using whole genome sequences, also showed close clustering of the selected reference strains and novel strains.Figure 1Maximum-parsimony (MP) tree reconstructed based on 16S rRNA gene sequences, showing the relationship between strains USB13T, AW1T, GTP1T, and HM2T and other closely related type strains. Bootstrap values based on 1000 replications are listed as percentages at branching points. Only bootstrap values exceeding 50% are shown. Bar, 50 substitutions per nucleotide position.Full size imageFigure 2Phylogenomic tree of strains USB13T, AW1T, GTP1T, and HM2T and their closely related taxa was reconstructed based on core genomes using UBCG version 3.0 pipeline42. NCBI GenBank accession numbers are shown in parentheses. Bootstrap analysis was carried out using 1000 replications. Percentage bootstrap values ( > 50%) are given at branching points. Bar, 0.050 substitution per position.Full size imageDraft genome sequences of the novel strains USB13T, AW1T, GTP1T, and HM2T were deposited in the GenBank database under the accession numbers JACORT000000000, JAEQNA000000000, JACORU000000000, and JADDIV000000000, respectively. In addition, the draft genome sequences of R. monticola KACC 19175T, R. alkalitolerans KACC 19305T, and R. ginsenosidimutans KACC 17527T were also deposited in GenBenk under the accession numbers JAEQNE000000000, JAEQND000000000, and JAEPWM000000000, respectively. The assembled genome size of the novel strains USB13T, AW1T, GTP1T, and HM2T was 5.53 Mbp, 5.11 Mbp, 6.15 Mbp, 4.31 Mbp, respectively. G + C content ranged from 67.9% to 69.9%, which was similar to those of the reference strains. The genomic features of the novel strains and their closely related Ramlibacter strains are presented in Table S2. CheckM analysis showed the following estimations for each strain: USB13T, had a 99.84% completeness and 0.68% contamination; AWIT, had a 99.84% completeness and 0.86% contamination; GTP1T, had a 99.38% completeness and 1.32% contamination; HM2T, had a 97.51% completeness and 0.16% contamination. These results indicated that the draft genome results for all strains were reliable. ANI values between the novel strains and reference strains ranged from 76.5–83.4% while dDDH values ranged from 20.7–26.7%, and AAI values ranged from 65.7–80.4%. All values were below the threshold for delineation of a new species54. ANI values between the novel strains and their reference strains are presented in Fig. 3, while a detailed comparison of GGDC and AAI values are shown in Table 2.Figure 3Heatmap of strains USB13T, AW1T, GTP1T, and HM2T and other closely related strains within the genus Ramlibacter, generated with OrthoANI values calculated using OAT software45. Bacterial strains and accession numbers are indentical to those of Fig. 2.Full size imageTable 2 Average amino acid identity (AAI) and digital DNA-DNA hybridization (dDDH) value comparisons between the closely related Ramlibacter type species and the novel strains, USB13T, AW1T, GTP1T, and HM2T. AAI values were calculated by two-way AAI, while dDDH values were calculated based on formula 246.Full size tableBased on NCBI PGAP annotation and CAZyme prediction results, strain USB13T, which was the only strain to show cellulolytic activity, possessed a total of four protein CDs encoding CAZymes, namely, two GH15 proteins, one glycosyl hydrolase protein, and one GH99-like domain-containing protein. Despite not showing any cellulolytic activity, strain AW1T possessed eight CAZyme CDs; the most amount among the novel strains. The enzymes include, two GH2 proteins, one GH5 protein, three GH15 proteins, one glycoside hydrolase protein, and one cellulase family glycosyl hydrolase. Strain GTP1T possessed two CDs encoding one GH15 protein and one GH16 protein; strain HM2T possessed three CDs encoding one GH2, one GH15, and one GH18 protein. All strains possessed GH15, which is known for its glucoamylase activity in fungi55. A detailed summary of the novel strains CAZymes are presented in Table S3 and a comparison of CAZyme numbers between strains USB13T, AW1T, GTP1T, and HM2T is summarized in Table S4. The presence of these genes may suggest the cellulolytic activity of strain USB13T, while it is uncertain why GH families responsible for endoglucanase (GH 5–8, 12, 16, 44, 45, 48, 51, 64, 71, 74, 81, 87, 124, and 128), exoglucanase (GH 5–7, and 48), and β-glucosidase (GH 1, 3, 4, 17, 30, and 116) were not present in the genome11.COG predictions (Fig. 4) revealed that the majority of the core genes of the four novel strains accounted for genes belonging to the functional categories C (energy production and conversion), E (amino acid transport and metabolism), I (lipid transport and metabolism), T (signal transduction mechanisms), and K (transcription). Meanwhile, the number of core genes belonging in category G, carbohydrate transport and metabolism, was the highest for strain USB13T (258), followed by GTP1T (230), HM2T (212), and AW1T (181). The high number of genes in strain USB13T may be a contributing factor in the strain’s cellulolytic activity. A comparison of COG gene count distribution of the novel strains is presented in Table S5.Figure 4Comparison of total number of matched genes of strains USB13T, AW1T, GTP1T, and HM2T according to functional classes based on Cluster of Orthologous Groups of proteins (COG) predictions48.Full size imageAntiSMASH analysis results showed four gene clusters within the genome of strain USB13T: ribosomally synthesized and post-translationally modified peptides (RIPP)-like cluster (989,516–1,000,916 nt; JACORT010000001), terpene synthesis (8,622–30,347 nt; JACORT010000003), RIPP precursor peptide recognition element (RRE)-containing cluster (311,469–333,619 nt; JACORT010000004), and redox-cofactor (281,860–303,948 nt; JACORT010000007). Among the clusters, the RRE-containing cluster showed 11% similarity to streptobactin, a tricatechol-type siderophore isolated from Streptomyces sp. YM5-79956. Strain AW1T had a total of eight gene clusters which encoded for: arylpolyene (165,946–207,130 nt), terpene (618,322–640,854 nt), RIPP-like proteins (804,411–819,137 nt), non-ribosomal peptide synthetase cluster (NRPS)-like (61,798–104,764 nt), betalactone (323,399–348,739 nt), N-acetylglutaminylglutamine amide (NAGGN; 106,834–121,648 nt), type I polyketide synthase (T1PKS; 56,584–107,578 nt), and heterocyst glycolipid synthase-like polyketide synthase (hglE-KS; 75,419–113,566 nt). Strain GTP1T possessed four gene clusters that encoded for RRE-containing cluster (175,155–199,102 nt), homoserine lactone (110,293–130,892 nt), a signaling molecule known for its involvement in bacterial quorum sensing, the RIPP-like cluster (38,002–48,856 nt), and terpene synthesis (47,942–69,701 nt). Strain HM2T had two gene clusters that encoded for resorcinol (403,967–445,901 nt), an organic compound known for its antiseptic properties, and terpene (697,660–721,242 nt), which showed 100% similarity for carotenoid synthesis. BRIG analysis results showed that a majority of the regions within the four analyzed genomes were conserved with at least 70% similarity (Figure S4).Cellulolytic potential and FE-SEM analysis of strain USB13T
A USB13T-inoculated basal salt medium sample containing degraded filter paper was examined under FE-SEM to observe the morphological interactions between cellulose fibers and USB13T cells. Images in Fig. 5 show individual rod cells of strain USB13T surrounding filter paper fibers, indicating bacterial adherence.Figure 5Field emission-scanning electron microscopy (FE-SEM) images of adhesion of strain USB13T to degraded filter paper fibers. Arrows indicate filter paper fibers. (A) low magnification (5000(times)) and (B), high magnification (20,000(times)) images of strain USB13T surrounding filter paper fibers.Full size imageThe enzymatic assay results showed endoglucanase, exoglucanase, β-glucosidase, and filter paper cellulase (FPCase) activities of strain USB13T, wherein activities for endoglucanase was the highest and β-glucosidase was the lowest in all experiments. As seen in Fig. 6A, enzyme activity for all cellulolytic enzymes increased along with its cultivation time. In addition, enzyme activities showed the highest results when tested on buffer solutions of pH 6.0 (Fig. 6B), indicating the enzymes’ resistance to moderately acidic conditions. The pH of the buffer solution seemed to be an important factor in enzyme activity, as activity of endoglucanase, exoglucanase, and FPCase drastically decreased when the pH was altered from pH 6.0 to pH 7.0. Meanwhile, β-glucosidase activity was relatively resistant to pH change as its activity decreased less than 50%. On day 7, enzyme activities were measured as 1.91 IU/mL for endoglucanase, 1.77 IU/mL for exoglucanase, 0.76 IU/mL for β-glucosidase, and 1.12 IU/mL for FPCase at pH 6.0. When measured at pH 8.0, where enzyme activity was the lowest, enzyme activities were measured as 0.51 IU/mL for endoglucanase, 0.25 IU/mL for exoglucanase, 0.45 IU/mL for β-glucosidase, and 0.23 IU/mL for FPCase; all values were less than half of the measured activity at pH 6.0. The results of strain USB13T are comparable to FPCase results of other species such as Mucilaginibacter polytrichastri RG4-7T (0.98 U/mL) isolated from the moss Polytrichastrum formosum14, Paenibacillus lautus BHU3 (2.9 U/mL) isolated from a landfill site57, and Serratia rubidaea DBT4 (0.5 U/mL) isolated from the gastrointestinal tract of a black Bengal goat58.Figure 6Cellulolytic enzyme activity of strain USB13T. Enzyme activity was defined in international units (IU); one unit of enzymatic activity was defined as the amount of enzyme that releases 1 μmol of glucose per mL per 1 min of reaction. (A) cellulase activity results under different cultivation time; (B) cellulase activity under different buffer solution pH. Values in the figure are mean values of triplicate data with standard deviation.Full size imageDespite the absence of the main three cellulolytic enzymes, endoglucanase, exoglucanase, and β-glucosidase, the cellulolytic activity of strain USB13T was confirmed through SEM images, CMC agar screening, and enzymatic assay results. However, because PGAP annotation results showed that other non-cellulolytic strains also possessed CAZymes, in some cases more than strain USB13T, further research is necessary to understand the mechanics of how CAZymes and other cellulases interact to degrade cellulose, and how these genes are expressed under certain conditions. Furthermore, the cellulolytic activity of strain USB13T can be further optimized for commercial use by adjusting growth conditions such as pH, temperature, and growth media.While cellulolytic bacteria are known to inhabit animal intestinal tracts, the rumen, and soil, they can be found almost everywhere, such as ocean floors, municipal landfills, and even extreme environments such as hot springs59. In these habitats, cellulolytic bacteria utilize cellulose while cohabiting with non-cellulolytic bacteria. There have been many studies suggesting the synergistic role non-cellulolytic bacteria play in cellulose degradation, where non-cellulolytic bacteria aid cellulose degradation by neutralizing pH or removing harmful metabolites60,61,62.Bacterial cellulases have shown immense value in various industries such as animal feed processing, food and brewery production, and agriculture, not to mention biofuel synthesis through biomass utilization11. Due to the versatile uses of bacterial cellulases, the cellulolytic strain USB13T has the potential to become an invaluable resource. However, further research of the novel strain’s cellulose-degradation mechanisms is necessary to develop and commercially make use of its bacterial cellulases in the future. In addition, research regarding co-culturing non-cellulolytic bacteria and strain USB13T may also help in developing effective methods to use an otherwise underutilized bioresource.Taxonomy of novel Ramlibacter speciesWhile phylogenetic analyses indicated that the novel strains USB13T, AW1T, GTP1T, and HM2T should be assigned to the genus Ramlibacter, differences in fatty acid compositions, polar lipid profiles, and physiological characteristics suggested that the four novel strains are noticeably distinct from other validly published species of the genus. Additionally, genomic characteristics such as ANI, dDDH, and AAI values further supported the novel strains’ position as a distinct species within the genus Ramlibacter. Therefore, we propose that the strains USB13T, AW1T, GTP1T, and HM2T represent novel species within the genus Ramlibacter.Description of the novel Ramlibacter speciesThe descriptions of the novel species are given according to the standards of the Judicial Commission of the International Committee on Systematic Bacteriology63.Description of Ramlibacter cellulosilyticus sp. nov
Ramlibacter cellulosilyticus (cel.lu.lo.si.ly’ti.cus. N.L. n. cellulosum, cellulose; N.L. adj. lyticus from Gr. lytikos, dissolving; N.L. masc. adj. cellulosilyticus, cellulose-dissolving).Cells of strain USB13T are Gram-negative, rod-shaped, non-flagellated and motile by gliding. The strain is positive for both oxidase and catalase activity, while cells have a width of 0.3–0.5 μm and length of 2.0–2.4 μm. When observed on R2A agar, colonies are reddish white, flat with entire margins, and have a diameter of 1–2 mm. Growth of strain USB13T is observed at 7–50 °C (optimum, 28–30 °C), at pH 5.0–10.0 (optimum, pH 6.0), and at NaCl concentrations of 0–7% (optimum, 0–3%). The strain is unable to grow in anaerobic conditions. Produces siderophores and hydrolyzes Tween 20, Tween 80, CMC, and esculin. According to the API ZYM results, the strain showed positive results for alkaline phosphatase, esterase lipase (C8), leucine arylamidase, acid phosphatase, β-galactosidase, α-glucosidase, and β-glucosidase. In the API 20NE assay, strain USB13T showed positive results only for β-galactosidase. The predominant respiratory quinone is ubiquinone 8 (Q-8). The major fatty acids are C16:0, C10:0 3-OH, and summed feature 3 (consisting of C16:1 ω7c and/or C16:1 ω6c). The polar lipid profile consists of diphosphatidylglycerol (DPG), phosphatidylglycerol (PG), phosphatidylethanolamine (PE), one unidentified phosphoaminolipid, two unidentified phosphoglycoaminolipids, and six unidentified polar lipids. The G + C content is 69.7%. The GenBank/EMBL/DDBJ accession numbers for the 16S rRNA gene sequence and the assembled genome sequence of strain USB13T are MN603953 and JACORT000000000, respectively.The type strain USB13T (= KACC 21656T = NBRC 114839T) was isolated from shallow coastal water at Haeundae Beach, Busan, Republic of Korea.Description of Ramlibacter aurantiacus sp. nov
Ramlibacter aurantiacus (au.ran.ti’a.cus. L. masc. adj. aurantiacus, orange-colored, referring to the orange colonies of the strain).Cells of strain AW1T are Gram-negative, coccoid to short rod-shaped, non-flagellated, and motile by gliding. The strain is negative for oxidase activity, and positive for catalase activity. When observed on R2A agar, colonies are orange, convex, with entire margins, and 0.5–1.0 mm in diameter. Under TEM cells have and approximate width of 0.3–0.5 μm and length of 0.6–0.8 μm. Growth of strain AW1T can be observed at 7–45 °C (optimum, 30 °C), at pH 7.0–10.0 (optimum, 7.0–8.0), and at NaCl concentrations of 0–3% (optimum, 0–1%). The strain does not grow under anaerobic conditions but is able to hydrolyze Tween 80. In addition, AW1T is not able to produce siderophores. In the API ZYM assay, positive for alkaline phosphatase, esterase (C4), esterase lipase (C8), leucine arylamidase, and β-glucosidase. In the API 20NE assay, positive for esculin hydrolysis. The predominant respiratory quinone is ubiquinone 8 (Q-8). The major fatty acids are C16:0, C17:0 cyclo, summed feature 3 (consisting of C16:1 ω7c and/or C16:1 ω6c), and summed feature 8 (consisting of C18:1 ω7c and/or C18:1 ω6c). The polar lipid profile consists of diphosphatidylglycerol (DPG), phosphatidylglycerol (PG), phosphatidylethanolamine (PE), one unidentified phosphoglycolipid, one unidentified lipid, and three unidentified glycolipids. The G + C content is 68.6%. The GenBank/EMBL/DDBJ accession numbers for the 16S rRNA gene sequence and the assembled genome sequence of strain AW1T are MN498045 and JAEQNA000000000, respectively.The type strain AW1T (= KACC 21544T = NBRC 114862T) was isolated from soil at Aewol, Jeju Island, Republic of Korea.Description of Ramlibacter albus sp. nov
Ramlibacter albus (al’bus. L. masc. adj. albus, white, referring to the white colonies of the strain).Strain GTP1T is non-motile, Gram-negative, strictly aerobic, positive for oxidase activity, and negative for catalase activity. When observed on R2A, colonies are white, convex, with entire margins, and 1–2 mm in diameter. Under TEM, cells lack flagella, are rod-shaped, and have a width of 0.7–0.8 μm and length of 1.6–1.9 μm. Growth of strain GTP1T can be observed at 10–45 °C (optimum, 30 °C), at pH 5.0–8.0 (optimum, pH 7.0), and at NaCl concentrations of 0–2% (optimum, 0%). The strain shows positive results for Tween 20 and Tween 80 hydrolysis. GTP1T does not produce siderophores when tested on CAS-blue agar. According to API ZYM results, strain GTP1T is positive for alkaline phosphatase, esterase (C4), esterase lipase (C8), and leucine arylamidase, while the API 20NE assay results show negative results for all substrates. The predominant respiratory quinone is ubiquinone 8 (Q-8). The major fatty acids are C16:0 and summed feature 3 (consisting of C16:1 ω7c and/or C16:1 ω6c). The polar lipid profile consists of diphosphatidylglycerol (DPG), phosphatidylglycerol (PG), phosphatidylethanolamine (PE), and two unidentified phosphoaminolipids. The predominant respiratory quinone is ubiquinone 8 (Q-8). The major fatty acids are C16:0, C17:0 cyclo, summed feature 3 (consisting of C16:1 ω7c and/or C16:1 ω6c), and summed feature 8 (consisting of C18:1 ω7c and/or C18:1 ω6c). The polar lipid profile consists of diphosphatidylglycerol (DPG), phosphatidylglycerol (PG), phosphatidylethanolamine (PE), one unidentified phosphoaminolipid, one unidentified phosphoglycolipid, one unidentified phosphoglycoaminolipid, and two unidentified polar lipids. The G + C content is 67.9%. The GenBank/EMBL/DDBJ accession numbers for the 16S rRNA gene sequence and the assembled genome sequence of strain GTP1T are MN498046 and JACORU000000000, respectively.The type strain GTP1T (= KACC 21702T = NBRC 114488T) was isolated from soil at Seogwipo, Jeju Island, Republic of Korea.Description of Ramlibacter pallidus sp. nov
Ramlibacter pallidus (pal’li.dus. L. masc. adj. pallidus, pale, referring to the color of the colonies).Cells of strain HM2T are Gram-negative, and positive for both oxidase and catalase activities. When observed on R2A agar, colonies are cream-colored, transparent, 1.0–2.5 mm in diameter, and flat with entire margins. Under TEM, monotrichous flagella are observed, and cells are rod-shaped with a width of 0.4–0.78 μm and length of 1.7–1.8 μm. The strain shows the fastest growth at a temperature range of 25–35 °C and at pH values between 8.0 and 9.0. When NaCl is present, growth is observed at concentrations of 0–3% (w/v), with optimal growth was observed at concentrations of 0–1% (w/v). The strain is not able to tolerate anaerobic conditions. Strain HM2T hydrolyzes Tween 80 and weakly hydrolyzes casein. However, siderophore production cannot be observed when tested on CAS-blue agar. According to API ZYM tests, strain HM2T shows positive results for alkaline phosphatase, esterase (C4), esterase lipase (C8), leucine arylamidase, valine arylamidase, acid phosphatase, and naphthol-AS-BI-phosphohydrolase. In addition, API 20NE tests show positive results for nitrate (NO3) to nitrite (NO2-) reduction and esculin hydrolysis. The G + C content is 69.9%. The GenBank/EMBL/DDBJ accession numbers for the 16S rRNA gene sequence and the assembled genome sequence of strain HM2T are MN498047 and JADDIV000000000, respectively.The type strain HM2T (= KCTC 82557T = NBRC 114489T) was isolated from soil at Seopjikoji, Jeju Island, Republic of Korea. More

Indonesia’s bleeding toad (Leptophryne cruentata) is critically endangered.Credit: Pepew Fegley/Shutterstock

One-quarter of all plant and animal species are threatened with extinction owing to factors such as climate change and pollution. Starting this week, negotiators and ministers from more than 190 countries are meeting at a United Nations biodiversity summit called COP15 in Montreal, Canada, to address the emergency.
10 startling images of nature in crisis — and the struggle to save it
From 7 to 19 December, they will be trying to seal a new deal to save Earth’s biodiversity. The treaty, known as the post-2020 Global Biodiversity Framework, is intended to establish precise targets for countries to protect and restore nature, including conserving 30% of the planet by 2030 and cutting nutrient pollution, such as reducing nitrogen fertilizer loss from farmland.Time is running out. “We’re driving species to extinction at a rate about 1,000 times faster than they are created through evolution,” says Stuart Pimm, an ecologist at Duke University in Durham, North Carolina, and head of Saving Nature, a non-profit conservation organization.As COP15 kicks off, researchers and policy experts are concerned that countries still disagree on too many issues to secure a deal that will protect species and ecosystems effectively. Here, Nature looks at the extent of the crisis, and what scientists say countries must do to succeed.Which species are most at risk, and what’s threatening them?Among the most at-risk groups are amphibians and reef-forming corals. A global assessment shows that more than 40% of amphibians are threatened with extinction1, including the critically endangered bleeding toad (Leptophryne cruentata), which lives in Mount Gede Pangrango National Park in Java, Indonesia.These toads were thought to be extinct until the year 2000, when some were spotted by a team led by Mirza Kusrini, a herpetologist at Bogor Agricultural University in Indonesia. But the researchers found that the amphibians were infected with chytrid (Chytridiomycota sp.), a fungus that has devastated global amphibian populations. Kusrini says that climate change is probably making life hard for the tiny toad, which got its common name from the crimson, splatter-like spots covering its body. Warm weather can stimulate fungal outbreaks and shift the timing of behaviours, such as the toads’ breeding season, making the amphibians vulnerable.

Source: Red List Index/IUCN

Global warming, which has been raising sea temperatures, is also responsible for harming coral reefs around the globe (see ‘Threat assessment’). Over a period of 9 years, up to 2018, 14% of the world’s coral died out — a massive problem, because today, coral reefs support one-quarter of all marine species.Research shows that climate change is quickly becoming a large threat to biodiversity2. But still, the most-destructive forces are the conversion of land and seas for agricultural uses and people exploiting natural resources through fishing, logging, hunting and the wildlife trade. About 75% of land and 66% of ocean areas have been significantly altered, usually for producing food.What might happen if species disappear?It’s difficult to predict, because doing so requires knowledge of which species are present in a particular ecosystem, such as a rainforest, and what functions they have, says Shahid Naeem, an ecologist at Columbia University in New York City. Much of that information is often unknown. However, scientists have shown3 that ecosystems with less biodiversity are not as good at capturing and converting resources into biomass, such as happens when plants capture nutrients or sunlight used for growth.
Why deforestation and extinctions make pandemics more likely
Neither are less-diverse ecosystems as good at decomposing and recycling biological materials and nutrients. For example, studies show that dead organisms are broken down, and their nutrients recycled, more quickly when a high variety of plant litter covers the forest floor4. Ecosystems with low biodiversity also have low resilience — they are not as able to bounce back after a perturbation or shock, such as a fire, as more-diverse systems are, Naeem says.“If we lose parts of our system, it simply won’t function very efficiently, and it won’t be very robust,” he adds. “The science behind that is rock solid.”Ecosystems also provide clean water and can sometimes prevent diseases from spreading to humans. When species are lost, these services deteriorate, Kusrini says. For example, most amphibians eat insects, many of which are considered pests, such as cockroaches, termites and mosquitoes. Studies have shown a rise in cases of malaria — spread by mosquitoes — in areas in Central America where amphibian populations have collapsed5. “You know when they disappear”, Kusrini says, because insect numbers rise and people start using more pesticides to kill them.What solutions do researchers say are needed to protect biodiversity?Protecting and conserving habitats is central to saving species. This idea is captured in the framework being negotiated at COP15. The draft includes the goal of conserving at least 30% of the world’s land and sea by 2030. But for protections to be most effective, they must include regions that are rich in biodiversity, such as tropical forests, Pimm says. Despite an increase in protected areas worldwide over the past ten years, species numbers have still declined, because these safeguards were not in the right places, studies show6.

Delegates at COP15 in Montreal show their support for a new agreement among nations to protect Earth’s biodiversity.Credit: UN Convention on Biological Diversity (CC BY 2.0)

“What we’re going to be looking for at COP15 is more quality, not just more quantity,” Pimm says.Eradicating invasive species is another important conservation strategy, and the framework’s draft currently calls for cutting the introduction of such species in half. Some estimates suggest that invasive predators, such as cats and rats, are responsible for more than half of all extinctions of birds, mammals and reptiles7.It’s important that nations agree on a framework with at least some quantifiable targets, so that progress can be measured, and so that countries can be held accountable if they fail to meet their targets, researchers say. “I’m afraid what will happen is, they will produce a long list of ‘waffle’,” Pimm says. “We need quantification.”Will nations manage to agree on a new deal to protect nature?As COP15 begins, the outlook is not good. The text of the draft is still littered with unresolved issues. At a press conference on 6 December, Elizabeth Mrema, executive secretary of the Convention on Biological Diversity — the global treaty that underpins the new biodiversity deal — said that national negotiators had made insufficient progress in a final round of discussions before the start of the summit. She urged countries to compromise, otherwise they will fail to reach a deal. “The state of the planet is in crisis,” Mrema said. “This is our last chance to act.”
Troubled biodiversity plan gets billion-dollar funding boost
One key contentious issue is how to finance biodiversity conservation, particularly in low- and middle-income countries, which are home to much of the world’s biodiversity. These nations, including Brazil and Gabon, would like a new fund to be established with US$100 billion added per year in aid. So far, that proposal has not gained traction with wealthier countries. “They really need to have the financial commitments, because things don’t get done without the money,” Naeem says.Despite the pessimism, Naeem is certain that scientists and advocates will keep pushing for a deal. “There would be real change” if countries were able to achieve a universal decrease in biodiversity loss, he says. More

Field parameters and major ionsThe results of the hydrochemistry and environmental isotopes for the 40 lakes are presented spatially in Figs. S1–S11 and are located in Tables S1 and S2.The lake waters are oxic (8.6–12.6 mg l−1) and range from slightly acidic (pH 6.0) to slightly alkaline (pH 9.2). Lake water temperatures are generally highest for lakes along the west coast (greater than 10 °C, Table S2). Phosphate concentrations are below detection level (0.1 mg l−1) for all lakes and nitrate was low ranging from below detection limit ( More

Aagaard-Tillery KM, Grove K, Bishop J, Ke X, Fu Q, McKnight R et al. (2008) Developmental origins of disease and determinants of chromatin structure: maternal diet modifies the primate fetal epigenome. J Mol Endocrinol 41:91–102Article 
CAS 

Google Scholar 
Alekseev V, Lampert W (2001) Maternal control of resting – egg production in Daphnia. Nature 414:899–901Article 
CAS 

Google Scholar 
Anderson OS, Sant KE, Dolinoy DC (2012) Nutrition and epigenetics: an interplay of dietary methyl donors, one-carbon metabolism and DNA methylation. J Nutritional Biochem 23:853–859Article 
CAS 

Google Scholar 
Agrawal AA, Laforsch C, Tollrian R (1999) Transgenerational induction of defenses in animals and plants. Nature 401:60–63Article 
CAS 

Google Scholar 
Asselman J, De Coninck DIM, Vandegehuchte MB, Jansen M, Decaestecker E, De Meester L, Busshe JV, Vanhaecke L, Janssen CR, De Schamphelaere KAC (2015) Global cytosine methylation in Daphnia magna depends on genotype, environment, and their interaction. Environ Toxicol 34:5
Google Scholar 
Bernardo J (1996) Maternal effects in animal ecology. Am Zool 36:83–105Article 

Google Scholar 
Berger SL (2007) The complex language of chromatin regulation during transcription. Nature 447:407–412Article 
CAS 

Google Scholar 
Bewick AJ, Vogel KJ, Moore AJ, Schmitz RJ (2017) Evolution of DNA methylation across insects. Mol Biol Evol 34:654–665CAS 

Google Scholar 
Bird A (2007) Perceptions of epigenetics. Nature 447:396–398Article 
CAS 

Google Scholar 
Boersma M (1995) The allocation of resources to reproduction in Daphnia galeata: against the odds? Ecology 76(4):121–1261Article 

Google Scholar 
Boersma M (1997) Offspring size in Daphnia: does it pay to be overweight? Hydrobiologia 360:79–88Article 

Google Scholar 
Boycott AE, Diver C (1923) On the inheritance of the sinistrality in Limnea peregra. Proc R Soc Lond B 95:207–213Article 

Google Scholar 
Brett MT (1993) Resource quality effects on Daphnia longispina offspring fitness. J Plankton Res 15(4):403–412Article 

Google Scholar 
Burns CW (1995) Effects of crowding and different food levels on growth and reproductive investment of Daphnia. Oecologia 101:234–244Article 

Google Scholar 
Cameron NM, Shahrokh D, Del Corpo A, Dhir SK, Szyf M, Champagne FA, Meaney MJ (2008) Epigenetic programming of phenotypic variations in reproductive strategies in the rat through maternal care. J Neuroendocrinol 20:795–801Article 
CAS 

Google Scholar 
Champagne FA (2012) Epigenetics and developmental plasticity across species. Dev Psychobiol 55:33–41Article 

Google Scholar 
Chan SY, Vasilopoulou E, Kilby MD (2009) The role of the placenta in thyroid hormone delivery to the fetus. Nat Clin Pract Endocrinol Metab 5:45–54Article 
CAS 

Google Scholar 
Chong S, Whitelaw E (2004) Epigenetic germline inheritance. Curr Opin Genet Dev 14(6):692–696Article 
CAS 

Google Scholar 
Clark J, Garbutt JS, McNally L, Little TJ (2017) Disease spread in age structured populations with maternal age effects. Ecol Lett 20:445–451Article 

Google Scholar 
Colbourne JK, Herbert PDN, Taylor DJ (1997) Evolutionary origins of phenotypic diversity. In: Givnish TJ, Systma KJ (eds) Daphnia in molecular evolution and adaptive radiation. Cambridge University Press. p 163–188Colbourne JK, Pfrender ME, Gilbert D, Thomas WK, Tucker A, Oakley TH et al. (2011) The ecoresponsive genome of Daphnia pulex. Science 331(6017):555–561Article 
CAS 

Google Scholar 
Desmarais KH (1997) Keeping Daphnia out of the surface film with cetyl alcohol. J Plankton Res 19(1):149–154Article 

Google Scholar 
Dorts J, Falisse E, Schoofs E, Flamion E, Kestermont P, Silvestre F (2016) DNA methyltransferases and stress-related genes expression in zebrafish larvae after exposure to heat and copper during reprogramming of DNA methylation. Sci Rep 6:34254Article 
CAS 

Google Scholar 
Ducker GS, Rabinowitz JD (2016) One-carbon metabolism in health and disease. Cell Metab 25:27–42. https://doi.org/10.1016/j.cmet.2016.08Dudycha JL, Brandon CS, Deitz KC (2012) Population genomics of resource exploitation: insights from gene expression profiles of two Daphnia ecotypes fed alternate resources. Ecol Evol 2:329–340Dzialowski EM, Reed WL, Sotherland PR (2009) Effects of egg size on double-crested cormorant (Phalacrocorax auritus) egg composition and hatchling phenotype. Comp Biochem Physiol A Mol Integr Physiol 152:262–267Article 

Google Scholar 
Frost PC, Ebert D, Larson JH, Marcus MA, Wagner ND, Zalewski A (2010) Transgenerational effects of poor elemental food quality on Daphnia magna. Oecologia 162(4):865–872Article 

Google Scholar 
Gabsi F, Glazier DS, Hammers-Wirtz M, Ratte HT, Preuss TG (2014) How to interactive maternal traits and environmental factors determine offspring size in Daphnia magna?. Ann Limnol 50:9–18Article 

Google Scholar 
Garbutt JS, Little TJ (2016) Bigger is better: changes in body size explain a maternal effect of food on offspring disease resistance. Ecol Evolution 7:1403–1409Article 

Google Scholar 
Gibney ER, Nolan CM (2010) Epigenetics and gene expression. Heredity 105:4–13Article 
CAS 

Google Scholar 
Gillis MK, Walsh MR (2019) Individual variation in plasticity dulls transgenerational responses to stress. Funct Ecol 33:1993–2002Glazier DS (1992) Effects of food, genotype, and maternal size and age on offspring investment in Daphnia magna. Ecology 73(3):910–926Article 

Google Scholar 
Gliwicz ZM, Guisande C (1992) Family planning in Daphnia: resistance to starvation in offspring born to mothers grown at different food levels. Oceologia 91:463–467Article 

Google Scholar 
Goos JM, Swain CJ, Munch SB, Walsh MR (2018) Maternal diet and age alter direct and indirect relationships between lifer-history traits across multiple generations. Funct Ecol 33:491–502Article 

Google Scholar 
Goulden CE, Horning LL (1980) Population oscillations and energy reserves in planktonic cladocera and their consequences to competition. Proc Natl Acad Sci USA 77:1716–1720Article 
CAS 

Google Scholar 
Groothuis TG, Schwabl H (2008) Hormone-mediated maternal effects in birds: mechanisms matter but what do we know of them? Philos Trans R Soc Lond B Biol Sci 363:1647–1661Article 
CAS 

Google Scholar 
Guisande C, Gliwicz ZM (1992) Egg size and clutch size in two Daphnia species at different food levels. J Plankton Res 14(7):997–1007Article 

Google Scholar 
Hearn J, Chow FW-N, Barton H, Tung M, Wilson P, Blaxter M et al. (2018) Daphnia magna microRNAs respond to nutritional stress and ageing but are not transgenerational. Mol Ecol 27:1402–1412Article 
CAS 

Google Scholar 
Hearn J, Pearson M, Blaxter M, Wilson PJ, Little TJ (2019) Genome-wide methylation is modified by caloric restriction in Daphnia magna. BCM Genetics 20:197Hearn J, Plenderleith F, Little TJ (2021) DNA methylation differs extensively between strains of the same geographical origin and changes with age in Daphnia magna. Epigenetics Chromatin 14:4. https://doi.org/10.1186/s13072-020-00379-zHead JA (2014) Patterns of DNA methylation in animals: an ecotoxicological perspective. Integr Comp Biol 54:77–86Article 
CAS 

Google Scholar 
Hebert PDN (1981) Obligate asexuality in Daphnia. Am Nate 117:784–789Article 

Google Scholar 
Herman JJ, Sultan SE (2016) DNA methylation mediates genetic variation for adaptive transgenerational plasticity. Proc Biol Sci 283(1838):20160988. https://doi.org/10.1098/rspb.2016.0988Article 
CAS 

Google Scholar 
Hiruta C, Nishida C, Tochinai S (2010) Abortive meiosis in the oogenesis of parthenogenetic Daphnia pulex. Chromosome Res 18:833–840Article 
CAS 

Google Scholar 
Ho DH, Burggren WW (2010) Epigenetics and transgenerational transfer: a physiological perspective. J Exp Biol 213:3–16Article 
CAS 

Google Scholar 
Ho DH (2008) Morphological and physiological developmental consequences of parental effects in the chicken embryo (Gallus gallus domesticus) and the zebrafish larva (Danio rerio). Diss: University of North TexasInnes DJ, Fox CJ, Winsor GL (2000) Avoiding the cost of males in obligately asexual Daphnia pulex (Leydig). Proc: Biol Sci 267(1447):991–997CAS 

Google Scholar 
Jeremias G, Barbosa J, Marques SM, De Schamphelaere KAC, Van Nieuwerburgh F, Deforce D, Gonçalves FJM, Pereira JL, Asselman J (2018) Transgenerational inheritance of dna hypomethylation in Daphnia magna in response to salinity stress. Environ Sci Technol 52(17):10114–10123Article 
CAS 

Google Scholar 
Jian X, Yang W, Zhao S, Liang H, Zhao Y, Chen L et al. (2013) Maternal effects of inducible tolerance against the toxic cyanobacterium Microcystis aeruginosa in the grazer Daphnia carinata. Environ Pollut 178:142–146Article 

Google Scholar 
Keating KI (1985) The influence of vitamin-B12 deficiency on the reproduction of Daphnia-Pulex Leydig (Cladocera). J Crustacean Biol 5:30–136Article 

Google Scholar 
Kleiven OT, Larsson P, Hobaek A (1992) Sexual reproduction in Daphnia magna requires three stimulie. Oikos 65:197–206Article 

Google Scholar 
Kusari F, O’Doherty AM, Hodges NJ, Wojewodzic MW (2017) Bi-directional effects of vitamin B12 and methotrexate on Daphnia magna fitness and genomic methylation. Sci Rep 7:11872Article 

Google Scholar 
Kvist J, Athanasio CG, Solari OS, Brown JB, Colbourne JK, Pfrender ME, Mirbahai L (2018) Pattern of DNA methylation in Daphnia: evolutionary perspective. Genome Biol Evolution 10(8):1988–2007Article 
CAS 

Google Scholar 
Lamka GF, Harder AM, Sundaram M, Schwartz TS, Christie MR, DeWoody JA, Willoughby JR (2022) Epigenetics in ecology, evolution, and conservation. Front Ecol Evol 10. https://doi.org/10.3389/fevo.2022.871791LaMontagne JM, McCauley E (2001) Maternal effects in Daphnia: what mothers are telling their offspring and do they listen. Ecol Lett 4:64–71Article 

Google Scholar 
Li Q, Jiang X (2014) Offspring tolerance to toxic Microcystis aeruginosa in Daphnia pulex shaped by maternal food availability and age. Fundam Appl Limnol 185:315–319Article 

Google Scholar 
Mastorci F, Vicentini M, Viltart O, Manghi M, Graiani G, Quaini Fet al. (2009) Long-term effects of prenatal stress: changes in adult cardiovascular regulation and sensitivity to stress. Neurosci Biobehav Rev 33:191–203Article 

Google Scholar 
Mkee D, Ebert D (1996) The interactive effects of temperature, food level and maternal phenotype on offspring size in Daphnia magna. Oecologia 107(2):189–196Article 

Google Scholar 
Mousseau TA, Fox CW (1998) The adaptive significance of maternal effects. Trends Ecol Evol 13:403–407Article 
CAS 

Google Scholar 
Nguyen ND, Matsuura T, Kato Y, Watanabe H (2020) Caloric restriction upregulates the expression of DNMT3.1, lacking the conserved catalytic domain, in Daphnia magna. Genesis 58:12Article 

Google Scholar 
Nguyen ND, Matsuura T, Kato Y, Watanabe H (2021) DNMT3.1 controls trade-offs between growth, reproduction, and life span under starved conditions in Daphnia magna. Sci Rep 11:7326Article 
CAS 

Google Scholar 
Nusslein-Volhard C, Frohnhofer HG, Lehmann R (1987) Determination of anteroposterior polarity in Drosophila. Science 238:1675–1681Article 
CAS 

Google Scholar 
Pieters BJ, Liess M (2006) Maternal nutritional state determines the sensitivity of Daphnia magna offspring to short-term fenvalerate exposure. Aquat Toxicol 76:286–277Article 

Google Scholar 
R Core Team (2021) R: a language and environmental for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. https://www.R-project.org/Richards EJ (2006) Inherited epigenetic variation – revisiting soft inheritance. Nat Rev Genet 7:395–401Article 
CAS 

Google Scholar 
Stollewerk A (2010) The water flea Daphnia – a new model system for ecology and evolution? J Biol 9(2):21Article 

Google Scholar 
Sturtevant AH (1923) Inheritance of direction of coiling in Limnea. Science 58:269Article 
CAS 

Google Scholar 
Suter MA, Chen A, Burdine MS, Choudhury M, Harris RA, Lane RH et al. (2012) A maternal high-fat diet modulates fetal SIRT1 histone and protein deacetylase activity in nonhuman primates. FASEB J 26:5106–5114Article 
CAS 

Google Scholar 
Tessier AJ, Consolatti NL (1989) Variation in offspring size in Daphnia and consequences for individual fitness. Oikos 56:269–276Article 

Google Scholar 
Tessier AJ, Consolatti NL (1991) Resource quantity and offspring quality in Daphnia. Ecology 72(2):468–478Article 

Google Scholar 
Trerotola M, Relli V, Simeone P, Alberti S (2015) Epigenetic inheritance and the missing heritability. Hum Genomics 9(1):17. https://doi.org/10.1186/s40246-015-0041-3Article 
CAS 

Google Scholar 
Trijau M, Asselman J, Armant O, Adam-Guillermin C, De Schamphelaere KAC, Alonzo F (2018) Transgenerational DNA methylation changes in Daphnia magna exposed to chronic γ irradiation. Environ Sci Technol 52(7):4331–4339Article 
CAS 

Google Scholar 
Urabe J, Sterner RW (2001) Contrasting effects of different types of resource depletion on life-history traits in Daphnia. Funct Ecol 15:165–174Article 

Google Scholar 
Vandegehuchte MB, Kyndt T, Vanholme B, Haegeman A, Gheysen G, Janssen CR (2009a) Occurrence of DNA methylation in Daphnia magna and influence of multigeneration Cd exposure. Environ Int 35(4):700–706Article 
CAS 

Google Scholar 
Vandegehuchte MB, Lemiere F, Janssen CR (2009b) Quantitative DNA-methylation in Daphnia magna and effects of multigeneration Zn exposure. Comp Biochem Physiol C Toxicol Pharmacol 150:343–348Article 
CAS 

Google Scholar 
Vandegehuchte MB, Janssen CR (2011) Epigenetics and its implications for ecotoxicology. Ecotoxicology 20:607–624Article 
CAS 

Google Scholar 
Vandegehuchte MB, Janssen CR (2014) Epigenetics in an ecotoxicological context. Mutat Res Genet Toxicol Environ Mutagen 764–765:36–45Article 

Google Scholar 
Walsh MR, La Pierre KJ, Post DM (2014) Phytoplankton composition modifies predator-driven life history evolution in Daphnia. Evol Ecol 28:397–411Article 

Google Scholar 
Walsh MR, Cooley F, Biles K, Munch SB (2015) Predator-induced phenotypic plasticity within- and across generations: a challenge for theory? Proc R Soc B Biol Sci 282:20142205Article 

Google Scholar 
Wickham H (2016) ggplot2: elegant graphics for data analysis. Springer-Verlag New York. ISBN 978-3-319-24277-4. https://ggplot2.tidyverse.orgWolf JB, Wade MJ (2009) What are maternal effects (and what are they not)? Philos Trans R Soc Lond B Biol Sci 364(1520):1107–1115Article 

Google Scholar 
Zaffagnini F (1987) Reproduction in Daphnia. Mem Ist Ital Idrobiol 45:245–284
Google Scholar 
Zhao S, Fernald RD (2005) Comprehensive algorithm for quantitative real-time polymerase chain reaction. J Comput Biol 12(8):1045–1062Article 

Google Scholar  More

The isotopic spread for each study siteThe overall linear relationship between δ2H and δ18O values for hair (n = 81) and toenails (n = 39), respectively, were (Fig. 2):$$delta^{2} {text{H}}_{{text{hair(VSMOW)}}} = , 0.89 times delta^{18} {text{O}}_{{text{hair(VSMOW)}}} {-} , 86.16,;{text{R}}^{2} = , 0.19,;p , < , 0.01$$ (1) $$delta^{2} {text{H}}_{{text{toenail(VSMOW)}}} = , 0.15 times delta^{18} {text{O}}_{{text{toenail(VSMOW)}}} {-} , 91.69,;{text{R}}^{2} = , 0.00,;p , = , 0.69$$ (2) Figure 2δ2H and δ18O values (‰) of all samples for both hair (δ2H: n = 81, δ18O: n = 82) and toenails (δ2H and δ18O: n = 39). The solid black line represents the Global Meteoric Water Line (GMWL) [δ2H = 8 (times) δ18O + 10] and is included in the graph for comparison purposes. The regression lines between oxygen and hydrogen values for hair [δ2Hhair(VSMOW) = 0.89 × δ18Ohair(VSMOW) − 86.16, R2 = 0.19, p  − 82‰ were then split further where any samples with δ2Hhair values less than − 73‰ were initially classified as Site 2. These samples were then split again to either Site 2 (δ2Hhair ≥ 76‰) or Site 4 (δ2Hhair  − 73‰ were classified as Site 4. No samples could be classified as originating from Site 3. The second CART model was built for stable hydrogen and oxygen isotopes of toenails (Model 2) (Fig. 5b). The model included only two decision nodes in which the first predictor variable was δ2Htoenail value, where samples with values less than − 93‰ were predicted to be from Site 1. For toenail samples with hydrogen values greater than − 93‰, oxygen values were used to determine whether they could be classified as Site 2 or Site 4. Those samples with δ18Otoenail values less than 9.6‰ were classified as Site 2 and those with values greater than 9.6‰ were predicted as Site 4. No samples were predicted to be from Site 3 purely from stable hydrogen and oxygen isotopes in toenails. Finally, the third model consisted of stable hydrogen and oxygen isotope values in both hair and toenail samples (Model 3) (Fig. 5c). Model 3 selected toenails as the best attribute for classification, which indicates that toenail isotope values are the better predictor when both hair and toenail samples are present for analysis from Sites 1–4. The model was similar to that of Model 2.Figure 5Decision trees developed from both δ2H and δ18O values of (a) hair [Model 1, trained with n = 65], (b) toenails [Model 2, trained with n = 32] and (c) of both hair and toenails [Model 3, trained with n = 28]. The predicted study site numbers are shown on the first row within each bubble. The proportions of samples in each node are shown as decimals for Sites 1, 2, 3, 4, respectively. The percentages indicate the proportion of samples within each sub-partition.Full size imageConfusion matrices (Table 1) were constructed for all three models to evaluate the performance of the classification models. Of the three models, Model 3 proved to be the most accurate model with an overall accuracy of 71.4% (see Supplementary Fig S2. online). The performance evaluation summary, including measures for sensitivity, specificity, positive predictive value, and negative predictive value for all three models, is provided in (see Supplementary Table S3. online).Intra-individual differencesBoth hair and toenail samples were retrieved from 35 of the 86 individuals. The paired difference between δ2H values in hair and toenails of the same individual was tested using the Wilcoxon Signed Rank's test for non-normal data as the dataset failed the Shapiro–Wilk's normality test at the α = 0.05 significance level. Significant differences were found between δ2H values of hair (n = 35, mean = − 78.0‰, s.d. = 3.06) and toenails (n = 35, mean = − 90.9‰, s.d. = 3.27) from the same individual (p  0.05. Overall, the isotopic values of δ2H in hair were higher than those of toenail from the same individual by 13.0‰, on average, with a standard deviation of 8.4‰. For δ18O, the average was 1.5‰ with a standard deviation of 4.6‰ (Fig. 6).Figure 6(a) δ2H and (b) δ18O values in hair and toenails for all individuals that provided both tissue types (n = 35). Study site information are also shown by shapes. The standard deviations of each sample, ran in either duplicates or triplicates, are shown by error bars. Note that error bars cannot be seen for some samples due to small standard deviations. The average difference between the isotopic values of hair and toenail from the same individual were 13.0‰ with a standard deviation of 8.4‰ for δ2H and 1.5‰ with a standard deviation of 4.6‰ for δ18O.Full size imageThe linear relationships between δ2H in hair and toenails for all individuals were (see Supplementary Fig S3. online):$$delta^{2} {text{H}}_{{{text{hair}}}} = , 0.48 times delta^{2} {text{H}}_{{text{toenail }}} {-} , 34.72,;{text{R}}^{2} = , 0.16,;p , < , 0.05$$ (19) and for δ18O:$$delta^{18} {text{O}}_{{{text{hair}}}} = , 0.55 times delta^{18} {text{O}}_{{{text{toenail}}}} + , 5.16,;{text{R}}^{2} = , 0.13,;p , < 0.05$$ (20) Overall, both equations showed a weak relationship, as seen by the small R2 values. More

Barbosa, P., Krischik, V. A. & Jones, C. G. Microbial mediation of plant-herbivore interactions (John Wiley & Sons, 1991).
Google Scholar 
Berenbaum, M. R. Allelochemicals in insect–microbe–plant interactions; agents provocateurs in the
coevolutionary arms race. In Nov. Asp. Insect-Plant Interact. (eds Barbosa, P. & Letourneau, D. K.) 97–123 (1988).Mason, C. J., Jones, A. G. & Felton, G. W. Co-option of microbial associates by insects and their impact on plant–folivore interactions. Plant Cell Environ. 42, 1078–1086 (2019).Article 
CAS 

Google Scholar 
Sugio, A., Dubreuil, G., Giron, D. & Simon, J.-C. Plant–insect interactions under bacterial influence: Ecological implications and underlying mechanisms. J. Exp. Bot. 66, 467–478 (2015).Article 
CAS 

Google Scholar 
Hansen, A. K. & Moran, N. A. The impact of microbial symbionts on host plant utilization by herbivorous insects. Mol. Ecol. 23, 1473–1496 (2014).Article 

Google Scholar 
Mendes, R., Garbeva, P. & Raaijmakers, J. M. The rhizosphere microbiome: Significance of plant beneficial, plant pathogenic, and human pathogenic microorganisms. FEMS Microbiol. Rev. 37, 634–663 (2013).Article 
CAS 

Google Scholar 
Pineda, A. et al. Helping plants to deal with insects: The role of beneficial soil-borne microbes. Trends Plant Sci. 15, 507–514 (2010).Article 
CAS 

Google Scholar 
Hammer, T. J. & Bowers, M. D. Gut microbes may facilitate insect herbivory of chemically defended plants. Oecologia 179, 1–14 (2015).Article 
ADS 

Google Scholar 
Liu, H. et al. An ecological loop: Host microbiomes across multitrophic interactions. Trends Ecol. Evol. 34, 1118–1130 (2019).Article 

Google Scholar 
Grunseich, J. M., Thompson, M. N., Aguirre, N. M. & Helms, A. M. The role of plant-associated microbes in mediating host-plant selection by insect herbivores. Plants 9, 6 (2020).Article 
CAS 

Google Scholar 
Ferrari, J. et al. Linking the bacterial community in pea aphids with host-plant use and natural enemy resistance. Ecol. Entomol. 29, 60–65 (2004).Article 

Google Scholar 
McLean, A. H. et al. Insect symbionts in food webs. Philos. Trans. R. Soc. B Biol. Sci. 371, 20150325 (2016).Article 

Google Scholar 
Giron, D., Dedeine, F., Dubreuil, G. et al. Influence of microbial symbionts on plant–insect interactions. In: Advances in botanical research. Elsevier, pp 225–257 (2017).Jones, A. G., Mason, C. J., Felton, G. W. & Hoover, K. Host plant and population source drive diversity of microbial gut communities in two polyphagous insects. Sci. Rep. 9, 1–11 (2019).Article 

Google Scholar 
Xu, T.-T., Jiang, L.-Y., Chen, J. & Qiao, G.-X. Host plants influence the symbiont diversity of Eriosomatinae (Hemiptera: Aphididae). Insects 11, 217. https://doi.org/10.3390/insects11040217 (2020).Article 

Google Scholar 
Qin, M. et al. Microbiota associated with Mollitrichosiphum aphids (Hemiptera: Aphididae: Greenideinae): Diversity, host species specificity and phylosymbiosis. Environ. Microbiol. 23(4), 2184–2198. https://doi.org/10.1111/1462-2920.15391 (2021).Article 
CAS 

Google Scholar 
Douglas, A. E. Microbial brokers of insect-plant interactions revisited. J. Chem. Ecol. 39, 952–961 (2013).Article 
CAS 

Google Scholar 
Engel, P. & Moran, N. A. The gut microbiota of insects–diversity in structure and function. FEMS Microbiol. Rev. 37, 699–735 (2013).Article 
CAS 

Google Scholar 
Chung, S. H. et al. Host plant species determines symbiotic bacterial community mediating suppression of plant defenses. Sci. Rep. 7, 1–13 (2017).
Google Scholar 
Holt, J. R. et al. Differences in microbiota between two multilocus lineages of the sugarcane aphid (Melanaphis sacchari) in the continental United States. Ann. Entomol. Soc. Am. 113(4), 257–265 (2020).Article 
CAS 

Google Scholar 
McLean, A. H., Godfray, H. C. J., Ellers, J. & Henry, L. M. Host relatedness influences the composition of aphid microbiomes. Environ. Microbiol. Rep. 11, 808–816 (2019).Article 

Google Scholar 
Jones, R. T., Sanchez, L. G. & Fierer, N. A cross-taxon analysis of insect-associated bacterial diversity. PLoS ONE 8, e61218 (2013).Article 
ADS 
CAS 

Google Scholar 
Najar-Rodríguez, A. J. et al. The microbial flora of Aphis gossypii: Patterns across host plants and geographical space. J. Invertebr. Pathol. 100, 123–126. https://doi.org/10.1016/j.jip.2008.10.005 (2009).Article 

Google Scholar 
Blankenchip, C. L., Michels, D. E., Braker, H. E. & Goffredi, S. K. Diet breadth and exploitation of exotic plants shift the core microbiome of tropical herbivorous beetles. PeerJ. Prepr. 6, e26692v1 (2018).
Google Scholar 
Gauthier, J.-P., Outreman, Y., Mieuzet, L. & Simon, J.-C. Bacterial communities associated with host-adapted populations of pea aphids revealed by deep sequencing of 16S ribosomal DNA. PLoS ONE 10, e0120664 (2015).Article 

Google Scholar 
Wagner, S. M. et al. Facultative endosymbionts mediate dietary breadth in a polyphagous herbivore. Funct. Ecol. 29, 1402–1410 (2015).Article 

Google Scholar 
Guidolin, A. S. & Cônsoli, F. L. Symbiont diversity of Aphis (Toxoptera) citricidus (Hemiptera: Aphididae) as influenced by host plants. Microb. Ecol. 73, 201–210 (2017).Article 

Google Scholar 
Leonardo, T. E. & Muiru, G. T. Facultative symbionts are associated with host plant specialization in pea aphid populations. Proc. R. Soc. Lond. B Biol. Sci. 270, S209–S212 (2003).Article 

Google Scholar 
Xu, S., Jiang, L., Qiao, G. & Chen, J. The bacterial flora associated with the polyphagous aphid Aphis gossypii Glover (Hemiptera: Aphididae) is strongly affected by host plants. Microb. Ecol. 79, 971–984. https://doi.org/10.1007/s00248-019-01435-2 (2020).Article 
CAS 

Google Scholar 
Ferrari, J., West, J. A., Via, S. & Godfray, H. C. J. Population genetic structure and secondary symbionts in host-associated populations of the pea aphid complex. Evolution 66, 375–390. https://doi.org/10.1111/j.1558-5646.2011.01436.x (2012).Article 

Google Scholar 
Brady, C. M. et al. Worldwide populations of the aphid Aphis craccivora are infected with diverse facultative bacterial symbionts. Microb. Ecol. 67, 195–204. https://doi.org/10.1007/s00248-013-0314-0 (2014).Article 

Google Scholar 
Henry, L. M., Maiden, M. C., Ferrari, J. & Godfray, H. C. J. Insect life history and the evolution of bacterial mutualism. Ecol. Lett. 18, 516–525 (2015).Article 

Google Scholar 
Simon, J.-C. et al. Host–based divergence in populations of the pea aphid: Insights from nuclear markers and the prevalence of facultative symbionts. Proc. R. Soc. Lond. B Biol. Sci. 270, 1703–1712. https://doi.org/10.1098/rspb.2003.2430 (2003).Article 

Google Scholar 
Brady, C. M. & White, J. A. Cowpea aphid (Aphis craccivora) associated with different host plants has different facultative endosymbionts. Ecol. Entomol. 38, 433–437. https://doi.org/10.1111/een.12020 (2013).Article 

Google Scholar 
Blackman, R. L. & Eastop, V. F. Aphids on the world’s herbaceous plants and shrubs, 2 Vol. set (John Wiley & Sons, 2008).
Google Scholar 
Züst, T. & Agrawal, A. A. Population growth and sequestration of plant toxins along a gradient of specialization in four aphid species on the common milkweed Asclepias syriaca. Funct. Ecol. 30, 547–556 (2016).Article 

Google Scholar 
Zytynska, S. E. & Weisser, W. W. The natural occurrence of secondary bacterial symbionts in aphids. Ecol. Entomol. 41, 13–26 (2016).Article 

Google Scholar 
Harrison, J. S. & Mondor, E. B. Evidence for an invasive aphid “Superclone”: Extremely low genetic diversity in Oleander aphid (Aphis nerii) populations in the Southern United States. PLoS ONE 6, e17524. https://doi.org/10.1371/journal.pone.0017524 (2011).Article 
ADS 
CAS 

Google Scholar 
Mooney, K., Jones, P. & Agrawal, A. Coexisting congeners: Demography, competition, and interactions with cardenolides for two milkweed-feeding aphids. Oikos 117, 450–458 (2008).Article 
CAS 

Google Scholar 
Groeters, F. R. Geographic and clonal variation in the milkweed-oleander aphid, Aphis nerii (Homoptera: Aphididae), for winged morph production, life history, and morphology in relation to host plant permanence. Evol. Ecol. 3, 327–341 (1989).Article 

Google Scholar 
Dolan, R. W., Moore, M. E. Indiana Plant Atlas. [S.M. Landry and K.N. Campbell (original application development), USF Water Institute. University of South Florida]. Butler University Friesner Herbarium, Indianapolis, Indiana (2022).McMartin, K. A., Malcolm, S. B. Defense expression in the aphid Myzocallis asclepiadis. Final Report. Pierce Cedar Creek Institute, Hastings, MI (2008).Zaya, D. N., Pearse, I. S. & Spyreas, G. Long-term trends in Midwestern Milkweed abundances and their relevance to monarch butterfly declines. Bioscience 67, 343–356. https://doi.org/10.1093/biosci/biw186 (2017).Article 

Google Scholar 
Binetruy, F., Dupraz, M., Buysse, M. & Duron, O. Surface sterilization methods impact measures of internal microbial diversity in ticks. Parasit. Vectors 12, 268 (2019).Article 

Google Scholar 
Gohl, D. M. et al. Systematic improvement of amplicon marker gene methods for increased accuracy in microbiome studies. Nat. Biotechnol. 34, 942–949 (2016).Article 
CAS 

Google Scholar 
Caporaso, J. G. et al. Global patterns of 16S rRNA diversity at a depth of millions of sequences per sample. Proc. Natl. Acad. Sci. 108, 4516–4522. https://doi.org/10.1073/pnas.1000080107 (2011).Article 
ADS 

Google Scholar 
Bolger, A. M., Lohse, M. & Usadel, B. Trimmomatic: A flexible trimmer for Illumina sequence data. Bioinformatics 30, 2114–2120 (2014).Article 
CAS 

Google Scholar 
Martin, M. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet J. 17, 10–12 (2011).Article 

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

Google Scholar 
Jousselin, E. et al. Assessment of a 16S rRNA amplicon Illumina sequencing procedure for studying the microbiome of a symbiont-rich aphid genus. Mol. Ecol. Resour. 16, 628–640. https://doi.org/10.1111/1755-0998.12478 (2016).Article 
CAS 

Google Scholar 
McMurdie, P. J. & Holmes, S. phyloseq: an R package for reproducible interactive analysis and graphics of microbiome census data. PLoS ONE 8, e61217 (2013).Article 
ADS 
CAS 

Google Scholar 
Dixon, P. VEGAN, a package of R functions for community ecology. J. Veg. Sci. 14, 927–930 (2003).Article 

Google Scholar 
Love, M. I., Huber, W. & Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 15, 550 (2014).Article 

Google Scholar 
Wright, E. S. Using DECIPHER v2. 0 to analyze big biological sequence data in R. R J. 8(1), 352 (2016).Article 

Google Scholar 
Schliep, K., Potts, A. A., Morrison, D. A. & Grimm, G. W. Intertwining phylogenetic trees and networks (No. e2054v1). PeerJ Preprints (2016).Hannula, S. E., Zhu, F., Heinen, R. & Bezemer, T. M. Foliar-feeding insects acquire microbiomes from the soil rather than the host plant. Nat. Commun. 10, 1–9 (2019).Article 
CAS 

Google Scholar 
Gomes, S. I. et al. Microbiomes of a specialist caterpillar are consistent across different habitats but also resemble the local soil microbial communities. Anim. Microbiome 2, 1–12 (2020).Article 
CAS 

Google Scholar 
Malacrinò, A. Host species identity shapes the diversity and structure of insect microbiota. Mol. Ecol. 31, 723–735. https://doi.org/10.1111/mec.16285 (2022).Article 

Google Scholar 
Colman, D. R., Toolson, E. C. & Takacs-Vesbach, C. D. Do diet and taxonomy influence insect gut bacterial communities?. Mol. Ecol. 21, 5124–5137 (2012).Article 
CAS 

Google Scholar 
Pons, I., Renoz, F., Noël, C. & Hance, T. Circulation of the cultivable symbiont Serratia symbiotica in aphids is mediated by plants. Front. Microbiol. 10, 764. https://doi.org/10.3389/fmicb.2019.00764 (2019).Article 

Google Scholar 
Li, Q. et al. Plant-mediated horizontal transmission of Hamiltonella defensa in the wheat aphid Sitobion miscanthi. J. Agric. Food Chem. 66, 13367–13377. https://doi.org/10.1021/acs.jafc.8b04828 (2018).Article 
CAS 

Google Scholar 
Jousselin, E., Cø eur d’Acier, A., Vanlerberghe-Masutti, F. & Duron, O. Evolution and diversity of A rsenophonus endosymbionts in aphids. Mol. Ecol. 22, 260–270 (2013).Article 

Google Scholar 
Nováková, E., Hypša, V. & Moran, N. A. Arsenophonus, an emerging clade of intracellular symbionts with a broad host distribution. BMC Microbiol. 9, 143 (2009).Article 

Google Scholar 
Chong, R. A. & Moran, N. A. Evolutionary loss and replacement of Buchnera, the obligate endosymbiont of aphids. ISME J. 12, 898–908 (2018).Article 
CAS 

Google Scholar 
Wulff, J. A. & White, J. A. The endosymbiont Arsenophonus provides a general benefit to soybean aphid (Hemiptera: Aphididae) regardless of host plant resistance (Rag). Environ. Entomol. 44, 574–581 (2015).Article 
CAS 

Google Scholar 
Ivens, A. B., Gadau, A., Kiers, E. T. & Kronauer, D. J. Can social partnerships influence the microbiome? Insights from ant farmers and their trophobiont mutualists. Mol. Ecol. 27, 1898–1914 (2018).Article 

Google Scholar 
Fischer, C. Y. et al. Bacteria may enhance species association in an ant–aphid mutualistic relationship. Chemoecology 25, 223–232 (2015).Article 
CAS 

Google Scholar 
Smith, R. A., Mooney, K. A. & Agrawal, A. A. Coexistence of three specialist aphids on common Milkweed, Asclepias syriaca. Ecology 89, 2187–2196 (2009).Article 

Google Scholar 
Katayama, N., Tsuchida, T., Hojo, M. K. & Ohgushi, T. aphid genotype determines intensity of ant attendance: Do endosymbionts and honeydew composition matter?. Ann. Entomol. Soc. Am. 106, 761–770 (2013).Article 
CAS 

Google Scholar 
Hansen, T. E. & Enders, L. S. Host Plant species influences the composition of milkweed and Monarch microbiomes. Front. Microbiol. https://doi.org/10.3389/fmicb.2022.840078 (2022).Article 

Google Scholar  More

Friedlingstein, P. et al. Earth Syst. Sci. Data 12, 3269–3340 (2020).Article 

Google Scholar 
LeBauer, D. S. & Treseder, K. K. Ecology 89, 371–379 (2008).Article 
PubMed 

Google Scholar 
Vitousek, P. M. & Howarth, R. W. Biogeochemistry 13, 87–115 (1991).Article 

Google Scholar 
Barker, W. et al. Nature https://doi.org/10.1038/s41586-022-05502-6 (2022).Article 

Google Scholar 
Sprent, J. I. Legume Nodulation: A Global Perspective (Wiley-Blackwell, 2009).
Google Scholar 
Gei, M. et al. Nature Ecol. Evol. 2, 1104–1111 (2018).Article 
PubMed 

Google Scholar 
Peng, J. et al. Glob. Biogeochem. Cycles 34, e2019GB006296 (2020).Article 

Google Scholar 
Batterman, S. A. et al. Nature 502, 224–227 (2013).Article 
PubMed 

Google Scholar 
Taylor, B. N. & Menge, D. N. L. Nature Plants 4, 655–661 (2018).
Google Scholar 
McCulloch, L. A. & Porder, S. New Phytol. 231, 1734–1745 (2021).Article 
PubMed 

Google Scholar 
Sheffer, E., Batterman, S. A., Levin, S. A. & Hedin, L. O. Nature Plants 1, 15182 (2015).Article 

Google Scholar 
Barron, A. R., Purves, D. W. & Hedin, L. O. Oecologia 165, 511–520 (2011).Article 
PubMed 

Google Scholar 
Adams, M. A., Turnbull, T. L., Sprent, J. I. & Buchmann, N. Proc. Natl Acad. Sci. USA 113, 4098–4103 (2016).Article 
PubMed 

Google Scholar 
Taylor, B. N. & Ostrowsky, L. R. J. Trop. Ecol. 35, 270–279 (2019).Article 

Google Scholar 
Sprent, J. I. New Phytol. 174, 11–25 (2007).Article 
PubMed 

Google Scholar  More

SamplingSediment samples were obtained from the Kap København Formation in North Greenland (82° 24′ 00″ N 22° 12′ 00″ W) in the summers of 2006, 2012 and 2016 (see Supplementary Table 3.1.1). Sampled material consisted of organic-rich permafrost and dry permafrost. Prior to sampling, profiles were cleaned to expose fresh material. Samples were hereafter collected vertically from the slope of the hills either using a 10 cm diameter diamond headed drill bit or cutting out ~40 × 40 × 40 cm blocks. Sediments were kept frozen in the field and during transportation to the lab facility in Copenhagen. Disposable gloves and scalpels were used and changed between each sample to avoid cross-contamination. In a controlled laboratory environment, the cores and blocks were further sub-sampled for material taking only the inner part of sediment cores, leaving 1.5–2 cm between the inner core and the surface that provided a subsample of approximately 6–10 g. Subsequently, all samples were stored at temperatures below −22 °C.We sampled organic-rich sediment by taking samples and biological replicates across the three stratigraphic units B1, B2 and B3, spanning 5 different sites, site: 50 (B3), 69 (B2), 74a (B1), 74b (B1) and 119 (B3). Each biological replicate from each unit at each site was further sampled in different sublayers (numbered L0–L4, Source Data 1, sheet 1).Absolute age datingIn 2014, Be and Al oxide targets from 8× 1 kg quartz-rich sand samples collected at modern depths ranging from 3 to 21 m below stream cut terraces were analysed by accelerator mass spectrometry and the cosmogenic isotope concentrations interpreted as maximum ages using a simple burial dating approach1 (26Al:10Be versus normalized 10Be). The 26Al and 10Be isotopes were produced by cosmic ray interactions with exposed quartz in regolith and bedrock surfaces in the mountains above Kap København prior to deposition. We assume that the 26Al:10Be was uniform and steady for long time periods in the upper few metres of these gradually eroding palaeo-surfaces. Once eroded by streams and hillslope processes, the quartz sand was deposited in sandy braided stream sediment, deltaic distributary systems, or the near-shore environment and remained effectively shielded from cosmic ray nucleons buried (many tens of metres) under sediment, intermittent ice shelf or ice sheet cover, and—at least during interglacials—the marine water column until final emergence. The simple burial dating approach assumes that the sand grains experienced only one burial event. If multiple burial events separated by periods of re-exposure occurred, then the starting 26Al:10Be before the last burial event would be less than the initial production ratio (6.75 to 7.42, see discussion below) owing to the relatively faster decay of 26Al during burial, and therefore the calculated burial age would be a maximum limiting age. Multiple burial events can be caused by shielding by thick glacier ice in the source area, or by sediment storage in the catchment prior to final deposition. These shielding events mean that the 26Al:10Be is lower, and therefore a calculated burial age assuming the initial production ratio would overestimate the final burial duration. We also consider that once buried, the sand grains may have been exposed to secondary cosmogenic muons (their depth would be too great for submarine nucleonic production). As sedimentation rates in these glaciated near-shore environments are relatively rapid, we show that even the muonic production would be negligible (see Supplemental Information). However, once the marine sediments emerged above sea level, in-situ production by both nucleogenic and muogenic production could alter the 26Al:10Be. The 26Al versus 10Be isochron plot reveals this complex burial history (Supplementary Information, section 3) and the concentration versus depth composite profiles for both 26Al and 10Be reveal that the shallowest samples may have been exposed during a period of time (~15,000 years ago) that is consistent with deglaciation in the area (Supplemental Information). While we interpret the individual simple burial age of all samples as a maximum limiting age of deposition of the Kap København Formation Member B, we recommend using the three most deeply shielded samples in a single depth profile to minimize the effect of post-depositional production. We then calculate a convolved probability distribution age for these three samples (KK06A, B and C). However, this calculation depends on the 26Al:10Be production ratio we use (that is, between 6.75 and 7.42) and on whether we adjust for erosion in the catchment. So, we repeat the convolved probability distribution function age for the lowest and highest production ratio and zero to maximum possible erosion rate, to obtain the minimum and maximum limiting age range at 1σ confidence (Supplementary Information, section 3). Taking the midpoint between the negative and positive 3σ confidence limits, we obtain a maximum burial age of 2.70 ± 0.46 Myr. This age is also supported by the position of those three samples on the isochron plot, which suggests the true age may not be significantly different that this maximum limiting age.Thermal ageThe extent of thermal degradation of the Kap København DNA was compared to the DNA from the Krestovka Mammoth molar. Published kinetic parameters for DNA degradation64 were used to calculate the relative rate difference over a given interval of the long-term temperature record and to quantify the offset from the reference temperature of 10 °C, thus estimating the thermal age in years at 10 °C for each sample (Supplementary Information, section 4). The mean annual air temperature (MAT) for the the Kap København sediment was taken from Funder et al. (2001)6 and for the Krestovka Mammoth the MAT was calculated using temperature data from the Cerskij Weather Station (WMO no. 251230) 68.80° N 161.28° E, 32 m from the International Research Institute Data Library (https://iri.columbia.edu/) (Supplementary Table 4.4.1).We did not correct for seasonal fluctuation for the thermal age calculation of the Kap København sediments or from the Krestovka Mammoth. We do provide theoretical average fragment length for four different thermal scenarios for the DNA in the Kap København sediments (Supplementary Table 4.4.2). A correction in the thermal age calculation was applied for altitude using the environmental lapse rate (6.49 °C km−1). We scaled the long-term temperature model of Hansen et al. (2013)65 to local estimates of current MATs by a scaling factor sufficient to account for the estimates of the local temperature decline at the last glacial maximum and then estimated the integrated rate using an activation energy (Ea) of 127 kJ mol−1 (ref. 64).Mineralogic compositionThe minerals in each of the Kap København sediment samples were identified using X-ray diffraction and their proportions were quantified using Rietveld refinement. The samples were homogenized by grinding ~1 g of sediment with ethanol for 10 min in a McCrone Mill. The samples were dried at 60 °C and added corundum (CR-1, Baikowski) as the internal standard to a final concentration of 20.0 wt%. Diffractograms were collected using a Bruker D8 Advance (Θ–Θ geometry) and the LynxEye detector (opening 2.71°), with Cu Kα1,2 radiation (1.54 Å; 40 kV, 40 mA) using a Ni-filter with thickness of 0.2 mm on the diffracted beam and a beam knife set at 3 mm. We scanned from 5–90° 2θ with a step size of 0.1° and a step time of 4 s while the sample was spun at 20 rpm. The opening of the divergence slit was 0.3° and of the antiscatter slit 3°. Primary and secondary Soller slits had an opening of 2.5° and the opening of the detector window was 2.71°. For the Rietveld analysis, we used the Profex interface for the BGMN software66,67. The instrumental parameters and peak broadening were determined by the fundamental parameters ray-tracing procedure68. A detailed description of identification of clay minerals can be found in the supporting information.AdsorptionWe used pure or purified minerals for adsorption studies. The minerals used and treatments for purifying them are listed in Supplementary Table 4.2.6. The purity of minerals was checked using X-ray diffraction with the same instrumental parameters and procedures as listed in the above section i.e., mineralogical composition. Notes on the origin, purification and impurities can be found in the Supplementary Information section 4. We used artificial seawater69 and salmon sperm DNA (low molecular weight, lyophilized powder, Sigma Aldrich) as a model for eDNA adsorption. A known amount of mineral powder was mixed with seawater and sonicated in an ultrasonic bath for 15 min. The DNA stock was then added to the suspension to reach a final concentration between 20–800 μg ml−1. The suspensions were equilibrated on a rotary shaker for 4 h. The samples were then centrifuged and the DNA concentration in the supernatant determined with UV spectrometry (Biophotometer, Eppendorf), with both positive and negative controls. All measurements were done in triplicates, and we made five to eight DNA concentrations per mineral. We used Langmuir and Freundlich equations to fit the model to the experimental isotherm and to obtain adsorption capacity of a mineral at a given equilibrium concentration.PollenThe pollen samples were extracted using the modified Grischuk protocol adopted in the Geological Institute of the Russian Academy of Science which utilizes sodium pyrophosphate and hydrofluoric acid70. Slides prepared from 6 samples were scanned at 400× magnification with a Motic BA 400 compound microscope and photographed using a Moticam 2300 camera. Pollen percentages were calculated as a proportion of the total palynomorphs including the unidentified grains. Only 4 of the 6 samples yielded terrestrial pollen counts ≥50. In these, the total palynomorphs identified ranged from 225 to 71 (mean = 170.25; median = 192.5). Identifications were made using several published keys71,72. The pollen diagram was initially compiled using Tilia version 1.5.1273 but replotted for this study using Psimpoll 4.1074.DNA recoveryFor recovery calculation, we saturated mineral surfaces with DNA. For this, we used the same protocol as for the determination of adsorption isotherms with an added step to remove DNA not adsorbed but only trapped in the interstitial pores of wet paste. This step was important because interstitial DNA would increase the amount of apparently adsorbed DNA and overestimate the recovery. To remove trapped DNA after adsorption, we redispersed the minerals in seawater. The process of redispersing the wet paste in seawater, ultracentrifugation and removal of supernatant lasted less than 2.5 min. After the second centrifugation, the wet pastes were kept frozen until extraction. We used the same extraction protocol as for the Kap København sediments. After the extraction, the DNA concentration was again determined using UV spectrometry.MetagenomesA total of 41 samples were extracted for DNA75 and converted to 65 dual-indexed Illumina sequencing libraries (including 13 negative extraction- and library controls)30. 34 libraries were thereafter subjected to ddPCR using a QX200 AutoDG Droplet Digital PCR System (Bio-Rad) following manufacturer’s protocol. Assays for ddPCR include a P7 index primer (5′-AGCAGAAGACGGCATAC-3′) (900nM), gene-targeting primer (900 nM), and a gene-targeting probe (250nM). We screened for Viridiplantae psbD (primer: 5′-TCATAATTGGACGTTGAACC-3′, probe: 5′-(FAM)ACTCCCATCATATGAAA(BHQ1)-3′) and Poaceae psbA (primer: 5′-CTCACAACTTCCCTCTAGAC-3′, probe 5′-(HEX)AGCTGCTGTTGAAGTTC(BHQ1)-3′). Additionally, 34 of the 65 libraries were enriched using targeted capture enrichment, for mammalian mitochondrial DNA using the PaleoChip Arctic1.0 bait-set31 and all libraries were hereafter sequenced on an Illumina HiSeq 4000 80 bp PE or a NovaSeq 6000 100 bp PE. We sequenced a total of 16,882,114,068 reads which, after low complexity filtering (Dust = 1), quality trimming (q ≥ 25), duplicate removal and filtering for reads longer than 29 bp (only paired read mates for NovaSeq data) resulted in 2,873,998,429 reads that were parsed for further downstream analysis. We next estimated kmer similarity between all samples using simka32 (setting heuristic count for max number of reads (-max-reads 0) and a kmer size of 31 (-kmer-size 31)), and performed a principal component analysis (PCA) on the obtained distance matrix (see Supplementary Information, ‘DNA’). We hereafter parsed all QC reads through HOLI33 for taxonomic assignment. To increase resolution and sensitivity of our taxonomic assignment, we supplemented the RefSeq (92 excluding bacteria) and the nucleotide database (NCBI) with a recently published Arctic-boreal plant database (PhyloNorway) and Arctic animal database34 as well as searched the NCBI SRA for 139 genomes of boreal animal taxa (March 2020) of which 16 partial-full genomes were found and added (Source Data 1, sheet 4) and used the GTDB microbial database version 95 as decoy. All alignments were hereafter merged using samtools and sorted using gz-sort (v. 1). Cytosine deamination frequencies were then estimated using the newly developed metaDMG, by first finding the lowest common ancestor across all possible alignments for each read and then calculating damage patterns for each taxonomic level36 (Supplementary Information, section 6). In parallel, we computed the mean read length as well as number of reads per taxonomic node (Supplementary Information, section 6). Our analysis of the DNA damage across all taxonomic levels pointed to a minimum filter for all samples at all taxonomic levels with a D-max ≥ 25% and a likelihood ratio (λ-LR) ≥ 1.5. This ensured that only taxa showing ancient DNA characteristics were parsed for downstream profiling and analysis and resulted in no taxa within any controls being found (Supplementary Information, section 6).Marine eukaryotic metagenomeWe sought to identify marine eukaryotes by first taxonomically labelling all quality-controlled reads as Eukaryota, Archaea, Bacteria or Virus using Kraken 276 with the parameters ‘–confidence 0.5 –minimum-hit-groups 3’ combined with an extra filtering step that only kept those reads with root-to-leaf score >0.25. For the initial Kraken 2 search, we used a coarse database created by the taxdb-integration workflow (https://github.com/aMG-tk/taxdb-integration) covering all domains of life and including a genomic database of marine planktonic eukaryotes63 that contain 683 metagenome-assembled genomes (MAGs) and 30 single-cell genomes (SAGs) from Tara Oceans77, following the naming convention in Delmont et al.63, we will refer to them as SMAGs. Reads labelled as root, unclassified, archaea, bacteria and virus were refined through a second Kraken 2 labelling step using a high-resolution database containing archaea, bacteria and virus created by the taxdb-integration workflow. We used the same Kraken 2 parameters and filtering thresholds as the initial search. Both Kraken 2 databases were built with parameters optimized for the study read length (–kmer-len 25 –minimizer-len 23 –minimizer-spaces 4).Reads labelled as eukaryota, root and unclassified were hereafter mapped with Bowtie278 against the SMAGs. We used MarkDuplicates from Picard (https://github.com/broadinstitute/picard) to remove duplicates and then we calculated the mapping statistics for each SMAG in the BAM files with the filterBAM program (https://github.com/aMG-tk/bam-filter). We furthermore estimated the postmortem damage of the filtered BAM files with the Bayesian methods in metaDMG and selected those SMAGs with a D-max ≥ 0.25 and a fit quality (λ-LR) higher than 1.5. The SMAGs with fewer than 500 reads mapped, a mean read average nucleotide identity (ANI) of less than than 93% and a breadth of coverage ratio and coverage evenness of less than 0.75 were removed. We followed a data-driven approach to select the mean read ANI threshold, where we explored the variation of mapped reads as a function of the mean read ANI values from 90% to 100% and identified the elbow point in the curve (Supplementary Fig. 6.11.1). We used anvi’o79 in manual mode to plot the mapping and damage results using the SMAGs phylogenomic tree inferred by Delmont et al. as reference. We used the oceanic signal of Delmont et al. as a proxy to the contemporary distribution of the SMAGs in each ocean and sea (Fig. 5 and Supplementary Information, section 6).Comparison of DNA, macrofossil and pollenTo allow comparison between records in DNA, macrofossil and pollen, the taxonomy was harmonized following the Pan Arctic Flora checklist43 and NCBI. For example, since Bennike (1990)18, Potamogeton has been split into Potamogeton and Stuckenia, Polygonym has been split to Polygonum and Bistorta, and Saxifraga was split to Saxifraga and Micranthes, whereas others have been merged, such as Melandrium with Silene40. Plant families have changed names—for instance, Gramineae is now called Poaceae and Scrophulariaceae has been re-circumscribed to exclude Plantaginaceae and Orobancheae80. We then classified the taxa into the following: category 1 all identical genus recorded by DNA and macrofossils or pollen, category 2 genera recorded by DNA also found by macrofossils or pollen including genus contained within family level classifications, category 3 taxa only recorded by DNA, category 4 taxa only recorded by macrofossils or pollen (Source Data 1).Phylogenetic placementWe sought to phylogenetically place the set of ancient taxa with the most abundant number of reads assigned, and with a sufficient number of reference sequences to build a phylogeny. These taxa include reads mapped to the chloroplast genomes of the flora genera Salix, Populus and Betula, and to the mitochondrial genomes of the fauna families Elephantidae, Cricetidae, Leporidae, as well as the subfamilies Capreolinae and Anserinae. Although the evolution of the chloroplast genome is somewhat less stable than that of the plant mitochondrial genome, it has a faster rate of evolution, and is non-recombining, and hence is more likely to contain more informative sites for our analysis than the plant mitochondria81. Like the mitochondrial genome, the chloroplast genome also has a high copy number, so that we would expect a high number of sedimentary reads mapping to it.For each of these taxa, we downloaded a representative set of either whole chloroplast or whole mitochondrial genome fasta sequences from NCBI Genbank82, including a single representative sequence from a recently diverged outgroup. For the Betula genus, we also included three chloroplast genomes from the PhyloNorway database34,83. We changed all ambiguous bases in the fasta files to N. We used MAFFT84 to align each of these sets of reference sequences, and inspected multiple sequence alignments in NCBI MSAViewer to confirm quality85. We trimmed mitochondrial alignments with insufficient quality due to highly variable control regions for Leporidae, Cricetidae and Anserinae by removing the d-loop in MegaX86.The BEAST suite49 was used with default parameters to create ultrametric phylogenetic trees for each of the five sets of taxa from the multiple sequence alignments (MSAs) of reference sequences, which were converted from Nexus to Newick format in Figtree (https://github.com/rambaut/figtree). We then passed the multiple sequence alignments to the Python module AlignIO from BioPython87 to create a reference consensus fasta sequence for each set of taxa. Furthermore, we used SNPSites88 to create a vcf file from each of the MSAs. Since SNPSites outputs a slightly different format for missing data than needed for downstream analysis, we used a custom R script to modify the vcf format appropriately. We also filtered out non-biallelic SNPs.From the damage filtered ngsLCA output, we extracted all readIDs uniquely classified to reference sequences within these respective taxa or assigned to any common ancestor inside the taxonomic group and converted these back to fastq files using seqtk (https://github.com/lh3/seqtk). We merged reads from all sites and layers to create a single read set for each respective taxon. Next, since these extracted reads were mapped against a reference database including multiple sequences from each taxon, the output files were not on the same coordinate system. To circumvent this issue and avoid mapping bias, we re-mapped each read set to the consensus sequence generated above for that taxon using bwa89 with ancient DNA parameters (bwa aln -n 0.001). We converted these reads to bam files, removed unmapped reads, and filtered for mapping quality  > 25 using samtools90. This produced 103,042, 39,306, 91,272, 182 and 129 reads for Salix, Populus, Betula, Elephantidae and Capreolinae, respectively.We next used pathPhynder62, a phylogenetic placement algorithm that identifies informative markers on a phylogeny from a reference panel, evaluates SNPs in the ancient sample overlapping these markers, and traverses the tree to place the ancient sample according to its derived and ancestral SNPs on each branch. We used the transversions-only filter to avoid errors due to deamination, except for Betula, Salix and Populus in which we used no filter due to sufficiently high coverage. Last, we investigated the pathPhynder output in each taxon set to determine the phylogenetic placement of our ancient samples (see Supplementary Information for discussion on phylogenetic placement).Based on the analysis described above we further investigated the phylogenetic placement within the genus Mammut, or mastodons. To avoid mapping reference biases in the downstream results, we first built a consensus sequence from all comparative mitochondrial genomes used in said analysis and mapped the reads identified in ngsLCA as Elephantidae to the consensus sequence. Consensus sequences were constructed by first aligning all sequences of interest using MAFFT84 and taking a majority rule consensus base in Geneious v2020.0.5 (https://www.geneious.com). We performed three analyses for phylogenetic placement of our sequence: (1) Comparison against a single representative from each Elephantidae species including the sea cow (Dugong dugon) as outgroup, (2) Comparison against a single representative from each Elephantidae species, and (3) Comparison against all published mastodon mitochondrial genomes including the Asian elephant as outgroup.For each of these analyses we first built a new reference tree using BEAST v1.10.4 (ref. 47) and repeated the previously described pathPhynder steps, with the exception that the pathPhynder tree path analysis for the Mammut SNPs was based on transitions and transversions, not restricting to only transversions due to low coverage.
Mammut americanum
We confirmed the phylogenetic placement of our sequence using a selection of Elephantidae mitochondrial reference sequences, GTR+G, strict clock, a birth-death substitution model, and ran the MCMC chain for 20,000,000 runs, sampling every 20,000 steps. Convergence was assessed using Tracer91 v1.7.2 and an effective sample size (ESS)  > 200. To determine the approximate age of our recovered mastodon mitogenome we performed a molecular dating analysis with BEAST47 v1.10.4. We used two separate approaches when dating our mastodon mitogenome, as demonstrated in a recent publication92. First, we determined the age of our sequence by comparing against a dataset of radiocarbon-dated specimens (n = 13) only. Secondly, we estimated the age of our sequence including both molecularly (n = 22) and radiocarbon-dated (n = 13) specimens using the molecular dates previously determined92. We utilized the same BEAST parameters as Karpinski et al.92 and set the age of our sample with a gamma distribution (5% quantile: 8.72 × 104, Median: 1.178 × 106; 95% quantile: 5.093 × 106; initial value: 74,900; shape: 1; scale: 1,700,000). In short, we specified a substitution model of GTR+G4, a strict clock, constant population size, and ran the Markov Chain Monte Carlo chain for 50,000,000 runs, sampling every 50,000 steps. Convergence of the run was again determined using Tracer.Molecular dating methodsIn this section, we describe molecular dating of the ancient birch (Betula) chloroplast genome using BEAST v1.10.4 (ref. 47). In principle, the genera Betula, Populus and Salix had both sufficiently high chloroplast genome coverage (with mean depth 24.16×, 57.06× and 27.04×, respectively, although this coverage is highly uneven across the chloroplast genome) and enough reference sequences to attempt molecular dating on these samples. Notably, this is one of the reasons we included a recently diverged outgroup with a divergence time estimate in each of these phylogenetic trees. However, our Populus sample clearly contained a mixture of different species, as seen from its inconsistent placement in the pathPhynder output. In particular, there were multiple supporting SNPs to both Populus balsamifera and Populus trichocarpa, and both supporting and conflicting SNPs on branches above. Furthermore, upon inspection, our Salix sample contained a surprisingly high number of private SNPs which is inconsistent with any ancient or even modern age, especially considering the number of SNPs assigned to the edges of the phylogenetic tree leading to other Salix sequences. We are unsure what causes this inconsistency but hypothesize that our Salix sample is also a mixed sample, containing multiple Salix species that diverged from the same placement branch on the phylogenetic tree at different time periods. This is supported by looking at all the reads that cover these private SNP sites, which generally appear to be from a mixed sample, with reads containing both alternate and reference alleles present at a high proportion in many cases. Alternatively, or potentially jointly in parallel, this could be a consequence of the high number of nuclear plastid DNA sequences (NUPTs) in Salix93. Because of this, we continued with only Betula.First, we downloaded 27 complete reference Betula chloroplast genome sequences and a single Alnus chloroplast genome sequence to use as an outgroup from the NCBI Genbank repository, and supplemented this with three Betula chloroplast sequences from the PhyloNorway database generated in a recent study29, for a total of 31 reference sequences. Since chloroplast sequences are circular, downloaded sequences may not always be in the same orientation or at the same starting point as is necessary for alignment, so we used custom code (https://github.com/miwipe/KapCopenhagen) that uses an anchor string to rotate the reference sequences to the same orientation and start them all from the same point. We created a MSA of these transformed reference sequences with Mafft84 and checked the quality of our alignment by eye in Seqotron94 and NCBI MsaViewer. Next, we called a consensus sequence from this MSA using the BioAlign consensus function87 in Python, which is a majority rule consensus caller. We will use this consensus sequence to map the ancient Betula reads to, both to avoid reference bias and to get the ancient Betula sample on the same coordinates as the reference MSA.From the last common ancestor output in metaDMG36, we extracted read sets for all units, sites and levels that were uniquely classified to the taxonomic level of Betula or lower, with at a minimum sequence similarity of 90% or higher to any Betula sequence, using Seqtk95. We mapped these read sets against the consensus Betula chloroplast genome using BWA89 with ancient DNA parameters (-o 2 -n 0.001 -t 20), then removed unmapped reads, quality filtered for read quality ≥25, and sorted the resulting bam files using samtools89. For the purpose of molecular dating, it is appropriate to consider these read sets as a single sample, and so we merged the resulting bam files into one sample using samtools. We used bcftools89 to make an mpileup and call a vcf file, using options for haploidy and disabling the default calling algorithm, which can slightly biases the calls towards the reference sequence, in favour of a majority call on bases that passed the default base quality cut-off of 13. We included the default option using base alignment qualities96, which we found greatly reduced the read depths of some bases and removed spurious SNPs around indel regions. Lastly, we filtered the vcf file to include only single nucleotide variants, because we do not believe other variants such as insertions or deletions in an ancient environmental sample of this type to be of sufficiently high confidence to include in molecular dating.We downloaded the gff3 annotation file for the longest Betula reference sequence, MG386368.1, from NCBI. Using custom R code97, we parsed this file and the associated fasta to label individual sites as protein-coding regions (in which we labelled the base with its position in the codon according to the phase and strand noted in the gff3 file), RNA, or neither coding nor RNA. We extracted the coding regions and checked in Seqotron94 and R that they translated to a protein alignment well (for example, no premature stop codons), both in the reference sequence and the associated positions in the ancient sequence. Though the modern reference sequence’s coding regions translated to a high-quality protein alignment, translating the associated positions in the ancient sequence with no depth cut-off leads to premature stop codons and an overall poor quality protein alignment. On the other hand, when using a depth cut-off of 20 and replacing sites in the ancient sequence which did not meet this filter with N, we see a high-quality protein alignment (except for the N sites). We also interrogated any positions in the ancient sequence which differed from the consensus, and found that any suspicious regions (for example, with multiple SNPs clustered closely together spatially in the genome) were removed with a depth cut-off of 20. Because of this, we moved forward only with sites in both the ancient and modern samples which met a depth cut-off of at least 20 in the ancient sample, which consisted of about 30% of the total sites.Next, we parsed this annotation through the multiple sequence alignment to create partitions for BEAST47. After checking how many polymorphic and total sites were in each, we decided to use four partitions: (1) sites belonging to protein-coding positions 1 and 2, (2) coding position 3, (3) RNA, or (4) non-coding and non-RNA. To ensure that these were high confidence sites, each partition also only included those positions which had at least depth 20 in the ancient sequence and had less than 3 total gaps in the multiple sequence alignment. This gave partitions which had 11,668, 5,828, 2,690 and 29,538 sites, respectively. We used these four partitions to run BEAST47 v1.10.4, with unlinked substitution models for each partition and a strict clock, with a different relative rate for each partition. (There was insufficient information in these data to infer between-lineage rate variation from a single calibration). We assigned an age of 0 to all of the reference sequences, and used a normal distribution prior with mean 61.1 Myr and standard deviation 1.633 Myr for the root height48; standard deviation was obtained by conservatively converting the 95% HPD to z-scores. For the overall tree prior, we selected the coalescent model. The age of the ancient sequence was estimated following the overall procedures of Shapiro et al. (2011)98. To assess sensitivity to prior choice for this unknown date, we used two different priors, namely a gamma distribution metric towards a younger age (shape = 1, scale = 1.7); and a uniform prior on the range (0, 10 Myr). We also compared two different models of rate variation among sites and substitution types within each partition, namely a GTR+G with four rate categories, and base frequencies estimated from the data, and the much simpler Jukes Cantor model, which assumed no variation between substitution types nor sites within each partition. All other priors were set at their defaults. Neither rate model nor prior choice had a qualitative effect on results (Extended Data Fig. 10). We also ran the coding regions alone, since they translated correctly and are therefore highly reliable sites and found that they gave the same median and a much larger confidence interval, as expected when using fewer sites (Extended Data Fig. 10). We ran each Markov chain Monte Carlo for a total of 100 million iterations. After removing a burn-in of the first 10%, we verified convergence in Tracer91 v1.7.2 (apparent stationarity of traces, and all parameters having an Effective Sample Size  > 100). We also verified that the resulting MCC tree from TreeAnnotator47 had placed the ancient sequence phylogenetically identically to pathPhynder62 placement, which is shown in Extended Data Fig. 9. For our major results, we report the uniform ancient age prior, and the GTR+G4 model applied to each of the four partitions. The associated XML is given in Source Data 3. The 95% HPD was (2.0172,0.6786) for the age of the ancient Betula chloroplast sequence, with a median estimate of 1.323 Myr, as shown in Fig. 2.Reporting summaryFurther information on research design is available in the Nature Portfolio Reporting Summary linked to this article. More