Archived natural DNA samplers reveal four decades of biodiversity change across the tree of life

doi:10.21203/rs.3.rs-5139547/v1

Download PDF

Biological Sciences - Article

Archived natural DNA samplers reveal four decades of biodiversity change across the tree of life

https://doi.org/10.21203/rs.3.rs-5139547/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Detecting the imprints of global environmental change on biological communities is a paramount task for ecological research1,2,3,4. But due to a lack of standardized long-term biomonitoring data, community assembly in the Anthropocene remains poorly understood5,6,7. Novel sources of data for analyzing biodiversity change across time and space are urgently needed8,9. By metabarcoding highly standardized biota samples from a long-term pollution monitoring archive in Germany, we here analyze four decades of community diversity for tens of thousands of species across the tree of life. The samples – tree leaves, marine macroalgae, and marine and limnic mussels – represent natural community DNA samplers10,11, preserving a taxonomically diverse imprint of their associated biodiversity at the time of collection9,12. We find no evidence for localized diversity declines1, but a strong, gradual compositional turnover as a universal pattern of temporal biodiversity change in Germany’s terrestrial and aquatic ecosystems5,13,14. This turnover results in biotic homogenization in terrestrial ecosystems, indicative of country-wide diversity decline. Our work uncovers a massive cryptic biodiversity loss across ecosystems, resulting from out-of-equilibrium ecological dynamics. This highlights the immense promise of alternative sample sources to provide standardized time series data of biodiversity change in the Anthropocene.

Biological sciences/Ecology/Community ecology

Biological sciences/Ecology/Biodiversity

Earth and environmental sciences/Ecology/Ecological genetics

Global ecosystems have undergone unprecedented change in the past decades. To understand the consequences of this transformation for biodiversity, standardized, long-term and taxonomically broad biomonitoring data are essential. Such data, however, are lacking for most taxa and ecosystems^15,16. Available time series are often short or incomplete and limited to a few target taxa and study locations^7,17.

A promising solution to this problem is offered by environmental specimen banks (ESBs), long-term pollution monitoring archives, which collect indicator organisms from various terrestrial and aquatic ecosystems. ESB samples are collected according to highly standardized protocols and stored at low temperatures¹⁸, ensuring an excellent preservation of the sample-associated nucleic acids. The indicator species collected by ESBs comprise metaorganisms representing diverse communities of interacting taxa, each of which has left a trace of its DNA in the sample. Recent work has shown that such metaorganisms are excellently suited “natural samplers” for studying the surrounding biota via DNA metabarcoding^9,10. Examples include sponges and mussels, which filter DNA particles from the water column, hence providing a snapshot of their aquatic communities, or plant material which retains DNA traces of interacting arthropods^11,12. The long-term archives of ESBs can thus serve as community DNA samplers, providing the standardized time series data so urgently needed to understand biodiversity change.

Here, we use metabarcoding of samples from the German ESB, one of the largest, technologically advanced and longest-operating ESBs, to reconstruct biodiversity change across broad taxonomic, spatial and temporal scales. We focus on four sample types from terrestrial, limnic and marine ecosystems and recover sample associated communities from across the tree of life. Using archived leaves from tree canopies, we characterize communities of canopy-associated fungi, bacteria and arthropods. Samples of a dominant marine macroalgal species along the European coastline reveal coastal bacterial and animal communities associated with the alga. Finally, marine and limnic mussels provide an imprint of the surrounding bacterial and eukaryotic communities in oceans and rivers. Our analysis yields highly standardized community diversity data for tens of thousands of species in parallel, many of which represent cryptic and largely understudied taxonomic groups.

Using these data, we explore common patterns of biodiversity change across the tree of life in terrestrial, marine and limnic habitats in Germany over the past four decades. Recent work has highlighted different responses of biota to changing environmental conditions in the Anthropocene. We explore the generality of these responses by testing three hypotheses for the temporal variation of biodiversity: 1) Stressful conditions have led to extinctions of species at individual sites, i.e. local declines of α-diversity¹. 2) The losses of species are countered by the immigration of novel taxa, leading to a pattern of biotic turnover (𝛽-diversity) without an overall loss of ɑ-diversity. This turnover could occur 2a) gradually with the changing environment⁵, or 2b) rapidly, when the community reaches a tipping point^2,20. 3) The immigration of species across broad geographic scales, for example of widespread invasive taxa, leads to a pattern of biotic homogenization, i.e. a loss of 𝛽-diversity across space¹⁴. To test these hypotheses, we developed a dynamic model for community assembly, based on the equilibrium theory of island biogeography¹⁹, which generates null expectations of diversity trends in the absence of disturbance (Extended Data Figure 3A).

Figure 1. Overview of sampling sites and recovered biodiversity for ESB samples from terrestrial, limnic and marine ecosystems spanning a time window from 1985 to 2022. A) ESB sampling sites of tree leaves, zebra mussels, bladderwrack and blue mussels across Germany. B) Total OTU and order level richness (log scale) of different taxonomic groups across the tree of life associated with tree leaves, zebra mussels, bladderwrack and blue mussels. Taxonomic groups are represented by different colors. Icons refer to natural DNA sampler species.

We metabarcoded 550 samples of archived natural sampler organisms from three marine, nine limnic and nine terrestrial sites in Germany (Figure 1B; Supplementary Table 1). Using time series of leaves from trees (1985-2022), marine macroalgae (1985-2021), and marine (1985-2020) as well as limnic mussels (1994-2020), we reconstructed communities associated with these organisms at the time of sampling. Rarefaction and bootstrapping analyses indicated sufficient sequencing depth and sampling size for biodiversity estimations at both local and regional scales (Extended Data Figure 1 + Extended Data Table 1). Our analysis recovered highly diverse prokaryote and eukaryote communities (Figure 1B), a total of 66,184 operational taxonomic units (OTUs) in 751 orders and 102 phyla. Tree leaves recovered 5,183 OTUs of bacteria in 94 orders, 6,250 fungal OTUs in 113 orders and 3,271 metazoan (mainly arthropod) OTUs in 24 orders. We found 5,474 bacterial OTUs in 101 orders and 787 metazoan OTUs in 78 orders in marine macroalgae. 21,266 OTUs of bacteria in 180 orders and 3,551 OTUs of microeukaryotes (mainly algae and protozoa) in 160 orders were found in marine blue mussels. In limnic zebra mussels, we found 14,292 bacterial OTUs in 184 orders, 5,587 microeukaryote (mainly algae and protozoa) OTUs in 173 orders and 523 metazoan OTUs in 71 orders (Figure 1B, Supplementary Table 2).

