Local biodiversity patterns are expected to strongly reflect variation in topography, land use, dispersal boundaries, nutrient supplies, contaminant spread, management practices and other anthropogenic influences. In contrast, studies focusing on specific taxa revealed a biodiversity homogenization effect in areas subjected to long-term intensive industrial agriculture. We investigated whether land use affects biodiversity levels and community composition (α & β diversity) in 67 kettle holes (KH) representing small aquatic islands embedded in the patchwork matrix of a largely agricultural landscape comprising grassland, forest, and arable fields. These KH, similar to millions of standing water bodies of glacial origin, spread across northern Europe, Asia, and North America, are physico-chemically diverse, differ in the degree of coupling with their surroundings. We assessed biodiversity patterns of eukaryotes, Bacteria and Archaea in relation to environmental features of the KH, using deep-amplicon-sequencing of environmental DNA (eDNA). First, we asked whether deep sequencing of eDNA provides a representative picture of KH biodiversity across the Bacteria, Archaea, and Eukaryotes. Second, we investigated if and to what extent KH biodiversity is influenced by the surrounding land-use. Our data shows that deep eDNA amplicon sequencing is useful for in-depth assessments of cross-domain biodiversity comprising both micro- and macro-organisms, but, has limitations with respect to single-taxa conservation studies. Using this broad method, we show that sediment eDNA, integrating several years to decades, depicts the history of agricultural land-use intensification. The latter, coupled with landscape wide nutrient enrichment (including by atmospheric deposition), groundwater connectivity between KH and organismal (active and passive) dispersal in the tight network of ponds, resulted in a biodiversity homogenization in the KH water, levelling off today’s detectable differences in KH biodiversity between land-use types.

Study sites and sampling

Samples for eDNA analysis were collected during 5 sampling campaigns of 2-3 days each in December 2016, and March, May, June and October 2017. All samples were taken in a set of 67 kettle holes located in northeastern Germany (Fig. 1). The area, one of the least populated in Germany, has a long history of farming, with >90% of the land now being covered by arable fields (Kalettka and Rudat, 2006), although some of that land has been reconverted to grassland nearly two decades ago (Serrano et al., 2017). Routine monitoring of the KH water and riparian vegetation in the area started in 1993, shortly after the reunification of Germany (Kalettka and Rudat, 2006). Each KH was categorized based on land-use type within a perimeter of ca. 50 m around the kettle holes, i.e. distinguishing KH in arable fields, grasslands, and forest patches (Fig. 1).

Water samples were collected whenever enough water was present in the KH. Some occasionally fell dry, however (Nitzsche et al., 2017), and thus could not be sampled at all times, particularly in October 2017. To obtain representative samples, total volumes of ca. 20 L were collected at 5-15 locations selected within each KH, with the number of individual samples depending on KH size. The samples were pooled in cleaned buckets and 2 L were subsampled in the field, placed in ice chests containing a mixture of ice and table salt to lower the temperature during transport, and subsequently frozen at -80 °C in the laboratory for later eDNA analysis.

Sediment cores were taken at three time points (Table S1). In March 2017, sediment samples were collected from 54 of the 67 KH, both wet and dry. In some instances, a dense mat of belowground plant parts prevented sediment coring. Subsequently, sediment cores were only collected from wet KH that had recently dried out, or from previously dry KH that had refilled. Between 3-7 cores were taken per KH, depending on KH size, covering both littoral and central areas. The cores were sectioned into surface (upper 5 cm) and lower (5-20 cm depth) sediment layers to try to separate current benthic communities from older resting stages and preserved eDNA. The sections were separately transferred into plastic bags and subsampled (1 g wet weight) for eDNA extraction. Both the complete samples and subsamples were stored at -80 °C for further processing. DNA extractions from multiple cores representing surface or lower sediment layers of a given KH at each sampling date were pooled. A compilation of the collected samples is given in Table S1.

Analysis of water physico-chemical properties

Temperature, conductivity, pH, redox potential, and oxygen concentration and saturation were measured on site during sampling using a multiparameter field probe (HI98194, Hanna Instruments, Vöhringen, Germany). Additional water (1 L) was collected to determine concentrations of nutrients and major ions. These samples were immediately frozen in an ice chest containing crushed ice mixed with table salt (NaCl) and analyzed within 48 h. Water analysis followed German standard methods (DIN 38405, 2018). Ca²⁺, Mg²⁺, K⁺, Na⁺, and total Fe were analyzed by inductively coupled plasma optical emission spectrometry (ICP-iCAP 6300 DUO, ThermoFisher Scientific GmbH, Dreieich, Germany). Br^-, Cl^-, NO₃^-, NO₂^- and SO₄^2- were analyzed using ion chromatography (882 Compact IC plus, Deutsche Metrohm GmbH & Co. KG, Filderstadt, Germany). Ammonium (NH₄⁺) and soluble reactive phosphorus (ortho-phosphate; o-PO₄^3--P) were measured spectrophotometrically (SPECORD 210 plus, Analytik Jena AG, Jena, Germany). Total phosphorus (TP) was measured as soluble reactive phosphorus after microwave digestion (Gallery™ Plus, Microgenics GmbH, Hennigsdorf, Germany). Dissolved organic carbon (DOC), total organic carbon (TOC) and total nitrogen (TN) were determined using an elemental analyzer (TOC-VCPH, Shimadzu Deutschland GmbH, Duisburg, Germany) with chemiluminescence detection. The specific absorption coefficient (SAC) was measured on a spectrophotometer (SPECORD 210 plus, Analytik Jena AG, Germany) as a proxy of dissolved aromatic carbon content (Weishaar et al., 2003). Finally, the SAC:DOC ratio was used as a rough measure of DOC quality.

