Environmental DNA from archived leaves reveals widespread temporal turnover and biotic homogenization in forest arthropod communities

  1. Henrik Krehenwinkel  Is a corresponding author
  2. Sven Weber
  3. Rieke Broekmann
  4. Anja Melcher
  5. Julian Hans
  6. Rüdiger Wolf
  7. Axel Hochkirch
  8. Susan Rachel Kennedy
  9. Jan Koschorreck
  10. Sven Künzel
  11. Christoph Müller
  12. Rebecca Retzlaff
  13. Diana Teubner
  14. Sonja Schanzer
  15. Roland Klein
  16. Martin Paulus
  17. Thomas Udelhoven
  18. Michael Veith
  1. University of Trier, Germany
  2. German Federal Environment Agency, Germany
  3. Max Planck Institute for Evolutionary Biology, Germany
  4. Ludwig Maximilians University, Germany

Abstract

A major limitation of current reports on insect declines is the lack of standardized, long-term, and taxonomically broad time series. Here, we demonstrate the utility of environmental DNA from archived leaf material to characterize plant-associated arthropod communities. We base our work on several multi-decadal leaf time series from tree canopies in four land use types, which were sampled as part of a long-term environmental monitoring program across Germany. Using these highly standardized and well-preserved samples, we analyze temporal changes in communities of several thousand arthropod species belonging to 23 orders using metabarcoding and quantitative PCR. Our data do not support widespread declines of α-diversity or genetic variation within sites. Instead, we find a gradual community turnover, which results in temporal and spatial biotic homogenization, across all land use types and all arthropod orders. Our results suggest that insect decline is more complex than mere α-diversity loss, but can be driven by β-diversity decay across space and time.

Editor's evaluation

This landmark study reveals novel temporal arthropod biodiversity insights that can be leveraged from environmental DNA traces, that have been cryopreserved on leaf tissue as part of a long-term monitoring scheme. The strength of the evidence underlying the major conclusions is convincing and limitations in the quantitative aspects of the data synthesis are acknowledged appropriately. The work will be of interest to a breadth of ecological practitioners.

https://doi.org/10.7554/eLife.78521.sa0

eLife digest

Insects are a barometer of environmental health. Ecosystems around the world are being subjected to unprecedented man-made stresses, ranging from climate change to pollution and intensive land use. These stresses have been associated with several recent, dramatic declines in insect populations, particularly in areas with heavily industrialised farming practices.

Despite this, the links between insect decline, environmental stress, and ecosystem health are still poorly-understood. A decline in one area might look catastrophic, but could simply be part of normal, longer-term variations. Often, we do not know whether insect decline is a local phenomenon or reflects wider environmental trends. Additionally, most studies do not go far back enough in time or cover a wide enough geographical range to make these distinctions.

To understand and combat insect decline, we therefore need reliable methods to monitor insect populations over long periods of time. To solve this problem, Krehenwinkel, Weber et al. gathered data on insect communities from a new source: tree leaves. Originally, these samples were collected to study air pollution, but they also happen to contain the DNA of insects that interacted with them before they were collected – for example, DNA deposited in chew marks where the insects had nibbled on the leaves. This is called environmental DNA, or eDNA for short.

To survey the insect communities that lived in these trees, Krehenwinkel, Weber et al. first extracted eDNA from the leaves and sequenced it. Analysis of the different DNA sequences from the leaf samples revealed not only the number of insect species, but also the abundance (or rarity) of each species within each community. Importantly, the leaves had been collected and stored in stable conditions over several decades, allowing changes in these insect populations to be tracked over time.

eDNA analysis revealed subtle changes in the make-up of forest insect communities. In the forests where the leaves were collected, the total number of insect species remained much the same over time. However, many individual species still declined, only to be replaced by newcomer species. These ‘colonisers’ are also widespread, which will likely lead to an overall pattern of fewer species that are more widely distributed – in other words, more homogeneity.

The approach of Krehenwinkel, Weber et al. provides a reliable method to study insect populations in detail, over multiple decades, using archived samples from environmental studies. The information gained from this has real-world significance for environmental issues with enormous social impact, ranging from conservation, to agriculture and even public health.

Introduction

Dramatic declines of terrestrial insects have been reported in recent years, particularly in areas of intensified land use (Hallmann et al., 2017; van Klink et al., 2020; Seibold et al., 2019; Sánchez-Bayo and Wyckhuys, 2019). However, some authors have urged caution in generalizing these results (Didham et al., 2020; Thomas et al., 2019; Cardoso et al., 2019), suggesting that reported patterns of decline may be more localized than currently assumed or reflect long-term natural abundance fluctuations (Macgregor et al., 2019; Crossley et al., 2020). Most studies on insect decline suffer from a lack of long-term time series data and are limited in geographic and taxonomic breadth, often using biomass as a proxy for diversity estimates (Hallmann et al., 2017; Daskalova et al., 2021). Hence, what is needed are methods and sample types that yield standardized long-term time series data for the diversity of arthropod communities across broad taxonomic and geographic scales (Forister et al., 2021).

In recent years, environmental DNA metabarcoding has offered a promising new approach to monitor biological communities (Taberlet et al., 2018; Tautz et al., 2002; Thomsen and Sigsgaard, 2019). This includes terrestrial arthropods, whose eDNA can be recovered from various substrates, for example plant material (Nakamura et al., 2017). Here, we develop a DNA metabarcoding and quantitative PCR (qPCR) protocol to simultaneously recover diversity and relative DNA copy number of arthropod community DNA from powdered leaf samples. We then analyze 30-year time series data of arthropod communities from canopy leaf material of four tree species from 24 sites across Germany. These sites represent four land use types with different degrees of anthropogenic disturbance: urban parks, agricultural areas, timber forests, and national parks (Figure 1). The samples were collected using a highly standardized protocol by the German Environmental Specimen Bank (ESB), a large biomonitoring effort for Germany’s ecosystems, and stored at below −150°C. By basing our analyses on DNA sequences, we can measure diversity from haplotype variation within species to taxonomic diversity of the whole community.

Overview of different sampling sites and samples in this study.

(A) Sampling sites, land use types, and the four tree species (European Beech Fagus sylvatica (N = 98), Lombardy Poplar Populus nigraitalica’ (N = 65), Norway Spruce Picea abies (N = 123), and Scots Pine Pinus sylvestris (N = 26)) sampled by the ESB in Germany. (B) Representative images of sampling sites for the four land use types and average number of detected pesticides (and 95 % confidence interval) across these land use types in three time periods (before 2000, 2000–2009, and 2010–2018). Gas chromatography–tandem mass spectrometry (GC–MS/MS) analysis (Löbbert et al., 2021) shows that the detected pesticide load distinguishes the different land use types, with agricultural sites continuously showing the highest number of pesticides.

Current studies on insect decline primarily focus on site-based assessments of α-diversity and biomass (Hallmann et al., 2017; Seibold et al., 2019). Yet, these metrics alone are insufficient for characterizing ongoing biodiversity change (Marta et al., 2021; Kortz and Magurran, 2019; Magurran and Henderson, 2010). Significant temporal community change and declines can also occur at the scale of β- or γ-diversity, without affecting local richness. This may be driven by community turnover and spatial biotic homogenization (Karp et al., 2012). But diversity may also vary temporally. Fluctuating occurrences of transient species considerably increase diversity within single sites over time (D’Souza and Hebert, 2018). The loss of such transient species in favor of taxa with a temporally stable occurrence results in an increasingly predictable community and hence a temporal diversity decline. Biotic turnover may occur gradually following changing environmental conditions, but can also occur abruptly, when ecosystems reach tipping points (Barnosky et al., 2012). In the latter case, a rapid and considerable biotic remodeling of the ecosystem may be found. The relevance of these spatial and temporal factors in insect decline remains largely elusive.

Here, we use our high-resolution data on arthropods from canopy leaf samples to test the hypotheses that (1) α-diversity and biomass of canopy-associated arthropod communities have declined in the last 30 years (Hallmann et al., 2017; Seibold et al., 2019), or (2) community change has occurred in the form of turnover and possibly homogenization of communities across space and time (Kortz and Magurran, 2019; Dornelas et al., 2014; Thomsen et al., 2016). Last, we hypothesize that (3) biodiversity declines will be particularly pronounced in areas of intensified land use (Outhwaite et al., 2022).

Results

A standardized protocol to characterize plant-associated arthropod communities

We developed a standardized and robust protocol to reproducibly recover plant-associated arthropod communities from powdered leaf material. We controlled for effects of the amount of leaf material per sample, rainfall before sampling, amount of leaf homogenate used for DNA extraction, extraction replication, and primer choice on the recovered diversity (see Methods and Figure 2—figure supplements 15 for details on standardization).

Using our optimized protocol, we analyzed 312 ESB leaf samples. We recovered 2054 OTUs from our samples, with different tree species having significantly different OTU numbers on average (Figure 2A, linear mixed model [LMM], p < 0.05). Nevertheless, all tree species showed a balanced and relatively similar taxonomic composition at the order level (Figure 2A, Figure 2—figure supplement 5E). We identified 23 orders, 218 families, and 413 genera. The richest order was Diptera (600 OTUs in 48 families), followed by Hymenoptera (369 OTUs in 21 families), Acari (293 OTUs in 21 families), Lepidoptera (233 OTUs in 32 families), Hemiptera (152 OTUs in 19 families), Coleoptera (133 OTUs in 29 families), and Araneae (99 OTUs in 15 families). The recovered species assemblages were ecologically diverse, including herbivores, detritivores, predators, parasites, and parasitoids (Figure 2—figure supplement 6). Each tree species harbored a unique arthropod community (Figure 2B, Figure 2—figure supplement 5C, D), with typical monophagous taxa exclusively recovered from their respective host trees. The arthropod communities from different sites and land use types were also differentiated within tree species (Figure 2—figure supplement 7, PERMANOVA, p < 0.05). In addition to arthropod–host plant associations, we were able to detect interactions between arthropods. For example, abundances of the spruce gall midge Piceacecis abietiperda and its parasitoid, the chalcid wasp Torymus sp., were well correlated across all analyzed spruce sites (LM, p < 0.05). Both underwent coupled abundance cycles, with similar maxima every 6–8 years (Figure 2C).

Figure 2 with 7 supplements see all
Recovery of diversity and interactions in canopy-associated arthropod communities from leaf material.

(A) Barplot showing recovered order composition of OTUs across the four tree species (N = 312). In addition to ESB samples, results from freeze-dried leaves stored at room temperature for 6–8 years (N = 5) and hand-collected bulk insect samples (N = 5) from European Beech are shown. Orders amounting to less than 1% of the total OTU number are merged as ‘Other’. Numbers below each barplot show the total number of arthropod OTUs followed by the mean OTU number per sample. (B) Non-metric multidimensional scaling (NMDS) plot showing tree-specific composition of arthropod communities for the same samples. (C) Temporal changes in abundance of the spruce gall midge Piceacecis abietiperda and its parasitoid Torymus sp. between 1987 and 2018 in a spruce forest in the Saarland. Inset (D) shows correlation of relative abundance between the two species across all ESB spruce samples.