The detected species accurately represented their respective ecosystems and natural sampler organisms (Extended Data Figure 2; Supplementary Table 3). For example, various typical coastal metazoans were found in bladderwrack samples and numerous species of eukaryotic algae reflect the phytoplankton community surrounding mussels. Typical canopy-dwelling arthropods and leaf-associated fungi and bacteria were recovered from leaves (Supplementary Table 3). The species detected in tree canopies showed a high specificity for their respective host tree, with each tree species harboring a unique community of leaf-associated taxa (Extended Data Figure 2). Typical monophagous arthropods or host-specific microorganisms were exclusively found on their respective host tree (Supplementary Table 3). The communities recovered from all sample types were also site-specific, with many taxa limited to single sites, likely indicating their specific ecological requirements or biogeographic affinities. For example, highly disparate communities were associated with bladderwrack and blue mussels from the Baltic versus the Northern Sea; these two seas are distinguished by pronounced salinity differences. Additionally, our two Northern Sea sampling sites, which are separated by about 200 km, harbored different sets of taxa. The same held true for the leaf-associated communities at different forest sites across Germany and the zebra mussel-associated communities in different river systems (Extended Data Figure 2), which have different biogeographic affinities. These results underline the accuracy and efficiency of natural sampler time series for the targeted enrichment of biological communities.

Figure 2. Multidecadal trends of α-diversity, as well as temporal- and spatial 𝛽-diversity in communities across the tree of life natural DNA sampler organisms from the ESB. A) Trends in OTU richness (α-diversity) of the associated communities between 1991 and 2021. B) Temporal OTU turnover (temporal 𝛽-diversity) of the associated communities calculated from the turnover component of Jaccard distance is shown as a function of the distance in years between samples of the same sampling site. C) Spatial OTU turnover (spatial 𝛽-diversity) of the associated communities between 1991 and 2021. Trends describe the degree of dissimilarity in community composition between different sampling locations over time. OTU richness and turnover values were summarized as mean with standard error bars across sampling locations and time windows. D)-F) Diversity trends from A)-C) reduced to their respective slope. Filled circles indicate significant departures from the null expectations generated through the dynamic model for community assembly, suggesting an out-of-equilibrium dynamic.

The communities recovered by the natural DNA samplers provide highly standardized and well resolved temporal biodiversity data across ecosystems and the tree of life. This offers an unprecedented opportunity to test for general patterns of biodiversity change. We analyzed temporal patterns of α-diversity (OTU richness), compositional OTU turnover (temporal 𝛽-diversity) and biotic homogenization (spatial 𝛽-diversity). To evaluate the significance of the observed trends, we compared them to null expectations of community changes in the absence of anthropogenic disturbances. These expectations were based on a new dynamic model of community ecology we developed, built upon the equilibrium theory of island biogeography^5,19 (Extended Data Figure 3A; see Methods).

Recent and often prominently featured work on insect decline has highlighted losses of α-diversity (richness) as a driving feature of biodiversity change in the Anthropocene^1,17,21. Our data did not support this assumption, refuting our first hypothesis. No universal trend for α-diversity was identified. Irrespective of ecosystems and taxonomic groups, richness either remained stable, increased, or declined (Figure 2A+D, Extended Data Figure 4A). A particularly pronounced drop in richness was found in marine prokaryotes. In contrast, richness strongly increased in limnic prokaryotes (Figure 2A+D) across the studied river systems. Interestingly, a reverse pattern was found for limnic and marine micro-eukaryotes, which showed α-diversity losses of in limnic habitats, but increases in marine sites (Figure 2A+D). A more common pattern was found in terrestrial ecosystems, where most studied taxa showed a slight, increase of α-diversity at individual. Surprisingly, we even found a slight albeit significant increase in richness for terrestrial arthropod communities, supporting recent work suggesting that forest canopy arthropods are not affected by pronounced insect decline ⁹.

In contrast to α-diversity, a clear universal trend was found for compositional turnover (temporal 𝛽-diversity). This trend significantly deviated from the null expectations in all studied communities, supporting our second hypothesis⁵. We found considerable losses of OTUs in all communities, which however, were countered by the immigration of novel, likely better adapted taxa (Figure 2B+E, Extended Data Figure 4B). This led to an out-of-equilibrium dynamic with a stronger-than-expected temporal turnover of community composition (Extended Data Table 2). Slightly different temporal dynamics of compositional turnover were observed across taxa (Extended Data Figure 3B). The higher rates of local extinction and immigration observed in metazoans and algae led to OTU-poor, rapidly-changing communities. Bacterial communities were also characterized by higher immigration rates, but combined with lower local extinction rates. This resulted in more diverse, dynamic communities. In contrast, lower extinction and immigration rates observed in fungi generated slower taxonomic turnover compared with other communities, but still faster than expected. The observed turnovers were gradual in all communities: no abrupt breaks of community composition were detected (Figure 2B+E). The observation of a gradual compositional turnover as a predominant pattern of biodiversity change in the Anthropocene is well supported by other recent work^4,5,22. We found no indication of rapid state shifts in communities reaching tipping points². This is also supported by the patterns at higher taxonomic levels. While we observed a significant compositional OTU turnover in all the ecosystems and taxonomic groups, the turnover did not affect higher taxonomic ranks: the temporal composition of phyla, classes or orders remained remarkably stable in all the natural samplers (Extended Data Figure 5). Hence, those taxa that disappeared were primarily replaced by closely related and likely functionally similar ones, which are better adapted to changing environmental conditions.

