Ecology

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

Keeling, P. J. & Burki, F. Progress towards the tree of eukaryotes. Curr. Biol. 29, R808–R817 (2019).Article 
CAS 

Google Scholar 
Gawryluk, R. M. R. et al. Non-photosynthetic predators are sister to red algae. Nature 572, 240–243 (2019).Article 
CAS 

Google Scholar 
Janouškovec, J. et al. A new lineage of eukaryotes illuminates early mitochondrial genome reduction. Curr. Biol. 27, 3717–3724 (2017).Article 

Google Scholar 
Lax, G. et al. Hemimastigophora is a novel supra-kingdom-level lineage of eukaryotes. Nature 564, 410–414 (2018).Article 
ADS 
CAS 

Google Scholar 
Oren, A. Prokaryote diversity and taxonomy: current status and future challenges. Philos. Trans. R. Soc. Lond. B 359, 623–638 (2004).Article 
CAS 

Google Scholar 
Shu, W. S. & Huang, L. N. Microbial diversity in extreme environments. Nat. Rev. Microbiol. 20, 219–235 (2022).Article 
CAS 

Google Scholar 
Massana, R., del Campo, J., Sieracki, M. E., Audic, S. & Logares, R. Exploring the uncultured microeukaryote majority in the oceans: reevaluation of ribogroups within stramenopiles. ISME J. 8, 854–866 (2014).Article 

Google Scholar 
de Vargas, C. et al. Eukaryotic plankton diversity in the sunlit ocean. Science 348, 1261605 (2015).Article 

Google Scholar 
Flegontova, O. et al. Extreme diversity of diplonemid eukaryotes in the ocean. Curr. Biol. 26, 3060–3065 (2016).Article 
CAS 

Google Scholar 
Ahlering, M. A. & Carrel, J. E. Predators are rare even when they are small. Oikos 95, 471–475 (2001).Article 

Google Scholar 
Hehenberger, E. et al. Novel predators reshape holozoan phylogeny and reveal the presence of a two-component signaling system in the ancestor of animals. Curr. Biol. 27, 2043–2050 (2017).Article 
CAS 

Google Scholar 
Tikhonenkov, D. V. et al. Description of Colponema vietnamica sp. n. and Acavomonas peruviana n. gen. n. sp., two new alveolate phyla (Colponemidia nom. nov. and Acavomonidia nom. nov.) and their contributions to reconstructing the ancestral state of alveolates and eukaryotes. PLoS ONE 9, e95467 (2014).Article 
ADS 

Google Scholar 
Tikhonenkov, D. V. et al. New lineage of microbial predators adds complexity to reconstructing the evolutionary origin of animals. Curr. Biol. 30, 4500–4509 (2020).Article 
CAS 

Google Scholar 
Mylnikov, A. P. & Tikhonenkov, D. V. The new alveolate carnivorous flagellate Colponema marisrubri sp. n. (Colponemida, Alveolata) from the Red Sea. Zool. Zh. 88, 1163–1169 (2009).
Google Scholar 
Strassert, J. F. H., Irisarri, I., Williams, T. A. & Burki, F. A molecular timescale for eukaryote evolution with implications for the origin of red algal-derived plastids. Nat. Commun. 12, 1879 (2021).Article 
ADS 
CAS 

Google Scholar 
Rodriguez-Ezpeleta, N. et al. Detecting and overcoming systematic errors in genome-scale phylogenies. Syst. Biol. 56, 389–399 (2007).Article 
CAS 

Google Scholar 
Strassert, J. F. H., Jamy, M., Mylnikov, A. P., Tikhonenkov, D. V. & Burki, F. New phylogenomic analysis of the enigmatic phylum Telonemia further resolves the eukaryote tree of life. Mol. Biol. Evol. 36, 757–765 (2019).Article 
CAS 

Google Scholar 
Lanfear, R., Kokko, H. & Eyre-Walker, A. Population size and the rate of evolution. Trends Ecol. Evol. 29, 33–41 (2014).Article 

Google Scholar 
Bahler, M. & Rhoads, A. Calmodulin signaling via the IQ motif. FEBS Lett. 513, 107–113 (2002).Article 
CAS 

Google Scholar 
Schaffer, D. E., Iyer, L. M., Burroughs, A. M. & Aravind, L. Functional innovation in the evolution of the calcium-dependent system of the eukaryotic endoplasmic reticulum. Front. Genet. 11, 34 (2020).Article 

Google Scholar 
Morita-Yamamuro, C. et al. The Arabidopsis gene CAD1 controls programmed cell death in the plant immune system and encodes a protein containing a MACPF domain. Plant Cell Physiol. 46, 902–912 (2005).Article 
CAS 

Google Scholar 
Rosado, C. J. et al. The MACPF/CDC family of pore-forming toxins. Cell. Microbiol. 10, 1765–1774 (2008).Article 
CAS 

Google Scholar 
Ishino, T., Chinzei, Y. & Yuda, M. A Plasmodium sporozoite protein with a membrane attack complex domain is required for breaching the liver sinusoidal cell layer prior to hepatocyte infection. Cell. Microbiol. 7, 199–208 (2005).Article 
CAS 

Google Scholar 
Satoh, H., Oshiro, N., Iwanaga, S., Namikoshi, M. & Nagai, H. Characterization of PsTX-60B, a new membrane-attack complex/perforin (MACPF) family toxin, from the venomous sea anemone Phyllodiscus semoni. Toxicon 49, 1208–1210 (2007).Article 
CAS 

Google Scholar 
Tikhonenkov, D. V., Mazei, Y. A. & Embulaeva, E. A. Degradation succession of heterotrophic flagellate communities in microcosms. Zh. Obs. Biol. 69, 57–64 (2008).CAS 

Google Scholar 
Tikhonenkov, D. V. et al. On the origin of TSAR: morphology, diversity and phylogeny of Telonemia. Open Biol. 12, 210325 (2022).Article 
CAS 

Google Scholar 
Picelli, S. et al. Full-length RNA-seq from single cells using Smart-seq2. Nat. Protoc. 9, 171–181 (2014).Article 
CAS 

Google Scholar 
Keeling, P. J., Poulson, N. & McFadden, G. I. Phylogenetic diversity of parabasalian symbionts from termites, including the phylogenetic position of Pseudotrypanosoma and Trichonympha. J. Eukaryot. Microbiol. 45, 643–650 (1998).Article 
CAS 

Google Scholar 
Medlin, L., Elwood, H. J., Stickel, S. & Sogin, M. L. The characterization of enzymatically amplified eukaryotic 16S-like rRNA-coding regions. Gene 71, 491–499 (1988).Article 
CAS 

Google Scholar 
Tikhonenkov, D. V., Janouškovec, J., Keeling, P. J. & Mylnikov, A. P. The morphology, ultrastructure and SSU rRNA gene sequence of a new freshwater flagellate, Neobodo borokensis n. sp. (Kinetoplastea, Excavata). J. Eukaryot. Microbiol. 63, 220–232 (2016).Article 
CAS 

Google Scholar 
Andrews, S. FastQC: a quality control tool for high throughput sequence data (Babraham Bioinformatics, 2010); https://www.bioinformatics.babraham.ac.uk/projects/fastqc/.Zhang, J., Kobert, K., Flouri, T. & Stamatakis, A. PEAR: a fast and accurate Illumina Paired-End reAd mergeR. Bioinformatics 30, 614–620 (2013).Article 

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 
Grabherr, M. G. et al. Full-length transcriptome assembly from RNA-seq data without a reference genome. Nat. Biotechnol. 29, 644–652 (2011).Article 
CAS 

