Philippa Kaur

Coyne, J. A. & Orr, H. A. Speciation 83–178 (Sinauer, 2004).Dieckmann, U., Doebeli, M., Metz, J. A. & Tautz, D. Adaptive Speciation (Cambridge University Press, 2004).Butlin, R. K., Galindo, J. & Grahame, J. W. Sympatric, parapatric or allopatric: the most important way to classify speciation? Philos. T. Roy. Soc. B 363, 2997–3007 (2008).Article 

Google Scholar 
Smadja, C. M. & Butlin, R. K. A framework for comparing processes of speciation in the presence of gene flow. Mol. Ecol. 20, 5123–5140 (2011).Article 

Google Scholar 
Seehausen, O. et al. Genomics and the origin of species. Nat. Rev. Genet. 15, 176–192 (2014).Article 
CAS 

Google Scholar 
Kulmuni, J., Butlin, R. K., Lucek, K., Savolainen, V. & Westram, A. M. Towards the completion of speciation: the evolution of reproductive isolation beyond the first barriers. Philos. T. Roy. Soc. B 375, 20190528 (2020).Article 

Google Scholar 
Hofreiter, M. & Stewart, J. Ecological change, range fluctuations and population dynamics during the Pleistocene. Curr. Biol. 19, R584–R594 (2009).Article 
CAS 

Google Scholar 
Longman, J., Mills, B. J. W., Manners, H. R., Gernon, T. M. & Palmer, M. R. Late Ordovician climate change and extinctions driven by elevated volcanic nutrient supply. Nat. Geosci. 14, 924–929 (2021).Article 
ADS 
CAS 

Google Scholar 
Thomson, R. C., Spink, P. Q. & Shaffer, H. B. A global phylogeny of turtles reveals a burst of climate-associated diversification on continental margins. Proc. Natl Acad. Sci. USA 118, e2012215118 (2021).Article 
CAS 

Google Scholar 
Chaboureau, A. C., Sepulchre, P., Donnadieu, Y. & Franc, A. Tectonic-driven climate change and the diversification of angiosperms. Proc. Natl Acad. Sci. USA 111, 14066–14070 (2014).Article 
ADS 
CAS 

Google Scholar 
Hewitt, G. The genetic legacy of the Quaternary ice ages. Nature 405, 907–913 (2000).Article 
ADS 
CAS 

Google Scholar 
Schmitt, T. Molecular biogeography of Europe: Pleistocene cycles and postglacial trends. Front. Zool. 4, 11 (2007).Article 

Google Scholar 
Haffer, J. Speciation in Amazonian forest birds. Science 165, 131–137 (1969).Article 
ADS 
CAS 

Google Scholar 
Ebdon, S. et al. The Pleistocene species pump past its prime: evidence from European butterfly sister species. Mol. Ecol. 30, 3575–3589 (2021).Article 

Google Scholar 
Excoffier, L., Foll, M. & Petit, R. J. Genetic consequences of range expansions. Annu. Rev. Ecol. Evol. Syst. 40, 481–501 (2009).Article 

Google Scholar 
Baker, H. G. Self-compatibility and establishment after ‘long-distance’ dispersal. Evolution 9, 347–349 (1955).
Google Scholar 
Fisher, R. The Genetical Theory of Natural Selection 125–129 (Oxford University Press, 1930).Endler, J. A. Geographic Variation, Speciation, and Clines. Monographs in Population Biology Vol. 10, 53–65, 142–150 (Princeton University Press, 1977).Doebeli, M. & Dieckmann, U. Speciation along environmental gradients. Nature 421, 259–264 (2003).Article 
ADS 
CAS 

Google Scholar 
Ispolatov, J. & Doebeli, M. Diversification along environmental gradients in spatially structured populations. Evol. Ecol. Res. 11, 295–304 (2009).
Google Scholar 
Rettelbach, A., Servedio, M. R. & Hermisson, J. Speciation in peripheral populations: effects of drift load and mating systems. J. Evol. Biol. 29, 1073–1090 (2016).Article 
CAS 

Google Scholar 
Wright, S. I., Kalisz, S. & Slotte, T. Evolutionary consequences of self-fertilization in plants. Proc. R. Soc. Lond. Ser. B 280, 20130133 (2013).
Google Scholar 
Hu, X.-S. Mating system as a barrier to gene flow. Evolution 69, 1158–1177 (2015).Article 
CAS 

Google Scholar 
Glémin, S. How are deleterious mutations purged? Drift versus nonrandom mating. Evolution 57, 2678–2687 (2003).
Google Scholar 
Warwick, S. I., Francis, A. & Al-Shehbaz, I. A. Brassicaceae: species checklist and database on CD-Rom. Plant Syst. Evol. 259, 249–258 (2006).Article 

Google Scholar 
Warwick, S. I., Al-Shehbaz, I. A. & Sauder, C. A. Phylogenetic position of Arabis arenicola and generic limits of Aphragmus and Eutrema (Brassicaceae) based on sequences of nuclear ribosomal DNA. Can. J. Bot. 84, 269–281 (2006).Article 
CAS 