The observed temporal change of 𝛽-diversity, i.e. a compositional turnover of communities, can also be seen at the level of thousands of individual OTUs within our data. Each sampler type and data set showed replacements of various OTUs, with gains and losses approximately balanced (Extended Data Figure 6). Various novel colonizers were detected, among them typical invasive species. For example, the occurrence of the invasive Pacific oyster (Crassostrea gigas) in coastal bladderwrack in our 2004 data reflects the timing of its actual introduction into the German coastline²³. Similarly, the colonization of various plant pathogens can be traced over time. For example, Norway spruce and European beech at the national park Harz have been colonized by a Taphrina species and spruces at Lake Belau were infested by Coniothyrium (Supplementary Table 4). At the same time, declines in the occurrence of several taxa are evident, for example the common periwinkle (Littorina littorea) in marine ecosystems and the green silver-lines (Pseudoips prasinana) in forest ecosystems. Both are considered common species and their decline was previously overlooked. ESB samples thus can serve as an important early warning system for both the decline of local species and the emergence of problematic invaders ^9,10.

We last explored patterns of spatial 𝛽-diversity, e.g. biotic homogenization. Similar to α-diversity, no universal trend across ecosystems was evident from our data (Figure 2C+F). Most aquatic sampler organisms did not show a clear trend of spatial homogenization over time. Some of the communities became even more heterogeneous across space, for example in marine microeukaryotes. The increasing spatial heterogeneity was especially evident in prokaryotes and microeukaryotes in limnic sites (Figure 2C+F). Different river systems in Germany have become more heterogenous over time, possible due to colonization by novel taxa from different biogeographic regions, for example between Danube and Rhine system. Interestingly, a general pattern of homogenization across space was evident in terrestrial canopy ecosystems. All the taxonomic groups, from prokaryotes to fungi and arthropods, showed a significant spatial homogenization of communities over time in the terrestrial samples (Figure 2C+F, Extended Data Figure 4C). As in the temporal 𝛽-diversity, the observed change was gradual rather than abrupt. This finding is in line with recent observations that widespread generalist species benefit from environmental change and replace more specialized ones^14,24,25,26. Terrestrial ecosystems did not lose diversity at the level of individual sites, in fact many sites even showed an increasing α-diversity. Despite the local increase of diversity, a clear loss of diversity was evident at larger spatial scales. This indicates the immigration of widespread taxa across the tree of life in Germany’s canopy ecosystems.

In summary, our work closes critical gaps in understanding of temporal and spatial biodiversity change in the Anthropocene. Archived natural samplers provide time series data of unprecedented taxonomic resolution and standardization and enable the study of patterns of biodiversity across the tree of life¹⁸, from prokaryotes to multicellular eukaryotes. We report an out-of-equilibrium dynamic in both aquatic and terrestrial ecosystems in Germany, potentially driven by anthropogenic disturbances²⁷. The widespread taxonomic replacement and biotic homogenization we observe across deeply divergent taxonomic groups constitute serious environmental issues, which are hitherto insufficiently characterized. Our data indicate that thousands of species have disappeared from Germany's aquatic and terrestrial ecosystems in the past decades. The majority of the increasing or declining taxa we observed are highly cryptic and rarely detected by monitoring programs and biodiversity research, which focus on prominent taxa like plants²⁸, vertebrates²⁹ or insects¹. At the same time, many of the taxa disappearing from our dataset represent critical elements in food webs and are essential to ecosystem stability, for example phytoplankton in aquatic ecosystems or arthropods in tree canopies³⁰. Our data highlight a hitherto poorly understood facet of environmental change: a massive cryptic loss of biodiversity by turnover and biotic homogenization.

Data availability

Raw sequencing data will be available in the Science Data Bank upon publication with the identifier DOI 10.57760/sciencedb.13553.

The R code for the dynamic model will be uploaded to a public repository upon publication.

Acknowledgements

We thank Karin Fischer for assistance with laboratory work. This work was funded by departmental funds of the department of Biogeography at Trier University and by the TrendDNA project of the German Environment Agency (Umweltbundesamt). We want to especially thank the ESB project group at Trier University and the German Environment Agency for providing the time series samples.

Author contributions

IJ, JH, AM, SW, MS, EG, CS, AS, and NS performed the laboratory work and analyzed the data. HK, IJ, JH, BP, and HM conceptualized the study, analyzed the data, and wrote the manuscript. BP and HM developed the dynamic model for community ecology. SK performed the sequencing. MP, RK, and DT led the sampling. JK enabled access to the samples. All authors contributed to revising the final version of the manuscript.

The authors declare no competing interests.

Supplementary Information is available for this paper.

Correspondence and requests for materials should be addressed to Henrik Krehenwinkel.

