Changes in land use and agricultural intensification threaten biodiversity and ecosystem functioning of small water bodies. We studied 67 kettle holes (KH) in an agricultural landscape in northeastern Germany using landscape-scale metatranscriptomics, to understand the responses of active bacterial, archaeal, and eukaryotic communities, to land-use type. These KH are proxies of the millions of small standing water bodies of glacial origin spread across the northern hemisphere. Like other landscapes in Europe, the study area has been used for intensive agriculture since the 1950s. In contrast to a parallel eDNA study which revealed the homogenization of biodiversity across KH conceivably resulting from long-lasting intensive agriculture, land-use type affected the structure of the active KH communities during spring crop fertilization, but not a month later. This effect was more pronounced in eukaryotes than in bacteria. In contrast, gene expression patterns did not differ between months or across land-use type, suggesting a high degree of functional redundancy across the KH communities. Variability in gene expression was best explained by active bacterial and eukaryotic community structures, suggesting that these changes in functioning are primarily driven by interactions between organisms. Our results show that influences of the surrounding landscape result in temporary changes in the activity of different community members. Thus, even in KH where biodiversity has been homogenized, communities continue to respond to land management. This needs to be considered when developing sustainable management options for restoration purposes and for successful mitigation of further biodiversity loss in agricultural landscapes.

The sampling focused on 67 kettle holes (KH) in northeastern Germany (Uckermark district, State of Brandenburg; Fig. S1), 52 of which were sampled in May and 43 five weeks later in June. No samples were taken in dried-up KH, resulting in a total of 41 KH sampled on both occasions. Of the samples KH 36, 7, 9, and 28, 6, 9 were in arable fields, grasslands and forest in May and June, respectively. The area is among the least populated regions in Germany. The study area has long been used for extensive agriculture, with >90 % of the land used as arable fields (Kalettka & Rudat, 2006). This includes areas where land use was changed from arable fields to grasslands nearly two decades ago (Serrano et al., 2017). Since the 1950s, agriculture in the area was industrialized, which included increased fertilizer and pesticide use.

KH were categorized according to the predominant land-use type within a perimeter of ca. 50 m. Accordingly, all KH in crop fields (rapeseed, corn, wheat, barley, rye, triticale), are referred to as “arable field KH,” both those directly adjacent to the fields and those surrounded by natural vegetation. KH in grasslands are referred to as “grassland KH”. “Forest KH”, located in the Kiecker nature reserve (Nordwestuckermark, Brandenburg), comprised KH in vast mixed forests (beech and oak) as well as in forest patches (> 100 m in diameter) surrounded by arable fields (Fig. S1). However, the last category was treated as “arable fields” in analyses where we applied a stricter definition of forests.

Sampling

Water Samples for RNA analysis were collected during two sampling campaigns (each 2-3 days) in late spring and early summer 2017, together with samples collected for eDNA analysis (Ionescu et al. submitted). Water samples were taken whenever water was available. To obtain a representative sample from each water body, total volumes of ca. 20 L were collected from 5-15 different locations in each KH, with the number of individual samples varying with KH size. The water was combined in prewashed buckets and mixed, before 1.7 L were resampled for RNA analysis into plastic canisters containing 800 mL RNA-stabilizing solution (15 mM EDTA, 18.5 mM sodium citrate, 4 M ammonium sulfate). Samples were placed in iceboxes containing a mixture of ice and table salt to lower the freezing point. Upon arrival in the laboratory, the samples were frozen at -80 °C until further analyses.

RNA extraction and processing

Before RNA extraction, standard volumes of water (2.3 L: sample + fixative) were sequentially filtered on a Nalgene filtration tower (ThermoFisher Scientific, Dreieich, Germany). Polycarbonate filters with pore sizes of 10 and 5 µm (Millipore TCTP04700, TMTP04700, Merck, Darmstadt, Germany) were used, as well as combusted GF/F and polycarbonate filters with 0.2 µm pore size (Whatman WHA1825047, Millipore GTTP04700, Merck, Darmstadt, Germany). All filter diameters were 47 mm. The entire water volume was passed through all filters. The filters were rinsed twice with 50 mL autoclaved MQ water to remove salts and subsequently flash frozen.