The recovered community composition from ESB leaf homogenate samples was similar to hand-collected branch clipping samples (Figure 2A, B). Branch clipping recovered a larger diversity of spiders, a taxon which is exclusively found on leaf surfaces. In contrast, about 25% of the recovered taxa from ESB leaf powder likely inhabited the insides of leaves, e.g., gallers and miners (Figure 2—figure supplement 6B). Overall, arthropod DNA in leaf homogenates appears temporally very stable: Even freeze-dried leaf material that had been stored at room temperature for 8 years yielded surprisingly similar arthropod communities to ESB samples (Figure 2A, B). Besides analyzing diversity, we generated information on relative arthropod 18S rDNA copy number in relation to the corresponding plant 18S rDNA copy number by qPCR. Relative eDNA copy number should be a predictor for relative biomass. We designed a standardized qPCR assay based on lineage-specific blocking SNPs in the 18SrDNA gene (Figure 3A). We tested two primer combinations with 3′-blocking SNPs in (1) only the forward or (2) both forward and reverse primer sequences. The primer combination with mismatches in both forward and reverse primers led to a near complete suppression of non-arthropod amplification in all tested samples (5.41% vs. 44.49% on average) and was hence chosen for the qPCR experiment (Figure 3A, B). The qPCR assay accurately predicted changes in relative copy number of arthropod DNA on plant material across a 10,000-fold dilution series. Even when comparing taxonomically heterogeneous mock communities, very similar CT values were recovered (Figure 3C, Figure 3—figure supplement 1), highlighting the accuracy and wide applicability of our approach. Overall, we found a significant positive correlation of OTU richness and relative arthropod copy number in the ESB samples (LMM, p < 0.05; see Methods: ‘Statistical analysis’), supporting recent work suggesting a biomass–diversity relationship (Hallmann et al., 2021).

Figure 3 with 1 supplement see all
Quantitative PCR (qPCR) experiment to detect relative rDNA copy number of arthropods in plant homogenates.

(A) Schematic overview of the blocking approach to amplify homologous 18SrDNA fragments for either arthropod or plant DNA, based on lineage-specific priming mismatches. (B) Effect of primer mismatches on the recovery of arthropod sequences. Barplots show recovered read proportion of different higher taxa from 15 ESB leaf samples. The left plot shows the effect of a diagnostic mismatch in forward and reverse primers, while the right plot shows the effect of only a forward primer mismatch. (C) Boxplot showing CT values recovered from the seven mock communities of arthropod species from 13 different orders (see Figure 3—figure supplement 1 for community composition) and across a 1.10000 dilution series. Separate CT values for each community are indicated by the dots (Figure 3—figure supplement 1).

Temporal changes of diversity, copy number, and species composition in canopy arthropod communities

Based on our time series data of archived ESB leaf samples, we tested the hypothesis that α-diversity (including intraspecific genetic diversity) and biomass (relative rDNA copy number) have undergone widespread temporal declines, particularly in areas of intensive land use. Our statistical analysis does not support previously reported widespread temporal α-diversity declines (Figure 4A, Figure 2—figure supplement 5H, I, LMM, p > 0.05), even when different land use types are analyzed separately (Figure 4—figure supplement 1, LMM, p > 0.05). Instead, warm summers and cold winters were negatively associated with richness (LMM, p < 0.05). The temporal pattern of diversity was also largely independent of taxonomy: most orders did not show temporal trends when analyzed separately (Figure 4C; Figure 4—figure supplement 2). Exceptions include a significant loss of lepidopteran diversity, which is primarily driven by OTU loss at urban sites, and an overall increasing diversity of mites (LMM, p < 0.05). The overall temporally stable diversity is also visible at separate sites; a diversity decline across all orders was observed at only a single site (Figure 4—figure supplement 3). Similar to α-diversity, we did not find widespread temporal declines of genetic diversity. Neither community-level zero radius OTU (zOTU) richness (which was well correlated to OTU richness, R2 = 0.89) nor within-OTU haplotype richness declined significantly over time (Figure 4B; Figure 4—figure supplement 1C, D).

Figure 4 with 3 supplements see all
Temporal changes of diversity and copy number across all sites and samples (N = 312).

(A) Arthropod OTU richness (representing α-diversity). (B) Haplotype richness within OTUs (representing genetic variation). (C) Relative OTU richness per order. (D) Relative copy number of arthropod DNA (representing biomass).

In contrast to the stable α-diversity, relative copy number showed an overall decrease over time (Figure 4D, LMM, p < 0.05), suggesting that arthropod biomass may indeed be declining in woodlands (Seibold et al., 2019). The effect appears to be particularly driven by urban sites, coinciding with a loss of lepidopteran diversity (Figure 4—figure supplements 13). However, declines of arthropod DNA copy number are also visible in several agricultural and timber forest sites, particularly in the last 10 years of our time series (Figure 4—figure supplement 1, Figure 4—figure supplement 3).

We next explored temporal changes in abundance for 413 separate OTUs from a total of 19 sites. In line with our hypothesis (1), we predicted a majority of declining species. However, we found no significant difference between the average number of declining (6.94%) and increasing (10.04%) OTUs (t-test, p > 0.05, Figure 5A). With the exception of Acari, which showed an overrepresentation of increasing OTUs, declines and increases in OTU read abundance were independent of arthropod order and land use type (Figure 5A, B, Fisher’s exact test p > 0.05). The observed replacement of about 15% of OTUs within sites translates into a significant temporal change of taxonomic β-diversity (Figure 5C, Figure 2—figure supplement 5J and K). We found a strong positive correlation of temporal distance and Jaccard dissimilarity for most sites (PERMANOVA, p < 0.05). Thus, species are continuously replaced in all land use types (Figure 4—figure supplement 1F). In the majority of analyzed sites, β-diversity did not show a correlation with differences in copy number (Figure 5—figure supplement 2).

Figure 5 with 2 supplements see all
Temporal changes of species composition and β-diversity within and between sampling sites.

(A) Boxplots of the proportion of 413 separate OTUs that significantly increased, decreased, or did not show temporal abundance changes. The colored symbols overlaying the boxplots represent land use type for each site. (B) Stacked barplot showing the recovered order composition of the same three categories of OTUs. Orders amounting to less than 1% of the total OTU number are merged as ‘Other’. (C) Within-site Jaccard dissimilarity as a function of number of years between sampling events, showing a pronounced taxonomic turnover (N = 312). (D) Jaccard dissimilarity between consecutive sampling years in each decade of sampling (N = 15, 44 & 39), calculated within sites for beech forests. We find a loss of temporal β-diversity over time. Points indicate means and error bars show 95% confidence interval. (E) Jaccard dissimilarity between sites for the same samples as in D. We find a loss of average between-site β-diversity, that is spatial homogenization. (F) NMDS plot showing community dissimilarity within and between Bavarian Forest and Harz Mountains national parks over three decades (<2000: sampled before 2000; <2010: sampled between 2000 and 2009; <2020: sampled between 2010 and 2018).

While the community turnover did not affect local α-diversity, we still observed associated losses of overall diversity. The first noteworthy pattern concerns a loss of temporal β-diversity within sites. β-Diversity between consecutive sampling years dropped significantly in many sites, particularly in beech forests (Figure 5D). Thus, diversity within sites is increasingly homogenized over time. We also found a significant decrease of β-diversity between sites for beech forests (Figure 5E). Our data suggest a loss of site-specific species and a gain of more widespread generalists, irrespective of land use. This pattern also emerges at the level of individual OTUs. Several novel colonizers spread rapidly in woodlands and showed similar abundance trends across various sites in parallel (Figure 5—figure supplement 1). The spatial and temporal change of β-diversity is illustrated by an NMDS plot of arthropod communities from two beech forests in National Parks, the Harz and the Bavarian Forest (Figure 5F). While the two sites are well separated by the first NMDS axis, the second axis shows a pronounced temporal turnover of communities. In the past decades, this turnover has led to temporally more predictable communities within sites and increasingly similar communities between the two national parks.

Discussion

Here, we show that DNA from archived leaf material provides a robust source of data to reconstruct temporal community change across the arthropod tree of life. Leaf samples should also cover a broad phenological window: Adults of many insect species are only active during a short time period of the year, but their larvae spend the whole year on their host plant (Gagné and Graney, 2014). A Malaise trap will miss these taxa most times of the year, while an eDNA approach should detect the larvae throughout. As sampled leaves make up the habitat and often food source of arthropods, it is also possible to infer the exposure of arthropod communities to chemical pollution by analyzing chemicals in the leaf material. This is of critical importance, as pesticide use has often been invoked as a driver of insect decline (Goulson et al., 2015; Siviter et al., 2021). The high temporal stability of arthropod DNA in 8-year-old dried leaf samples also suggests the utility of other plant archives, for example, herbaria, for arthropod eDNA. However, for such less standardized and less well-preserved sample types, careful evaluation of cross-contamination or chemical DNA modification (Orlando et al., 2021), which could inflate the recovered diversity, may be warranted.

We here analyze a unique leaf archive and provide unprecedented insights into arthropod community change in the tree canopy, an ecosystem known for its high and often cryptic diversity (Nakamura et al., 2017). Our results do not confirm the hypothesis of widespread losses of α-diversity. Initial reports of insect declines originate mainly from grassland ecosystems that have undergone massive changes in land use (Hallmann et al., 2017). Central European canopy communities may be less affected by such change. Interestingly, the only sites that showed declines of richness were agricultural and urban sites, suggesting that land use may at least locally affect neighboring canopy communities (Hallmann et al., 2014). Our data also suggest negative effects of warm summers on richness. Climate warming has recently been suggested to act in conjunction with land use change to drive insect declines (Outhwaite et al., 2022).

Instead of declining richness, we detected DNA copy number declines in all land use types, suggesting that overall biomass may indeed be in decline. Extinction is the endpoint of a long trajectory of decline and increasingly affects common species (Collen et al., 2011), hence biomass decline could foreshadow future biodiversity loss. Alternatively, the dropping copy number may reflect a taxonomic turnover of species with either different eDNA shedding rates or different rDNA copy numbers in the different communities. However, if the latter were true, an association of copy number changes and turnover would be expected, which we did not find in our data.

Instead of widespread losses of α-diversity, we indeed found pronounced taxonomic turnover in nearly all communities (Dornelas et al., 2014), supporting our second hypothesis. While resident species are continuously lost, they are mostly replaced by novel colonizers. This turnover can result in biotic homogenization across space and time. Less interannual variation of the occurrence of taxa within sites and reduced spatial variation of their occurrence between sites both cause a decline in overall β-diversity. Biotic homogenization is often associated with intensification of land use (Karp et al., 2012; Gossner et al., 2016) and landscape simplification (Holland et al., 2005). However, the pattern we observed affected all land use types equally. The universality of these changes suggests that neither site- nor taxon-specific factors are responsible. Possible explanations include factors that act at a larger scale, such as climate change-induced range shifts (Marta et al., 2021), nitrogen deposition (Gámez-Virués et al., 2015), and the introduction of invasive species (Soroye et al., 2020; Lister and Garcia, 2018). Given that our leaf samples recover a fairly broad phenological window, the alternative explanation that the observed community-wide turnover pattern may have resulted from shifting phenologies (Cohen et al., 2018) is unlikely. The gradual replacement of species also suggests that we are not yet observing ecosystems reaching tipping points (Barnosky et al., 2012).

In summary, our work shows the great importance of standardized time series data to accurately reveal biodiversity change in the Anthropocene (Thomsen et al., 2016) across space and time, beyond the decline of α-diversity and biomass (Dornelas et al., 2014). Taxonomic replacement and biotic homogenization, even in seemingly pristine habitats such as national parks, signify an important and hitherto insufficiently recognized facet of the current insect crisis.

Materials and methods

Samples and metadata used in this study

Tree samples of the German Environmental Specimen Bank – standardized time series samples stored at ultra-low temperatures

Request a detailed protocol