Seibold, S. et al. Arthropod decline in grasslands and forests is associated with landscape-level drivers. Nature, 574, 671-674 (2019). https://doi.org/10.1038/s41586-019-1684-3.
Barnosky, A. D. et al. Approaching a state shift in Earth’s biosphere. Nature, 486, 52-58 (2012). https://doi.org/10.1038/nature11018.
Jandt, U. et al. More losses than gains during one century of plant biodiversity change in Germany. Nature, 611, 512-518 (2022). https://doi.org/10.1038/s41586-022-05320-w.
Blowes, S. A. et al. The geography of biodiversity change in marine and terrestrial assemblages. Science, 366, 339-345 (2019). DOI: 10.1126/science.aaw1620
Dornelas, M. et al. Assemblage time series reveal biodiversity change but not systematic loss. Science, 344, 296-299 (2014). DOI: 10.1126/science.1248484.
Johnson, T. F. et al. Revealing uncertainty in the status of biodiversity change. Nature, 628, 788-794 (2024). https://doi.org/10.1038/s41586-024-07236-z
Crossley, M. S. et al. No net insect abundance and diversity declines across US Long Term Ecological Research sites. Nature Ecology & Evolution, 4, 1368-1376 (2020). https://doi.org/10.1038/s41559-020-1269-4.
Outhwaite, C. L., Gregory, R. D., Chandler, R. E., Collen, B. & Isaac, N. J. B. Complex long-term biodiversity change among invertebrates, bryophytes and lichens. Nature Ecology & Evolution, 4, 384-392 (2020). https://doi.org/10.1038/s41559-020-1111-z.
Krehenwinkel, H. et al. Environmental DNA from archived leaves reveals widespread temporal turnover and biotic homogenization in forest arthropod communities. eLife, 11, e78521 (2022a). https://doi.org/10.7554/eLife.78521.
Junk, I., Schmitt, N. & Krehenwinkel, H. Tracking climate-change-induced biological invasions by metabarcoding archived natural eDNA samplers. Current Biology, 33, 943-944 (2023). https://doi.org/10.1016/j.cub.2023.07.035.
Mariani, S., Baillie, C., Colosimo, G. & Riesgo, A. Sponges as natural environmental DNA samplers. Current Biology, 29, 401-402 (2019). https://doi.org/10.1016/j.cub.2019.04.031.
Weber, S. et al. Molecular diet analysis in mussels and other metazoan filter feeders and an assessment of their utility as natural eDNA samplers. Molecular Ecology Resources, 23, 471-485 (2023). https://doi.org/10.1111/1755-0998.13710.
Pilotto, F. et al. Meta-analysis of multidecadal biodiversity trends in Europe. Nature Communications, 11, 3486 (2020). https://doi.org/10.1038/s41467-020-17171-y.
Blowes, S. A. et al. Synthesis reveals approximately balanced biotic differentiation and homogenization. Science Advances, 10, eadj9395. (2024). DOI: 10.1126/sciadv.adj9395.
Bálint, M. et al. Environmental DNA Time Series in Ecology. Trends in Ecology & Evolution, 33, 945-957 (2018). DOI: 10.1016/j.tree.2018.09.003.
Magurran, A. E. et al. Long-term datasets in biodiversity research and monitoring: assessing change in ecological communities through time. Trends in Ecology & Evolution, 25, 574-582 (2010). https://doi.org/10.1016/j.tree.2010.06.016.
Wagner, D., Grames, E. M., Forister, M. L., Stopak, D. Insect Decline in the Anthropocene: Death by a thousand cuts. PNAS, 118, e2023989118 (2021a). https://doi.org/10.1073/pnas.2023989118.
Fliedner, A. et al. Environmental specimen banks and the European Green Deal. Science of the Total Environment, 852, 158430 (2022). https://doi.org/10.1016/j.scitotenv.2022.158430.
MacArthur, R.H., & Wilson, E.O. An equilibrium theory of insular zoogeography. Evolution, 17, 373-387 (1963). https://doi.org/10.2307/2407089.
Maucieri, D.G., Starko, S. & Baum, J.K. Tipping points and interactive effects of chronic human disturbance and acute heat stress on coral diversity. Proceedings of the Royal Society B, 290, 20230209 (2023). https://doi.org/10.1098/rspb.2023.0209.
Wagner, D., Fox, R, Salcido, D. M. & Dyer, L. A. A window to the world of global insect declines: Moth biodiversity trends are complex and heterogeneous. PNAS, 118, e2002549117 (2021b). https://doi.org/10.1073/pnas.2002549117.
Magurran, A. E. et al. Divergent biodiversity change within ecosystems. PNAS, 11, 1843-1847 (2018). https://doi.org/10.1073/pnas.1712594115.
Brandt, G., Wehrmann, A. & Wirtz, K.W. Rapid invasion of Crassostrea gigas into the German Wadden Sea by larval supply. Journal of Sea Research, 59, 279-296 (2008). https://doi.org/10.1016/j.seares.2008.03.004.
Stuart-Smith, R. D., Mellin, C., Bates, A. E. & Edgar, E. J. Habitat loss and range shifts contribute to ecological generalization among reef fishes. Nature Ecology & Evolution, 5, 656–662 (2021). https://doi.org/10.1038/s41559-020-01342-7.
Dharmarajam, G., Gupta, P., Vishnudas, C. K. & Robin, V.V. Anthropogenic disturbance favours generalist over specialist parasites in bird communities: Implications for risk of disease emergence. Ecology Letters, 24, 1859-1868 (2021). https://doi.org/10.1111/ele.13818.
Xu, W. B. et al. Regional occupancy increases for widespread species but decreases for narrowly distributed species in metacommunity time series. Nature Communications, 14, 1463 (2023). https://doi.org/10.1038/s41467-023-37127-2.
Marta, S., Brunetti, M., Manenti, R., Provenzale, A. & Ficetola, G.F. Climate and land-use changes drive biodiversity turnover in arthropod assemblages over 150 years. Nature Ecology & Evolution, 5, 1291-1300 (2021). https://doi.org/10.1038/s41559-021-01513-0.
Jandt, U. et al. More losses than gains during one century of plant biodiversity change in Germany. Nature, 611, 512-518 (2022). https://doi.org/10.1038/s41586-022-05320-w.
Semenchuk, P. et al. Relative effects of land conversion and land-use intensity on terrestrial vertebrate diversity. Nature Communications, 13, 615 (2022). https://doi.org/10.1038/s41467-022-28245-4.
Field, C. B., Behrenfeld, M. J., Randerson, J. T. & Falkowski, P. Primary production of the biosphere: Integrating terrestrial and oceanic components. Science, 281, 237-240 (1998). DOI: 10.1126/science.281.5374.237.

Specimen Bank Data

The German Environmental Specimen Bank (ESB) has been in operation since the early 1980s. The ESB collects samples of indicator species from various terrestrial and aquatic ecosystems. These species serve as accumulators of environmental chemicals and provide a detailed image of pollution and ecosystem health. The long operational history of ESBs results in exceptionally long time series of environmental pollution. To make pollution analysis comparable between years, ESB samples are collected according to highly standardized protocols. Samples are taken at the same time of the year, at identical sites and using identical protocols. ESB collection is done using sterile equipment to avoid carryover of even trace amounts of pollutants between samples. To ensure preservation of unstable chemical compounds, the samples are stored over liquid nitrogen after collection and for the long term, halting all chemical and biological degradation. To acquire an integrative view of pollution in an ecosystem, ESB samples are large, each one including hundreds to thousands of specimens or tissue compartments (in case of trees)^18,31-36. Each sample is cryo-milled to a fine powder of a grain size of 200 µm, thoroughly homogenizing all traces of chemicals³⁷.