To avoid introducing batch effects (Bálint, Márton, Schatz, Düring, & Grossart, 2018), Eppendorf tubes containing the filters representing sample-fractions were shuffled and randomly allocated to separate batches. RNA was extracted following a phenol/chloroform procedure modified from Nercessian et al. (2005). In brief, a CTAB extraction buffer containing SDS and N-lauryl sarcosine was added to the samples together with an equal volume of phenol/chloroform/isoamylalcohol (25:24:1) solution. The samples underwent a bead-beating treatment, followed by centrifugation, cleaning with chloroform, and precipitation with PEG-6000 (Sigma-Aldrich, Taufkirchen, Germany). The precipitated DNA/RNA mix was rinsed with 1 mL 70 % ethanol, dried, and dissolved in water. Finally, all extractions belonging to a given sample were pooled.

DNA was removed by two sequential treatments with the TurboDNAfree Kit (Invitrogen ThermoFisher Scientific, Dreieich, Germany), after which the samples were transferred to an RNAstable 96-well plate (Sigma-Aldrich, Taufkirchen, Germany) for shipment. A total of 98 samples were sequenced at MrDNA (Molecular Research, Shallowater, Texas, USA) according to the following procedure: The RNA samples were resuspended in 30 µL of nuclease-free water and cleaned using the RNeasy PowerClean Pro Cleanup Kit (Qiagen, Germantown, MD, USA). The concentration of total RNA was determined using the Qubit® RNA Assay Kit (Life Technologies, Thermofisher, Grand Island, NY, USA). Next, 750 ng of total RNA were used to remove the remaining DNA contamination using Baseline-ZERO™ DNase (Epicentre, Lucigen, Middleton, WI, USA) according to the manufacturer's instructions, followed by a purification step with RNA Clean & Concentrator-5 columns (Zymo Research, Irvine, CA, USA). DNA-free RNA samples were used for library preparation using the TruSeq™ RNA LT Sample Preparation Kit (Illumina, Hayward, CA, USA) according to the manufacturer’s instructions. Following library preparation, the final concentration of all the libraries were measured using the Qubit® dsDNA HS Assay Kit (Life Technologies, Thermofisher), and the average library size was determined using the Agilent 2100 Bioanalyzer (Agilent Technologies, Cedar Creek, TX, USA). The libraries were then pooled in equimolar ratios of 2 nM, and 6 pM of the library pools was clustered using the cBot (Illumina, Hayward, CA, USA) and sequenced 2x125 paired end reads on 20 lanes for 250 cycles using the HiSeq 2500 system (Illumina, Hayward, CA, USA). The sequenced data was submitted to the NCBI short read archive under project number PRJNA640812 (https://www.ncbi.nlm.nih.gov/sra/PRJNA640812).

Raw files of paired end reads were quality-trimmed using Trimommatic (V 0.39) (Bolger, Lohse, & Usadel, 2014). Ribosomal RNA reads were removed by stringent mapping to a database of SSU, LSU and 5S rRNA assembled manually from the SSU and LSU Silva databases (V132) (Quast et al., 2013). Subsequently the SSU rRNA was annotated using PhyloFlash (Gruber-Vodicka, Seah, & Pruesse, 2020) and Kraken2 (Wood, Lu, & Langmead, 2019). The non-rRNA sequences were further checked using BARNAP (V 0.9). The clean non-rRNA reads of each sample were individually processed according to the Trinotate (https://github.com/Trinotate/Trinotate.github.io/wiki) pipeline, including assembly with Trinity V 2.6.5 (Grabherr et al., 2011), protein prediction using TransDecoder (https://github.com/TransDecoder/TransDecoder), and annotation with Diamond BlastP and BlastX (Buchfink, Xie, & Huson, 2015) against the Uniprot database. Sequences were also annotated with hmmsearch (Gough, Karplus, Hughey, & Chothia, 2001) and the pFam (Finn et al., 2014) database. Kallisto (V 0.44) (Bray, Pimentel, Melsted, & Pachter, 2016) was used to map the reads from each sample against the samples’ assembled transcripts resulting in TPM-normalized counts. The data was merged to generate abundance matrices for statistical analysis. BlastP, BlastX, EC-number and Subsystems’ matrices were obtained and separately analyzed. The presented results stem from the Subsystem annotation of the data. More information on SEED subsystems is available at: https://www.theseed.org/wiki/SEED_Viewer_Manual.