We used a total of 312 leaf samples of four common German tree species: the European Beech Fagus sylvatica (98 samples), the Lombardy Poplar Populus nigraitalica’ (65 samples), the Norway Spruce Picea abies (123 samples), and the Scots Pine Pinus sylvestris (26 samples). The samples have been collected annually or biannually by the German Environmental Specimen Bank (ESB) since the 1980s and serve as indicators for aerial pollutants (Schulze et al., 2007). A total of 24 sampling sites were included, covering sampling periods of up to 31 years and representing four land use types of varying degrees of anthropogenic disturbance (Figure 1). These include natural climax forest ecosystems in core zones of national parks (six sites, National Park), forests commercially used for timber (six sites, Timber Forest), tree stands in close proximity to agricultural fields (six sites, Agricultural), and trees in urban parks (three sites, Urban). The sites were initially chosen to represent their land use type permanently for long-term monitoring, and the corresponding land use categories have mostly remained temporally stable.

ESB samples are collected and processed according to a highly standardized protocol at the same time every year. Sampling events between different years of the time series usually do not differ by more than 2 weeks. All used equipment is sterilized before field work by several washes and heat treatments (Tarricone et al., 2018a; Tarricone et al., 2018b; Klein et al., 2018). A defined amount of leaf material (>1.100 g) is collected from a defined number of trees (15 at most sites) from each site and from 4 branches from each tree. The branches are distributed equally spaced in the outer crown area of the tree. The amount sampled translates to several thousand leaves from each site, which should suffice to saturate the recovered arthropod diversity. For a subset of samples, biometric analysis is performed, for example the weight of individual leaves and general condition of the tree are noted. Leaf weight has not changed over the time series at most sites. The sampled leaves are intended to represent the exact natural state of the tree. They are not washed or altered in any way before processing, and eDNA traces, and small arthropods on the leaves’ surfaces, as well as from galls and leaf mines, are included in the sample. Each sample is stored on liquid nitrogen immediately after collection and ground to a powder with an average diameter of 200 µm using a cryomill. The cryomill is thoroughly cleaned between samples to prevent cross-contamination. The resulting homogenates are then stored for long term on liquid nitrogen (Rüdel et al., 2009; Rüdel et al., 2015). The cold chain is not interrupted after collection and during processing, ensuring optimal preservation of nucleic acids in the samples. The homogenization of the sample also guarantees a thorough mixing, resulting in equal distribution of environmental chemicals and probably eDNA in the sample. Previous work suggests that very small subsamples of the homogenate suffice to detect even trace amounts of environmental chemicals in the sample (Gámez-Virués et al., 2015). eDNA from leaf surfaces may be affected by weather conditions before sampling, for example, washed away by heavy rain or damaged by strong UV exposure (Valentin et al., 2020). However, leaves are only collected dry by the ESB, that is, not immediately after rain. The date of the most recent rainfall is noted for each ESB sampling event, allowing us to explore the effects of recent weather conditions on the recovered arthropod communities.

The utility of plant material stored at room temperature to recover arthropod DNA

Request a detailed protocol

ESB samples are stored under optimal conditions for nucleic acid preservation. By contrast, most archived leaf samples are stored at room temperature, for example dried leaf material in herbaria. To test the general suitability of archived leaf material for arthropod community analysis, we included 25 additional beech leaf samples, each consisting of 100 leaves from a total of 5 sites. The samples were freeze-dried and then ground to a fine powder by bead beating. The resulting powder was comparable to our ESB homogenates, but unlike them, it was stored at room temperature for 6–8 years.

Hand-collected branch clipping samples to explore the accuracy of leaf-derived arthropod DNA

Request a detailed protocol

To evaluate the performance of our leaf DNA-based protocol in comparison to commonly used sampling methods, we generated a branch clipping dataset from five beech stands close to Trier University. Branch clipping is a widely used method to collect arthropods residing on leaf surfaces in trees (Delvare, 1997) and thus the best comparable traditional methodology to our protocol. We sampled five trees per site and collected ten branches of about 40 cm length from each tree. The branches were clipped off, stored in plastic bags and then brought to the laboratory. Here, arthropods were manually collected from each sample. All collected arthropod specimens were pooled by tree and stored for later DNA extraction in 99% ethanol.

Climate data

Request a detailed protocol

We downloaded monthly climate data for all study sites from the German Climate Center distributed as a raster dataset interpolated from the surrounding weather stations by the German Meteorological Service (Deutscher Wetterdienst – DWD). We collected data for average annual temperature and rainfall as well as summer and winter temperature and rainfall.

Measurement of pesticide content from archived leaf material

Request a detailed protocol

The ESB sampling was historically set up as a tool for pollution assessment (Schulze et al., 2007), and the samples are therefore stored to preserve any possible pollutant. Because these leaves serve as a habitat for the associated arthropods, such well-preserved material allows us to explore levels of chemical pollution occurring directly within the arthropods’ environment. Our samples were screened for pesticides and persistent organic pollutants with a modified QuEChERS approach (Löbbert et al., 2021). 2.0 g of sample material were extracted with acetonitrile (10 ml) and ultrapure water (10 ml), followed by a salting-out step using magnesium sulfate, sodium chloride and a citrate buffer (6.5 g; 8:2:3 [wt/wt]). After a dispersive solid-phase extraction cleanup step with magnesium sulfate, PSA (primary-secondary amine), and GCB (graphitized carbon black) (182.5 mg; 300:50:15 [wt/wt]), the supernatant was analyzed with a sensitive gas chromatography–tandem mass spectrometry (GC–MS/MS) instrument. All samples were analyzed for 208 GC-amenable compounds of different pesticide and pollutant classes, including pyrethroids, organochlorine and organophosphate pesticides, and polychlorinated biphenyls.

Molecular processing

DNA isolation

Request a detailed protocol

We developed a highly standardized protocol for the analysis of leaf-associated arthropod community DNA. We optimized various protocol steps to ensure the reliability and reproducibility of our data. We first explored the effect of DNA extraction on recovered diversity.

As mentioned above, the cryo-homogenization of ESB samples ensures a very homogeneous distribution of even trace amounts of chemicals in the sample. This should also hold true for DNA. Hence small subsamples of large homogenate samples should suffice for analysis. To test this hypothesis, we first performed a weight series extraction from 50, 100, 200, 400, 800, and 1600 mg of homogenate with several replicates for each weight. Additional extraction replicates of 16 beech samples at 200 mg were also included. This analysis confirms the pronounced homogenization of the samples, with 200 mg of homogenate sufficing to accurately recover α- and β-diversity (Figure 2—figure supplement 1A, B).

A single DNA extract was made from each ESB and freeze-dried sample, using the Puregene Tissue Kit according to the manufacturer’s protocols (Qiagen, Hilden, Germany). All samples were processed under a clean bench and kept over liquid nitrogen during processing to prevent thawing. Samples were transferred using a 1000-µl pipette with cutoff tips. The resulting wide bore tips were used to drill out cores of defined sizes from the leaf powder. To remove undesired coprecipitates, we performed another round of purification for each sample, using the Puregene kit following the manufacturer’s protocol. The hand-collected arthropod specimens from our branch clipping were pulverized in a Qiagen Tissuelyzer at 200 Hz for 2 min using new 5-mm stainless steel beads, and DNA was extracted from the pulverized samples using the Puregene kit as described in de Kerdrel et al., 2020. Branch clipping and leaf samples were processed in separate batches and using separate reagents to avoid possible carryover between sample types.

Primer choice, PCR amplification, and sequencing

Request a detailed protocol

As the standard DNA barcode marker for arthropods (Andújar et al., 2018), the mitochondrial COI gene offers the best taxonomic identification of German arthropod species. We thus selected COI for our metabarcoding analysis. We tested several primer pairs to optimize recovery of arthropod DNA from the leaf homogenates (Gibson et al., 2014; Leray et al., 2013; Jusino et al., 2019). Leaf homogenates are dominated by plant DNA, with arthropod eDNA only present in trace amounts. The majority of commonly used arthropod metabarcoding primers are very degenerate and will readily amplify mitochondrial COI of plants. Thus, for such degenerate primers, the vast majority of recovered reads will belong to the plant. We found an ideal tradeoff between suppressing plant amplification while still recovering a taxonomically broad arthropod community in the primer pair ZBJ-ArtF1c/ZBJ-ArtR2c (Zeale et al., 2011). While this primer pair is known to have taxonomic biases for certain arthropod groups (Piñol et al., 2015), it is still widely used as an efficient and reliable marker for community analysis (Thomsen and Sigsgaard, 2019; Eitzinger et al., 2021). Recently, we designed a novel and highly degenerate primer pair by modifying two degenerate metabarcoding primers (Gibson et al., 2014; Leray et al., 2013), which allows the suppression of plant amplification (NoPlantF_270/mICOIintR_W; Supplementary file 1; Krehenwinkel et al., 2022). To ensure the reproducibility of the diversity patterns recovered from our original ZBJ-ArtF1c/ZBJ-ArtR2c dataset, we additionally processed eleven complete ESB time series (174 samples) for this novel primer pair and compared results for species composition, α- and β-diversity. This analysis supports very similar patterns for temporal species abundance trends, as well as α- or β-diversity for both primer pairs (Figure 2—figure supplement 5).

All PCRs were run with 1 µl of DNA in 10 µl volumes, using the Qiagen Multiplex PCR kit according to the manufacturer’s protocol and with 35 cycles and an annealing temperature of 46°C. A subsequent indexing PCR of five cycles at an annealing temperature of 55°C served to attach sequencing adapters and 8-bp dual indexes (all with a minimum 2 bp difference) to each sample (using the layout described in Lange et al., 2014). We had previously tested the effect of DNA extraction and PCR replicates, which showed well correlated and reproducible OTU composition (R2extract1vs2 = 0.90, R2PCR1vs2 = 0.97, LM p < 0.05), as well as α-diversity patterns (R2extract1vs2 = 0.93, R2PCR1vs2 = 0.90, LM p < 0.05). PCR and extraction replicates also recovered a significantly lower β-diversity than within- and between-site comparisons (βPCR1vs2 = 0.14; βextract1vs2 = 0.19; βwithin site = 0.62; Pairwise Wilcoxon Test, p < 0.05) (Figure 2—figure supplement 1). Considering the similar recovered diversity patterns for PCR and extraction replicates and the observed saturation of diversity with single 200 mg homogenate extractions, we did not perform extraction replicates, but ran all PCRs in duplicate as technical replicates. To assure that dual replicates suffice to recover diversity patterns, we also added a PCR triplicate for two time series of beech samples (42 samples in total) and sequenced each of these samples to 78,000 reads on average. Patterns of diversity were highly correlated between duplicate and triplicate datasets (R2 = 0.95). The final libraries were quantified on a 1.5% agarose gel and pooled in approximately equal abundances based on gel band intensity. The final pooled sample was cleaned using 1× Ampure Beads XP (Beckmann-Coulter, Brea, CA, USA) and then sequenced on an Illumina MiSeq (Illumina, San Diego, CA, USA) using several V2 kits with 300 cycles at the Max Planck Institute for Evolutionary Biology in Plön, Germany. Branch clipping samples were amplified and sequenced separately using the above protocol. Negative control PCRs and blank extraction PCRs were run alongside all experiments and sequenced as well, to explore the effect of possible cross-contamination or index carryover between samples.

Test for DNA carryover in the cryomill

Request a detailed protocol

The sample processing pipeline of the ESB is laid out to be sterile and entirely avoid cross-contamination between samples. To test the efficiency of these protocols for eDNA sampling, we included a test on the possibility of carryover in the cryomill. Using the milling schedule of tree samples from 2015 to 2018, we compared the β-diversity between tree samples that were processed in the cryomill consecutively. Assuming an eDNA carryover takes place, the β-diversity should be significantly reduced compared to samples which are processed in different years. We did not find an effect of processing order in the cryomill on beta diversity for 18 within- and between-tree species comparisons (Figure 2—figure supplement 3). We also explored the effect of single species carryover in the cryomill. This was done using samples of different tree species, which were processed consecutively in the cryomill. We compared the read abundances of the 10 most abundant monophagous species with that found in the consecutively processed sample of a different tree species. The comparisons were done for one poplar and beech sample as well as one pine and spruce sample. To ensure even minor carryover would be detected, we sequenced all samples to a high depth of 78,000 reads on average. Yet, no signal of carryover was observed (Supplementary file 3).