Google Scholar 
Hohmann, N. et al. Taming the wild: resolving the gene pools of non-model Arabidopsis lineages. BMC Evol. Biol. 14, e224 (2014).Article 

Google Scholar 
Novikova, P. Y. et al. Sequencing of the genus Arabidopsis identifies a complex history of nonbifurcating speciation and abundant trans-specific polymorphism. Nat. Genet. 48, 1077–1082 (2016).Article 
CAS 

Google Scholar 
Perrier, A. & Willi, Y. Intraspecific variation in reproductive barriers between two recently-diverged, allopatric Arabidopsis species. J. Evol. Biol. https://doi.org/10.1111/jeb.14122 (2022). (in press).Griffin, P. C. & Willi, Y. Evolutionary shifts to self-fertilisation restricted to geographic range margins in North American Arabidopsis lyrata. Ecol. Lett. 17, 484–490 (2014).Article 
CAS 

Google Scholar 
Willi, Y., Fracassetti, M., Zoller, S. & Van Buskirk, J. Accumulation of mutational load at the edges of a species range. Mol. Biol. Evol. 35, 781–791 (2018).Article 
CAS 

Google Scholar 
Schmickl, R., Jørgensen, M. H., Brysting, A. K. & Koch, M. A. The evolutionary history of the Arabidopsis lyrata complex: a hybrid in the Amphi-Beringian area closes a large distribution gap and builds up a genetic barrier. BMC Evol. Biol. 10, e98 (2010).Article 

Google Scholar 
Pyhäjärvi, T., Aalto, E. & Savolainen, O. Time scales of divergence and speciation among natural populations and subspecies of Arabidopsis lyrata (Brassicaceae). Am. J. Bot. 99, 1314–1322 (2012).Article 

Google Scholar 
Dyke, A. S. in Quaternary Glaciations – Extent and Chronology, Part II: North America (Elsevier, Amsterdam, 2004).Kirkpatrick, M. & Ravigné, V. Speciation by natural and sexual selection: models and experiments. Am. Nat. 159, S22–S35 (2002).Article 

Google Scholar 
Igic, B., Lande, R. & Kohn, J. R. Loss of self‐incompatibility and its evolutionary consequences. Int. J. Plant Sci. 169, 93–104 (2008).Article 

Google Scholar 
Willi, Y. & Määttänen, K. Evolutionary dynamics of mating system shifts in Arabidopsis lyrata. J. Evol. Biol. 23, 2123–2131 (2010).Article 
CAS 

Google Scholar 
Lucek, K. & Willi, Y. Drivers of linkage disequilibrium across a species’ geographic range. PLoS Genet. 17, e1009477 (2021).Article 
CAS 

Google Scholar 
Pironon, S. et al. Geographic variation in genetic and demographic performance: new insights from an old biogeographical paradigm: the centre-periphery hypothesis. Biol. Rev. 92, 1877–1909 (2017).Article 

Google Scholar 
Encinas-Viso, F., Young, A. G. & Pannell, J. R. The loss of self-incompatibility in a range expansion. J. Evol. Biol. 33, 1235–1244 (2020).Article 

Google Scholar 
Jarne, P. & Auld, J. R. Animals mix it up too: the distribution of self-fertilization among hermaphroditic animals. Evolution 60, 1816–1824 (2006).
Google Scholar 
Foxe, J. P. et al. Reconstructing origins of loss of self-incompatibility and selfing in North American Arabidopsis lyrata: a population genetic context. Evolution 64, 3495–3510 (2010).Article 

Google Scholar 
Koski, M. H., Layman, N. C., Prior, C. J., Busch, J. W. & Galloway, L. F. Selfing ability and drift load evolve with range expansion. Evol. Lett. 3, 500–512 (2019).Article 

Google Scholar 
Prior, C. J. & Busch, J. W. Selfing rate variation within species is unrelated to life‐history traits or geographic range position. Am. J. Bot. 108, 2294–2308 (2021).Article 

Google Scholar 
Skeels, A. & Cardillo, M. Reconstructing the geography of speciation from contemporary biodiversity data. Am. Nat. 193, 240–254 (2019).Article 

Google Scholar 
Sánchez-Castro, D., Perrier, A. & Willi, Y. Reduced climate adaptation at range edges in North American Arabidopsis lyrata. Glob. Ecol. Biogeogr. 31, 1066–1077 (2022).Article 

Google Scholar 
Roessler, K. et al. The genome-wide dynamics of purging during selfing in maize. Nat. Plants 5, 980–990 (2019).Article 
CAS 

Google Scholar 
Li, H. Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM. Preprint at https://arxiv.org/abs/1303.3997 (2013).Hu, T. T. et al. The Arabidopsis lyrata genome sequence and the basis of rapid genome size change. Nat. Genet. 43, 476–481 (2011).Article 

Google Scholar 
Li, H. et al. The sequence alignment/map format and SAMtools. Bioinformatics 25, 2078–2079 (2009).Article 

Google Scholar 
McKenna, A. et al. The genome analysis toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data. Genome Res. 20, 1297–1303 (2010).Article 
CAS 