Google Scholar 
Altschul, S. F., Gish, W., Miller, W., Myers, E. W. & Lipman, D. J. Basic local alignment search tool. J. Mol. Biol. 215, 403–410 (1990).Article 
CAS 

Google Scholar 
Laetsch, D. R. & Blaxter, M. L. BlobTools: interrogation of genome assemblies. F1000Research 6, 1287 (2017).Article 

Google Scholar 
Haas, B. J. et al. Denovo transcript sequence reconstruction from RNA-seq using the Trinity platform for reference generation and analysis. Nat. Protoc. 8, 1494–1512 (2013).Article 
CAS 

Google Scholar 
Li, W. & Godzik, A. Cd-hit: a fast program for clustering and comparing large sets of protein or nucleotide sequences. Bioinformatics 22, 1658–1659 (2006).Article 
CAS 

Google Scholar 
Buchfink, B., Xie, C. & Huson, D. H. Fast and sensitive protein alignment using DIAMOND. Nat. Methods 12, 59–60 (2015).Article 
CAS 

Google Scholar 
Shen, W. & Ren, H. TaxonKit: a practical and efficient NCBI taxonomy toolkit. J. Genet. Genomics 48, 844–850 (2021).Richter, D. J. et al. EukProt: a database of genome-scale predicted proteins across the diversity of eukaryotes. Peer Community Journal 2, e56 (2022).Simao, F. A., Waterhouse, R. M., Ioannidis, P., Kriventseva, E. V. & Zdobnov, E. M. BUSCO: assessing genome assembly and annotation completeness with single-copy orthologs. Bioinformatics 31, 3210–3212 (2015).Article 
CAS 

Google Scholar 
Kanehisa, M., Furumichi, M., Sato, Y., Ishiguro-Watanabe, M. & Tanabe, M. KEGG: integrating viruses and cellular organisms. Nucleic Acids Res. 49, D545–D551 (2021).Article 
CAS 

Google Scholar 
Moriya, Y., Itoh, M., Okuda, S., Yoshizawa, A. C. & Kanehisa, M. KAAS: an automatic genome annotation and pathway reconstruction server. Nucleic Acids Res. 35, W182–W185 (2007).Article 

Google Scholar 
Burki, F. The eukaryotic tree of life from a global phylogenomic perspective. Cold Spring Harb. Perspect. Biol. 6, a016147 (2014).Article 

Google Scholar 
Waskom, M. et al. mwaskom/Seaborn: v0.8.1 (September 2017). Zenodo https://doi.org/10.5281/zenodo.883859 (2017).Eddy, S. R. Accelerated profile HMM searches. PLoS Comput. Biol. 7, e1002195 (2011).Article 
ADS 
MathSciNet 
CAS 

Google Scholar 
Finn, R. D. et al. The Pfam protein families database: towards a more sustainable future. Nucleic Acids Res. 44, D279–D285 (2016).Article 
CAS 

Google Scholar 
Letunic, I. & Bork, P. 20 years of the SMART protein domain annotation resource. Nucleic Acids Res. 46, D493–D496 (2018).Article 
CAS 

Google Scholar 
Almagro Armenteros, J. J. et al. SignalP 5.0 improves signal peptide predictions using deep neural networks. Nat. Biotechnol. 37, 420–423 (2019).Article 
CAS 

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

Google Scholar 
Burns, J. A., Pittis, A. A. & Kim, E. Gene-based predictive models of trophic modes suggest Asgard archaea are not phagocytotic. Nat. Ecol. Evol. 2, 697–704 (2018).Article 

Google Scholar 
Emms, D. M. & Kelly, S. OrthoFinder: phylogenetic orthology inference for comparative genomics. Genome Biol. 20, 238 (2019).Article 

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

Google Scholar 
Hall, T. A. BioEdit: a user-friendly biological sequence alignment editor and analysis program for Windows 95/98/NT. Nucleic Acids Symp. Ser. 41, 95–98 (1999).CAS 

Google Scholar 
Minh, B. Q. et al. IQ-TREE 2: new models and efficient methods for phylogenetic inference in the genomic era. Mol. Biol. Evol. 37, 1530–1534 (2020).Article 
CAS 

Google Scholar 
Whelan, S., Irisarri, I. & Burki, F. PREQUAL: detecting non-homologous characters in sets of unaligned homologous sequences. Bioinformatics 34, 3929–3930 (2018).CAS 

Google Scholar 
Capella-Gutierrez, S., Silla-Martinez, J. M. & Gabaldon, T. trimAl: a tool for automated alignment trimming in large-scale phylogenetic analyses. Bioinformatics 25, 1972–1973 (2009).Article 
CAS 

Google Scholar 
Roure, B., Rodriguez-Ezpeleta, N. & Philippe, H. SCaFoS: a tool for selection, concatenation and fusion of sequences for phylogenomics. BMC Evol. Biol. 7, S2 (2007).Article 

Google Scholar 
Lartillot, N., Rodrigue, N., Stubbs, D. & Richer, J. PhyloBayes MPI: phylogenetic reconstruction with infinite mixtures of profiles in a parallel environment. Syst. Biol. 62, 611–615 (2013).Article 
CAS 

Google Scholar 
Dayhoff, M., Schwartz, R. & Orcutt, B. in Atlas of Protein Sequence and Structure (ed. Dayhoff, M.) 345–352 (National Biomedical Research Foundation, 1978).Susko, E. & Roger, A. J. On reduced amino acid alphabets for phylogenetic inference. Mol. Biol. Evol. 24, 2139–2150 (2007).Article 
CAS 

Google Scholar 
Lartillot, N. & Philippe, H. A Bayesian mixture model for across-site heterogeneities in the amino-acid replacement process. Mol. Biol. Evol. 21, 1095–1109 (2004).Article 
CAS 

Google Scholar 
Quang le, S., Gascuel, O. & Lartillot, N. Empirical profile mixture models for phylogenetic reconstruction. Bioinformatics 24, 2317–2323 (2008).Article 

Google Scholar 
Wang, H. C., Minh, B. Q., Susko, E. & Roger, A. J. Modeling site heterogeneity with posterior mean site frequency profiles accelerates accurate phylogenomic estimation. Syst. Biol. 67, 216–235 (2018).Article 
CAS 

Google Scholar 
Kück, P. & Struck, T. H. BaCoCa—a heuristic software tool for the parallel assessment of sequence biases in hundreds of gene and taxon partitions. Mol. Phylogenet. Evol. 70, 94–98 (2014).Article 

Google Scholar 
Shimodaira, H. An approximately unbiased test of phylogenetic tree selection. Syst. Biol. 51, 492–508 (2002).Article 

Google Scholar 
Kumar, S., Stecher, G. & Tamura, K. MEGA7: molecular evolutionary genetics analysis version 7.0 for bigger datasets. Mol. Biol. Evol. 33, 1870–1874 (2016).Article 
CAS 

Google Scholar 
Bankevich, A. et al. SPAdes: a new genome assembly algorithm and its applications to single-cell sequencing. J. Comput. Biol. 19, 455–477 (2012).Article 
MathSciNet 
CAS 