Sequence processing

Request a detailed protocol

Reads were demultiplexed by dual indexes using CASAVA (Illumina, San Diego, CA, USA) allowing no mismatches in indexes. Demultiplexed reads were merged using PEAR (Zhang et al., 2014) with a minimum overlap of 50 and a minimum quality of 20. The merged reads were then quality filtered for a minimum of 90% of bases >Q30 and transformed to fasta files using FastX Toolkit (Gordon and Hannon, 2010). Primer sequences were trimmed off using sed in UNIX, with degenerate sites allowed to vary and only retaining sequences beginning with the forward and ending with the reverse primer. The reads were then dereplicated using USEARCH (Edgar, 2010). The dereplicated sequences were clustered into zero radius OTUs (hereafter zOTUs) using the unoise3 command (Edgar, 2016) and 3% radius OTUs using the cluster_otus command in USEARCH with a minimum coverage of 8 and a minimum occurrence of three reads in a sample. Chimeras were removed de novo during OTU clustering. All resulting sequences were translated in MEGA (Takahara et al., 2012) and only those with intact reading frames were retained. To assign taxonomic identity to the zOTU sequences, we used BLASTn (Altschul et al., 1990) against the complete NCBI nucleotide database (downloaded February 2021) and kept the top 10 hits. Sequences were identified to the lowest possible taxonomic level, with a minimum of 98% similarity to classify them as species. All non-arthropod sequences were removed. We then built Maximum Likelihood phylogenies from alignments of the zOTU sequences for all recovered arthropod orders separately using RaxML (Stamatakis, 2014). These phylogenies were used to perform another clustering analysis using ptp (Zhang et al., 2013) to generate OTUs from the data. Due to the well-developed German Barcode of Life database (Geiger et al., 2017), actual species identity can be reliably inferred by database comparisons for many arthropod groups. 3% radius OTUs often oversplit species, for example several 3% radius OTUs comprised one actual species. The ptp clustering often merged several 3% radius OTUs, but came closest to the actual species assignments by BLAST. Moreover, the recovered diversity values for 3% OTUs, zOTUs and ptp-based OTUs were well correlated (R2 = 0.90). We thus proceeded to use ptp-based OTUs (hereafter referred to as OTUs) for subsequent analysis on taxonomic diversity, as it should best approximate actual species diversity. zOTUs represent individual haplotypes in the dataset and were used as an indicator of genetic diversity. Using the taxonomic assignments, we estimated which taxonomic groups were particularly well represented in our data. Each tree species likely harbors a unique arthropod community with numerous monophagous species, a majority of which should be recovered by a broadly applicable molecular method. Where possible, we performed a finer scale ecological assessment for the recovered taxa, classifying them by trophic ecology and expected position on the outside or inside of the leaf. For example, mining taxa would likely be recovered from the inside of the leaf, while other taxa likely reside on the leaf’s surface.

Detection of relative arthropod DNA copy number using qPCR

Request a detailed protocol

Initial reports on insect decline were entirely based on biomass (Hallmann et al., 2017). Biomass, however, does not necessarily predict diversity (Gough et al., 1994). We therefore aimed to generate information not only on diversity, but also on relative biomass of arthropods in tree canopies. Previous eDNA studies show that DNA copy number is correlated with the biomass of a target taxon (Takahara et al., 2012), making qPCR a possible approach for biomass estimation. We developed a qPCR protocol to detect relative abundance of arthropod DNA copy number in leaf samples, using the plant DNA copy number as an internal reference for quantification. We used the nuclear 18SrRNA gene (hereafter 18S). Although 18S can show interspecific copy number variation, it provides relatively good approximations of actual taxon abundances in amplicon assays (Krehenwinkel et al., 2017; Krehenwinkel et al., 2019b). Primer pairs targeting plants and arthropods were designed to meet the following criteria: (1) Identical PCR fragments should be amplified for plants and arthropods so that PCR for both taxa will perform similarly. (2) The arthropod-specific primer should not amplify plants, and vice versa. (3) Fungi should be excluded from amplification, as DNA of fungal endophytes is probably at least as abundant in leaf samples as arthropod eDNA. (4) The primers should target conserved regions in order to amplify a broad spectrum of plants or arthropods. We used diagnostic SNPs at each primer’s 3′-end to achieve the lineage specificity (Krehenwinkel et al., 2019a).

Two possible qPCR primer pairs were designed, one targeting a 172 bp and the other a 176 bp fragment of 18S (Supplementary file 1). The first primer pair contained a 3′-AA-mismatch discriminating arthropods from plants and fungi in the forward primer and a 3′-TT-mismatch discriminating plants from arthropods in the reverse primer, while the second pair had the same 3′-AA-mismatch in the forward primer but no reverse primer mismatch (Figure 3A). To test the lineage specificity of both primer pairs, we performed an amplicon sequencing experiment with the arthropod-specific primers. Four samples from each of the four tree species were amplified with both primer pairs, then indexed, pooled, sequenced, and processed as described above. All reads were clustered into 3% radius OTUs. Taxonomy was assigned to the OTUs, and an OTU table was built, as described above for the ESB sample metabarcoding experiment. The proportion of arthropod reads was then estimated for each sample and primer pair. The first primer combination (3′-mismatch in both forward and reverse primers) led to a near complete suppression of plant and fungal amplification in all tested samples and was therefore used for the qPCR (Figure 3A, B).

To account for the low quantity of arthropod DNA in relation to plant DNA, we used a nested qPCR assay, with a high accuracy for low DNA copy numbers (Tran et al., 2014). The sample was first amplified in a regular PCR with 15 cycles using the Qiagen Multiplex PCR kit. Two separate PCRs were run: one using the arthropod primers with an undiluted DNA extract as template, and the other using the plant primers with a 1:100 dilution of the DNA extract. The primers included a 33-bp forward and 34-bp reverse tail, based on Illumina TruSeq libraries, which were complementary to sequences in the qPCR primers. After being cleaned of residual primers with 1× AMPure beads XP, the products of the first PCR were used as template in the qPCR. qPCR was run with the Power SYBR Green Mastermix (Fisher Scientific, Waltham, MA, USA) on an ABI StepOnePlus Real-Time PCR System (Applied Biosystems, Foster City, CA, USA) according to the manufacturer’s protocol, using 35 cycles and an annealing temperature of 55°C. All reactions were run in triplicate and the average of the three CT values used for analysis. CT values showed high reproducibility between triplicate PCRs (SDplant = 0.13; SDarthropod = 0.10). Non-template controls were run alongside all qPCRs to rule out contamination.

The qPCR efficiency was estimated for the plant and arthropod-specific marker using two ESB tree samples. A 10,000-fold dilution series was used for efficiency estimation. This assay was very stringent, as it corresponds to a dilution of the naturally occurring arthropod eDNA in a plant sample by 10,000. Both assays showed a high efficiency across the dilution series (EPlant = 94.77%, EArthropod = 99.73%). The 1:10000 dilution (CTplant 1.1000 = 15.4; CTinsect 1.1000 = 27.2) is far less than the actual amount of insect DNA in an ESB sample (average CTinsect = 21.3; average CTplant = 15.7), supporting the reliability of our experiment.

To estimate the accuracy of relative arthropod DNA copy detection across diverse arthropod communities, we also performed a spike-in assay, in which a dilution series (10,000-fold) of arthropod mock community DNA was added to a leaf extract and analyzed using qPCR. Seven mock communities were prepared, each containing varying amounts of DNA from 13 arthropod species representing 13 different orders (Figure 3—figure supplement 1). The relative copy number of arthropod DNA in relation to plant DNA was estimated using the Delta CT Method (Schmittgen and Livak, 2008). The optimized qPCR protocol was then used to quantify the relative DNA copy number of arthropods in all 312 ESB leaf samples.

Statistical analysis

Request a detailed protocol

Using USEARCH, an OTU table was built including all samples with the taxonomically annotated zOTU sequences as reference. A species-level OTU table was then generated by merging the zOTUs into their respective ptp-based OTU clusters.

The negative control samples were mostly free of arthropod sequences. We found 1.88 arthropod reads per control on average (0–5 reads per control). The recovered reads belonged to taxa that were highly abundant in one of the analyzed tree species, suggesting minor carryover during PCR or sequencing. Based on the negative control samples, we removed all entries in the OTU table with fewer than three reads to counter this possible carryover.

We used two approaches to rarefy our OTU table. First, using rarefaction analysis in vegan (Oksanen et al., 2013) in R (v 4.1.0) (Team RS, 2015), we explored saturation of diversity. Based on this analysis, 5000 reads were randomly sampled for each of the two PCR duplicates using GUniFrac (Chen et al., 2018), and the duplicates were merged into a final sample of 10,000 reads (Figure 2—figure supplement 4) after the replicates were checked for reproducible patterns of species composition, α- and β-diversity. To ensure that undersampling did not affect our results, we performed an additional analysis with the unrarefied dataset, which yielded an average coverage of 21,676 reads per sample. A second rarefaction was informed by the relative copy number of arthropod DNA recovered from the tree samples with our qPCR assay. Assuming the copy number reflects biomass, we sampled read numbers proportional to the specific relative copy number for each sample. The copy number-informed and unrarefied datasets showed highly correlated diversity patterns with the dataset rarefied to 10,000 reads (R2 ≥ 0.91).

Taxonomic α- and β-diversity were calculated in vegan in R. Quantitative biodiversity assessments by metabarcoding at the community level are likely biased (Krehenwinkel et al., 2017). We therefore limited our assessments of α-diversity to richness and β-diversity to binary dissimilarity.

We also measured temporal abundance changes of single OTUs within sites. Within OTUs, temporal changes in read abundance at a site should reflect the relative abundance with reasonable accuracy (Krehenwinkel et al., 2017). Only sites spanning a minimum time series of 10 years were included, and we only used species that occurred in at least three sampling events for a particular site and for which at least 100 reads were recovered. This filtering served to exclude rare species, which imitate abundance increases or declines by randomly occurring early or late in the time series. To account for likely fluctuations in abundance, we used the log + 1 of read abundance. Significant increases or declines of abundance over time were estimated for each OTU and site using non-parametric Spearman correlation in R. Both the qPCR-informed and the rarefied datasets were used to calculate species abundance changes.

As mentioned above, zOTUs represent individual haplotypes and thus genetic variation within species. Using the zOTU data, we calculated the haplotype (zOTU) richness within each individual OTU as a complementary measure for genetic variation. A decline in biodiversity could manifest itself in an overall loss of species, which should be detectable at the OTU level. Alternatively, biodiversity decline could initially only affect genetic variation within species, for example, be the result of declining population sizes without actual extinctions. This should be detectable by losses of overall zOTU diversity and zOTU richness within single OTUs. To derive within-OTU genetic diversity, we identified OTUs that consisted of more than one zOTU. zOTU variation could be affected by low abundance sequence noise. Hence, we used the same filtering criteria as described above for the species abundance change. Moreover, we only included zOTUs in our calculations that were present in both technical replicates of a sample. The richness of the remaining zOTUs within each of these OTUs was then calculated.