Google Scholar 
Cingolani, P. et al. A program for annotating and predicting the effects of single nucleotide polymorphisms, SnpEff: SNPs in the genome of Drosophila melanogaster strain w1118; iso-2; iso-3. Fly 6, 80–92 (2012).Article 
CAS 

Google Scholar 
Alexander, D. H., Novembre, J. & Lange, K. Fast model-based estimation of ancestry in unrelated individuals. Genome Res. 19, 1655–1664 (2009).Article 
CAS 

Google Scholar 
Danecek, P. et al. The variant call format and VCFtools. Bioinformatics 27, 2156–2158 (2011).Article 
CAS 

Google Scholar 
Stamatakis, A. RAxML version 8: a tool for phylogenetic analysis and post-analysis of large phylogenies. Bioinformatics 30, 1312–1313 (2014).Article 
CAS 

Google Scholar 
Pickrell, J. K. & Pritchard, J. K. Inference of population splits and mixtures from genome-wide allele frequency data. PLoS Genet. 8, e1002967 (2012).Article 
CAS 

Google Scholar 
Excoffier, L., Dupanloup, I., Huerta-Sánchez, E., Sousa, V. C. & Foll, M. Robust demographic inference from genomic and SNP data. PLoS Genet. 9, e1003905 (2013).Article 

Google Scholar 
Marchi, N. et al. The genomic origins of the world’s first farmers. Cell 185, 1842–1859 (2022).Article 
CAS 

Google Scholar 
Li, H. & Durbin, R. Inference of human population history from individual whole-genome sequences. Nature 475, 493–496 (2011).Article 
CAS 

Google Scholar 
Genete, M., Castric, V. & Vekemans, X. Genotyping and de novo discovery of allelic variants at the Brassicaceae self-incompatibility locus from short-read sequencing data. Mol. Biol. Evol. 7, 1193–1201 (2020).Article 

Google Scholar 
Lynch, M. et al. Genome-wide linkage-disequilibrium profiles from single individuals. Genetics 198, 269–281 (2014).Article 

Google Scholar 
R Core Team. R: a language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria (2021).Paradis, E. & Schliep, K. ape 5.0: an environment for modern phylogenetics and evolutionary analyses in R. Bioinformatics 35, 526–528 (2019).Article 
CAS 

Google Scholar 
Nychka, D., Furrer, R., Paige, J. & Sain, S. fields: tools for spatial data. R package version 14.1 https://github.com/dnychka/fieldsRPackage (2021).Asquith, W. lmomco—L-moments, censored L-moments, trimmed L-moments, L-comoments, and many distributions. R package version 2.4.7 (2022).Revell, L. J. phytools: an R package for phylogenetic comparative biology (and other things). Methods Ecol. Evol. 3, 217–223 (2012).Article 

Google Scholar 
Lemon, J. Plotrix: a package in the red light district of R. R. N. 6, 8–12 (2006).
Google Scholar 
Pebesma, E. J. & Bivand, R. S. Classes and methods for spatial data in R. R. N. 5, 9–13 (2005).
Google Scholar 
Bivand, R. S., Pebesma, E. & Gomez-Rubio, V. Applied Spatial Data Analysis with R Second edition (Springer, 2013). More

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

Lenoir, J. et al. Species better track climate warming in the oceans than on land. Nat. Ecol. Evol. 4(8), 1044–1059 (2020).Article 

Google Scholar 
Hobbs, R. J., Valentine, L. E., Standish, R. J. & Jackson, S. T. Movers and stayers: Novel assemblages in changing environments. Trends Ecol. Evol. 33, 116–128 (2017).Article 

Google Scholar 
Gilman, S. E., Urban, M. C., Tewksbury, J., Gilchrist, G. W. & Holt, R. D. A framework for community interactions under climate change. Trends Ecol. Evol. 25, 325–331 (2010).Article 

Google Scholar 
Ockendon, N. et al. Mechanisms underpinning climatic impacts on natural populations: Altered species interactions are more important than direct effects. Glob. Change Biol. 20, 2221–2229 (2014).Article 
ADS 

Google Scholar 
Gómez-Aparicio, L., García-Valdés, R., Ruíz-Benito, P. & Zavala, M. A. Disentangling the relative importance of climate, size and competition on tree growth in Iberian forests: Implications for forest management under global change. Glob. Change Biol. 17, 2400–2414 (2011).Article 
ADS 

Google Scholar 
Pecl, G. T. et al. Biodiversity redistribution under climate change: Impacts on ecosystems and human well-being. Science https://doi.org/10.1126/science.aai9214 (2017).Article 

Google Scholar 
Scheffers, B. R. et al. The broad footprint of climate change from genes to biomes to people. Science 354, aaf7671. https://doi.org/10.1126/science.aaf7671 (2016).Article 
CAS 

Google Scholar 
Vergés, A. et al. The tropicalization of temperate marine ecosystems: Climate-mediated changes in herbivory and community phase shifts. Proc. R. Soc. B-Biol. Sci. 281, 20140846. https://doi.org/10.1098/rspb.2014.0846 (2014).Article 

Google Scholar 
Poore, A. G. B. et al. Global patterns in the impact of marine herbivores on benthic primary producers. Ecol. Lett. 15, 912–922. https://doi.org/10.1111/j.1461-0248.2012.01804.x (2012).Article 