Google Scholar 
Dierckxsens, N., Mardulyn, P. & Smits, G. NOVOPlasty: de novo assembly of organelle genomes from whole genome data. Nucleic Acids Res. 45, e18 (2017).
Google Scholar 
Kuznetsov, A. & Bollin, C. J. in Multiple Sequence Alignment (ed. Katoh, K.) 261–295 (Springer, 2021).Lohse, M., Drechsel, O., Kahlau, S. & Bock, R. OrganellarGenomeDRAW—a suite of tools for generating physical maps of plastid and mitochondrial genomes and visualizing expression data sets. Nucleic Acids Res. 41, W575–W581 (2013).Article 

Google Scholar 
Johnson, P. Z., Kasprzak, W. K., Shapiro, B. A. & Simon, A. E. RNA2Drawer: geometrically strict drawing of nucleic acid structures with graphical structure editing and highlighting of complementary subsequences. RNA Biol. 16, 1667–1671 (2019).Article 

Google Scholar 
Burger, G., Gray, M. W., Forget, L. & Lang, B. F. Strikingly bacteria-like and gene-rich mitochondrial genomes throughout jakobid protists. Genome Biol. Evol. 5, 418–438 (2013).Article 

Google Scholar 
Criscuolo, A. & Gribaldo, S. BMGE (Block Mapping and Gathering with Entropy): a new software for selection of phylogenetic informative regions from multiple sequence alignments. BMC Evol. Biol. 10, 210 (2010).Article 

Google Scholar 
Zhang, D. et al. PhyloSuite: an integrated and scalable desktop platform for streamlined molecular sequence data management and evolutionary phylogenetics studies. Mol. Ecol. Resour. 20, 348–355 (2020).Article 

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

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

Google Scholar 
Massana, R. et al. Marine protist diversity in European coastal waters and sediments as revealed by high-throughput sequencing. Environ. Microbiol. 17, 4035–4049 (2015).Article 
CAS 

Google Scholar 
Gendron, E. M. S., Darcy, J. L., Hell, K. & Schmidt, S. K. Structure of bacterial and eukaryote communities reflect in situ controls on community assembly in a high-alpine lake. J. Microbiol. 57, 852–864 (2019).Article 
CAS 

Google Scholar 
Minerovic, A. D. et al. 18S-V9 DNA metabarcoding detects the effect of water-quality impairment. Ecol. Indic. 113, 106225 (2020).Article 
CAS 

Google Scholar 
Pearman, J. K. et al. Cross-shelf investigation of coral reef cryptic benthic organisms reveals diversity patterns of the hidden majority. Sci. Rep. 8, 8090 (2018).Article 
ADS 
CAS 

Google Scholar 
Rodas, A. M. et al. Eukaryotic plankton communities across reef environments in Bocas del Toro Archipelago, Panamá. Coral Reefs 39, 1453–1467 (2020).Article 

Google Scholar 
Schoenle, A. et al. High and specific diversity of protists in the deep-sea basins dominated by diplonemids, kinetoplastids, ciliates and foraminiferans. Commun. Biol. 4, 501 (2021).Article 
CAS 

Google Scholar 
Schulhof, M. A. et al. Sierra Nevada mountain lake microbial communities are structured by temperature, resources and geographic location. Mol. Ecol. 29, 2080–2093 (2020).Article 
CAS 

Google Scholar 
Yi, Z. et al. High-throughput sequencing of microbial eukaryotes in Lake Baikal reveals ecologically differentiated communities and novel evolutionary radiations. FEMS Microbiol. Ecol. 93, fix073 (2017). More

Fernández-Martínez, M. et al. Nutrient availability as the key regulator of global forest carbon balance. Nat. Clim. Chang. 4, 471–476 (2014).Article 
ADS 

Google Scholar 
Wright, S. J. Plant responses to nutrient addition experiments conducted in tropical forests. Ecol. Monogr. 89, e01382 (2019).Article 

Google Scholar 
Levy-Varon, J. H. et al. Tropical carbon sink accelerated by symbiotic dinitrogen fixation. Nat. Commun. 10, 5637 (2019).Article 
ADS 
CAS 

Google Scholar 
Batterman, S. A. et al. Key role of symbiotic dinitrogen fixation in tropical forest secondary succession. Nature 502, 224–227 (2013).Article 
ADS 
CAS 

Google Scholar 
Ter Steege, H. et al. Continental-scale patterns of canopy tree composition and function across Amazonia. Nature 443, 444–447 (2006).Article 
ADS 

Google Scholar 
Hedin, L. O., Brookshire, E. N. J., Menge, D. N. L. & Barron, A. R. The nitrogen paradox in tropical forest ecosystems. Annu. Rev. Ecol. Evol. Syst. 40, 613–635 (2009).Article 

Google Scholar 
Menge, D. N. L. et al. Patterns of nitrogen-fixing tree abundance in forests across Asia and America. J. Ecol. 107, 2598–2610 (2019).Article 
CAS 

Google Scholar 
Matson, W. J.Jr Herbivory in relation to plant nitrogen content. Annu. Rev. Ecol. Syst. 11, 119–161 (1980).Article 

Google Scholar 
Coley, P. D., Bateman, M. L. & Kusar, T. A. The effects of plant quality on caterpillar growth and defense against natural enemies. Oikos 115, 219–228 (2006).Article 

Google Scholar 
Wieder, W. R., Cleveland, C. C., Smith, W. K. & Todd-Brown, K. Future productivity and carbon storage limited by terrestrial nutrient availability. Nat. Geosci. 8, 441–444 (2015).Article 
ADS 
CAS 

Google Scholar 
Barron, A. R., Purves, D. W. & Hedin, L. O. Facultative nitrogen fixation by canopy legumes in a lowland tropical forest. Oecologia 165, 511–520 (2011).Article 
ADS 

Google Scholar 
McCulloch, L. A. & Porder, S. Light fuels while nitrogen suppresses symbiotic nitrogen fixation hotspots in neotropical canopy gap seedlings. New Phytol. 231, 1734–1745 (2021).Article 
CAS 

Google Scholar 
Brookshire, E. N. J. et al. Symbiotic N fixation is sufficient to support net aboveground biomass accumulation in a humid tropical forest. Sci Rep. 9, 7571 (2019).Article 
ADS 
CAS 

Google Scholar 
Gei, M. et al. Legume abundance along successional and rainfall gradients in Neotropical forests. Nat. Ecol. Evol. 2, 1104–1111 (2018).Article 

Google Scholar 
Vance, C. P. in Nitrogen-fixing Leguminous Symbioses. Nitrogen Fixation: Origins, Applications, and Research Progress, Vol. 7 (eds Dilworth, M. J. et al.) (Springer, 2008).Vitousek, P. M. & Howarth, R. W. Nitrogen limitation on land and in the sea: how can it occur? Biogeochemistry 13, 87–115 (1991).Article 

Google Scholar 
Menge, D. N. L., Levin, S. A. & Hedin, L. O. Evolutionary tradeoffs can select against nitrogen fixation and thereby maintain nitrogen limitation. Proc. Natl Acad. Sci. USA 105, 1573–1578 (2008).Article 
ADS 
CAS 

Google Scholar 
Sheffer, E., Batterman, S. A., Levin, S. A. & Hedin, L. O. Biome-scale nitrogen fixation strategies selected by climatic constraints on nitrogen cycle. Nat. Plants 1, 15182 (2015).Article 
CAS 