Arthropods are an ecologically very heterogeneous taxon, with different groups showing very different life histories and possibly responses to ecosystem change. To account for this heterogeneity, we calculated α- and β-diversity metrics for the complete arthropod dataset, as well as the 10 most common arthropod orders in the dataset. Using NMDS (k = 2, 500 replications, Jaccard dissimilarity) in vegan, we visualized differentiation of the recovered arthropod communities. We then tested for effects of tree species, sampling year, sampling site, land use type, weather before sampling, amount of leaf material in a sample, climatic variables and detected pesticide load on α- and β-diversity. Also, hand-collected branch clipping samples, as well as the freeze-dried samples, were compared with the ESB samples for their recovered arthropod community composition and diversity.

Factors contributing to β-diversity were evaluated using a PERMANOVA in vegan. To evaluate an association of community turnover and copy number variation (e.g., biomass), we also explored patterns of association between these two variables for each site (Figure 5—figure supplement 2). Statistical analysis for temporal changes of α-diversity and relative copy number were performed using the nlme (v 3.1-159; 2022) (Bates et al., 2014) package in R. LMMs were applied to analyze the statistical importance of involved fixed and random effects. Temperature (annual, summer, and winter temperatures), corresponding rainfall data, year, and land use type were treated as fixed effects. Site ID was included as random effect. The Akaike information criterion was used in stepwise regression to identify the final models. A continuous-time first-order autocorrelation model term (corCAR1) was included in the lme function to account for serial autocorrelation. The dominant predictor variable is the tree species, which always contributes most to the marginal R2 in all models. Arthropod DNA copy number (marginal R2 = 0.26) showed negative associations with time (p < 0.001) and winter temperatures (p = 0.038). None of the other diversity metrics (OTU and zOTU richness, copy number-corrected richness, saturated richness, and genetic diversity) showed an association with time. However, all richness values were positively correlated to winter temperatures and negatively to summer temperatures (marginal R2 = 0.45–0.48, p < 0.05). Genetic diversity (marginal R2 = 0.34) was correlated to winter rainfall (p < 0.001).

Data availability

All raw reads are available in the Dryad Digital Repository (https://doi.org/10.5061/dryad.x0k6djhmp). The OTU table with metadata and qPCR results has been uploaded as supplementary files.

The following data sets were generated
    1. Krehenwinkel H
    (2022) Dryad Digital Repository
    eDNA from archived leaves reveals no losses of α-diversity, but widespread community turnover and biotic homogenization as drivers of forest insect decline.
    https://doi.org/10.5061/dryad.x0k6djhmp

References

    1. Collen B
    2. McRae L
    3. Deinet S
    4. De Palma A
    5. Carranza T
    6. Cooper N
    7. Loh J
    8. Baillie JEM
    (2011) Predicting how populations decline to extinction
    Philosophical Transactions of the Royal Society of London. Series B, Biological Sciences 366:2577–2586.
    https://doi.org/10.1098/rstb.2011.0015
    1. Delvare G
    (1997)
    A review of methods for sampling arthropods in tree canopies
    Canopy Arthropods 27:52.
    1. Geiger MF
    2. Rulik B
    3. Waegele WW
    (2017)
    Overview on the activities in the German barcode of life project phase II
    Genome 60:936–938.
  1. Software
    1. Gordon A
    2. Hannon GJ
    (2010)
    FASTX-TOOLKIT, version 0.14
    Computer program and documentation distributed by the author.
    1. Magurran AE
    2. Henderson PA
    (2010) Temporal turnover and the maintenance of diversity in ecological assemblages
    Philosophical Transactions of the Royal Society of London. Series B, Biological Sciences 365:3611–3620.
    https://doi.org/10.1098/rstb.2010.0285
  2. Book
    1. Taberlet P
    2. Bonin A
    3. Zinger L
    4. Coissac E
    (2018)
    Environmental DNA: For Biodiversity Research and Monitoring
    Oxford University Press.
  3. Website
    1. Team RS
    (2015) RStudio: Integrated development for R
    RStudio. Accessed March 12, 2015.

Decision letter

  1. Simon Creer
    Reviewing Editor; Bangor University, Bangor, United Kingdom
  2. Detlef Weigel
    Senior Editor; Max Planck Institute for Biology, Germany
  3. Rafael Valentin
    Reviewer; Princeton University, United States
  4. Thomas Gilbert
    Reviewer

Our editorial process produces two outputs: (i) public reviews designed to be posted alongside the preprint for the benefit of readers; (ii) feedback on the manuscript for the authors, including requests for revisions, shown below. We also include an acceptance summary that explains what the editors found interesting or important about the work.

Decision letter after peer review:

Thank you for submitting your article "Environmental DNA from archived leaves reveals widespread temporal turnover and biotic homogenization in forest arthropod communities" for consideration by eLife. Your article has been reviewed by 3 peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Detlef Weigel as the Senior Editor. The following individuals involved in the review of your submission have agreed to reveal their identity: Rafael Valentin (Reviewer #1); Thomas Gilbert (Reviewer #3).

The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this to help you prepare a revised submission.

My colleagues have performed an excellent assessment of the breakthroughs and importance of this novel study and I will not paraphrase their valuable insights and feedback but will also join in the discussion.

It can be frequently observed that low diversity environmental DNA samples can be under sequenced using high throughput sequencing, whereas high diversity samples normally yield higher levels of sequence rates. By rarefying their analysis to 5000 reads (an approach that is increasingly not used, c.f. the McMurdy and Holmes debate https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1003531), there is a risk that low diversity samples may be over-sampled and high diversity samples may be under-sampled. I appreciate that 5000 data points are likely to be appropriate when capturing leaf-based arthropod communities, but further evidence of the rarefaction, per sample coverage and rationale, would be valuable.

One of the conclusions of the research, according to qPCR of the 18S ribosomal DNA marker, is that (cellular) biomass appears to have decreased over time; a narrative that is coherent amongst the declining population trends that prevail in contemporary ecology. The challenge here is the use of the degenerate 18S marker, which is not linked to the interspecific copy number variation that will be prevalent amongst the different taxa in the study (e.g. taxon 1, 5 copies; taxon 2, 15 copies; taxon x, y copies…). A potentially equally plausible explanation for the declining copy number is that the copy number of the homogenised communities is lower than the older, non-homogenised communities. Without exhaustive per taxon calibration, the inability to compare interspecific abundance between taxa is one of the downsides of any metabarcoding study, that is focused on multicopy markers that differ either in copy number, or the amount of tandem repeats, as is the case with mitochondrial and nuclear ribosomal markers respectively. I suggest therefore that the investigators devise an additional test to measure, e.g., the level of association between qPCR copy number, time, and the level of homogenization (e.g. Jaccard dissimilarity) simultaneously in order to address this concern head-on.

In addition to my insights above, I will highlight the revisions that have been highlighted as essential amongst the review panel.

Essential revisions:

1) The biomass assay – Explore whether the decrease in qPCR (aka here as biomass) values is associated with homogenization, time, or other identifiable factors simultaneously (e.g. leaf characteristics) and discuss the results in an objective fashion.

2) The data have been rarefied to 5000 reads – is this enough and can the team clearly demonstrate their baseline data strategy? Would anything change by using either proportional data or using read depth as an offset on the total dataset (e.g. https://www.nature.com/articles/s42003-020-01562-4)?

3) Clarity on the nature of controls, decontamination procedures, and how these guided data filtering in the metabarcoding study.

4) Revisit the simplistic modelling and devise appropriate analyses that will take into consideration non-linear trends and different factors where appropriate.

5) Justify the replicability/robustness of the qPCR strategy.

6) Justify OTU picking strategy.

I would also invite the authors to address all the reviewer's comments via standard rebuttal letter and submission of track changed version of the original submitted text where possible. Where reviewers have different opinions, please discuss clearly your views in the context of your study and the evidence presented (cf. discussion of read numbers, qPCR, and abundance).

Reviewer #1 (Recommendations for the authors):

While reading through the methods I was at first deeply concerned that PCR errors remained within the zOTU table, but only later (in the statistical analysis section) did I see that OTUs with three or fewer reads were removed from the dataset. This should be moved up with the rest of the bioinformatic information to make it clear earlier and alleviate that concern.

I saw no mention of methods that would filter out reads, or samples, due to contamination levels. There is always some level of contamination that takes place in eDNA metabarcoding, and addressing this is paramount. Filtering out contaminants can range from removing contaminant reads in technical replicates to removing technical replicates entirely. This extends to the decision to use just two technical replicates for the metabarcoding portion of the study. While not bad, should technical replicates be filtered out it leaves just a single technical replicate to represent the sample, which isn't sufficient?

On line 185 only 413 OTUs were said to be used to assess temporal changes. Why were just these OTUs selected from the larger dataset? This should be explained and justified.

Were beads for the tissuelyzer and cryomil reused? If so, how were they decontaminated to ensure no contamination of subsequent samples took place? If they were not reused this should be made clear in the text.

What was the justification for a 3% OTU radius? Was this a precedent set from another paper, or was it a random selection? More detail here is needed.

Reviewer #2 (Recommendations for the authors):

I elaborate on my concerns regarding the sequence and statistical analyses below.

One of my greatest concerns is in the statistical analyses of the temporal trend. In most cases, the authors used linear models (LM) to test whether there is a temporal trend (i.e., increase or decrease). This is inappropriate because empirical time series often contain temporal autocorrelation structures and because the application of LM to such time series may often result in the false-positive detection of the trend. This means that LM can often detect a "significant" trend even in a random walk time series. This is a well-known issue, and the authors should apply a more appropriate statistical method, e.g., the autocorrelation model, state-space modeling, or some other methods, to judge whether there is a temporal trend.

Another concern is the use of read abundance as an explained variable in the statistical models. In L452-454, the authors claimed that the DNA copy numbers may be a proxy of biomass. I agree with this statement. However, sequence read abundance is not the DNA copy numbers and cannot be a proxy of abundance in most cases. If I understood correctly, the authors rarefied the sequence reads to 5000 reads for each sample, and the read abundance in the statistical model seems to be the relative abundance (e.g., if an OTU produced 500 reads in a sample, the relative abundance of the OTU in the sample is 10%). The authors quantified 18S DNA copy numbers of arthropods, so multiplying the relative abundance by the total 18S DNA copy numbers may produce a better proxy of the abundance of arthropods. I would recommend the authors reconsider the read abundance issue carefully.

Reviewer #3 (Recommendations for the authors):

Line 91. 312 ESB samples. Perhaps expand here on distribution through time etc? Could be part of Figure 1?

Line 125+ – stability of arthropod DNA comment. I'm actually surprised there is no higher diversity in the 8-year RT stored samples, as a result of post-mortem modifications to the DNA. Having said that does the PCR strategy negate that? If yes maybe worth stating this clearly. Actually in general given how critical the assumption that any zOTUs used represent true biological sequence variation (as opposed to PCR error, sequencing error, post mortem DNA damage) I advise the authors to early in the paper make it clear why they believe their zOTU estimates to be accurate.

Line 129. rDNA is typically used for many things. Ribosomal DNA. recombinant DNA. Even relative DNA. Maybe for clarity here state 'information on relative arthropod 18s rDNA copy' then there is no uncertainty.

Figure 5F's label is confusing (e.g. Does < 2020 mean 2010-2020 or all years before 2020 etc?). I suspect the former but that's not how it reads as written in the figure.

Can the authors elaborate on when the leaves were collected in the year? In line 287 the authors state it was the same time point every year. How is that shaped by the time of year the leaves come out every year? (Which must fluctuate a lot with annual weather variation). Also in line 257, the authors state the leaves represent a 'fairly broad phenological window'. Can this be elaborated on?

Line 287/288. It seems odd to write 'a defined amount' of both leaf and number of trees, and then give values prefixed by '>'. And then say 'defined number of branches from each tree' but then not give a value. I think I know what the intended meaning is, but perhaps consider rewording this sentence! Also, I assume the exact numbers are somewhere – maybe refer to where at this point?

