Sampling strategy
All sampling sites were located in the Eifel National park, situated in the south-western part of Germany close to the Belgian border (Supplementary Fig. 1, Supplementary Table 1).
In this study the sampling site comprised a forest conversion gradient from a Norway Spruce (Picea abies) monoculture to a European Beech forest (Fagus sylvatica). To reflect the different stages of conversion from spruce to beech, four forest types were defined: pure beech (PB), old beech (OB), young beech (YB) and pure spruce (PS) (Supplementary Table 1, Supplementary Fig. 2).
The forest types differed in tree species composition and tree age. The pure beech and pure spruce forest types were monoculture stands. The pure beech stands were approximately 180 years old and partly under special protection through North-Rhine Westphalia (Naturwaldzelle) (Sampling site 01). The spruce monocultures were substantially younger with ca. 60 years old. Spruces of the same age dominated the young beech sampling sites that had only recently been underplanted with beeches. At the old beech sampling sites, beeches had already reached a height of more than 3 m and actions to remove spruces from the forest were conducted.
A total of 12 Townsend Malaise traps (three per forest type) were set up in the Eifel National Park, North-Rhine Westphalia, Germany, during July 2016. To ensure that the orientation of the Malaise traps was consistent and to minimize potential biases caused by wind direction and position of sun, the highest point of each trap was set facing south. The traps were left in the field for the full duration of the experiment until April 2017, ensuring that insects were collected from exactly the same locations. In October 2016, two additional traps (Malaise Trap 13, pure spruce and Malaise Trap 14, old beech) were installed at two further sampling sites (Sample Site 13 and Sample Site 14). All traps were equipped with a bottle filled with approximately one litre of absolute ethanol (99,96%) over a 2-week period in July 2016 (13.07–27.07), October 2016 (13.10–27.10), January 2017 (11.01–25.01) and April 2017 (12.04–26.04) (Supplementary Table 2). The ethanol was replaced every week to ensure that the concentration of the preservative ethanol was stable and to avoid loss of insects a mesh filter was used (MICROFIL V Filter White Gridded 0.45 µm-diameter 47 mm & 100 ml Funnel Sterilized). Due to heavy snow during the winter period, new traps were set at the start of the new sampling season in January 2017.
Three soil samples were collected around each Malaise trap, from the organic horizon of the top 10 cm layer (excluding the litter layer). Soil sample sites were located 4 m and 5 m away from the trap, forming a triangle in the centre of which the Malaise trap was located (Supplementary Fig. 3). One corner of the sampling triangle was pointing south, while both remaining corners were pointing north west and north east, respectively.
Each sampling site was sampled four times in the course of a 1-year period. Soil sampling and Malaise trapping were synchronized and soil sampling was done on day 14 of each Malaise trapping period, when the last bottles were collected (Supplementary Table 3). Each soil sample consisted of approximately twenty 44 mm × 100 mm cores, taken 5 cm apart. A total of 162 soil samples were collected and stored at − 20 °C until further processing.
DNA extraction
Bulk samples from Malaise traps
Non-destructive DNA extraction was performed by overnight incubation in lysis buffer, using a modified protocol of Aljanabi and Martinez (1997). The arthropods were sieved from the collecting ethanol using a mesh filter (MICROFIL V Filter White Gridded 0.45 µm-diameter 47 mm and 100 ml Funnel Sterilized), which was processed with the specimens. The insects were dried for 10 min at room temperature. Depending on biomass, between 15 and 25 ml of extraction buffer (0.4 M NaCl, 10 mM Tris-HCl pH 8.0, 2 mM EDTA pH 8.0) and 2% Sodium dodecyl sulphate (SDS) were added to each bulk sample. Finally, 400 µg Proteinase K was added per ml of lysis buffer and samples were lysed overnight at 52 °C on an orbital shaker at 30 rpm. The next day, the lysate was poured out of the bottles using the MICROFIL V Filter (White Gridded 0.45 µm-Dia 47 mm and 100 ml Funnel Sterilized) and a 6 M NaCl solution was added to the lysate to a concentration of 4 mmol. The samples were vortexed for 30 s, centrifuged at 4700 rpm for 30 s and the supernatant was transferred to a falcon tube and an equal volume of isopropanol was added. After careful mixing by inversion, the tubes were left at − 20 °C for 1 h and subsequently centrifuged at 4700 rpm for 60 min. The supernatant was discarded and the resulting pellet was washed with 20 ml of ice cold 70% ethanol, by centrifuging at 4700 rpm for 15 min. The remaining ethanol was discarded and the pellet was left to dry at 20 °C overnight. The pellet was then resuspended in 1 ml of sterile H2O and stored at − 20 °C until further processing.
Soil samples
DNA extraction from the soil samples was conducted using two different extraction methods: a commercial (lysis-based) DNA extraction kit (Macherey-Nagel NucleoSpin Soil) and a (no lysis) phosphate buffer protocol from Taberlet et al. 201240. Each of the triplicated samples were processed individually. After defrosting the soil overnight at 4 °C, the samples were thoroughly mixed, DNA was extracted from 0.5 g of soil per sample using the Macherey-Nagel NucleoSpin Soil Kit following the manufacturer’s protocol.
The second DNA extraction method allowed extracellular DNA to be extracted from larger amounts of starting material using a phosphate buffer and did not include a lysis step. Each of the three samples taken per sample site and season were treated individually. Soil samples were removed from the − 20 °C chamber approximately 12 h before DNA extraction and stored at 4 °C overnight. The next morning, each sample was thoroughly mixed and an equal weight of saturated phosphate buffer solution (Na2HPO4; 0.12 M; pH 8)40 was added. Samples were placed in an orbital shaker at 120 rpm for 15 min. Thereafter, duplicates were processed, where two 2 ml Eppendorf tubes were filled with 1.7 ml of the resulting mixture and centrifuged for 10 min at 10,000 g. Four hundred microliters of the resulting supernatant were transferred to a new 2 ml collection tube and 200 μl of SB binding buffer from the Macherey-Nagel NucleoSpin Soil Kit was added. Duplicate lysates were merged by loading onto a single NucleoSpin Soil Column and centrifuged at 10,000 g for 1 min. From this step onwards, the standard manufacturer’s protocol for the Macherey-Nagel NucleoSpin Kit was followed from step 8. DNA was eluted with 50 μl of SE Buffer (Macherey-Nagel). Ten microliters of the resulting DNA eluate was diluted in 90 μl of pure H2O (Sigma), followed by DNA purification using the PowerClean Pro DNA Clean-Up Kit (MO Bio Laboratories, Inc.) following the manufacturer’s protocol. For the purposes of this study, results from the two types of soil extraction were merged.
Choice of primers and amplicon library preparation
For amplicon library preparation a primer pair targeting the 313 bp ‘mini barcode’ region of the mitochondrial Cytochrome c Oxidase subunit I gene (COI) was used41. The ‘mini barcode’ primer pair consisted of the forward primer mlCOIintF 5′–ACACTCTTTCCCTACACGACGCTCTTCCGATCTGGWACWGGWTGAACWGTWTAYCCYCC -3′41 and the reverse primer dgHCO2198 5′– GTGACTGGAGTTCAGACGTGTGCTCTTCCGATCTTAAACTTCAGGGTGACCAAARAAYCA-3′41 (Illumina overhang in regular font and primer in bold). Library preparation was carried out using a two-step PCR approach42,43 , whereby PCR1 amplifies the gene region of interest and PCR2 adds the sample index together with the Illumina overhang (indexed primers). For each sample, a unique combination of indexes was chosen.
DNA extracts were quantified using the Quantus Fluorometer (Promega). Ten nanograms of template DNA was used for PCR1. PCR1 consisted of 7.5 µl Q5 Hot Start High-Fidelity 2 × Master Mix (New England BioLabs), 1 μl Sigma H2O, 0.5 µl forward Primer (10 µM), 0.5 µl reverse primer (10 µM), 0.5 μl Bovine Serum Albumin (Thermo Scientific) and 1 µl template DNA, making up a total of 15 µl. PCR1 cycling conditions were as follows: 2 min at 98 °C (1 ×); 40 s at 98 °C, 40 s at 45 °C, 30 s at 72 °C (20 ×); 3 min at 72 °C (1 ×). PCR products were purified with 4 µl of HT ExoSAP-IT (Applied Biosystems) to each sample, following the manufacturers’ protocol. For PCR2 (index PCR) the purified PCR products were split into two PCR tubes.
For PCR2, each tube contained 12.5 µl Q5 Hot Start High-Fidelity 2X Master Mix (New England BioLabs), 3 µl Sigma H2O, 1.2 µl of index forward primer (10 µM) (AATGATACGGCGACCACCGAGATCTACAC NNNNNNNN ACACTCTTTCCCTACACGACGC TC), 1.2 µl of index reverse primer (10 µM) CAAGCAGAAGACGGCATACGAGAT NNNNNNNN GTGACTGGAGTTCAGACGTGTGCTC) and 8 µl of purified PCR1 product. PCR2 cycling conditions were as follows: 2 min at 98 °C (1 ×); 40 s at 98 °C, 30 s at 55 °C, 30 s at 72 °C (20 ×); 3 min at 72 °C (1 ×). All tagged PCR products were visualised by gel electrophoresis and PCR bands with the expected size were excised and purified using the QIAquick Gel Extraction Kit (Qiagen). Purified PCR products were quantified using the Quantus Fluorometer (Promega) and pooled in equal concentrations. The resulting purified amplicon library pool (3 ng/µl) was sequenced on an Illumina MiSeq (MiSeq Reagent Kit v3, 2 × 300 bp) sequencing platform at Liverpool University’s Centre for Genomic Research (Liverpool, UK). Raw sequence data were deposited in the GenBank short read archive (SRA) under accession numberPRJNA681091 and PRJNA706915.
Bioinformatics and data analysis
The raw fastq files were trimmed for the presence of Illumina adapter sequences using Cutadapt version 1.2.1. at the Centre for Genomic Research (Liverpool, UK). Sequences were trimmed using Sickle version 1.200 with a minimum window quality score of 20 and reads shorter than 20 bp were removed after trimming.
The fastq sequences were then checked for the presence of the COI primers with Cutadapt version 1.1844 using the following settings: maximum error rate (-e): 0.1, minimum overlap (-O): 20, minimum sequence length (-m): 50. Sequences lacking either the forward or reverse primer were removed and primer pairs were trimmed off from the remaining sequences. Subsequently, paired-end reads were merged with vsearch version 2.7.045. Merged sequences with a length of 293–333 bp were retained for further analysis and filtered with a maxEE threshold of 1.0 using vsearch (version 2.7.0)45 before demultiplexing the fastq sequences using the script split_libraries_fastq.py implemented in QIIME146 using a phred quality threshold of 19. Dereplication, size sorting, denovo chimera detection as well as Operational Taxonomic Unit (OTU) clustering with a 97% cutoff was conducted with vsearch 2.7.045. Finally, an OTU table was built by using the –usearch_global function in vsearch 2.7.045 followed by the python script “uc2otutab.py” (https://drive5.com/python/uc2otutab_py.html). For taxonomy assignment, representative sequences were blasted against the GBOL database (https://www.bolgermany.de/gbol1/identifications downloaded on 2nd of July 2019) using blastn 2.9.0+47.
The resulting OTU table was curated with LULU48. Curation started with an initial blasting of OTU representative sequences against each other using blastn (version 2.9.0). The following parameter settings were chosen: ‘query coverage high-scoring sequence pair percent’ (-qcov_hsp_perc) was set to 80, meaning that a sequence was reported as a match when 80% of the query formed an alignment with an entry of the reference file. Secondly, minimum percent identity (-perc_identity) was set to 84, requiring the reference and query sequence to match by at least 84% to be reported as a match. The format of the output file was customized using the –outfmt settings ‘6 qseqid sseqid pident’. The output file included the name of the query sequence and the name of the reference sequence next to the percentage match. The resulting OTU match list was uploaded into R (version 3.5)49 and the R-package ‘lulu’ (version 0.1.0)48 was used to perform post clustering curation using standard settings. The LULU algorithm filters the dataset for artificial OTUs and these were either classified as “daughter OTU” and merged with the corresponding “parent OTU” or were discarded from the dataset.
The resulting curated OTU table was loaded into Excel where data was formatted to upload into R (R studio running R version 3.5). Only OTUs with an assignment at species level (blastID ≥ 99%) were used for subsequent analysis. Furthermore, results from the two types of soil extraction were merged.
UpsetR plots were prepared using the R package UpSetR (version 1.4.0)50 for visualization of shared arthropod OTUs between sample types in each season. Differences in number of OTU proportions are shown in a Marioko plot prepared with the R package ggplot251. To analyze dissimilarities between communities depending on season and sample type, Permutational Multivariate Analysis of Variance (PERMANOVA) using Jaccard distance matrices for incidence data of detected arthropod species (blastID ≥ 99%) were performed using dplyr (version 0.8.3)52, betapart (version 1.5.1)53 and vegan (version 2.5–6)54. In order to analyse dissimilarity differences in arthropod community composition between the different forests and seasons the Jaccard similarity index (J) was used on a presence-absence matrix based on arthropod species. Calculated Jaccard indices were visualized on a heatmap using the R package ggplot251. Sample completeness curves and sample-size-based R/E curve with extrapolations of Hill numbers for incidence data based on the combined dataset for all forests and seasons were prepared using the R-package iNEXT55 at default settings (40 knots, 95% confidence intervals generated by the bootstrap procedure (50 bootstraps)).
To correlate community structure and diversity levels with the different seasons and forest types a Permutational Multivariate Analysis of Variance (PERMANOVA) based on the Jaccard similarity index for a presence-absence matrix of detected arthropod species (blastID ≥ 99%) was performed with the adonis function in R. Differences in arthropod community composition between the different forests and seasons was assessed using the Jaccard similarity index (J), where the higher the index, the more similar the communities.
Source: Ecology - nature.com