Google Scholar 
Vitousek, P. M. & Field, C. B. Ecosystem constraints to symbiotic nitrogen fixers: a simple model and its implications. Biogeochemistry 46, 179–202 (1999).Article 
CAS 

Google Scholar 
Coley, P. D. & Barone, J. A. Herbivory and plant defenses in tropical forests. Annu. Rev. Ecol. Syst. 27, 305–335 (1996).Article 

Google Scholar 
Fyllas, N. M. et al. Basin-wide variations in foliar properties of Amazonian forest: phylogeny, soils and climate. Biogeosciences 6, 2677–2708 (2009).Article 
ADS 

Google Scholar 
Batterman, S. A. et al. Phosphatase activity and nitrogen fixation reflect species differences, not nutrient trading or nutrient balance, across tropical rainforest trees. Ecol. Lett. 21, 1486–1495 (2018).Article 

Google Scholar 
Menge, D. N. L., Wolf, A. A. & Funk, J. L. Diversity of nitrogen fixation strategies in Mediterranean legumes. Nat. Plants 1, 15064 (2015).Article 
CAS 

Google Scholar 
Ritchie, M. E. & Tilman, D. Responses of legumes to herbivores and nutrients during succession on a nitrogen-poor soil. Ecol. Soc. Am. 76, 2648–2655 (1995).
Google Scholar 
Taylor, B. N. & Ostrowsky, L. R. Nitrogen-fixing and non-fixing trees differ in leaf chemistry and defence but not herbivory in a lowland Costa Rican rain forest. J. Trop. Ecol. 35, 270–279 (2019).Article 

Google Scholar 
Endara, M.-J. et al. Coevolutionary arms race versus host defense chase in a tropical herbivore–plant system. Proc. Natl Acad. Sci. USA 114, E7499–E7505 (2017).Article 
CAS 

Google Scholar 
Kursar, T. A. & Coley, P. D. Convergence in defense syndromes of young leaves in tropical rainforests. Biochem. Syst. Ecol. 31, 929–949 (2003).Article 
CAS 

Google Scholar 
Kursar, T. A. et al. The evolution of antiherbivore defenses and their contribution to species coexistence in the tropical tree genus Inga. Proc. Natl Acad. Sci. USA 106, 18073–18078 (2009).Article 
ADS 
CAS 

Google Scholar 
Taylor, B. N. & Menge, D. N. L. Light regulates tropical symbiotic nitrogen fixation more strongly than soil nitrogen. Nat. Plants 4, 655–661 (2018).Article 
CAS 

Google Scholar 
Adams, M., Turnbull, T., Sprent, J. & Buchmann, N. Legumes are different: leaf nitrogen, photosynthesis, and water use efficiency. Proc. Natl Acad. Sci. USA 113, 4098–4103 (2016).Article 
ADS 
CAS 

Google Scholar 
Coley, P. D. Effects of plant growth rate and leaf lifetime on the amount and type of anti-herbivore defense. Oecologia 74, 531–536 (1988).Article 
ADS 
CAS 

Google Scholar 
Batterman, S. A., Wurzburger, N. & Hedin, L. O. Nitrogen and phosphorus interact to control tropical symbiotic N2 fixation: a test in Inga punctata. J. Ecol. 101, 1400–1408 (2013).Article 
CAS 

Google Scholar 
Eichhorn, M. P., Nilus, R., Compton, S. G., Hartley, S. E. & Burslem, D. F. R. P. Herbivory of tropical rain forest tree seedlings correlates with future mortality. Ecology 91, 1092–1101 (2010).Article 

Google Scholar 
Wink, M. Evolution of secondary metabolites in legumes (Fabaceae). South African J. Bot. 89, 164–175 (2013).Article 
CAS 

Google Scholar 
Currano, E. D. & Jacobs, B. F. Bug-bitten leaves from the early Miocene of Ethiopia elucidate the impacts of plant nutrient concentrations and climate on insect herbivore communities. Glob. Planet. Change 207, 103655 (2021).Article 

Google Scholar 
Wieder, W. R., Cleveland, C. C., Lawrence, D. M. & Bonan, G. B. Effects of model structural uncertainty on carbon cycle projections: biological nitrogen fixation as a case study. Environ. Res. Lett. 10, 044016 (2015).Article 
ADS 

Google Scholar 
Sprent, J. I. Legume Nodulation: A Global Perspective (John Wiley, 2009).Leigh, E. G. Jr Tropical Forest Ecology: A View from Barro Colorado Island (Oxford Univ. Press, 1999).Comita, L. S., Muller-Landau, H. C., Aguilar, S. & Hubbell, S. P. Asymmetric density dependence shapes species abundances in a tropical tree community. Science 329, 330–332 (2010).Article 
ADS 
CAS 

Google Scholar 
Queenborough, S. A., Metz, M. R., Valencia, R. & Wright, S. J. Demographic consequences of chromatic leaf defence in tropical tree communities: do red young leaves increase growth and survival? Ann. Bot. 112, 677–684 (2013).Article 

Google Scholar 
Schneider, C. A., Rasband, W. S. & Eliceiri, K. W. NIH Image to ImageJ: 25 years of image analysis. Nat. Methods 9, 671 (2012).Article 
CAS 

Google Scholar 
Pasquini, S. C. & Santiago, L. S. Nutrients limit photosynthesis in seedlings of a lowland tropical forest tree species. Oecologia 168, 311–319 (2012).Article 
ADS 
CAS 

Google Scholar 
Collalti, A. & Prentice, I. C. Is NPP proportional to GPP? Waring’s hypothesis 20 years on. Tree Physiol. 39, 1473–1483 (2019).Article 
CAS 

Google Scholar 
Westbrook, J. W. et al. What makes a leaf tough? Patterns of correlated evolution between leaf toughness traits and demographic rates among 197 shade-tolerant woody species in a Neotropical forest. Am. Nat. 177, 800–811 (2011).Article 

Google Scholar 
Wright, S. J. et al. Functional traits and the growth–mortality trade‐off in tropical trees. Ecology 91, 3664–3674 (2010).Article 

Google Scholar 
Kitajima, K. et al. How cellulose-based leaf toughness and lamina density contribute to long leaf lifespans of shade-tolerant species. New Phytol. 195, 640–652 (2012).Article 

Google Scholar 
Kitajima, K., Wright, S. J. & Westbrook, J. W. Leaf cellulose density as the key determinant of inter- and intra-specific variation in leaf fracture toughness in a species-rich tropical forest. Interface Focus https://doi.org/10.1098/rsfs.2015.0100 (2016).Sedio, B. E., Echeverri, J. C. R., Boya, C. A. & Wright, S. J. Sources of variation in foliar secondary chemistry in a tropical forest tree community. Ecology 98, 616–623 (2017).Article 

Google Scholar 
Brooks, M. E. et al. glmmTMB balances speed and flexibility among packages for zero-inflated generalized linear mixed modeling. R J. 9, 378–400 (2017).Article 

Google Scholar 
Smithson, M. & Verkuilen, J. A better lemon squeezer? Maximum-likelihood regression with beta-distributed dependent variables. Psychol. Methods 11, 54–71 (2006).Article 

Google Scholar 
Murphy, S. J., Xu, K. & Comita, L. S. Tree seedling richness, but not neighborhood composition, influences insect herbivory in a temperate deciduous forest community. Ecol. Evol. 6, 6310–6319 (2016).Article 