Methods in general, it took me a long time to work out exactly how the samples were collected and processed, and I'm still not sure. If I understand it correct – leaves were originally picked from trees immediately into liquid nitrogen, and then shortly after ground to a powder? Thus all happening years/decades ago? If correct please consider re-reading the text assuming you know nothing about the history and add any small clarifications that might be needed so the reader can immediately jump to this conclusion. If however, I misunderstand…then again the text needs clarification.

Line 295 I assume the cryomill was somehow decontaminated between sample batches? Please clarify how.

Ca Line 365. Please state clearly if extraction replicates were made on the samples used in the full experiment, (or not). It's not clear to me as written now. I can see they were used in the sample mass trial and perhaps given the replicate dissimilarity results were ok for the PCR and extraction replicate the decision was taken to not do replicates on the large scale? If so please state it clearly.

https://doi.org/10.7554/eLife.78521.sa1

Author response

My colleagues have performed an excellent assessment of the breakthroughs and importance of this novel study and I will not paraphrase their valuable insights and feedback but will also join in the discussion.

It can be frequently observed that low diversity environmental DNA samples can be under sequenced using high throughput sequencing, whereas high diversity samples normally yield higher levels of sequence rates. By rarefying their analysis to 5000 reads (an approach that is increasingly not used, c.f. the McMurdy and Holmes debate https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1003531), there is a risk that low diversity samples may be over-sampled and high diversity samples may be under-sampled. I appreciate that 5000 data points are likely to be appropriate when capturing leaf-based arthropod communities, but further evidence of the rarefaction, per sample coverage and rationale, would be valuable.

We have now expanded on this in the Methods, added more data and performed some additional analyses on the read coverage question. 5000 reads from two replicates each per sample should indeed achieve saturation of diversity in this case (Figure 2—figure supplement 4). We went for an equal coverage approach for all samples, which means the coverage for each sample was based on the lowest covered sample. Hence, most samples have considerably more reads than 5000. To address potential biases of this approach, we repeated the analysis with unrarefied data, which led to 21,676 reads per replicate sample on average. The results essentially remain unchanged (see Methods: Statistical Analysis; Suppl. File 2), with diversity of the two datasets highly correlated (R2=0.95). This analysis can be repeated by all readers based on the raw read data, which allows for much deeper coverage. We have included this in the Methods now (lines 581-588).

In addition, we explored the effect of including triplicate samples and an even deeper coverage in two complete timeseries of 42 beech samples. Here, we sequenced triplicates of every sample, resulting in an average coverage of 78,667 reads per sample. Again, the result remained unchanged compared to the reduced coverage dataset (lines 435-439).

In addition, we generated a second dataset for most samples with another primer pair, which also confirmed our findings from the first dataset. The fact that the second dataset recovered essentially identical patterns of α and β diversity also suggest the validity of our dataset (Figure 2—figure supplement 5).

We liked the recommendation by Reviewer 1 to combine data from qPCR and metabarcoding for a copy number-based rarefaction. We have also included this analysis now and find no change to our results (lines 583-588).

Considering this background, we are confident that the result we aim to achieve, namely testing for a temporal change of richness or arthropod community turnover, is not affected by our read sampling strategy.

One of the conclusions of the research, according to qPCR of the 18S ribosomal DNA marker, is that (cellular) biomass appears to have decreased over time; a narrative that is coherent amongst the declining population trends that prevail in contemporary ecology. The challenge here is the use of the degenerate 18S marker, which is not linked to the interspecific copy number variation that will be prevalent amongst the different taxa in the study (e.g. taxon 1, 5 copies; taxon 2, 15 copies; taxon x, y copies…). A potentially equally plausible explanation for the declining copy number is that the copy number of the homogenised communities is lower than the older, non-homogenised communities. Without exhaustive per taxon calibration, the inability to compare interspecific abundance between taxa is one of the downsides of any metabarcoding study, that is focused on multicopy markers that differ either in copy number, or the amount of tandem repeats, as is the case with mitochondrial and nuclear ribosomal markers respectively. I suggest therefore that the investigators devise an additional test to measure, e.g., the level of association between qPCR copy number, time, and the level of homogenization (e.g. Jaccard dissimilarity) simultaneously in order to address this concern head-on.

We agree that we cannot fully rule out the possibility that a change in copy number is due to shifts in community composition, and we have already addressed this possibility in the discussion (lines 257-260).

The 18S we used here is indeed a multi-copy marker, and shifts in community composition could also lead to shifts in copy number. This could mimic the pattern of copy number loss over time. However, we have tested our assay quite extensively. The test included diverse mock communities in a very pronounced dilution series (1:10,000), which still suggested that it predicts copy number changes accurately and with a high efficiency (Figure 3). The 1:10,000 dilution we used to estimate qPCR efficiency is also far lower than arthropod copy number in any actual ESB sample, yet with this assay we accurately predicted the amount of arthropod DNA. We also have found 18S to reflect coarse changes to abundance quite well in two previous studies, which we cited in the manuscript (Krehenwinkel et al. 2017 and 2019).

A strong change in copy number due to change in taxonomic composition would likely be found if turnover affects divergent taxa with very different copy number. This is contrary to our findings. We find that usually closely related taxa replace each other. E.g., a mining moth gets replaced by another one (see Figure 5—figure supplement 1).

The test for an association of homogenization and/or turnover with copy number is a very good idea, and we have now implemented this in our analysis. We have explored the correlation between community dissimilarity and the difference in copy number between sampling years (Figure 5—figure supplement 2). With the exception of one urban site, we do not find a clear association. In particular the beech forest sites, which show increasing homogenization over time, do not show an association of copy number differences and turnover. We thus believe that the decrease in copy number is better explained by a loss of biomass than simple turnover. This is also in line with recent work on insect decline, suggesting that abundance losses of common species can drive insect decline. However, we still cannot rule it out entirely and hence left the discussion of the topic (lines 257-260).

On a side note: Just last week I saw a very interesting presentation at a conference supporting our findings. The authors (Samu et al.) found stable richness of forest spider communities over time, but a loss of biomass as a general pattern for the studied spider communities.

In addition to my insights above, I will highlight the revisions that have been highlighted as essential amongst the review panel.

Essential revisions:

1) The biomass assay – Explore whether the decrease in qPCR (aka here as biomass) values is associated with homogenization, time, or other identifiable factors simultaneously (e.g. leaf characteristics) and discuss the results in an objective fashion.

We have now addressed this in the Methods (lines 626-629), Results (lines 197-198) and Discussion (lines 259-260). We did not find an association between copy number and β diversity. Assuming that turnover is driving the changes in copy number, we would find an association of β diversity and dissimilarity in copy number. With the exception of one site, this was not found (Figure 5—figure supplement 2).

Also, the observed turnover mostly affects closely related taxa. No community tipping point or pronounced change to the higher-level taxonomic compositions of the community is found in our data. Pronounced changes in copy number would rather be expected with replacement of distantly related taxa. This also makes copy number changes due to turnover less likely. Our mock community assay also shows that the qPCR assay quite accurately recovers copy number changes even in quite divergent communities.

We can also rule out changes in leaf characteristics as drivers of changes to relative copy number. The specimen bank collects biometric data for most its samples, including leaf weight, which has not changed over time in most sampling sites (lines 310-311).

Considering this background, we consider an actual loss of biomass a more likely explanation. Though we cannot fully rule out the effect of turnover, an overall loss of biomass affecting many common species would provide a simple explanation for the stable richness but drop of biomass.

2) The data have been rarefied to 5000 reads – is this enough and can the team clearly demonstrate their baseline data strategy? Would anything change by using either proportional data or using read depth as an offset on the total dataset (e.g. https://www.nature.com/articles/s42003-020-01562-4)?

We have now added several additional analyses and new data to underline the reliability of our data, which we address in the “Statistical analysis” part of the Methods (lines 581-588). 5,000 reads per replicate appear to suffice to saturate recovered diversity, as seen in this rarefaction curve (Figure 2—figure supplement 4). To make sure we did not undersample the communities, we repeated the analysis with unrarefied data, which led to 21,676 reads per replicate sample on average. The results essentially remain unchanged, with diversity highly correlated (R2=0.95). This analysis can be repeated by all readers based on the raw read data, which allows for deeper coverage than 5,000.

In addition, we explored the effect of including triplicate samples (which we newly generated) and an even deeper coverage in two complete timeseries of 42 beech samples. Here we sequenced triplicates of every sample, resulting in an average coverage of 78,667 reads per sample. Again, the result remained unchanged compared to the reduced coverage dataset (lines 435-439).

In addition, we generated a second dataset for most samples using a second primer pair, which also confirmed our findings from the first dataset. Using a completely different primer set most likely significantly increased the recovered taxonomic diversity. The fact that the second dataset recovered essentially identical patterns of α and β diversity also suggests the validity of our dataset (Figure 2—figure supplement 5).

We liked the idea of proportional sampling and have also incorporated a novel analysis on this (lines 583-588). We have used the relative copy number as a proxy to rarefy our dataset; i.e., the read coverage was based on relative arthropod copy number in a sample. Even incorporating the overall declining copy number in our rarefaction did not result in a temporal loss of richness for the studied communities, additionally underling the robustness of our data.

3) Clarity on the nature of controls, decontamination procedures, and how these guided data filtering in the metabarcoding study.

We have expanded on this in the Methods section of the manuscript, particularly under “Tree samples of the German Environmental Specimen Bank – Standardized time series samples stored at ultra-low temperatures”, “Test for DNA carryover in the cryomill” and “Statistical analysis”. We added details on extraction and PCR controls and more information on decontamination protocols assuring the sterility of the sampling procedure of the ESB (lines 303-304, 316). The ESB sampling is only performed with sterile equipment, which is thoroughly cleaned between sampling, making cross-contamination unlikely. The thorough cleaning and decontamination are also done for the cryomill.

We have now also included a new test to explore the possibility of contamination by the cryomill: “Test for DNA carryover in the cryomill” in the Methods section (lines 448-464). We have collected information on the milling schedule of samples processed between 2016 and 2019. We then identified tree samples that were processed consecutively in the cryomill. Assuming carryover is an issue, tree samples that were processed consecutively should show a significantly reduced β diversity between each other than between samples from the same site processed in different years. This effect was not found (Figure 2—figure supplement 3). We also tested single species carryover in more detail. This was done with a poplar and beech sample as well as a spruce and pine sample processed on two consecutive days. Assuming there is carryover in the mill, tree species-specific arthropod species should show signs of cross-contamination between the samples. This should only affect consecutively processed samples, but not samples from the same site processed in different years. The according samples were sequenced in triplicates to deep coverage (78,667 reads per sample on average) and the read abundances of the most common tree-species-specific arthropods measured. No sign of cross-contamination was found. This can be seen in the new (Suppl. File 3).

4) Revisit the simplistic modelling and devise appropriate analyses that will take into consideration non-linear trends and different factors where appropriate.

We have reworked statistical analysis now accounting for autocorrelation as reported under “Statistical analysis” in the Methods (lines 629-644). As a reviewer rightfully pointed out that abundance changes should not follow linear trends only, we have now used a more appropriate non-parametric test to explore species with significant increase or decline for the single species abundance change analysis. We have also used the qPCR-corrected read count for this analysis (Figure 5 A and B). The new analyses resulted in no overall change of the results.

5) Justify the replicability/robustness of the qPCR strategy.

We have now explained this in more detail in the methods, results and discussion, particularly under “Detection of relative arthropod DNA copy number using quantitative PCR” (lines 504-564). We based our idea of using nested qPCR on a paper (Tran et al. 2014), which we cite in the manuscript. This works shows that nested qPCR has a considerably higher sensitivity than standard qPCR, while retaining a very high accuracy. This is also what we found with extensive testing in our assay.