Google Scholar 
Bennett, S., Wernberg, T., Harvey, E. S., Santana-Garcon, J. & Saunders, B. J. Tropical herbivores provide resilience to a climate-mediated phase shift on temperate reefs. Ecol. Lett. 18, 714–723 (2015).Article 

Google Scholar 
Vergés, A. et al. Long-term empirical evidence of ocean warming leading to tropicalization of fish communities, increased herbivory and loss of kelp. Proc. Natl. Acad. Sci. 113(48), 13791–13796 (2016).Article 
ADS 

Google Scholar 
Vergés, A. et al. Tropical rabbitfish and the deforestation of a warming temperate sea. J. Ecol. 102, 1518–1527. https://doi.org/10.1111/1365-2745.12324 (2014).Article 

Google Scholar 
Kumagai, N. H. et al. Ocean currents and herbivory drive macroalgae-to-coral community shift under climate warming. Proc. Natl. Acad. Sci. 115, 8990–8995 (2018).Article 
ADS 
CAS 

Google Scholar 
Demko, A. M. et al. Declines in plant palatability from polar to tropical latitudes depend on herbivore and plant identity. Ecology 98, 2312–2321. https://doi.org/10.1002/ecy.1918 (2017).Article 

Google Scholar 
Floeter, S. R., Behrens, M. D., Ferreira, C. E. L., Paddack, M. J. & Horn, M. H. Geographical gradients of marine herbivorous fishes: Patterns and processes. Mar Biol 147, 1435–1447 (2005).Article 

Google Scholar 
Longo, G. O., Hay, M. E., Ferreira, C. E. L. & Floeter, S. R. Trophic interactions across 61 degrees of latitude in the Western Atlantic. Glob. Ecol. Biogeogr. 28, 107–117. https://doi.org/10.1111/geb.12806 (2019).Article 

Google Scholar 
Bolser, R. & Hay, M. Are tropical plants better defended? Palatability and defenses of temperate versus tropical seaweeds. Ecology 77, 2269–2286 (1996).Article 

Google Scholar 
Borer, E. T. et al. Global biogeography of autotroph chemistry: is insolation a driving force?. Oikos 122, 1121–1130. https://doi.org/10.1111/j.1600-0706.2013.00465.x (2013).Article 
CAS 

Google Scholar 
Miranda, T. et al. Convictfish on the move: Variation in growth and trophic niche space along a latitudinal gradient. ICES J. Mar. Sci. https://doi.org/10.1093/icesjms/fsz098%JICESJournalofMarineScience (2019).Article 

Google Scholar 
Linton, S. M. The structure and function of cellulase (endo-β-1, 4-glucanase) and hemicellulase (β-1, 3-glucanase and endo-β-1, 4-mannase) enzymes in invertebrates that consume materials ranging from microbes, algae to leaf litter. Comp. Biochem. Physiol. B Biochem. Mol. Biol. 240, 110354 (2020).Article 
CAS 

Google Scholar 
Poloczanska, E. S. et al. Global imprint of climate change on marine life. Nat. Clim. Change 3, 919–925. https://doi.org/10.1038/nclimate1958 (2013).Article 
ADS 

Google Scholar 
Nakamura, Y., Feary, D. A., Kanda, M. & Yamaoka, K. Tropical fishes dominate temperate reef fish communities within western Japan. PLoS ONE 8, e81107 (2013).Article 
ADS 

Google Scholar 
Tanaka, K., Taino, S., Haraguchi, H., Prendergast, G. & Hiraoka, M. Warming off southwestern Japan linked to distributional shifts of subtidal canopy-forming seaweeds. Ecol. Evol. 2, 2854–2865. https://doi.org/10.1002/ece3.391 (2012).Article 

Google Scholar 
Pessarrodona, A. et al. Homogenization and miniaturization of habitat structure in temperate marine forests. Glob. Change Biol. 27, 5262–5275 (2021).Article 
CAS 

Google Scholar 
Yamano, H., Sugihara, K. & Nomura, K. Rapid poleward range expansion of tropical reef corals in response to rising sea surface temperatures. Geophys. Res. Lett. 38, L04601. https://doi.org/10.1029/2010gl046474 (2011).Article 
ADS 

Google Scholar 
Mezaki, T. & Kubota, S. Changes of hermatypic coral community in coastal sea area of Kochi, high-latitude Japan. Aquabiology 201, 332–337 (2012).
Google Scholar 
Serisawa, Y., Imoto, Z., Ishikawa, T. & Ohno, M. Decline of the Ecklonia cava population associated with increased seawater temperatures in Tosa Bay, southern Japan. Fish Sci 70, 189–191. https://doi.org/10.1111/j.0919-9268.2004.00788.x (2004).Article 
CAS 