Google Scholar 
Bates, D., Mächler, M., Bolker, B. & Walker, S. Fitting linear mixed-effects models using lme4. Preprint at https://arxiv.org/abs/1406.5823 (2014).Moles, A. T. & Westoby, M. Do small leaves expand faster than large leaves, and do shorter expansion times reduce herbivore damage? Oikos 90, 517–524 (2000).Article 

Google Scholar 
Bürkner, P. C. brms: An R package for Bayesian multilevel models using Stan. J. Stat. Softw. https://doi.org/10.18637/jss.v080.i01 (2017). More

RESEARCH BRIEFINGS
07 December 2022

Ancient environmental DNA from northern Greenland opens a new chapter in genetic research, demonstrating that it is possible to track the ecology and evolution of biological communities two million years ago. The record shows an open boreal-forest ecosystem inhabited by large animals such as mastodons and reindeer. More

Restoring the connection between people and the rest of nature hinges on whole-system science, actions and negotiations.
Those who think about and practise sustainability are constantly looking for holistic interpretations of the world and are trying to understand systemic relations, networks and connections. Biodiversity has all of these things. It shows how every species needs other species to exist and thrive. It shows that all living organisms are part of a sophisticated and fascinating system made up of myriads of links. And humans are undoubtedly a part of it.
Credit: Pulsar Imagens / Alamy Stock PhotoIn the realm of sustainability, experts also ponder about time: how can life exist and thrive over time? Indeed, the above mentioned fascinating system evolves over time. And, over time, it has to adapt to unexpected change. It does that well when it is healthy, and less well when it is ill and constantly disturbed.For a long time, man-made impacts kept accumulating almost completely unchecked by societies, until the consequences for human well-being became untenable. Nowadays, environmental crises make the headlines regularly. They are nothing but the result of a broken connection between people and the rest of nature.Climate change is one major outcome of the broken human–rest of nature connection and has wide ramifications for both people and the planet. We now face imminent disaster, unequally across the world, yet addressing climate change remains an incredibly thorny task. Country representatives from most nations around the world meet regularly at the Conference of the Parties (COP) to the United Nations Framework Convention on Climate Change (UNFCC) — most recently at COP27, which was held in Egypt — to continue the debate on what actions are needed to move the climate agenda forward, all while disasters continue to hit the most vulnerable populations. The world has seen 27 COP meetings to the UNFCC so far; one wonders how many more meetings will be needed to see real change happen.Interestingly, country representatives also meet regularly to discuss biodiversity protection; biodiversity decline — the other major consequence of the broken human–rest of nature connection — is just as worrying, with severe and ramified implications that are still largely underappreciated by decision-makers. These gatherings are the COP meetings to the Convention on Biological Diversity (CBD). Last year, we wrote about the then forthcoming COP15 to the CBD (Nat. Sustain. 4, 189; 2021), the meeting in which the new conservation targets to be met by 2030 were to be agreed. We highlighted the extent to which experts worried that those new targets might not go far enough. The meeting was postponed more than once due to the COVID-19 pandemic, and it is finally happening on 7 December 2022, in Montreal, Canada. The world has already seen 15 COP meetings to the CBD, how many more meetings will be needed for the biodiversity crisis to be averted?But let’s go back to thinking about sustainability. Experts look for holistic visions of the world. Here is an interesting example of what holism means. Biodiversity decline and climate change are both the result of the broken connection between people and the rest of nature, they ultimately have the same, deep roots. They are mutually reinforcing phenomena: unhealthy biodiversity contributes to climate change, and climate change makes biodiversity ill. All this is bad news for human and planetary well-being. The climate–biodiversity conundrum, at least to some degree, has been recognized at a higher level — during COP27, leaders dedicated one day to biodiversity.Yet, given that these issues are highly interconnected and have the same origin, why is the world insisting on discussing them as separate agendas? Why are we still holding two separate COPs? How are these meetings going to promote any fruitful synergy? How will they lead people to reconnect with the rest of nature? Country representatives should be breaking silos, embracing holism and bringing these intertwined issues, and their multiple ramifications, to the same negotiating table.Nature Sustainability welcomes the long-awaited COP15 to the CBD and hopes that countries will agree on feasible yet ambitious 2030 targets to protect and enhance biodiversity. But most of all, we hope that all of the experts and leaders involved in addressing the environmental crises embrace holism to promote meaningful actions across the world aimed at restoring people’s connection with the rest of nature. We are eager to see progress to this end. In the meantime, the collection we started in March 2021 with Nature Ecology & Evolution has been updated to renew our support to the biodiversity community. More