Recent work shows that ESBs, being set up as chemical pollution archives, also makes them ideal for studies on biodiversity change. ESB indicator species can be considered metaorganisms or natural community DNA samplers, which preserve an imprint of the surrounding biological community at the time of sampling. The highly standardized and contamination-free sampling and sample processing conditions, coupled with storage at ultra-low temperatures, make ESB samples perfectly suited for metabarcoding. The cryomilling of large ESB samples also guarantees an even distribution of community DNA traces among the sample and breaks open cell walls of various microorganisms, whose DNA is then uniformly released. Former studies have already extensively tested and highlighted the suitability of different ESB samples for retrospective biodiversity monitoring^9,10,12. Here, we use four different types of ESB samples from terrestrial, limnic and marine habitats as natural community DNA samplers to measure biodiversity change across four decades.

1) Tree leaves: The ESB collects leaves from three common tree species in Germany, namely European beech (Fagus sylvatica), Norway spruce (Picea abies) and Lombardy poplar (Populus nigra). The leaves are collected once annually or biannually and serve as samplers for aerial pollutants deposited on the leaves^31,32,33. ESB leaf samples are sampled from different forest ecosystem types, spanning a land use gradient from core zones of national parks to timber forests, forests neighboring agricultural sites, and urban parks. Sample series from nine sites were included in this study, starting from 1985. Each sample contains hundreds of leaves from at least 15 individual trees, milled to a fine powder. These samples contain DNA traces of all organisms that interacted with the tree canopy at the time of collection⁹. Here we characterize communities of canopy-associated arthropods, fungi and bacteria.

2) Bladderwrack (Fucus vesiculosus): This macroalgae is widespread along the European coastline and makes up a significant part of the biomass in Europe’s coastal ecosystems. Marine pollutants are enriched in the tissue of the algae, making it an ideal sentinel species for pollution³⁴. Three sites have been sampled annually or biannually for bladderwrack thalli beginning in 1985. ESB samples from two North Sea sites are collected at intervals of two months, six times a year, and then merged into a pooled annual sample. The third site at the Baltic Sea is sampled twice a year. Bladderwrack is a critical species in coastal ecosystems, providing a habitat for countless taxa. All these taxa leave detectable DNA traces in the ESB sample¹². Here we characterize communities of animals and bacteria that interacted with the bladderwrack.

3) Blue mussels (Mytilus edulis): This is the most common mussel in coastal ecosystems of Northern and Central Europe. Blue mussels constantly filter the water column for planktonic organisms and organic particles. In doing so, they enrich pollutants in their tissue, making them an excellent sentinel species for pollution monitoring³⁵. The ESB has collected blue mussels at three coastal sites in Germany since 1985. The mussel’s entire soft tissue including respiratory water is used for the sample. Annual or biannual samples of hundreds of mussels are compiled from six sampling events at the North Sea and two at the Baltic Sea. With each mussel filtering roughly 1 liter of water per hour, these samples contain a comprehensive imprint of the annual planktonic biodiversity at the sampling site¹². Here we characterize communities of eukaryotic plankton and mussel-associated bacteria.

4) Zebra mussels (Dreissena polymorpha): The limnic zebra mussel is an invasive species from the Black Sea region, which has colonized nearly all major rivers of Germany since the 1960s. Like the blue mussel, zebra mussels are highly efficient filter feeders. Since the 1990s, zebra mussels are reared by the ESB on special plate stacks, which are then placed in four major German rivers for about one year, allowing the mussels to accumulate pollutants in their tissue. The mussels are then collected from the plate stacks, immediately deep-frozen and a sample of soft tissue including respiratory water is compiled from thousands of mussels³⁶. The samples from nine sites used here provide an overview of limnic biodiversity from major rivers of Germany¹². Here we characterize communities of animals, eukaryotic plankton and mussel-associated bacteria.

DNA extraction, library preparation, sequencing and sequence processing

Samples were processed as described in Krehenwinkel et al. (2022a) and Junk et al. (2023). Work steps were performed under clean benches to avoid carryover and cross-contamination. We isolated DNA from 200 mg of homogenate from each sample, which was shown to be sufficient to recover the sample-associated diversity⁹. DNA was extracted in one or two replicates depending on the sample type (Supplementary Table 5) using a CTAB protocol (OPS Diagnostics, New Jersey, USA), which proved best suited to extract high purity DNA from these sample types. The DNA extracts were then amplified for different DNA barcode markers to enrich various taxonomic groups from the samples (for a list of barcode markers and PCR conditions see Supplementary Table 5). PCR was performed using the Qiagen Multiplex PCR kit (Hilden, Germany) in 10-µl volumes according to the manufacturer’s protocols. Primers were chosen to amplify the associated community, but not the ESB indicator species itself, whose DNA dominates the extract. To characterize bacterial communities (all sample types), we amplified the V1 or V5-7 region of 16SrDNA^39,40,41. For terrestrial arthropods (tree leaf samples), we used a mitochondrial Cytochrome Oxidase I (COI) marker³⁸. For fungi (tree leaf samples) we used the ITS1 region of the nuclear ribosomal cluster^42,43. The variable V9 region of nuclear 18SrDNA was targeted to characterize communities of aquatic animals and eukaryotic plankton (bladderwrack and mussel samples¹²; for primer details, see Supplementary Table 5). PCR success was checked on 1.5 % agarose gels, and the PCR products were then amplified in another round of PCR to add Illumina Truseq adapters and unique combinations of dual indexes⁴⁴ (Supplementary Table 5). All final libraries were pooled in approximately equimolar amounts, cleaned of leftover primers using 1X AMPure beads XP (Beckman-Coulter, California, USA), and then sequenced on an Illumina Miseq (California, USA) using paired-end sequencing with 500-cycle V2 and 600-cycle V3 kits. To ensure reproducibility of our data and to recover rare species, we ran several PCR replicates for every sample, which were indexed and sequenced separately. The number of PCR replicates was adapted based on sample type and marker and varied between three and six (Supplementary Table 5). Blank DNA extractions were included in every batch of extractions, and non-template control PCRs were run alongside all PCR reactions. All controls were sequenced along with the samples to provide a baseline for carryover or cross-contamination during processing.