Google Scholar 
Kiriyama, T., Mitsunaga, N., Yasumoto, S., Fujii, A. & Yotsui, T. Undergrown phenomenon of brown alga, Hizikia fusiformis, thought to be caused by grazing of herbivores at Tsutsuura, Tsushima Islands [Japan]. Bulletin of Nagasaki Prefectural Institute of Fisheries (Japan) (1999).Kiriyama, T., Fujii, A. & Fujita, Y. Feeding and characteristic bite marks on Sargassum fusiforme by several herbivorous fishes. Aquac. Sci. 53, 355–365 (2005).
Google Scholar 
Yatsuya, K., Kiriyama, T., Kiyomoto, S., Taneda, T. & Yoshimura, T. On the deterioration process of Ecklonia and Eisenia beds observed in 2013 at Gounoura, Iki Island, Nagasaki Prefecture, Japan.-Initiation of the bed degradation due to high water temperature in summer and subsequent cascading effect by the grazing of herbivorous fish in autumn. Algal Resour. 7, 79–94 (2014).
Google Scholar 
Noda, M., Ohara, H., Murase, N., Ikeda, I. & Yamamoto, K. The grazing of Eisenia bicyclis and several species of Sargassaceous and Cystoseiraceous seaweeds by Siganus fuscescens in relation to the differences of species composition of their seaweed beds. Nippon Suisan Gakkaishi 80, 201–213 (2014).Article 

Google Scholar 
Noda, M., Kinoshita, J., Tanada, N. & Murase, N. Characteristics of bite scars observed in kelp forests of Lessoniaceae denuded by short-term foraging damages of the herbivorous fish Siganus fuscecens. J. Natl. Fish. Univ. 66, 111–122 (2018).
Google Scholar 
Wernberg, T. et al. Seaweed communities in retreat from ocean warming. Curr. Biol. 21, 1828–1832. https://doi.org/10.1016/j.cub.2011.09.028 (2011).Article 
CAS 

Google Scholar 
Terazono, Y., Nakamura, Y., Imoto, Z. & Hiraoka, M. Fish response to expanding tropical Sargassum beds on the temperate coasts of Japan. Mar. Ecol. Prog. Ser. 464, 209–220. https://doi.org/10.3354/meps09873 (2012).Article 
ADS 

Google Scholar 
Duffy, J. E. & Hay, M. E. Seaweed adaptations to herbivory – chemical, structural, and morphological defenses are often adjusted to spatial or temporal patterns of attack. Bioscience 40, 368–375 (1990).Article 

Google Scholar 
Endo, H., Suehiro, K., Kinoshita, J. & Agatsuma, Y. Combined effects of temperature and nutrient enrichment on palatability of the brown alga Sargassum yezoense (Yamada) Yoshida & T. Konno. Am. J. Plant Sci. 6, 275 (2015).Article 
CAS 

Google Scholar 
Clements, K. D., German, D. P., Piché, J., Tribollet, A. & Choat, J. H. Integrating ecological roles and trophic diversification on coral reefs: Multiple lines of evidence identify parrotfishes as microphages. Biol. J. Linn. Soc. 120, 729–751. https://doi.org/10.1111/bij.12914 (2017).Article 

Google Scholar 
Wang, Y., Naumann, U., Wright, S. T. & Warton, D. I. mvabund–an R package for model-based analysis of multivariate abundance data. Methods Ecol. Evol. 3, 471–474 (2012).Article 

Google Scholar 
Wilson, S. K., Bellwood, D. R., Choat, J. H. & Furnas, M. J. Detritus in the epilithic algal matrix and its use by coral reef fishes. Oceanogr. Mar. Biol. Annu. Rev. 41, 279–309 (2003).
Google Scholar 
Helfman, G. S. in The Behaviour of Teleost Fishes 366–387 (Springer, 1986).Prince, J., LeBlanc, W. & Maciá, S. Design and analysis of multiple choice feeding preference data. Oecologia 138, 1–4 (2004).Article 
ADS 

Google Scholar 
Hartig, F. DHARMa: residual diagnostics for hierarchical (multi-level/mixed) regression models. R package version 0.3 3 (2020).Ohno, M. & Ishikawa, M. Physiological ecology of brown alga, Ecklonia on coast of Tosa Bay, southern Japan. I. Seasonal variation of Ecklonia bed. Rep. USA Marine Biol. Inst. Kochi Univ. 4, 59–73 (1982).
Google Scholar 
Agostini, S. et al. Simplification, not “tropicalization”, of temperate marine ecosystems under ocean warming and acidification. Glob. Change Biol. 27, 4771–4784 (2021).Article 
CAS 

Google Scholar 
Clements, K. & Choat, J. Influence of season, ontogeny and tide on the diet of the temperate marine herbivorous fish Odax pullus (Odacidae). Mar. Biol. 117, 213–220 (1993).Article 

Google Scholar 
Mizuta, H., Hayasaki, J. & Yamamoto, H. Relationship between nitrogen content and sorus formation in the brown alga Laminaria japonica cultivated in southern Hokkaido, Japan. Fish. Sci. 64, 909–913 (1998).Article 
CAS 

Google Scholar 
Kumura, T., Yasui, H. & Mizuta, H. Nutrient requirement for zoospore formation in two alariaceous plants Undaria pinnatifida (Harvey) Suringar and Alaria crassifolia Kjellman (Phaeophyceae: Laminariales). Fish. Sci. 72, 860–869 (2006).Article 
CAS 