Fossil record of insectsWe compiled all species-level fossil occurrences of insects using https://paleobiodb.org/ (PBDB) as a starting point (downloaded October 12, 2021). The dataset obtained from PBDB contained initially 5808 occurrences for a period ranging from the Asselian to the Rhaetian. The dataset was cleaned of synonyms, outdated combinations, nomina dubia, and other erroneous and doubtful records, based on revisions provided in the literature and/or on the expertise of the authors. After correction, including data addition from the literature, our dataset was composed of 3636 species (1784 genera, and 418 families) for 17,250 occurrences resulting from an in-depth study and curation of the entire bibliography of fossil insects, spanning from the Asselian (lowermost Permian) to the Rhaetian (uppermost Triassic). Although most of the taxa included in the datasets are nominal taxa (published and named), a few unnamed taxa (genera or species) that are considered separate from others were also included, although not formally named in the literature or not published yet. These unpublished taxa are identifiable by the notation ‘fam. nov.’ or ‘gen. nov.’ following their names.Occurrences used here are specimens originating from a given stratigraphic horizon assigned to a given taxon. The age of each occurrence is based on data from PBDB, corrected with a more precise age (generally stage, sometimes substage), and the age of each time bin boundaries relies on the stratigraphic framework proposed in the International Chronostratigraphic Chart (updated to correspond with the ICS 2022/0295). Similarly, the ages of some species assigned to the wrong stage were corrected. In fact, some species from the French Permian deposit of Lodève were initially considered to be of Artinskian age in PBDB but most species from this deposit originate from the Merifons member, which is of Kungurian age96.Our data compilation allows a robust integration of data before and after our period of interest (i.e. the lower Permian and all geologic stages after the Carnian) to encompass occurrences of genera that may survive until the Late Triassic and to generate a sufficient background for the model to correctly estimate the extinction events around the P/T boundary. Since we used different datasets, the differences between genus-level or family-level occurrence numbers are explained by the systematic placement of some specimens that can only be placed confidently in a family but not in a genus (Supplementary Table 1). Tentative species identifications originally placed with uncertainty (reported as ‘aff.’ or ‘?’) were always included at a higher taxonomic level. Uncertain generic attributions were integrated as occurrences at the family level (e.g. a fossil initially considered Tupus? is recorded as an occurrence of Meganeuridae). Our total dataset was subdivided into smaller datasets, which represent orders or other subclades of insects (e.g. Mecoptera, Holometabola and Polyneoptera). Note that all the ichnospecies—a species name assigned to trace fossils (e.g. resting trace, nest and leaf damage)— and insect eggs (e.g. Clavapartus latus, Furcapartus exilis and Monilipartus tenuis) were not included in the analyses97. To prevent potential issues regarding the diversification estimates for clades with poor delineation, we refrained from analysing several orders that serve as taxonomic ‘wastebaskets’ (e.g. Grylloblattodea). These groups are poorly defined, likely polyphyletic or paraphyletic, and not supported by apomorphic characters—e.g. the monophyly of the ‘Grylloblattodea’ (Grylloblattida Walker, 1914 plus numerous fossil families and genera of uncertain affinities) is not supported by any synapomorphy, nor the relationships within this group. The occurrences assigned to these orders were rather included in analyses conducted at a higher taxonomic level (at the Polyneoptera level in the case of the ‘Grylloblattodea’). The detail of the composition of all the datasets is given in Supplementary Table 14, and each dataset is available in Supplementary Data 1.Studying extinction should, when possible, rely on species-level diversity to better circumscribe extinction events at this taxonomic rank, which is primarily affected by extinction98,99,100. However, in palaeoentomology, species-level occurrence data may contain less information than genus-level data, mainly because species are most of the time only known from one deposit, resulting in reduced life span, and are also sometimes poorly defined. Insects are also less prone to long-lasting genera or species than other lineages, maybe because of the relatively short time between generations (allowing for rapid evolution) or because morphological characters are better preserved or more diagnostic than in other lineages (i.e. wing venation), allowing easier differentiation. Another argument for the use of genus-level datasets is the possibility to add occurrences represented by fossils that cannot be assigned at the species level because of poor preservation or an insufficient number of specimens/available characters. By extension, the genus life span provides clues as to survivor taxa and times of origination during periods of post-extinction or recovery. A genus encompassing extinction events indicates that at least one species of this genus crossed the extinction. To get the best signal and infer a robust pattern of insect dynamics around the P/T events, we have chosen to analyse our dataset at different taxonomic ranks (e.g. genus, family and order levels) to extract as much evidence as possible.To further support our choice to work at these different levels, most recent works aiming to decipher the diversification and extinction in insect lineages have worked using a combination of analyses21,22,26; this also applies to non-insect clades51,101,102. This multi-level approach should maximise our understanding of the Permo–Triassic events.Assessing optimal parameters and preliminary testsPrior to choosing the settings for the final analyses (see detail in Dynamics of origination and extinction), a series of tests were carried out to better evaluate the convergence of our analyses. First, we analysed our genus-level dataset with PyRate36 running for 10 million generations and sampling every 10,000 generations, on ten randomly replicated datasets using the reversible-jump Markov Chain Monte Carlo (RJMCMC) model37 and the parameters of PyRate set by default. As the convergence was too low, new settings were used, notably increasing the number of generations to 50 million generations and monitoring the MCMC mixing and effective sample size (ESS) each 10 million generations. We modified the minimal interval between two shifts (-min_dt option, testing 0.5, 1.5 and 2), and found no major difference in diversification patterns between our tests. We have opted for 50 million generations with a predefined time frame set for bins corresponding to the Permian and Triassic stages, and a minimum interval between two shifts of two Ma. These parameters allow for maintaining a short bin frame and high convergence values while correctly identifying periods of diversification and extinction. For each analysis, ten datasets were generated using the extract.ages function to randomly resample the age of fossil occurrences within their respective temporal ranges (i.e. resampled ages are randomly drawn between the minimum and the maximum ages of the geological stratum). We monitored chain mixing and ESS by examining the log files in Tracer 1.7.1103 after excluding the first 10% of the samples as a burn-in period. The parameters are considered convergent when their ESS are greater than 200.Dynamics of origination and extinctionWe carried out the analyses of the fossil datasets based on the Bayesian framework implemented in the programme PyRate36. We analysed the fossil datasets under two models: the birth–death model with constrained shifts (BDCS38) and the RJMCMC (-A 4 option37). These models allow for a simultaneous estimate for each taxon: (1) the parameters of the preservation process (Supplementary Fig. 17), (2) the times of origination (Ts) and extinction (Te) of each taxon, (3) the origination and extinction rates and their variation through time for each stage and (4) the number and magnitude of shifts in origination and extinction rates.All analyses were set with the best-fit preservation process after comparing (-PPmodeltest option) the homogeneous Poisson process (-mHPP option), the non-homogeneous Poisson process (default option), and the time-variable Poisson process (-qShift option). The preservation process infers the individual origination and extinction times of each taxon based on all fossil occurrences and on an estimated preservation rate, denoted q, expressed as expected occurrences per taxon per Ma. The time-variable Poisson process assumes that preservation rates are constant within a predefined time frame but may vary over time (here, set for bins corresponding to stages). This model is thus appropriate when rates over time are heterogeneous.We ran PyRate for 50 million MCMC generations and a sampling every 50,000 generations for the BDCS and RJMCMC models with time bins corresponding to Permian and Triassic stages (-fixShift option). All analyses were set with a time-variable Poisson process (-qShift option) of preservation and accounted for varying preservation rates across taxa using the Gamma model (-mG option), that is, with gamma-distributed rate heterogeneity with four rate categories36. As explained above, the minimal interval between two shifts (-min_dt option) was modified and a value of 2 was used. The default prior to the vector of preservation rates is a single gamma distribution with shape = 1.5 and rate = 1.5. We reduced the subjectivity of this parameter, and favoured a better adequation to the data, allowing PyRate to estimate the rate parameter of the prior from the data by setting the rate parameter to 0 (-pP option). Therefore, PyRate assigns a vague exponential hyper-prior to the rate and samples the rate along with all other model parameters. Similarly, because our dataset does not encompass the entire fossil record of insects, we assumed that a possible edge effect may interfere with our analyses, with a strong diversification during the lowermost Permian and, conversely a strong extinction during the uppermost Triassic. Because the RJMCMC and BDCS algorithms look for rate shifts, we constrained the algorithm to only search for shifts (-edgeShift option) within the following time range 295.0 to 204.5 Ma. We monitored chain mixing and ESS by examining the log files in Tracer 1.7.1103 after excluding the first 10% of the samples as a burn-in period. The parameters are considered convergent when their ESS are greater than 200.We then combined the posterior estimates of the origination and extinction rates across all replicates to generate rates through-time plots (origination, extinction, and net diversification). Shifts of diversification were considered significant when log Bayes factors were >6 in the RJMCMC model, while we considered shifts to be significant in the BDCS model when mean rates in a time bin did not overlap with the 95% credibility interval (CI) of the rates of adjacent time bins.We replicated all the analyses on ten randomly generated datasets of each clade and calculated estimates of the Ts and the Te as the average of the posterior samples from each replicate. Thus, we obtained ten posterior estimates of the Ts and Te for all taxa and we used these values to estimate the past diversity dynamics by calculating the number of living taxa at each time point. For all the subsequent analyses, we used the estimated Ts and Te of all taxa to test whether or not the origination and the extinction rate dynamics were correlated with particular abiotic factors, as suggested by the drastic changes in environmental conditions known during the Permo–Triassic. We used proxies for abiotic factors, such as global continental fragmentation or the dynamic of major clades of plants, and for biotic factors via species interaction within and between ecological guilds. This approach avoids re-modelling preservation and re-estimating times of origination and extinction, which reduces drastically the computational burden, while still allowing to account for the preservation process and the uncertainties associated with fossil ages. Similarly, the times of origination and extinction used in all the subsequent analyses were obtained while accounting for the heterogeneity of preservation, origination and extinction rates. To discuss the magnitude of the periods of extinction and diversification, we compared the magnitude of these events to the background origination and extinction rates (i.e. not during extinction or diversification peaks).The PyRate approach has proven to be robust following a series of tests and simulations that reflect commonly observed biases when modelling past diversity dynamics31,38. These simulations were based on datasets simulated under a range of potential biases (i.e. violations of the sampling assumptions, variable preservation rates, and incomplete taxon sampling) and reflecting the limitations of the fossil record. Simulation results showed that PyRate is able to correctly estimate the dynamics of origination and extinction rates, including sudden rate changes and mass extinction, even if the preservation levels are low (down to 1–3 fossil occurrences per species on average), the taxon sampling is partial (up to 80% missing) or if the datasets have a high proportion of singletons (exceeding 30% of the taxa in some cases). The strongest bias in birth–death rate estimates is caused by incomplete data (i.e. missing lineages) altering the distribution of taxa; a pervasive effect often mentioned for phylogeny-based models104,105,106. However, in the case of PyRate, the simulations confirm the absence of consistent biases due to an incomplete fossil record36. Finally, the recently implemented RJMCMC model was shown to be very accurate for estimating origination and extinction rates (i.e. more accurate than the BDCS model, the boundary-crossing and three-time methods) and is able to recover sudden extinction events regardless of the biases in the fossil dataset37.The severity of extinctions and survivorsFor each event—the Roadian–Wordian, the LPME, and the Ladinian–Carnian—we quantified the percentage of extinctions and survivors at the genus level. We used the Te and Ts from our RJMCMC analysis and computed the mean for the Te (Tem) and for the Ts (Tsm) of each genus. We then filtered our dataset to keep only the genera with a Tsm older than the upper boundary of the focal event, i.e., we only kept the genera that appeared before the end of the event. Then, we discarded the genera that have disappeared before the lower boundary of the focal event, i.e. Tem older that the lower boundary of the event. The remaining genera, which corresponds to all the genera (total) present during the crisis (Ttgen), can be classified into two categories, ‘survivor genera’ (Sgen), i.e. those that survived the crisis, and those that died: ‘extinct genera’ (Egen). The survivors have a Tem younger than the upper boundary of the focal event, while the ‘extinct genera’ died out during the event and have a Tem between the lower and upper boundaries of the event of interest. To obtain the percentage of survivors, we used the following formula: (Sgen/Ttgen) × 100. Similarly, the percentage of extinction is calculated as: (Egen/Ttgen) × 100.Age-dependent extinction modelWe assessed the effect of taxon age on the extinction probability by fitting the age-dependent extinction (ADE; -ADE 1 option) model50. This model estimates the probability for a lineage to become extinct as a function of its age, also named longevity, which is the elapsed time since its origination. It is recommended to run the ADE model over time windows with roughly constant origination and extinction rates, as convergence is difficult—but not impossible—to reach in extinction or diversification contexts50. We ran PyRate for 50 million MCMC generations with a sampling every 50,000 generations, with a time-variable Poisson process of preservation (-qShift option), while accounting for varying preservation rates across taxa using the Gamma model (-mG option). We replicated the analyses on ten randomised datasets and combined the posterior estimates across all replicates. We estimated the shape (Φ) and scale (Ψ) parameters of the Weibull distribution, and the taxon longevity in a million years. According to ref. 50, there is no evidence of age-dependent extinction rates if Φ = 1. However, the extinction rate is higher for young species and decreases with species age if Φ  1. Although ADE models are prone to high error rates when origination and extinction rates increase or decrease through time, simulations with PyRate have shown that fossil-based inferences are robust50. We investigated the effect of ADE during three different periods (-filter option) as follows: (1) between 264.28 Ma and 255 Ma (pre-decline), (2) between 254.5 Ma and 251.5 Ma (decline) and (3) between 234 Ma and 212 Ma (post-crisis). We monitored chain mixing and ESS by examining the log files in Tracer 1.7.1103 and considered the convergence of parameters sufficient when their ESS were greater than 200.Selection of abiotic and biotic variablesTo test correlations of insect diversification with environmental changes, we examined the link between a series of environmental variables and origination/extinction rates over a period encompassing the GEE, the LPME and the CPE but also for each extinction event. We focused on the role of nine variables, also called proxies, which have been demonstrated or assumed to be linked to extinctions and changes in insect diversity26,67.The variations in the atmospheric CO2 and O2 concentrations are thought to be correlated with the diversification of several insect lineages, including the charismatic giant Meganeuridae65,66,67. Because the increase of O2 concentration has likely driven the diversification of some insects, its diminution may have resulted in the extinction or decline of some lineages. Therefore, we investigated the potential correlation of the variations of this variable with insect dynamics using data from ref. 55. We extracted the data, with 1-million-year time intervals, spanning the Permo–Triassic.Similarly, the modification of CO2 concentration, notably its increase, is known to promote speciation in some modern insect groups107. Therefore, a similar effect may have occurred during the Permian and Triassic but remains to be tested. We based our analyses on the dataset of ref. 108. We used their cleaned dataset and extracted all verified values for the Permo–Triassic interval. Because the initial data (i.e. independent estimates) were made in various locations for the same age, different values of the CO2 concentration are provided. We incorporated all these values in our analysis, allowing PyRate to search for a correlation for each value of the CO2 concentration. We obtain a final correlation independent of the sampling location, in line with our large-scale analysis.The continental fragmentation, as approximated by plate tectonic change over time, has recently been proposed as a driver of Plecoptera dynamics26. Because the period studied encompasses a major geological event, the fragmentation of the supercontinent Pangea, we investigated the effect of continental fragmentation on insect diversification dynamics. We retrieved the index of continental fragmentation developed by ref. 69 using paleogeographic reconstructions for 1-million-year time intervals. This index approaches 1 when all plates are disjoined (complete plate fragmentation) and approaches 0 when the continental aggregation is maximal.Climate change (variations in warming and cooling periods) is a probable driver of diversification changes over the history of insects21,109. Temperature is likely directly linked with insect dynamics109 but also with their food sources, notably plants110. Because it was demonstrated that modification of temperature impacted floral assemblages110, we tested the correlation between temperature variations and the diversification dynamic of insects. Major trends in global climate change through time are typically estimated from relative proportions of different oxygen isotopes (δ18O) in samples of benthic foraminiferan shells111. We used the data from ref. 112, converted to absolute temperatures following the methodology described in Condamine et al.113 (see their section Global temperature variations through time). The resulting temperature data reflects planetary-scale climatic trends, with time intervals inferior to 1-million-year, which can be expected to have led to temporally coordinated diversification changes in several clades rather than local or seasonal fluctuations.The fluctuation in relative diversity of gymnosperms, non-Polypodiales ferns, Polypodiales ferns, spore-plants, and later the rise of angiosperms has likely driven the diversification of numerous insects57,60,61,114. Close interactions between insects and plants are well-recorded during the Permian and Triassic57,60,61. In fact, herbivorous insects are known to experience high selection pressure from bottom-up forces, resulting from interactions with their hosts or feeding plants30,72. Therefore, it appears crucial to investigate the effect of these modifications on the insects’ past dynamics. We used the data from ref. 38 for the different plant lineages (all with 1-million-year time intervals). All the datasets for these variables are available in the publications cited aside from each variable or in Supplementary Data 1.Multivariate birth–death modelWe used the multivariate birth–death (MBD) model to assess to what extent biotic and abiotic factors can explain temporal variation in origination and extinction rates55. The model is described in ref. 55, where origination and extinction rates can change through time in relation to environmental variables so that origination and extinction rates depend on the temporal variations of each factor. The strength and sign (positive or negative) of the correlations are jointly estimated for each variable. The sign of the correlation parameters indicates the sign of the resulting correlation. When their value is estimated around zero, no correlation is estimated. An MCMC algorithm combined with a horseshoe prior, controlling for over-parameterisation and for the potential effects of multiple testing, jointly estimates the baseline origination (λ0) and extinction (µ0) rates and all correlation parameters (Gλ and Gµ)55. The horseshoe prior is used to discriminate which correlation parameters should be treated as noise (shrunk around 0) and which represent a true signal (i.e. significantly different from 0). In the MBD model, a correlation parameter is estimated to quantify independently the role of each variable on the origination and the extinction.We ran the MBD model using 20 (for short intervals) or 50 million MCMC generations and sampling every 20,000 or 50,000 to approximate the posterior distribution of all parameters (λ0, µ0, nine Gλ, nine Gµ and the shrinkage weights of each correlation parameter, ωG). The MBD analyses used the Ts and the Te derived from our previous analyses under the RJMCMC model. The results of the MBD analyses were summarised by calculating the posterior mean and 95% CI of all correlation parameters and the mean of the respective shrinkage weights (across ten replicates), as well as the mean and 95% CI of the baseline origination and extinction rates. We carried out six analyses, over: (1) the Permo–Triassic (between 298.9 and 201.3 Ma); (2) the Roadian–Wordian (R/W) boundary (between 270 and 265 Ma), (3) the LPME (between 254.5 and 250 Ma), (4) the Ladinian–Carnian (L/C) boundary (between 240 to 234 Ma), (5) the Permian period (between 298.9 and 251.902 Ma) and (6) the Triassic period (between 251.902 and 201.3 Ma). We monitored chain mixing and ESS by examining the log files in Tracer 1.7.1103 and considered the convergence of parameters sufficient when their ESS were greater than 200.Multiple clade diversity-dependence modelTo assess the potential effect of diversity-dependence on the diversity dynamics of three or four insect guilds, we used the multiple clade diversity-dependence (MCDD) model in which origination and extinction rates are correlated with the diversity trajectory of other clades31. This model postulates that competitive interactions linked with an increase in diversity results in decreasing origination rates and/or increasing extinction rates. The MCDD model allows for testing diversity-dependence between genera of a given clade or between genera of distinct clades sharing a similar ecology.We estimated the past diversity dynamics for three (i.e. herbivores, predators, and a guild composed of generalists + detritivores/fungivores dubbed ‘others’) or four insect groups or guilds (i.e. herbivores, predators, generalists and detritivores/fungivores) by calculating the number of living species at every point in time based on the times of origination (Ts) and extinction (Te) estimated under the RJMCMC model (see above) (Supplementary Figs. 19–24). We defined our four insect groups with a cautious approach i.e. insect genera, families or orders for which nothing is known about the ecology or about the ecology of their close relatives were not considered for the analysis. For example, no diet was assigned to Diptera, Mecoptera or Glosselytrodea. The ecology of the Triassic Diptera and Permo–Triassic Mecoptera is difficult to establish because extant Diptera and Mecoptera have a wide diversity of ecology. Fossil Mecoptera are also putatively involved in numerous interactions with plants (species with elongated mouthparts), suggesting a placement in the herbivore group, while other species were likely predators. Therefore, we cannot decide to which group each species belongs. Similarly, nothing is known about the body and mouthparts of the Glosselytrodea, most of the time described based on isolated wings; we did not assign the order to any group. The definition and delineation of insect clades have also challenged the placement of several orders (e.g. ‘Grylloblattodea’) in one of our four groups. The order ‘Grylloblattodea’ is poorly delineated and mostly serves as a taxonomic ‘wastebasket’ to which it is impossible to assign a particular ecology. Finally, genera, species, or families not placed in a higher clade (e.g. Meshemipteron, Perielytridae) were not included in the analysis. Oppositely, the guilds ‘herbivores’ and ‘predators’ are well defined, and their ecology is evidenced by the morphology of their representatives and the principle of actualism. For example, the ecology of Meganeurites gracilipes (Meganeuridae) has been deeply studied, and its enlarged compound eyes, its sturdy mandibles with acute teeth, its tarsi and tibiae bearing strong spines, and the presence of a pronounced thoracic skewness are specialisations today found in dragonflies that capture their prey while in flight115. All Odonatoptera are well-known predator insects. The raptorial forelegs of the representatives of the order Titanoptera and their mouthparts with strong mandibles are linked with predatory habits81. The Palaeodictyopteroidea were herbivorous insects with long, beak-like, piercing mouthparts, and probably a sucking organ81,82. Most Hemiptera are confidently considered herbivorous insects by comparison with their extant representatives. For example, the Cicadomorpha or Sternorrhyncha are known to feed on plants and their fossil representatives likely possessed the same ecology because of similar morphologies116. Some hemipteran families (e.g. Nabidae) are predators and we cautiously distinguished herbivorous and carnivorous taxa among Hemiptera. The detail of the ecological assignations for the 1009 genera included in our analyses can be found in Supplementary Data 1 (Table MCCD).We calculated ten diversity trajectories from the ten replicated analyses under the RJMCMC model. The estimation of past species diversity might be biased by low preservation rates or taxonomic uncertainties. However, such trajectory curves are likely to provide a reasonably accurate representation of the past diversity changes in the studied clades, notably because the preservation during the Permian and Triassic period is relatively good for insects (i.e. no gaps).Our MCDD analyses comprise all the insect genera spanning from the lowermost Permian to the uppermost Triassic and were run and repeated on ten replicates (using the Te and Ts estimated under the RJMCMC model) with 50 million MCMC generations and a sampling frequency of 50,000. For each of the four insect groups, we computed the median and the 95% CI of the baseline origination and extinction rates (λi and µi), the within-group diversity-dependence parameters gλi and gµi, and the between-groups diversity dependence parameters gλij and gµij. The mean of the sampled diversity dependence parameters (e.g. gλij) was used as a measure of the intensity of the negative (if positive) or positive interactions (if negative) between each pair of groups. The interactions were considered significant when their median was different from 0 and the 95% CI did not overlap with 0. We monitored chain mixing and ESS by examining the log files in Tracer 1.7.1103 and considered the convergence of parameters sufficient when their ESS were greater than 200.We cross-validated the result of the MCDD model using the MBD model. The MBD model can be used to run a multiple clade diversity-dependence analysis by providing the diversity trajectories of insect guilds as a continuous variable. These data are directly generated by PyRate using the lineages-through-time generated by the RJMCMC analyses (-ltt option). We ran the MBD model using 50 million MCMC generations and sampling every 50,000 to approximate the posterior distribution of all parameters (λ0, µ0, four Gλ, four Gµ and the shrinkage weights of each correlation parameter, ωG). We carried out three analyses, over the period encompassing the three extinction events (between 275 and 230 Ma): (1) for herbivores; (2) for predators; and (3) for ‘others’. For each analysis, the lineages-through-time data of the two other guilds are used as continuous variables to investigate a diversity dependence effect. We monitored chain mixing and ESS by examining the log files in Tracer 1.7.1103 and considered the convergence of parameters sufficient when their ESS were greater than 200.Reporting summaryFurther information on research design is available in the Nature Portfolio Reporting Summary linked to this article. More