Forward and reverse reads were merged using PEAR⁴⁵ with a minimum quality of 20 and a minimum overlap of 50 bp. The merged reads were then quality-filtered by limiting the number of expected errors in a sequence to 1⁴⁶ and transformed to fasta format using USEARCH⁴⁷. Primer sequences were trimmed off using UNIX scripts. Long 18SrDNA amplicons (~350 bp) of aquaticmetazoan and microeukaryotic plankton generated from zebra mussels were trimmed to match the corresponding short amplicon of ~150 bp. Since both metazoan and phytoplankton amplicons span exactly the same nuclear 18SrDNA region, all sequences were combined into one file. Likewise, reads generated from the three tree species were saved to one file for each marker. After trimming, the resulting file for each marker and sample type was dereplicated and clustered into zero radius OTUs (hereafter OTUs) using the USEARCH pipeline. OTU tables were built for each sample type and marker, also using USEARCH. Taxonomy was annotated using blast2taxonomy script v1.4.2.⁴⁸ after BLAST searching⁴⁹ against the entire NCBI Genbank database for 18SrDNA and COI with a maximum number of 10 target sequences. The Silva database⁵⁰ was used for annotating 16SrDNA sequences, and the UNITE database⁵¹ for fungal ITS1. The FungalTraits database⁵² was used for the functional annotation of fungi. Taxonomic assignments based on BLAST hits with a base pair length of less than 80 % of the amplicon length and/or less than 85 % sequence identity were removed. Wed excluded all taxa except Bacteria, Algae/Protozoa, Metazoa or Fungi from the respective datasets.

Only OTUs with a minimum of 3 reads per sample replicate were retained. The OTU tables were checked for contamination using negative controls. PCR replicates were merged and all datasets were checked for sufficient sequencing depth and sampling size (Extended Data Figure 1; for resulting sampling and OTU count as well as number of phyla and orders see Supplementary Table 2 “cleaned dataset”).

Statistical model and analyses of community diversity

For each type of sample, we only (i) selected sites that were sampled for at least 5 years, (ii) only kept sampling years represented by at least 50 % of the sites and (iii) removed years that were isolated from the others (>2 years). We also removed samples with low read coverage (less than 50 % of the median number of reads). Finally, because OTU read abundances from metabarcoding datasets are subject to many biases⁵³, we converted OTU abundances into binary presence/absence data and only analyzed trends in terms of OTU occurrence. To limit cross-contamination, we considered an OTU as present in a sample if it represents at least 0.01 % of the total reads (for resulting sampling and OTU count as well as number of phyla and orders see Supplementary Table 2 “filtered dataset (model)”; Extended Data Table 1).

We measured community diversity trends in three different ways. First, in each community, in each year, we computed the ɑ-diversity using the OTU richness as a measure of local diversity at a given time. Second, we computed 𝛽-diversities between pairs of communities (1) sampled at the same site in different years or (2) sampled at different sites in the same year. (1) gives an idea of temporal turnover in community composition (temporal 𝛽-diversity), whereas (2) indicates the spatial turnover of the communities (spatial 𝛽-diversity). We measured 𝛽-diversities using the turnover component of the Jaccard distances (R-package betapart ⁵⁴).

To identify temporal trends in these diversity indices, temporal models were fitted using mixed linear models accounting for the temporal autocorrelation between sampling years and the effect of the different sampling sites. We used the lme function (R-package nlme⁵⁵) with the corAR1 temporal correlation and the different sites as random effects. We fitted these temporal models with either the ɑ-diversities or the temporal and spatial 𝛽-diversities as response variables.

The significance of the observed trends was evaluated by comparing them to null expectations of community changes. Changes in community composition occur as a result of immigration and extinction events, which are influenced by neutral and/or niche factors. This generates a dynamic equilibrium with ever-changing communities, even in the absence of any kind of disturbance. Following Dornelas et al. (2014), we built upon the equilibrium theory of island biogeography (ETIB)¹⁹ to set up a dynamic model for community assembly that generates null expectations of diversity trends in the absence of disturbance (Extended Data Figure 3A). The ETIB is a lineage-based model of species colonization of a local community (the island) from a metacommunity (the continent). In its simplest form¹⁹, at each time step, it assumes that each OTU has a probability to migrate from the metacommunity to the local community , and once settled in the community, each OTU has a probability to go extinct. The number of new immigration events (i.e. of species not already in the community) per time step is given by where is the total number of OTUs in the metacommunity and is the number of OTUs already present in community ; it declines as the number of OTUs in the community increases. The number of extinction events per time step, given by , increases with the number of OTUs in the community. An equilibrium is reached when the number of immigration events per time step equals the number of extinction events, i.e. . The equilibrium number of OTUs in the community is given by . This simple form of ETIB implies a linear decrease or increase of the number of new immigration events (resp. extinction events) per unit of time with the number of settled species in the local community. It assumes that all species have the same probabilities to migrate or go extinct (neutrality), and that these probabilities do not depend on the number of species in the community, implying that there is a negligible influence of interspecific competition on immigration and extinction. It thus applies best to communities that are far from carrying capacity. This model has the advantage of being very straightforward to simulate using a simple discrete-time Markov chain.

Given the incomplete sampling of communities typically achieved with metabarcoding techniques, we assumed that each OTU present in community is observed at each time step with a fixed probability . This extra parameter can be handled using hidden Markov models. We assumed that the rates , , and vary from one community to the other due to various extrinsic and intrinsic factors (e.g. distance to the metacommunity, community size, and environmental factors).