Google Scholar 
Qiu, Z. et al. Future climate change is predicted to affect the microbiome and condition of habitat-forming kelp. Proc. R. Soc. B 286, 20181887 (2019).Article 

Google Scholar 
Hoey, A. S. & Bellwood, D. R. Limited functional redundancy in a high diversity system: Single species dominates key ecological process on coral reefs. Ecosystems 12, 1316–1328. https://doi.org/10.1007/s10021-009-9291-z (2009).Article 

Google Scholar 
Streit, R. P., Hoey, A. S. & Bellwood, D. R. Feeding characteristics reveal functional distinctions among browsing herbivorous fishes on coral reefs. Coral Reefs 34, 1037–1047 (2015).Article 
ADS 

Google Scholar 
Van Alstyne, K. L. & Paul, V. J. The biogeography of polyphenolic compounds in marine macroalgae – Temperate brown algal defenses deter feeding by tropical herbivorous fishes. Oecologia 84, 158–163 (1990).Article 
ADS 

Google Scholar 
Targett, N. M., Boettcher, A. A., Targett, T. E. & Vrolijk, N. H. Tropical marine herbivore assimilation of phenolic-rich plants. Oecologia 103, 170–179 (1995).Article 
ADS 

Google Scholar 
Prado, P. & Heck, K. L. Seagrass selection by omnivorous and herbivorous consumers: Determining factors. Mar. Ecol. Prog. Ser. 429, 45–55. https://doi.org/10.3354/meps09076 (2011).Article 
ADS 

Google Scholar 
Montgomery, W. L. & Gerking, S. D. Marine macroalgae as foods for fishes: an evaluation of potential food quality. Environ. Biol. Fish. 5, 143–153 (1980).Article 

Google Scholar 
Duffy, J. & Paul & V.J.,. Prey nutritional quality and the effectiveness of chemical defenses against tropical reef fishes. Oecologia 90, 333–339 (1992).Article 
ADS 
CAS 

Google Scholar 
Michael, P. J., Hyndes, G. A., Vanderklift, M. A. & Vergés, A. Identity and behaviour of herbivorous fish influence large-scale spatial patterns of macroalgal herbivory in a coral reef. Mar. Ecol. Prog. Ser. 482, 227–240 (2013).Article 
ADS 

Google Scholar 
Bennett, S. & Bellwood, D. R. Latitudinal variation in macroalgal consumption by fishes on the Great Barrier Reef. Mar. Ecol. Prog. Ser. 426, 241–252 (2011).Article 
ADS 

Google Scholar 
Zarco-Perello, S., Wernberg, T., Langlois, T. J. & Vanderklift, M. A. Tropicalization strengthens consumer pressure on habitat-forming seaweeds. Sci. Rep. 7, 820. https://doi.org/10.1038/s41598-017-00991-2 (2017).Article 
ADS 
CAS 

Google Scholar 
Smith, S. M. et al. Tropicalisation and kelp loss shift trophic composition and lead to more winners than losers in fish communities. Glob. Change Biol. 27(11), 2537–2548 (2021).Article 
ADS 
CAS 

Google Scholar 
Zarco-Perello, S. et al. Range-extending tropical herbivores increase diversity, intensity and extent of herbivory functions in temperate marine ecosystems. Funct. Ecol. 34, 2411–2421. https://doi.org/10.1111/1365-2435.13662 (2020).Article 

Google Scholar  More