DNA extraction

The collected 2-L water samples were sequentially filtered (Nalgene filtration tower; ThermoFisher Scientific, Dreieich, Germany) to prevent clogging the filters used were polycarbonate membrane filters (pore size of 10 and 5 µm), combusted GF/F filters and finally polycarbonate filters with a pore size of 0.2 µm (47 mm diameter of all filters). The GF/F filter was included owing to its charge to capture naked eDNA and DNA released from cells lysed by freezing and thawing. All filters were rinsed twice with 50 mL autoclaved MilliQ water to remove salts, and subsequently flash frozen and stored at -80 °C.

Total (environmental) DNA was extracted from 329 samples consisting of 182 water samples, 75 surface sediment samples (< 5 cm), and 66 deeper sediment (5-20 cm) samples. To prevent analytical biases (Bálint et al., 2018), the different filtered subsamples were extracted in separate, randomly selected batches. DNA was extracted with phenol/chloroform according to a method modified by Nercessian et al. (2005). In brief, a CTAB extraction buffer containing SDS and N-laurylsarcosine was added to the samples together with an equal volume of phenol/chloroform/isoamylalcohol (25:24:1) solution. The samples were subject to bead-beating (FastPrep-24™ 5G Instrument, MP Biomedical, Eschwege, Germany), followed by centrifugation (14,000 g), a cleaning step with chloroform, and DNA precipitation with PEG-6000 (Sigma-Aldrich, Taufkirchen, Germany). The precipitated DNA was rinsed with 1 mL of 70% ethanol, dried and dissolved in water. Finally, all extracts from the same sample were pooled and kept at -80 °C till further processing.

Sequencing

Sequencing was conducted separately for the SSU rRNA gene of Archaea, Bacteria and eukaryotes at MrDNA (Shallowater, TX, USA) using the following primers: Arch2A519F (5’ CAG CMG CCG CGG TAA 3’) and Arch1071R (5’ – GGC CAT GCA CCW CCT CTC - 3’) for archaea (Fischer et al., 2016); 341F (5’ CCT ACG GGN GGC WGC AG 3’) and 785R (5’ GAC TAC HVG GGT ATC TAA TCC 3’) for bacteria (Thijs et al., 2017); and Euk1560F (5’ TGG TGC ATG GCC GTT CTT AGT 3’) and Euk2035R (5’ CAT CTA AGG GCA TCA CAG ACC 3’) for eukaryotes (Hardy et al., 2010). The primers were barcoded on the forward primer and used in a 30-cycle PCR using the HotStarTaq Plus Master Mix Kit (Qiagen, Hilden, Germany) under the following conditions: 94 °C for 3 min, followed by 30 cycles at 94 °C for 30 s, 53 °C for 40 s and 72 °C for 1 min, followed by a final elongation step at 72 °C for 5 min. The PCR products were checked in 2% agarose gel to determine success of the amplification and relative band intensity. To ensure high coverage of rare taxa, batches of 20 samples were pooled for each sequencing run in equal proportions based on their molecular weight and DNA concentrations. The PCR products were purified using calibrated Ampure XP beads and then used to prepare an Illumina DNA library. Paired end 2 x 300 bp sequencing was performed on a MiSeq sequencer (Illumina, Inc., San Diego, CA. USA) following the manufacturer’s instructions. Sequence data are available at the NCBI Short Read Archive under project number PRJNA641761.

Bioinformatic analysis

Paired end reads were merged using BBMerge from the BBMap package (part of JGI tools; https://sourceforge.net/projects/bbmap), after which the joined reads were quality trimmed and demultiplexed using cutadapt (V 1.16) to remove reads of low quality (q > 20) and shorter than 150 nt. Taxonomic annotation was performed for all reads from all samples without clustering based on the SILVA SSU NR99 data base (V132; Quast et al., 2013) . This was accomplished by using PhyloFlash (V 3.3 b1; https://github.com/HRGV/phyloFlash; Gruber-Vodicka et al., 2019) and Kraken 2 (Wood et al., 2019). To improve the annotation of eukaryotic taxa, a new database was created consisting of all eukaryotic sequences in the SILVA SSU Parc database (V138). The SILVA Parc database also includes eukaryotic sequences shorter than 900 nucleotides and hence covers a much broader range of species than the SILVA NR99 database. The eukaryotic sequences from all samples were annotated using both PhyloFlash and the SINA aligner (Pruesse et al., 2012; V 1.6; https://github.com/epruesse/SINA) requiring a minimum consensus of 3 sequences for last common ancestor assignments. The resulting annotations were merged according to taxonomic names and presence/absence matrices were generated to account for the qualitative nature of the eDNA method especially when merging data from separate assays (i.e. separately targeting Archaea, Bacteria, and eukaryotes). Statistical analyses (see next section) using matrices generated by different annotation tools resulted in identical patterns.

The functional potential of the bacterial community was derived from the taxonomic annotation using the FaProTax tool (Louca et al., 2016).