Assuming neutrality, i.e. that all OTUs have the same immigration and extinction rates, is a strong assumption often proven to be unrealistic⁵⁶. We therefore relaxed the assumption of neutrality. In our non-neutral model, we assumed that immigration (resp. extinction) rates for each OTU are sampled from a beta distribution with parameters and (resp. ). The rates are therefore sampled around (resp. ) with a variance that is inversely proportional to : a large corresponds to scenarios of neutrality whereas closer to 0 indicates that immigration and extinction rates are very different across OTUs. While each OTU has specific immigration and extinction rates in each community, we assumed that the ranks of the OTUs in terms of immigration and extinction rates are conserved across communities (i.e. an OTU with a low extinction rate in community compared with other OTUs will also have a low extinction rate in community ). We thus obtained a non-neutral model derived from the equilibrium theory of island biogeography, which assumes that the presence of an OTU in a community results from the balance between immigration and extinction, and that each OTU is characterized by specific rates of immigration and extinction centered around the average rates. At equilibrium, some OTUs are more likely to be present in the community (e.g. OTUs with higher immigration and lower extinction rates).

Instead of testing different parameters values chosen a priori⁵, we implemented a new inference strategy to adjust the model parameters to the empirical data using a sequential technique. First, in each site, the sampling fraction was inferred using ACE (R-package vegan⁵⁷). Second, we estimated the average rates and by fitting the neutral model to each community using a hidden Markov model (R-package seqHMM⁵⁸). We used these estimates as community-specific estimates of the average rates of immigration and extinction of the non-neutral model. Third, given , , and , we used simulation-based inferences using artificial neural networks to estimate the parameters and . We generated a large number of simulated datasets by sampling and from uniform prior distributions and simulating the corresponding non-neutral model of community assembly, and for each of these simulations, we recorded ɑ-, spatial and temporal 𝛽-, and γ-diversities through time across all the sampled sites. We specifically incorporated subsampling into the simulations (with probability ), such that not all OTUs present in the local communities are observed, mimicking the detection bias of metabarcoding: simulated ɑ-, spatial and temporal 𝛽-, and γ-diversities are therefore directly comparable to empirical diversites. For , we used a uniform prior distribution between the number of observed OTUs and three times the ACE estimate of γ-diversity; for , we used a uniform prior between 1 and 5. We started the simulations at year 1500 with a random community composition at each site (each OTU has a probability to be initially present in the community, where is the theoretical number of species at equilibrium in the neutral model, given , and : ). Next, we simulated community composition over time in each site until 2023, sampled the communities according to , and recorded community composition through time for the years of sampling. We then computed for each simulation the ɑ- and 𝛽-diversities across all the sampled sites. We used the same methods and sampling scheme as for empirical data (the number of sampling sites varied through time) to obtain measures of ɑ-, spatial and temporal 𝛽-diversities. We trained an artificial neural network to estimate and from time series of ɑ-, 𝛽-, and γ-diversities using the Python library Keras⁵⁹, with 100,000 simulations per dataset until reaching a sufficient predictive power. Once trained, the artificial neural network takes as input ɑ- and 𝛽-diversities through time and outputs estimates of and . We used a neural network with 3 intermediate layers, containing 132, 64, and 32 neurons respectively, and with ELU activation functions. We prevented overfitting by using a dropout of 0.5 at each intermediate layer. Input and output data were scaled between 0 and 1 before fitting and the simulations were split between the training set (90 %) and the test set (10 %). Once validated (Extended Data Figure 7), we finally applied the trained neural network to our empirical data, and obtained corresponding and values.

Given the estimated and values, we simulated our non-neutral model 1000 times with these parameters to generate time series of ɑ- and 𝛽-diversities under an equilibrium model of community assembly. We then compared empirical and simulated temporal trends and considered an empirical trend to be significant if it was above or below 95 % of the simulated trends. We interpreted significant deviations from the simulated equilibrium model of community assembly as indicative of out-of-equilibrium dynamics in the empirical data, potentially driven by anthropogenic disturbances.

We then investigated the temporal trends of OTU occurrences, as some OTUs may be recent invaders and others may have gone extinct at the regional scale. In each ecosystem, we only looked at OTUs present in at least 10 % of the samples and across at least 2 sites. For each OTU, we first tested whether its occurrence in the communities tended to vary through time. To do so, we fitted a generalized linear mixed model with a binomial response (presence/absence of the OTU in a community) and considered the sampling site as random effects.

References Methods

Tarricone, K., Klein, R., Paulus, M. & Teubner, D. Guideline for Sampling and Sample Processing Red beech (Fagus sylvatica). Umweltbundesamt (2018b).
https://www.umweltprobenbank.de/upb_static/fck/download/SOP_ESB_Red_Beech_V2.0.3_2018_en.pdf.

Klein, R., Tarricone, K., Teubner, D. & Paulus, M. Guideline for Sampling and Sample Processing Norway Spruce (Picea abies)/ Scots Pine (Pinus sylvestris). Umweltbundesamt (2018).
https://www.umweltprobenbank.de/upb_static/fck/download/SOP_ESB_Spruce_Pine_V2.0.2_2018_en.pdf.

Tarricone, K., Klein, R., Paulus, M. & Teubner, D. Guideline for Sampling and Sample Processing Lombardy Poplar (Populus nigra ‘Italica’). Umweltbundesamt (2018a).
https://www.umweltprobenbank.de/upb_static/fck/download/SOP_ESB_Lombardy_Poplar_V2.0.3_2018_en.pdf.

Quack, M. et al. Guideline for Sampling and Sample Processing Bladderwrack (Fucus vesiculosus). Umweltbundesamt (2010).
https://www.umweltprobenbank.de/upb_static/fck/download/SOP_Blasentang.pdf.

Paulus, M., Klein., R. & Teubner, D. Guideline for Sampling and Sample Processing Blue Mussel (Mytilus eduliscomplex). Umweltbundesamt (2018).
https://www.umweltprobenbank.de/upb_static/fck/download/SOP_ESB_Blue_Mussel_V2.1.0_2018_en.pdf.