We tested the efficiency of the qPCR assay in a 1:10000 dilution series of two ESB tree DNA extracts. Both showed a very high efficiency from 1X to 1:10000 dilution. Changes in copy number of plants and insects are thus very reliably recovered by our assay. The CT values at a 1:10000 dilution are considerably higher (CT1:10000~30) than what we found in our experiments with the actual ESB extracts (average CT plant assay = 15.7; average CT insect assay = 21.3). Hence our assay is accurate and efficient at the whole range of arthropod DNA abundance found in our samples.

The accuracy of the assay is also shown in our mock community experiment, where spiked-in community DNA of seven highly distinct mock communities showed very comparable changes in copy number (Figure 3).

The technical replicates of our qPCR are very comparable with each other, showing a CT SD of ~0.10 on average. Also, all qPCRs were run along with a negative control to make sure contamination did not affect the outcome. Technical replicate CT values are now provided in the Suppl. File 2.

6) Justify OTU picking strategy.

We have explained our strategy in more detail in the Methods under “Sequence processing” (lines 488-494). We used three different approaches for OTU picking. The first two are implemented in USEARCH. Using cluster_otus, we generated 3 % OTUs; using unoise3, we generated zero radius OTUs corresponding to haplotypes. The first should approximate species richness, while the latter will allow us to explore genetic variation within species. To make the assessment of genetic variation more stringent, we have also included additional quality filtering criteria to call zOTUs, which we report on in the methods.

The German Barcode of Life for arthropods is quite complete, hence, we could ground-truth the species assignment by 3 % OTUs quite well. We found the 3% OTUs to over-split many species, e.g., a single species consists of various 3 % OTUs. We hence included ptp clustering (Zhang et al. 2013) as third OTU picking strategy. PTP clusters approached actual species assignments by BOLD much better than 3 % radius OTUs; hence, we used only zOTUs and PTP OTUs for further analysis. However, patterns of diversity were highly correlated for all three OTU picking strategies as we report in the Methods and in Suppl. File 2.

I would also invite the authors to address all the reviewer's comments via standard rebuttal letter and submission of track changed version of the original submitted text where possible. Where reviewers have different opinions, please discuss clearly your views in the context of your study and the evidence presented (cf. discussion of read numbers, qPCR, and abundance).

Reviewer #1 (Recommendations for the authors):

While reading through the methods I was at first deeply concerned that PCR errors remained within the zOTU table, but only later (in the statistical analysis section) did I see that OTUs with three or fewer reads were removed from the dataset. This should be moved up with the rest of the bioinformatic information to make it clear earlier and alleviate that concern.

This information has now been added to the “Sequence processing” section (lines 573-575).

I saw no mention of methods that would filter out reads, or samples, due to contamination levels. There is always some level of contamination that takes place in eDNA metabarcoding, and addressing this is paramount. Filtering out contaminants can range from removing contaminant reads in technical replicates to removing technical replicates entirely. This extends to the decision to use just two technical replicates for the metabarcoding portion of the study. While not bad, should technical replicates be filtered out it leaves just a single technical replicate to represent the sample, which isn't sufficient?

As mentioned above, we have expanded the section on data cleaning. We used both technical replicates after careful evaluation of contamination with the help of controls. We made sure both technical replicates showed predictable and comparable taxon composition and were not affected by contamination. We have now also included a dataset of triplicates (lines 435-439) and an unrarefied dataset (lines 581-588) to show that additional replication/deeper sequencing does not change the recovered biodiversity patterns (as reported under “Statistical analysis”). The fact that we used a second primer pair for our analysis, which also supports the same biodiversity patterns as our first primer pair, additionally gives us high confidence in the reliability of our data (Figure 2—figure supplement 5).

On line 185 only 413 OTUs were said to be used to assess temporal changes. Why were just these OTUs selected from the larger dataset? This should be explained and justified.

We address this issue in Methods: “Statistical analysis” (lines 594-598). As mentioned above, we excluded OTUs to make sure to only include reproducibly occurring taxa from the timeseries. Using rare OTUs, which for example only occur with a few reads or only a single time in our dataset, makes it more likely that we base our analysis on artefacts. A taxon occurring a single time in the time series may show a pattern of decline if it occurs only in the first datapoint of the timeseries. Also, a taxon represented by very few reads is more likely to originate from contamination. This is why we chose a minimum of 3 occurrences and 100 reads to score a taxon. We only included time series of 10 years or longer as changes in shorter time series may also be part of natural abundance fluctuations instead of an actual decline.

Were beads for the tissuelyzer and cryomil reused? If so, how were they decontaminated to ensure no contamination of subsequent samples took place? If they were not reused this should be made clear in the text.

The tissuelyzer beads were new and not reused, which we report now in the Methods (lines 394-396). The cryomill is a very large device, which uses titanium cylinders to grind samples. These cylinders are reused after thorough cleaning and decontamination, which is performed on all parts of the mill after each run. As mentioned above, we have now also included an analysis that shows there is no signal of between-sample carryover between milling events. We have added details to the text to clarify this (“Test for DNA carryover in the cryomill”, lines 448-464; Figure 2—figure supplement 3, Suppl. File 3).

What was the justification for a 3% OTU radius? Was this a precedent set from another paper, or was it a random selection? More detail here is needed.

We have now expanded on our OTU picking strategy in the text and explained why different OTUs were used (see ”Sequence processing”, lines 488-494). We used three different OTU picking strategies:

1. Zero radius OTUs or zOTUs, which represent haplotypic variation within species. We have now added a stringent filter to filter these zOTUs to calculate within-OTU diversity as a measure of genetic diversity.

2. 3 % radius OTUs are very commonly used in metabarcoding studies and usually assumed to represent species-level diversity. We initially also used 3 % radius OTUs for this purpose. The barcode reference library for German arthropods is very complete for many groups, hence we can reliably assign species status to many OTUs. But when comparing the 3% OTUs to the barcode reference library, we found that 3 % OTUs very commonly over-split species. E.g., a single species consisted of several OTUs. We have hence used a third clustering approach:

3. Ptp clustering commonly merged several OTUs into a single cluster. These clusters hence much better represented actual species diversity than 3 % OTUs. We therefore used ptp-based OTUs to represent species diversity.

However, all three OTU clustering approaches resulted in highly correlated diversity metrics (Suppl. File 2). So, for the general message on diversity change in German ecosystems, the OTU clustering approach does not seem to have a detectable effect.

Reviewer #2 (Recommendations for the authors):

I elaborate on my concerns regarding the sequence and statistical analyses below.

One of my greatest concerns is in the statistical analyses of the temporal trend. In most cases, the authors used linear models (LM) to test whether there is a temporal trend (i.e., increase or decrease). This is inappropriate because empirical time series often contain temporal autocorrelation structures and because the application of LM to such time series may often result in the false-positive detection of the trend. This means that LM can often detect a "significant" trend even in a random walk time series. This is a well-known issue, and the authors should apply a more appropriate statistical method, e.g., the autocorrelation model, state-space modeling, or some other methods, to judge whether there is a temporal trend.

As mentioned above, we have reanalyzed our data accordingly and report this in the Methods and Results. The reported patterns essentially remain unchanged after the new analysis.

Another concern is the use of read abundance as an explained variable in the statistical models. In L452-454, the authors claimed that the DNA copy numbers may be a proxy of biomass. I agree with this statement. However, sequence read abundance is not the DNA copy numbers and cannot be a proxy of abundance in most cases. If I understood correctly, the authors rarefied the sequence reads to 5000 reads for each sample, and the read abundance in the statistical model seems to be the relative abundance (e.g., if an OTU produced 500 reads in a sample, the relative abundance of the OTU in the sample is 10%). The authors quantified 18S DNA copy numbers of arthropods, so multiplying the relative abundance by the total 18S DNA copy numbers may produce a better proxy of the abundance of arthropods. I would recommend the authors reconsider the read abundance issue carefully.

As mentioned above, we did not use read abundance, but qPCR-based copy number as a measure for abundance. We developed a very robust qPCR assay for this purpose. As suggested by the reviewer, we have rarefied our OTU tables using qPCR-based copy number instead of equal coverage (see “Statistical analysis”). Using this new dataset did not change the recovered diversity patterns. We have also repeated this exercise for the single species abundances. Here we also find no general change of the observed patterns (Figure 5).

Reviewer #3 (Recommendations for the authors):

Line 91. 312 ESB samples. Perhaps expand here on distribution through time etc? Could be part of Figure 1?

We feel this may overload the figure a bit.

Line 125+ – stability of arthropod DNA comment. I'm actually surprised there is no higher diversity in the 8-year RT stored samples, as a result of post-mortem modifications to the DNA. Having said that does the PCR strategy negate that? If yes maybe worth stating this clearly. Actually in general given how critical the assumption that any zOTUs used represent true biological sequence variation (as opposed to PCR error, sequencing error, post mortem DNA damage) I advise the authors to early in the paper make it clear why they believe their zOTU estimates to be accurate.

Thank you for this interesting comment; we had not thought about the possibility of chemical degradation here. We are aware of studies on fossil DNA, were this is certainly a big issue. Do you think this could affect samples only a few decades old? The dried homogenates we are using so far were only ~ 10 years old, so we would be surprised if they showed chemical degradation. But we are currently exploring even older dried homogenate samples and indeed find a higher diversity in the oldest samples (> 20 years) than the younger ones. We will now carefully explore whether chemical degradation has an effect here. We have added a sentence and citation in the Discussion urging for care in analysis of dried samples (lines 241-243). Luckily, we are basing our work here on samples that have been stored on liquid nitrogen, hence there is no possibility for chemical modification.

When it comes to richness and β diversity, zOTUs and OTUs show highly replicated patterns. There does not seem to be an inflation of diversity here for zOTUs. However, we have now added more details on the picking of zOTUs for the analysis of genetic diversity in the methods (“Statistical analysis”, lines 611-615). We have included several selection steps to use them, as elaborated above.

Line 129. rDNA is typically used for many things. Ribosomal DNA. recombinant DNA. Even relative DNA. Maybe for clarity here state 'information on relative arthropod 18s rDNA copy' then there is no uncertainty.

Has been changed (line 129).

Figure 5F's label is confusing (e.g. Does < 2020 mean 2010-2020 or all years before 2020 etc?). I suspect the former but that's not how it reads as written in the figure.

We have now clarified this in the figure caption (lines 227-228).

Can the authors elaborate on when the leaves were collected in the year? In line 287 the authors state it was the same time point every year. How is that shaped by the time of year the leaves come out every year? (Which must fluctuate a lot with annual weather variation). Also in line 257, the authors state the leaves represent a 'fairly broad phenological window'. Can this be elaborated on?

We have added details in the Methods (lines 301-303). The aim of the ESB is to span the complete growing period of the leaf. That means beech leaves, for example, are collected between late August and early September. We have included the month of sampling in our statistical analysis. Over the whole time series, there is no trend of shifts in sampling events. The events are randomized on a time window between mid-August and mid-September. No effect of sampling month on diversity patterns was found.

The broad phenological window is based on the fact that many arthropods that live on leaves will only be outside of the leaf to mate for a few weeks in the year. During this time, you could catch them in a Malaise trap. The rest of the year, they spend in or on the leaves of their host plants. A good example is the spruce gall midge, which only flies for a few weeks in May, but can be detected in the leaves throughout the year. We have explained this a little more in the text (lines 234-235).