Study areaPenang is situated in the northern part of Peninsular Malaysia and lies within the latitudes 5°12’N to 5°30′ N and longitudes 100°09’E to 100°26’E (Fig. 7). Penang with a land area of 295 Km2, has an estimated population of 720,000 and is regarded as the most populated island in Malaysia. Penang shares the same border on the north and east with Kedah State and the south with Perak State. There are two main parts of Penang State: Penang Island and the mainland which is also regarded as Seberang Perai. These two parts of the State are connected by the two bridges. The eastern part of Penang Island is the most urbanized area comprising industries, commercial centres and residential buildings. However, the western part is less developed comprising mainly hilly terrain and forests22. This study is focused on the Island part of Penang. This island is endowed with a yearly equatorial climate (hot and humid). It has a mean annual temperature ranging between 27 and 30 °C while the mean annual relative humidity ranges between 70 and 90%. Also, the mean annual rainfall is about 267–624 cm.Figure 7The map of Penang State showing the Penang Island (created by the authors using ArcMap 10.8 software).Full size imageData acquisitionThe flow chart of the methodology is presented in Fig. 8. Landsat satellite images were used for the assessment of changes in land use covering a period of 2010–2021 (11 years).Figure 8The flow chart of the methodology.Full size imageThese images were gotten from the website of the United State Geological Survey (https://earthexplorer.usgs.gov). The Landsat images include the Landsat 5 TM (thematic mapper) and Landsat 8 OLI / TIRS (operational land imager / thermal infrared sensor). These were downloaded from the Landsat level 1 dataset (Table 6) with additional criteria which reduced the.Table 6 The characteristics of the satellite data used.Full size tableDetermination of LST and NDVI for Landsat 5 and 8Band 6 of Landsat 5 and band 10 of Landsat 8 were used for the determination of the land surface temperature (LST). The LST and normalized difference vegetation index were determined using the following steps:Conversion of top of atmosphere (TOA) radianceUsing the radiance rescaling factor, thermal infra-red digital numbers were converted to TOA spectral radiance using the equation below29: (frac{Red – NIR}{{Red + NIR}}) (frac{Red – NIR}{{Red + NIR}}) For Landsat 8,$$ {text{L}}lambda = left( {{text{ML}} times {text{ Qcal}}} right) + left( {{text{AL}} – {text{Oi}}} right) $$
(1)
For Landsat 5,$$ {text{L}}lambda = left( {{ }frac{{{text{LMax}}lambda – {text{LMin}}lambda }}{{{text{QcalMax}} – {text{QcalMin}}}}} right) times left( {left( {{text{Qcal }} – {text{QcalMin}}} right) + {text{LMin}}lambda } right) $$
(2)
where Lλ is TOA spectral radiance, ML is radiance multiplicative band Number, AL is radiance add band number, Qcal is quantized and calibrated standard product pixel values (DN for band 6 or band 10), Oi is the correction value for the respective bands, LMaxλ is spectral radiance scaled to QcalMax, LMinλ is spectral radiance scaled to QcalMin, QcalMax is maximum quantized calibrated pixel value, and QcalMin is minimum quantized calibrated pixel value.Conversion to TOA brightness temperature (BT)Spectral radiance data were converted to TOA brightness temperature using the thermal constant values in the Metadata file29.Kelvin (K) to Celcius (°C) degrees$$ BT = {raise0.5exhbox{$scriptstyle {K2}$} kern-0.1em/kern-0.15em lower0.25exhbox{$scriptstyle {ln left( {frac{K1}{{{text{L}}lambda { } + { }1}}} right)}$}} – 273.15 $$
(3)
where BT is the Top of atmosphere brightness temperature (°C), Lλ is TOA spectral radiance (W.m−2 .sr−1 .µm−1)), K1 is the K1 constant band number, and K2 is the K2 constant band number. For Landsat 5, K1 is 607.76, and K2 is 1260.56.Normalized difference vegetation index (NDVI)The Normalized Difference Vegetation Index (NDVI) is a standardized vegetation index which reveals the intensity of greenness and surface radiant temperature of the area30,31. The index value of NDVI usually ranges from − 1 to 1. The higher NDVI value indicates that the vegetation of the area is denser and healthier. This shows that the NDVI values of normal healthy vegetation range from 0.1– 0.75, while it is almost zero for rock and soil, and negative value for water bodies24. The NDVI is calculated using the followings:$$ {text{NDVI }} = frac{{left( {{text{NIR }}{-}{text{ RED}}} right){ }}}{{left( {{text{NIR }} + {text{ RED}}} right)}} $$
(4)
In Landsat 4–7$$ {text{NDVI }} = , left( {{text{Band 4 }}{-}{text{ Band 3}}} right) , / , left( {{text{Band 4 }} + {text{ Band 3}}} right) $$In Landsat 8$$ {text{NDVI }} = , left( {{text{Band 5 }}{-}{text{ Band 4}}} right) , / , left( {{text{Band 5 }} + {text{ Band 4}}} right) $$where: RED = DN values from the RED band, and NIR = DN values from the Near Infra-red band.Land Surface Emissivity (LSE)Land Surface Emissivity is the average emissivity of an element on the surface of the earth calculated from NDVI values.$$ {text{PV }} = left{ {frac{{left( {{text{NDVI }} – {text{ NDVImin}}} right)}}{{left( {{text{NDVImax }} – {text{ NDVImin}}} right)}}} right}^{2} $$
(5)
where PV is the Proportion of vegetation, NDVI is the DN value from the NDVI image, NDVImin is the minimum DN value from the NDVI image, and NDVImax is the maximum DN value from the NDVI image.$$ {text{E }} = left( {0.004{ } times {text{PV}}} right) + 0.986 $$
(6)
where E is land surface emissivity, PV is the Proportion of vegetation, 0.986 corresponds to a correction value of the equation.Land Surface Temperature (LST)Land Surface Temperature (LST) is the radiative temperature which is calculated using top of atmosphere brightness temperature, the wavelength of emitted radiance and land surface emissivity.$$ {text{LST}} = {raise0.5exhbox{$scriptstyle {BT}$} kern-0.1em/kern-0.15em lower0.25exhbox{$scriptstyle {left( {1 + left( {{lambda } times { }frac{{{text{BT}}}}{{{text{c}}2}}} right) times {text{ln}}left( {text{E}} right)} right)}$}} $$
(7)
Here c2 is 14388. The value of λ for Landsat 5 (Band 6) is 11.5 µm and Landsat 8 (Band 10) is 10.8 µm.Where BT is the top of atmosphere brightness temperature, λ is the wavelength of emitted radiance, and E is land surface emissivity.c2 = h*c/s (1.4388*10–2 mK = 14388 mK), h is Planck’s constant (6.626*1034 Js), s is Boltzmann constant (1.38*1023 JK), c is velocity of light (2.998*108 m/s).Determination of land use and land cover (LULC) of the study areaThe Landsat images were pre-screened and subjected to clipping and classification32. The boundary shape file of Penang was used to clip out the area of study.Image classificationThe unsupervised method involving a random assignment of sample training points and supervised methods of satellite image classification was employed in this study for determining the LULC types. This mixture of image classification methods has been reported as vital in achieving a high accuracy level33. Bands 5, 4 and 3 were used to classify Landsat 8 while bands 4, 3 and 2 were used for classifying Landsat 5. We used the extraction by mask in the spatial analyst tool of ArcMap 10.2.1 software to extract the study area from the selected bands of the Landsat satellite images. A widely used supervised image classification method was adopted for classifying the Landsat bands in this study32,34. The principle of operation of this method involves the identification of known sample training points which are then used to classify other unknown points with related spectral signatures35. The three monochromatic satellite bands were combined to produce the false colour composite (FCC) using the data management tool36. This involves drawing polygons on the LULC type to select the training points. The LULC types adopted for this study include urbanized areas, agricultural land, rocks, forests, bare surfaces, and water bodies. These were modified LULC types from IPOC Good Practice Guidance37. To achieve this, a minimum of 40 sample points were selected randomly for each category of LULC type36. Having prior knowledge of the study area assisted in the selection of the training points38.The multivariate maximum likelihood classification (MLC) technique was used for transforming the images. Other image transformation techniques have been used by researchers. These include the fuzzy set classifier, neural networks (NN) classifier, extraction and classification of homogenous objects (ECHO) classifier, per-field classifier, sub-pixel classifier, decision trees (DTs), support vector machines (SVMs), minimum distance classifier (MDC) and so-on39. The adoption of any of these techniques is dependent on the knowledge of the area of study, band selection, accessibility of data, the complexity of the landscape, the classification algorithm, and the proficiency of the analyst39. We preferred MLC to other techniques in this study due to its reported high level of accuracy in tropical regions32,34. Another reason for choosing MLC is that it is readily incorporated in many widely used GIS software packages. This MLC algorithm operates based on assigning pixels to the highest probability class and establishing the class ownership of such pixels. It is also regarded as a parametric classifier whose data follows almost a normal distribution39. We ensured the accuracy of this classifier by assigning a large number of training sample points using our prior knowledge of the study area.Description of the LULC categoriesThe urbanized area is the developed part of the study area. This includes houses, roads, railways, and industries. This is also known to be a settlement in other literature40. Agricultural land is the part of the study area dominated by agricultural activities and herbaceous plants and grasses. Agricultural land is generally a product of deforestation36. Rocks are part of the study area comprising solid mineral materials (rocks). Bare land is the bare soil which is either made open by natural or human activities.Forests are parts of the study area dominated by trees. They can be primary or secondary forests depending on the rate of disturbances. According to41, forest land is an area having more than 0.5 ha of flora comprising trees (height is above 5 m) with a canopy greater than 10%. The forests in Penang are generally both primary and secondary42. Water bodies are parts of the study area covered by water seasonally or permanently. These include seas, rivers, lakes, ponds, streams, or reservoirs40.Determination of change in the LULCThe rate and extent of change in the LULC of Penang within the periods under consideration were determined following the formula below43:$$ {text{Changed area }}left( {{text{C}}_{{text{a}}} } right) , = {text{ T}}_{{text{a}}} left( {text{year 2}} right) , {-}{text{ T}}_{{text{a}}} left( {text{year 1}} right) $$
(8)
$$ {text{Changed extent }}left( {{text{C}}_{{text{e}}} } right) , = {text{ C}}_{{text{a}}} /{text{ T}}_{{text{a}}} left( {text{year 1}} right) $$
(9)
$$ {text{Percentage of change }} = {text{ C}}_{{text{e}}} {text{x 1}}00 $$
(10)
where Ta means the total area.Determination of relationship between LST and NDVIThe values of LST and NDVI at 20 random points of each LULC class were used. The relationship between the LST and NDVI across all the LULC classes in each year was determined using the bivariate linear regression analysis. This was done in Paleontological Statistical (PAST) package 3.0.Classification accuracy assessmentThe classification accuracy was assessed by taking ground truth coordinate data of the LULC of the study area using a geographical positioning system (GPS) device (Garmin Etrex 10). These data were compared with the LULC classified in this study32. Consequently, an error matrix was generated. This normally uses ground truth data to explain the accuracy of the classified LULC. The error matrix comprises the user’s accuracy, the producer’s accuracy, overall accuracy and the Kappa index32.The producer’s accuracy (omission error) represents the probability of the correctly classified reference pixel and it is determined using this formula below:$${text{Producer’s accuracy }}left( % right) , = { 1}00% , – {text{ error of omission}} $$
(11)
Also, the user’s accuracy (commission error) represents the probability that the classified pixel matches the one on the ground36 and it is determined using the formula below:$$ {text{User’s accuracy }}left( % right) , = { 1}00% , – {text{ error of commission}} $$
(12)
The statistical accuracy of the matrix was determined using the Kappa coefficient44. This Kappa coefficient ranges from − 1 to + 145. Therefore, the overall accuracy of the classification was determined by dividing the total number of correctly classified pixels by the total number of sampled ground truth data40. 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

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

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