Teubner, D., Klein, R., Tarricone, K. & Paulus, M. Guideline for Sampling and Sample Processing Zebra Mussel (Dreissena polymorpha). Umweltbundesamt (2018).
https://www.umweltprobenbank.de/upb_static/fck/download/SOP_ESB_Zebra_Mussel_V2.1.1_2018_en.pdf.

Rüdel, H., Uhlig, S. & Weingärtner, M. Zerkleinerung und Homogenisierung von Umweltproben durch Cryomahlung. Fraunhofer-Institut für Molekularbiologie und Angewandte Oekologie (2008).
https://www.umweltprobenbank.de/upb_static/fck/download/IME_SOP_Probenvorbereitung_Dez2008_V200.pdf

Krehenwinkel, H., Weber, S., Künzel, S. & Kennedy, S. R. The bug in a teacup-monitoring arthropod-plant associations with environmental DNA from dried plant material. Biology Letters, 18, 20220091 (2022b).
https://doi.org/10.1098/rsbl.2022.0091.

Donia, M. S. Fricke, W. F., Partensky, F. & Schmidt, E. W. Complex microbiome underlying secondary and primary metabolism in the tunicate Prochloron symbiosis. PNAS, 108, 1423-1432 (2011). https://doi.org/10.1073/pnas.1111712108.

Bodenhausen, N., Horton, M. W. & Bergelson, J. Bacterial communities associated with the leaves and the roots of Arabidopsis thaliana. PloS One, 8, e56329 (2013). https://doi.org/10.1371/journal.pone.0056329.

Chelius, M. K. & Triplett, E. W. The Diversity of Archaea and Bacteria in Association with the Roots of Zea mays L. Microbial Ecology, 41, 252-263 (2001). DOI: 10.1007/s002480000087.

42. Gardes, M. & Bruns, T. D. ITS primers with enhanced specificity for basidiomycetes--application to the identification of mycorrhizae and rusts. Molecular Ecology,2, 113-118 (1993). DOI: 10.1111/j.1365-294x.1993.tb00005.x.

White, T. J., Bruns, T., Lee, S.& Taylor, J. Amplification and direct sequencing of fungal ribosomal RNA genes for phylogenetics. In : T. J. White, coord.: PCR Protocols: Elsevier, 315-322 (1990). https://doi.org/10.1016/B978-0-12-372180-8.50042-1.

Lange, V. et al. Data from: Cost-efficient high-throughput HLA typing by MiSeq amplicon sequencing. BMC Genomics15, 63 (2014). https://doi.org/10.1186/1471-2164-15-63.

Zhang, J., Kobert, K., Flouri, T. & Stamatakis, A. PEAR: a fast and accurate Illumina Paired-End reAd mergeR. Bioinformatics, 30, 614-620 (2014). https://doi.org/10.1093/bioinformatics/btt593.

Edgar, R. C. & Flyvbjerg, H. Error filtering, pair assembly and error correction for next-generation sequencing reads. Bioinformatics, 31, 3476-3482 (2015). https://doi.org/10.1093/bioinformatics/btv401.

Edgar, R. Usearch. Lawrence Berkeley National Lab (LBNL), Berkeley, CA (United States) (2010).

Schöneberg, Y. „yschoeneberg/blast2taxonomy“. (2023). Zenodo. doi: 10.5281/zenodo.10009721.

Altschul S. F., Gish W., Miller W., Myers E. W. & Lipman D. J. Basic local alignment search tool. Journal of Molecular Biology, 215, 403-410 (1990). https://doi.org/10.1016/S0022-2836(05)80360-2.

Quast, C. et al. The SILVA ribosomal RNA gene database project: improved data processing and web-based tools. Nucleic Acids Research, 41, 590-596 (2012). https://doi.org/10.1093/nar/gks1219.

Abarenkov, K. et al. The UNITE database for molecular identification and taxonomic communication of fungi and other eukaryotes: sequences, taxa and classifications reconsidered. Nucleic Acids Research, 52, 791-797 (2024). https://doi.org/10.1093/nar/gkad1039.

Põlme, S. et al. FungalTraits: a user-friendly traits database of fungi and fungus-like stramenopiles. Fungal Diversity105, 1-16 (2020). https://doi.org/10.1007/s13225-020-00466-2.

Piñol J., Mir G., Gomez-Polo P. & Agustí N. Universal and blocking primer mismatches limit the use of high-throughput DNA sequencing for the quantitative metabarcoding of arthropods. Molecular Ecology Resources, 15, 819-830 (2015). https://doi.org/10.1111/1755-0998.12355.

Baselga, A. et al. Package ‘betapart’: Partitioning Beta Diversity into Turnover and Nestedness Components (2023). https://cran.r-project.org/web/packages/betapart/index.html.

Pinheiro, J. et al. Package ‘nlme’: Linear and Nonlinear Mixed Effects Models (2024). https://cran.r-project.org/web/packages/nlme/index.html.

Lomolino M. V. A call for a new paradigm of island biogeography. Global Ecology and Biogeography, 9, 1-6 (2001). https://doi.org/10.1046/j.1365-2699.2000.00185.x.

Oksanen, J. et al. Package ‘vegan’: Community Ecology Package (2024). https://cran.r-project.org/web/packages/vegan/vegan.pdf.

Helske, S. & Helske, J. Mixture Hidden Markov Models for Sequence Data : The seqHMM Package in R. Journal of Statistical Software, 88, Article 3 (2019). https://doi.org/10.18637/jss.v088.i03.

Chollet, F. et al. Keras. (2015). https://github.com/fchollet/keras.

There is NO Competing Interest.

TableS1S5.xlsx
Supplementary Table 1-5
ExtendedDataFigures.docx

Download PDF

Version 1

posted

You are reading this latest preprint version

Archived natural DNA samplers reveal four decades of biodiversity change across the tree of life

Status:

Version 1

Abstract

Figures

Main

Declarations

References

Methods

Additional Declarations

Supplementary Files

Status:

Version 1