Line 287/288. It seems odd to write 'a defined amount' of both leaf and number of trees, and then give values prefixed by '>'. And then say 'defined number of branches from each tree' but then not give a value. I think I know what the intended meaning is, but perhaps consider rewording this sentence! Also, I assume the exact numbers are somewhere – maybe refer to where at this point?

We have expanded on this and provided more details in the Methods (“Tree samples of the German Environmental Specimen Bank – Standardized time series samples stored at ultra-low temperatures”, lines 301-308). The protocol is always exactly the same for each site. All variation of the protocol has been incorporated into the statistical model.

Methods in general, it took me a long time to work out exactly how the samples were collected and processed, and I'm still not sure. If I understand it correct – leaves were originally picked from trees immediately into liquid nitrogen, and then shortly after ground to a powder? Thus all happening years/decades ago? If correct please consider re-reading the text assuming you know nothing about the history and add any small clarifications that might be needed so the reader can immediately jump to this conclusion. If however, I misunderstand…then again the text needs clarification.

We have explained the sampling in more detail in the Methods (“Tree samples of the German Environmental Specimen Bank – Standardized time series samples stored at ultra-low temperatures”, lines 314-319). The ESB is unique in the whole world in that samples really go into the nitrogen in the field and never leave the nitrogen again until we analyze them.

Line 295 I assume the cryomill was somehow decontaminated between sample batches? Please clarify how.

We have added more information on decontamination protocols implemented by the ESB in the Methods (“Tree samples of the German Environmental Specimen Bank – Standardized time series samples stored at ultra-low temperatures”, line 316; and especially “Test for DNA carryover in the cryomill”, lines 448-464). The ESB sampling is only performed with sterile equipment, which is thoroughly cleaned between sampling, making cross-contamination unlikely. The thorough cleaning and decontamination are also done for the cryomill.

However, we have now included a new test to explore the possibility of contamination by the cryomill (Figure 2—figure supplement 3 and Suppl. File 3). We have collected information on the milling schedule of samples processed between 2016 and 2019. We have then identified tree samples that were processed consecutively in the cryomill. Assuming carryover is an issue, tree samples that were processed consecutively should show a significantly reduced β diversity between each other than between samples from the same site processed in different years. This effect was not found. We also tested single species carryover in more detail. This was done with a poplar and beech sample as well as a spruce and pine sample processed on two consecutive days. Assuming there is carryover in the mill, tree species-specific arthropod species should show signs of cross-contamination between the samples. This should only affect consecutively processed samples, but not samples from the same site processed at different years. The according samples were sequenced in triplicate to deep coverage (~80,000 reads per sample on average) and the read abundances of the most common tree-species-specific arthropods measured. No sign of cross-contamination was found, as seen in the new Figure 2—figure supplement 3 and Suppl. File 3.

Ca Line 365. Please state clearly if extraction replicates were made on the samples used in the full experiment, (or not). It's not clear to me as written now. I can see they were used in the sample mass trial and perhaps given the replicate dissimilarity results were ok for the PCR and extraction replicate the decision was taken to not do replicates on the large scale? If so please state it clearly.

We now state this explicitly in the Methods (lines 432-435). Considering the similar recovery of taxa from PCR and extraction replicates, we decided to run PCR replicates as technical replicates and no extraction replicates. We recently performed several replication experiments with different ESB homogenates, underlining the high efficiency of the cryomill in homogenizing the samples. For the well homogenized samples, only a small subset of homogenate is needed to saturate diversity. We find this across different ESB samples and taxa. PCR and/or extraction replication will essentially have a similar effect: more rare taxa will be picked up with replication. However, even after 2 replicates, we find the samples to show acceptable saturation, as we show that additional replication does not change our results.

https://doi.org/10.7554/eLife.78521.sa2

Article and author information

Author details

  1. Henrik Krehenwinkel

    University of Trier, Trier, Germany
    Contribution
    Conceptualization, Formal analysis, Supervision, Investigation, Methodology, Writing – original draft, Project administration, Writing – review and editing
    Contributed equally with
    Sven Weber
    For correspondence
    krehenwinkel@uni-trier.de
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-5069-8601
  2. Sven Weber

    University of Trier, Trier, Germany
    Contribution
    Formal analysis, Visualization, Methodology, Writing – original draft, Writing – review and editing
    Contributed equally with
    Henrik Krehenwinkel
    Competing interests
    No competing interests declared
  3. Rieke Broekmann

    University of Trier, Trier, Germany
    Contribution
    Investigation
    Competing interests
    No competing interests declared
  4. Anja Melcher

    University of Trier, Trier, Germany
    Contribution
    Investigation
    Competing interests
    No competing interests declared
  5. Julian Hans

    University of Trier, Trier, Germany
    Contribution
    Investigation, Methodology, Writing – review and editing
    Competing interests
    No competing interests declared
  6. Rüdiger Wolf

    University of Trier, Trier, Germany
    Contribution
    Data curation, Methodology
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-8144-5954
  7. Axel Hochkirch

    University of Trier, Trier, Germany
    Contribution
    Investigation, Writing – review and editing
    Competing interests
    No competing interests declared
  8. Susan Rachel Kennedy

    University of Trier, Trier, Germany
    Contribution
    Supervision, Investigation, Visualization, Methodology, Writing – original draft, Writing – review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-1616-3985
  9. Jan Koschorreck

    German Federal Environment Agency, Berlin, Germany
    Contribution
    Investigation, Writing – review and editing
    Competing interests
    No competing interests declared
  10. Sven Künzel

    Max Planck Institute for Evolutionary Biology, Plön, Germany
    Contribution
    Resources, Methodology
    Competing interests
    No competing interests declared
  11. Christoph Müller

    Ludwig Maximilians University, Munich, Germany
    Contribution
    Investigation, Methodology
    Competing interests
    No competing interests declared
  12. Rebecca Retzlaff

    University of Trier, Trier, Germany
    Contribution
    Investigation, Methodology
    Competing interests
    No competing interests declared
  13. Diana Teubner

    University of Trier, Trier, Germany
    Contribution
    Resources, Methodology
    Competing interests
    No competing interests declared
  14. Sonja Schanzer

    Ludwig Maximilians University, Munich, Germany
    Contribution
    Investigation, Methodology
    Competing interests
    No competing interests declared
  15. Roland Klein

    University of Trier, Trier, Germany
    Contribution
    Conceptualization, Resources, Investigation
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-8735-0393
  16. Martin Paulus

    University of Trier, Trier, Germany
    Contribution
    Conceptualization, Resources, Funding acquisition, Investigation
    Competing interests
    No competing interests declared
  17. Thomas Udelhoven

    University of Trier, Trier, Germany
    Contribution
    Formal analysis, Investigation, Methodology
    Competing interests
    No competing interests declared
  18. Michael Veith

    University of Trier, Trier, Germany
    Contribution
    Conceptualization, Resources, Funding acquisition, Writing – review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-7530-4856

Funding

Deutsche Bundesstiftung Umwelt

  • Henrik Krehenwinkel

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

We thank Karin Fischer for assistance with lab work. The German Environment Agency made the leaf samples available. Thanks to Diethard Tautz, Natalie Graham, and Rosemary Gillespie for critical revision of an earlier version of the manuscript. Frank Thomas, Dorothee Krieger, and Bernhard Backes provided access to the dried leaf homogenate we used here. Andrea Koerner helped in acquiring the leaf samples from the ESB. SW was funded by a PhD fellowship of the German Federal Environmental Foundation (DBU).

Senior Editor

  1. Detlef Weigel, Max Planck Institute for Biology, Germany

Reviewing Editor

  1. Simon Creer, Bangor University, Bangor, United Kingdom

Reviewers

  1. Rafael Valentin, Princeton University, United States
  2. Thomas Gilbert

Version history

  1. Received: March 10, 2022
  2. Preprint posted: April 29, 2022 (view preprint)
  3. Accepted: November 6, 2022
  4. Accepted Manuscript published: November 10, 2022 (version 1)
  5. Version of Record published: December 20, 2022 (version 2)

Copyright

© 2022, Krehenwinkel, Weber et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 1,795
    Page views
  • 283
    Downloads
  • 6
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

  1. Henrik Krehenwinkel
  2. Sven Weber
  3. Rieke Broekmann
  4. Anja Melcher
  5. Julian Hans
  6. Rüdiger Wolf
  7. Axel Hochkirch
  8. Susan Rachel Kennedy
  9. Jan Koschorreck
  10. Sven Künzel
  11. Christoph Müller
  12. Rebecca Retzlaff
  13. Diana Teubner
  14. Sonja Schanzer
  15. Roland Klein
  16. Martin Paulus
  17. Thomas Udelhoven
  18. Michael Veith
(2022)
Environmental DNA from archived leaves reveals widespread temporal turnover and biotic homogenization in forest arthropod communities
eLife 11:e78521.
https://doi.org/10.7554/eLife.78521

Share this article

https://doi.org/10.7554/eLife.78521

Further reading

    1. Ecology
    2. Plant Biology
    Jamie Mitchel Waterman, Tristan Michael Cofer ... Matthias Erb
    Research Article

    Volatiles emitted by herbivore-attacked plants (senders) can enhance defenses in neighboring plants (receivers), however, the temporal dynamics of this phenomenon remain poorly studied. Using a custom-built, high-throughput proton transfer reaction time-of-flight mass spectrometry (PTR-ToF-MS) system, we explored temporal patterns of volatile transfer and responses between herbivore-attacked and undamaged maize plants. We found that continuous exposure to natural blends of herbivore-induced volatiles results in clocked temporal response patterns in neighboring plants, characterized by an induced terpene burst at the onset of the second day of exposure. This delayed burst is not explained by terpene accumulation during the night, but coincides with delayed jasmonate accumulation in receiver plants. The delayed burst occurs independent of day:night light transitions and cannot be fully explained by sender volatile dynamics. Instead, it is the result of a stress memory from volatile exposure during the first day and secondary exposure to bioactive volatiles on the second day. Our study reveals that prolonged exposure to natural blends of stress-induced volatiles results in a response that integrates priming and direct induction into a distinct and predictable temporal response pattern. This provides an answer to the long-standing question of whether stress volatiles predominantly induce or prime plant defenses in neighboring plants, by revealing that they can do both in sequence.

    1. Ecology
    Congnan Sun, Yoel Hassin ... Yossi Yovel
    Research Article

    Covid-19 lockdowns provided ecologists with a rare opportunity to examine how animals behave when humans are absent. Indeed many studies reported various effects of lockdowns on animal activity, especially in urban areas and other human-dominated habitats. We explored how Covid-19 lockdowns in Israel have influenced bird activity in an urban environment by using continuous acoustic recordings to monitor three common bird species that differ in their level of adaptation to the urban ecosystem: (1) the hooded crow, an urban exploiter, which depends heavily on anthropogenic resources; (2) the rose-ringed parakeet, an invasive alien species that has adapted to exploit human resources; and (3) the graceful prinia, an urban adapter, which is relatively shy of humans and can be found in urban habitats with shrubs and prairies. Acoustic recordings provided continuous monitoring of bird activity without an effect of the observer on the animal. We performed dense sampling of a 1.3 square km area in northern Tel-Aviv by placing 17 recorders for more than a month in different micro-habitats within this region including roads, residential areas and urban parks. We monitored both lockdown and no-lockdown periods. We portray a complex dynamic system where the activity of specific bird species depended on many environmental parameters and decreases or increases in a habitat-dependent manner during lockdown. Specifically, urban exploiter species decreased their activity in most urban habitats during lockdown, while human adapter species increased their activity during lockdown especially in parks where humans were absent. Our results also demonstrate the value of different habitats within urban environments for animal activity, specifically highlighting the importance of urban parks. These species- and habitat-specific changes in activity might explain the contradicting results reported by others who have not performed a habitat specific analysis.