Extracting DNA from soil or directly from isolated nematodes indicate dissimilar community structure for Europe-wide forest soils

Extracting DNA

Nematodes are numerous in soils and play a crucial role in soil food-webs.DNA metabarcoding offers a timeeffective alternative to morphology-based assessments of nematode diversity.However, it is unclear how different DNA extraction methods prior to metabarcoding could affect community analysis.We used soils with woody vegetation from a European latitudinal gradient (29 sites, 39 to 79 • N, ~4500 km, covering six biomes) to systematically evaluate the effect of two sources of nematode DNA either directly extracted from soils vs. extracted from nematodes previously isolated from soils hypothesizing that the DNA source material may produce different diversities, community structures and abundances of feeding types.Nematode-sample DNA exhibited a higher richness, while no difference in Shannon diversity was found between the approaches.The DNA sources also created significantly different community structures, with greater differences observed across soil-extracted DNA than nematode-sample DNA.The most overrepresented species in nematode-sample DNA were Heterocephalobus elongatus, Eucephalobus striatus and Hexatylus sp., whereas Phasmarhabditis sp. and

Introduction
Among soil biota, nematodes are important components of soil food webs (Bardgett and van der Putten, 2014;Bar-On et al., 2018;van den Hoogen et al., 2019) and play a critical role in the below-ground mineralization of plant-derived organic C and N by feeding on various food sources and releasing nutrients via their excrements (Wardle, 2006;Osler and Sommerkorn, 2007;Briones, 2014;Andriuzzi et al., 2016).Nematodes cover many of the main trophic groups or feeding types of soil fauna, i.e., bacterivores, fungivores, omnivores, predators and root herbivores (Yeates et al., 1993;Neher, 2001;Yeates, 2003) and are valuable indicators of the complexity of soil food web structure and soil quality (Bongers and Ferris, 1999;Ferris et al., 2001).In recent years, they have received increasing attention, particularly in agricultural soils (Stone et al., 2016b;Treonis et al., 2018;Bongiorno et al., 2019).Nematodes respond quickly and in a taxon-specific manner to environmental changes, and are therefore widely used as bioindicators to reflect changes in the soil environment (Bongers and Ferris 1999;Neher, 2001;Yeates, 2003;Eisenhauer et al., 2012).
In order to develop nematode-based guidelines for policy makers and sustainable land use at the international scale, it is pivotal to take into account the biogeography of nematode diversity and community composition to locally calibrate the nematode indicators (Stone et al., 2016a(Stone et al., , 2016b)).Survey-based sampling along latitudinal or altitudinal gradients allows exploring potential shifts in nematode community diversity and composition in response to environmental gradients, as has been shown for other soil organisms, such as bacteria, fungi and protists (Bates et al., 2013;Tedersoo et al., 2014;Bahram et al., 2018;Delgado-Baquerizo et al., 2018).To date, however, few assessments of nematode communities at large geographical scales are available and most of them are of relatively limited sample size.Existing biogeographical studies indicate that climate, vegetation and soil abiotic conditions were the main drivers shaping soil nematode communities (Nielsen et al., 2014;Sylvain et al., 2014;Chen et al., 2015;Song et al., 2017;van den Hoogen et al., 2019;Wilschut et al., 2019;Li et al., 2020).In a comprehensive global study, including 6759 samples across all terrestrial biomes and continents, van den Hoogen et al. ( 2019) found a peak of nematode abundances at high latitudes, associated with large soil carbon stocks.The abundance of feeding types also followed the spatial patterns in total abundances, although some distinct drivers for single feeding types were revealed: herbivore-dominated communities were most influenced by vegetation, whereas bacterivores were most influenced by soil edaphic factors (i.e.sand content and pH).
Traditionally, intact nematodes are isolated from soil and then the bulk nematode samples are classified based on morphological traits.However, it is known that methods to isolate intact nematodes from soil create a bias, with small nematodes as well as large nematodes with a long body being discriminated by these methods (Viglierchio and Schmitt, 1983;Verschoor and de Goede, 2000;EPPO, 2013;Cesarz et al., 2019).Recently, DNA metabarcoding has increasingly been used to assess the diversity of nematodes previously isolated from soils, representing a less laborious alternative to morphological quantification (Powers et al., 2011;Darby et al., 2013;Treonis et al., 2018;Bell et al., 2021).Currently, the 18S rRNA gene is the most widely used molecular marker for nematode identification (Porazinska et al., 2009;Sapkota and Nicolaisen, 2015;Avó et al., 2017;Treonis et al., 2018;Ahmed et al., 2019), although other markers, such as the cytochrome c oxidase subunit I (COI), have also successfully been used (Derycke et al., 2010;Macheriotou et al., 2019).Studying nematode communities using DNA directly extracted from soils (Griffiths et al., 2018), instead of extracting DNA from animals that were first isolated from soil, is less investigated, despite this method being a very promising tool that has successfully been implemented in several ongoing soil biodiversity monitoring programmes that target soil bacteria, fungi, protists and multicellular animals using metabarcoding (Zinger et al., 2016;George et al., 2019;Gschwend et al., 2021aGschwend et al., , 2021b)).
To date, the effect of the source material for DNA extraction (i.e., direct extraction of DNA from soil versus extraction of DNA from previously isolated nematode sample) on the assessment of nematode diversity and community composition has only been investigated in one study of limited sample size (Griffiths et al., 2018).In this study, DNA extracted from a nematode sample (nematode-sample DNA) was shown to contain a higher fraction of nematode sequences, and therefore yielded a higher diversity of nematode taxa (Griffiths et al., 2018).However, as a disadvantage, this method is laborious and may introduce biases, such as a size-based discrimination of some taxa, which can be circumvented by directly extracting DNA from soils (Griffiths et al., 2018).It therefore remains unclear how differences between these sources of DNA could affect the quantification of nematode diversity and community structure at larger spatial scales.This is relevant because global patterns revealed using nematode samples isolated from soil samples, as in van den Hoogen et al. (2019), might not match those inferred from DNA directly extracted from soil.To our knowledge, no studies have investigated how differences in the source material of DNA extraction could affect the characterization of nematode communities across biomes.Calibrating nematode information derived from the two approaches across large spatial scales is however essential to reliably use soil nematodes as bioindicators for soil food web structure and soil quality (Geisen et al., 2018;Griffiths et al., 2018).
In the present study, we used DNA metabarcoding on a comprehensive set of soil samples collected along a European latitudinal gradient (from 37 • N to 79 • N, ~4500 km) including 18 countries, covering six biomes and a wide range of climatic and edaphic conditions.We characterized nematode communities based on DNA extracted from bulk nematode samples previously isolated from soil ("nematodesample DNA") or DNA directly extracted from soil ("soil-extracted DNA").Our study has three main goals: (1) to compare two sources of nematode DNA obtained by two different extraction approaches prior to metabarcoding of 18S rRNA genes at biome scale, (2) to explore if nematode richness and community composition show geographical patterns in European forests, and (3) to evaluate whether the richness and composition of communities can be predicted using environmental data such as biome, climate and soil type.Specifically, we want to reveal whether the different approaches in DNA extraction could result in significant differences in the estimates of nematode diversity and community structure, and thus affect the ecological conclusions of largescale field studies.Finally, besides looking at the diversity and composition of communities, we also examined the abundances of nematode feeding types to test whether they were affected by the source material of the DNA extraction.

Sampling sites
We conducted a continent-scale field study of soils containing woody vegetation (forests) analysing soil nematode communities in six biomes (alpine, Arctic, Atlantic, boreal, continental and Mediterranean) across Europe (Fig. 1).We define forest here as terrain containing woody vegetation.For the Arctic sites in Svalbard we selected spots covered by Salix arctica."Alpine" refers to mountain regions below the treeline.At least three sites were included in each biome (Table 1), yielding a total of 29 forest sites.The sites cover a wide range of environments across a latitudinal gradient from 39 • N to 79 • N, a distance of ~4500 km.Most of these forest sites belong to existing long-term national monitoring programmes (e.g.LTER, ICP Forests, ICP IM), and detailed information is available on their soil characteristics (pH, texture and organic matter content), vegetation (species diversity, standing biomass), local climate, and forest management (Cools et al., 2014;De Vos et al., 2015;Martins da Silva et al., 2016;Stone et al., 2016a).The soils from the sites have been classified according to well-accepted classification schemes (Zanella et al., 2011).For each site, we extracted the mean values of mean annual temperature ( • C) and total annual precipitation (mm) from the nearest weather station.The mean annual temperatures range from − 5.2 to 16.5 • C across the sites, the mean annual precipitation from 375 to 1357 mm, and the elevation from 20 to 1900 m above sea level (Table 1).

Field sampling
Soil sampling was carried out in 2017 according to a standardized protocol (Stone et al., 2016a).At each site, four plots (1 m 2 each) at a distance of 10 m from each other were selected in areas with no understorey vegetation and similar soil type to minimize the heterogeneity.In each plot, five soil cores (5 cm Ø × 10 cm deep) were taken at short distances from each othertwo cores for the isolation of nematodes and three cores pooled and homogenized for direct soil DNA extraction and soil physico-chemical analyses.This resulted in 8 cores for nematode isolation and 12 cores for soil DNA extraction and soil physico-chemical analyses per site.It is a common approach to use separate soil cores for different analyses that require different processing in nation-wide soil monitoring programmes (George et al., 2017(George et al., , 2019)).Litter was also collected from the four plots at each site.Continental and Mediterranean sites, but also the Atlantic sites that experience drought, were sampled in late spring (May), while the sites regularly exposed to winter frost (Arctic, Atlantic, alpine and boreal) were sampled in early summer (June and July).In each country, samples were transported to a local laboratory in cooler boxes containing ice packs and then immediately shipped overnight to the laboratory at WSL.

Abiotic soil variables
Immediately upon arrival in the laboratory, the three cores for soil DNA extraction and soil physico-chemical analyses were homogenized by sieving them using stainless steel sieves (mesh size 4 mm).Sieves were cleaned between samples by rinsing them under tap water and spraying them with ethanol.An aliquot (>25 g) of the sieved soils was frozen at − 20 • C for direct soil DNA extraction.The basic soil properties such as texture, pH and organic matter content, as well as total C and N concentrations in both the soil and litter, were determined after drying the samples for 48 h at 60 • C. Soil particles were fractionated into sand, silt and clay using the pipette method (Gee and Bauder, 1986).Total C and N concentrations were analysed using fine-ground samples by dry combustion using a CN analyser NC 2500 (CE Instruments, Wigan, United Kingdom).Inorganic carbon in the samples with a pH > 6.5 was removed with acid vapour prior to the analysis of organic carbon (Frey et al., 2021).Bulk density was calculated by dividing the weight of the dried fine fraction (<2 mm) by the volume of the full soil sample.The pH of dry soil samples was measured potentiometrically in 0.01 M CaCl 2 (soil:solution ratio = 1:2; 30 min equilibration time).

Nematode isolation from soil cores
The soil cores for nematode isolation and DNA extractions were processed immediately upon arrival.Nematodes were isolated from the soil cores (8 cores per site = 4 replicated plots x 2 cores per plot) using Oostenbrink dish method (Oostenbrink, 1960;Verschoor and de Goede, 2000).The Oostenbrink dish method is a standard method proposed from EPPO (European and Mediterranean Plant Protection Organization) to isolate nematodes from soils prior to DNA metabarcoding for community analysis (EPPO, 2013).The soil of one core was mixed and split in 2 parts (approximately of 200 g of soil each) before these were placed on a double cotton milk filter (FT25 Sana Vliesstoff-Filter, Zeltner Systemtechnik AG, Fulenbach, Switzerland) on a sieve in a dish.Afterwards the dish was gently filled with tap water from the side until the bottom of the sieve just reached the water.The nematodes were allowed to migrate through the filter into the water at 23 • C.After 72 h the nematodes of both aliquots were collected from the dishes on a 20-μm sieve and transferred to centrifugation tubes.This resulted in relatively clean suspensions for nematode counting.All nematodes in 3 x 1-ml subsample of the 10 ml suspension were counted at 40 × magnification under a microscope.After counting, the subsamples were put back complementing the 10 ml suspension which was used for DNA extraction.Counts of all nematode taxa were extrapolated to the entire sample and expressed per 100 g dry soil for further analyses.

DNA extraction and amplicon sequencing
"Nematode-sample DNA" was extracted from Oostenbrink-extracted nematodes using the DNeasy Blood & Tissue Kit (Qiagen, Hilden, Germany), as previously described (Resch et al., 2022)."Soil-extracted DNA" was extracted directly from 10 g of soil using the PowerMax Soil DNA extraction kit (Qiagen), as previously described (Frey et al., 2016;Resch et al., 2021).DNA was quantified using high-sensitivity Qubit assay (Thermo Fisher Scientific, Reinach, Switzerland) and diluted to 2.69 ng uL − 1 for PCRs.Subsequently, we amplified a 360-bp fragment of the 18S rRNA gene (Primers NF-1 and 18Sr2b; V6-V8 region; Porazinska et al., 2009).These primers pairs give good coverage of soil nematodes and have been widely used, but are not nematode specific and also amplify other eukaryotes.The nuclear small subunit ribosomal RNA gene (18S) was chosen as the most commonly used marker in nematode barcoding due to the availability of extensive database information (Macheriotou et al., 2019;Ahmed et al., 2019).The 5′ ends of the primers were tagged with the CS1 (forward) and CS2 (reverse) adapters required for multiplexing samples using the Fluidigm Access Array System (Fluidigm, South San Francisco, CA, USA).The PCR conditions to amplify the 18S RNA gene fragments consisted of an initial denaturation at 95 • C for 2 min, 35 cycles of denaturation at 95 • C for 40 s, annealing at 58 • C for 40 s, and elongation at 72 • C for 1 min, followed by a final elongation at 72 • C for 10 min (Resch et al., 2022).Each sample was amplified in triplicates and pooled prior to purification with Agencourt AMPure XP beads (Beckman Coulter, Berea, CA, USA) and quantification with the Qubit 2.0 fluorometric system (Life Technologies, Paisley, UK).Amplicon pools were sent to the Génome Québec Innovation Centre at McGill University (Montréal, Canada) for library preparation, and paired-end sequencing (2 × 250 bp) was performed on the Illumina MiSeq v3 platform (Illumina Inc., San Diego, CA, USA).

Bioinformatics
Bioinformatics analyses were carried out using the Quantitative Insights Into Microbial Ecology 2 program (QIIME 2, ver. 2021.4, Bolyen et al., 2019).The raw paired-end FASTQ files were imported into the QIIME2 program and demultiplexed using a native plugin.Thereafter, the Cutadapt plugin was processed primer trimmed.The Divisive Amplicon Denoising Algorithm 2 (DADA2) plugin in QIIME2 was used to determine amplicon sequence variants (ASVs; Callahan et al., 2016;Callahan et al., 2017).To remove low quality bases and reads, we truncated forward reads at 240 bp and reverse reads at 215 bp and discarded reads not fulfilling quality criteria (maxN = 0, maxEE = 2).Taxonomic assignment was accomplished using the Naive Bayes q2-feature-classifier (Bokulich et al., 2018) in QIIME2 against the PR 2 database v4.14.0 (Guillou et al., 2013) with 99% as the species delineation clustering threshold.Raw sequences were deposited in the NCBI Sequence Read Archive under the BioProject accession identifier PRJNA562119 and PRJNA793632.Nematode taxa were assigned to feeding types (herbivores, fungivores, bacterivores, omnivores and predators) and C-p/P-p (colonizerpersister) classes using NINJA (Nematode INdicator Joint Analysis; https://shiny.wur.nl/ninja/,Sieriebriennikov et al., 2014) based on genus level taxonomy (from ASV table).Genera that could not be assigned were designated as unclassified.C-p/P-p classes divide nematodes in five groups based on their growth rate and therefore ability to colonize disturbed habitats for free living (C-p) and plant associated (P-p) nematodes, respectively (Bongers, 1990).These range from r-strategists (C-p/P-p 1) to K-strategists (C-p/P-p 5).

Statistical analyses
All statistical analyses were conducted in R version 4.1.2(R Core Team, 2021).All figures were created using the ggplot2 package (Wickham, 2009).The R script is provided in the Supplementary material.P-values <0.05 were considered significant.The effects of site, biome and latitude on the abundance of nematodes were calculated using a generalized linear (mixed) model (GLM or GLMM) with a negative binomial distribution.The best fitting distribution was chosen using the 'fitdist' function in the fitdistrplus package (Delignette-Muller and Dutang, 2015).For testing the effects of biome and latitude, the site was included as a random factor and P-values were calculated using a likelihood ratio test that compared the models against a null-model with the random factor only.
Using datasets that included nematode sequences only, we compared the alpha-diversity indices (species richness and Shannon index at the ASV level) of the two DNA extraction methods (DNA sources) across the biomes.For the calculation of α-diversity (ASV) indices, samples were rarefied to an even number of sequences using the phyloseq package (McMurdie and Holmes, 2013).Rarefying of data was only performed for the analysis of alpha-diversity indices.The effect of site on the observed ASV richness and Shannon diversity was estimated using a one-way analysis of variance (ANOVA).A Bartlett test and a Shapiro test were used to check the homoscedasticity and normality assumptions of the ANOVA.The effects of biome and DNA extraction method (DNA source material) were tested using a linear mixed model with site as a random factor.For these models, homoscedasticity was tested by plotting residuals against means and normality using qqnorm plot, respectively.Rarefaction curves were calculated using the vegan package (Oksanen et al., 2019).
Community structure was analysed based on Bray-Curtis distances using relative read abundances.Principal coordinate analyses (PCoAs) were performed using the 'ordinate' function in the phyloseq package (McMurdie and Holmes, 2013).A permutational multivariate analysis of variance (PERMANOVA) was conducted to assess differences in community structure using the 'adonis2' function in the vegan package.The effects of site and biome were calculated for data subsets corresponding to each DNA extraction method (DNA source).For site, unrestricted permutations were used, whereas for biome, samples from each site were kept together to account for the hierarchical structure of the data (site nested within biome).The effect of the DNA extraction method (DNA source) was calculated together with the effect of site for the full dataset.To determine the P-value for the effect of DNA extraction, samples were only allowed to permute within sites without permuting the sites.Environmental parameters were regressed against the axes of the PCoA using the 'envfit' function in the vegan package.The 'procrustes' function in vegan was used to compare differences in community structure across the dataset between the two DNA sources.
To reveal if particular taxa were over-or underrepresented for either of the extraction methods (DNA sources), we calculated log 2 -fold changes in the abundance of the species in soil-extracted DNA over the abundance of species in the nematode-sample DNA.For the calculation of log 2 -fold changes between the two extraction methods (DNA sources) we used the DESeq2 package (Love et al., 2014;Donhauser et al., 2021).Normalization (i.e.calculation of size factors) was done at the ASV level.We calculated log 2 -fold changes between extraction methods on counts aggregated at the species level, controlling for the effect of site.
As the relative abundances of the nematode feeding types derived from the ASV taxonomic assignment were strongly heteroscedastic and not normally distributed, GLMs/GLMMs were used to test their responses to site, biome and DNA extraction method (DNA source).The effect of site was tested using separate data subsets for each extraction method (DNA source), the effect of biome using site as a random factor in data subsets of each extraction method (DNA source), and the effect of extraction method (DNA source) using site as a random factor over the full dataset.To this end, count tables were normalized using the normalization from the DESeq2 package at the ASV level rounded to integers.The function 'fitdist' from the fitdistrplus package (Delignette-Muller and Dutang, 2015) was then used to decide the best fitting distribution, which was the negative binomial distribution for all feeding types.The relationship between environmental variables and feeding type abundances was tested with a separate GLMM model for each variable with site as random factor in a subset for each DNA extraction method (DNA source).P-values were calculated using a likelihood ratio test that compared each of the models against a null-model with the random factor only.

Soil properties and total nematode abundance
The sampled forest soils showed a wide range of properties (Table 1), with pH ranging from 3.2 to 7.4, total C from 2% to 33.7%, total N from 0.07% to 1.5%, sand content from 13.9% to 95.7%, silt content from 2.4% to 62.3%, and clay content from 1.9% to 42.7%.Soil pH was highest in Mediterranean samples and lowest in Atlantic and boreal samples.Total C and N were both highest in alpine samples and C was lowest in Mediterranean samples, whereas N was lowest in the Arctic samples.Sand content was lowest in the continental and highest in the Arctic soils, whereas silt and clay contents showed the inverse trends.Accordingly, the total abundance of nematodes ranged from 33 to individuals per 100 g of dry soil (Fig. 2) and varied significantly across sites (P < 0.001; GLM).Differences among biomes were not significant (P = 0.30; GLMM), but were lowest in the Atlantic sites and highest in the alpine sites.Abundances increased slightly with latitude (Fig. 2), but this effect was neither significant (P = 0.13, GLMM).

Effect of the DNA source on α-diversity
When the datasets of the two methods were rarefied together, ASV richness of the nematode-sample DNA (ranging from 16 to 34) was higher than the ASV richness of the soil-extracted DNA (ranging from to 25) (Fig. 3a, Table 2).Shannon diversity of the nematode-sample DNA (ranging from 1.6 to 3.1) did not significantly differ from the Shannon diversity of the soil-extracted DNA (ranging from 1.4 to 2.8) (Fig. 3b, Table 2).For the nematode-sample DNA, richness and Shannon diversity differed across sites, but not across biomes, whereas for the soilextracted DNA, richness and Shannon diversity differed across both sites and biomes, although site had a stronger effect (Table 2).We found similar patterns for diversity indices at the species level (Figs.S1a and  b).Species richness ranged from 4 to 14 for soil-extracted DNA and from 7 to 17 for nematode-sample DNA and was not significantly different between the two extraction methods opposed to ASV level.Shannon diversity ranged from 0.85 to 2.3 for soil-extracted DNA and from 0.97 to 2.3 for nematode-sample DNA.
Although much smaller numbers of nematode sequences were retained in the soil-extracted DNA than in the nematode-sample DNA, in most samples for both methods the number of sequences was sufficient to recover the full diversity of nematodes in each sample (Fig. 3c).However, the rarefaction curves reached a plateau at much lower ASV numbers (Richness) for the soil-extracted DNA than for the nematodesample DNA, indicating a much lower richness in these samples (no additional ASVs would be recovered with a higher sequencing depth) in these samples.The same pattern was found calculating rarefaction curves at the species level (Fig. S1c).

Effect of the DNA source and soil properties on community structure
Community structure differed across sites in both DNA extraction methods (Fig. 4a and b, Table 3), with a slightly stronger effect for soilextracted DNA (Table 3).Variation in community structure across sites (average distance from the median) was greater in the soil-extracted DNA (Fig. 4c, Table 3) compared to the nematode-sample DNA (0.63 compared to 0.61, F 1,217 = 5.4,P = 0.02) In the soil-extracted DNA, biome also had a significant effect on the community structure, whereas this was not the case in the nematode-sample DNA (Table 3).Of the soil properties, only total C concentration showed a significant relationship with community structure in the nematode-sample DNA (Fig. 4a), whereas in the soil-extracted DNA all nine soil properties had a significant relationship with community structure along the first two PCoA axes (Fig. 4b).Among these variables, pH, total soil C concentration, total soil N concentration, organic matter content and soil moisture showed the strongest relationships.We also assessed the correlations among environmental variables and found strong positive relationships among total C, total N, OM and soil moisture, all of which showed a strong negative relationship with bulk density.Moreover, we found strong negative relationships between sand and silt as well as sand and clay (Fig. S2).

Effect of extraction method (DNA source) on species abundance
We found 72 species whose abundance differed (P adjusted < 0.05) between the two DNA extraction methods (DNA source): of these, five were overrepresented and 67 were underrepresented in the soilextracted DNA.Among the 50 most abundant species showing a significant log 2 -fold change, Phasmarhabditis sp., Eumonhystera filiformis, Eumonhystera sp. and Prodesmodora circulata in the phylum Chromadorea were more abundant in the soil-extracted DNA, while Heterocephalobus elongatus, Eucephalobus striatus and Hexatylus sp.(synonym Neotylenchus) in the phylum Chromadorea and Prionchulus muscorum in the phylum Enoplea were the most overrepresented species/genera in the nematode-sample DNA (Fig. 5a).As nematode extraction from soil may exclude taxa of certain body size, we tested the relationship between log 2 -fold change and size (weight derived from the NINJA tool but not length was used as a proxy for size), but found no evidence of body size affecting the log 2 -fold change (Fig. 5b).Moreover, we did not find a relationship between persister-colonizer groups (C-p/P-p) and differential abundance between the extraction methods.Thus, underrepresented species for soil-extracted DNA or overrepresented species for nematode-sample DNA comprised all persister-colonizer groups, such as Rhabditis brassicae (C-p 1), Heterocephalobus elongatus (C-p 2), Prodesmodora (C-p 3), Prionchulus muscorum (C-p 4) and Xiphinema diversicaudatum (P-p 5).
We also tested if there was a relationship between the rank abundance of nematode species and relative representation between the two extraction methods, which could arise as a result of different sampling depths between the methods.We found that species with a lower rank abundance were more strongly underrepresented in the soil-extracted DNA compared to the nematode-sample DNA (Fig. S3).This effect was strongest for rank abundances in the soil-extracted dataset (R 2 = 0.11, F 1,70 = 8.9, P = 0.004), but was also significant for rank abundances in the nematode-sample dataset (R 2 = 0.67, F 1,70 = 139.8,P < 0.0001) and in both datasets together (R 2 = 0.16, F 1,70 = 13.24,P = 0.0005).
When we compared the read abundance of the seven most abundant species in the two DNA extraction methods (DNA sources) (Fig. S4), we found that the relative abundance of Eumonhystera sp. was consistently lower in soil-extracted DNA across all sites and that the relative abundance of Rhabditis terricola was higher in these samples in almost all sites (Fig. S4).For the other five most abundant species, the over-or underrepresentation of either extraction method was less clear across the sites.Relative read abundances of species varied among sites within biomes, and few differences were found between biomes in either extraction method.Only Tylolaimophorus typicus was clearly more abundant in the Atlantic and continental sites than in the alpine, Arctic and Mediterranean sites for both extraction methods (DNA sources).Aphelenchoides saprophilus showed the highest relative abundances in the Atlantic and continental sites in the nematode-sample DNA, but the highest values in alpine sites in the soil-extracted DNA, thus highlighting the potentially significant role of the DNA source material when examining species read abundances across habitats.

Effects of the DNA source and soil properties on the abundance of feeding types
In the nematode-sample DNA and soil-extracted DNA, 80% and 83% of the sequences, respectively, could be classified into feeding types using the NINJA tool.Bacterivores constituted the majority of the sequences in both the nematode-sample DNA and soil-extracted DNA (45% in both), followed by fungivores (19% and 15%, respectively).The extraction method had a highly significant effect on the abundance of all feeding types (Table 4).Site had a significant effect on the abundance of all feeding types for both methods (Table 4).Biome had a significant effect on the abundance of ectoparasitic and endoparasitic herbivores in the soil-extracted DNA, but an effect on the abundance of ectoparasitic nematodes only in the nematode-sample DNA (Table 4).The abundance of bacterivores increased with latitude in the DNA of both extraction methods, but this trend was significant for soil-extracted DNA only (Fig. 6, Table 5).In contrast, omnivores decreased significantly with latitude in the nematode-sample DNA.Overall, the variation in feeding type abundances among replicates within a site was large (Fig. 6).
When assessing the role of soil properties, we found 13 significant associations between feeding types and environmental variables in the nematode-sample DNA and 16 significant associations in the soilextracted DNA (Table 5).For nematode-sample DNA, most of the significant relationships appeared between the abundance of ectoparasitic herbivores and all variables except latitude, and between the abundance of omnivores and latitude.For soil-extracted DNA most of the associations were observed for bacterivores (with C soil concentration, C litter concentration, soil moisture and latitude) and for fungivores (with N litter concentration).Furthermore, significant relationships were found for ectoparasitic herbivores (with C soil concentration, C litter concentration, N soil concentration, N litter concentration and soil moisture) and for omnivores (with pH, C litter concentration, N soil concentration, N litter concentration, bulk density and soil moisture) (Table 5).

Total nematode abundance
We observed large variation in nematode abundances across sites, but no significant difference among biomes.Overall, our abundances were of a similar order of magnitude as observed in earlier studies (Song et al., 2017;van den Hoogen et al., 2019).These studies reported a peak in nematode abundances either at approx.50 • N (Song et al., 2017) or in Arctic regions (van den Hoogen et al., 2019), whereas we did not find a clear latitudinal trend.One possibility for the lack of a latitudinal trend is that we did not include sites south of 30 • N. Also, while in our study the highest nematode abundance was found in the Subarctic site (70 • N), the two high-Arctic sites (78 • N) showed the lowest abundances in the whole dataset.This contrast of abundances in the three northernmost sites suggests that latitude may not be the main driver of nematode abundances.Nematode communities are also known to differ among land-cover types (Neher et al., 2005;Li et al., 2020), and one reason for our results differing from the trends reported earlier might be that we entirely focused on forest soils.

Effect of DNA source on α-diversity
Using universal 18S rRNA primers, the proportion of nematode sequences in the DNA extracted directly from soil (soil-extracted DNA) was 5-10% while the proportion was 85% in the DNA extracted from nematode samples that were previously isolated using the Oostenbrink dish method (nematode-sample DNA).These results confirm those reported by Griffiths et al. (2018) who used the same primers.With direct soil DNA extraction a large amount of sequences were of non-nematode origin, which highlights the need for the development of nematode-specific primers (Sapkota and Nicolaisen, 2015;Peham et al., 2017;Sikder et al., 2020;Kawanobe et al., 2021).
After rarefying to an even sequencing depth, we observed a slightly higher richness of nematode ASVs in the nematode-sample DNA than in the soil-extracted DNA, which may result from the much larger amount of soil in the former method (see discussion below).The richness of ASVs was comparable to values reported by Griffiths et al. (2018) for mid-latitude agricultural soils based on operational taxonomic units (OTUs).We found a significant effect of biome on nematode species richness in the soil-extracted DNA, with higher values recorded in Arctic sites.This finding is in contrast to the distribution of richness, with a minimum at mid latitudes, observed by Song et al. (2017), and to the much lower richness in Antarctic than temperate sites found by Nielsen et al. (2014).Unlike species richness, Shannon diversity did not differ between the extraction methods (DNA source), indicating that a higher evenness in the soil-extracted DNA compensated for the lower richness detected.These results suggest that the extraction method (DNA source) modifies the rank abundance distributions of nematode ASVs, which could be due to the Oostenbrink nematode extraction selectively enriching the sample with particular ASVs, for instance by discriminating against taxa with long bodies or relatively immotile nematodes, as reported previously (Verschoor and de Goede, 2000;EPPO, 2013).We did not find a relationship between weight as a proxy for size and overrepresentation of taxa (see below), suggesting that activity and life status may be more important.Alternatively, lower amounts of soil used for the soil-extracted DNA may lead to less rare species in the sample resulting in more even rank abundance distribution among more frequent species.
When samples were not rarefied to an even sequencing depth, the nematode-sample DNA had a much higher richness than the soilextracted DNA although the rarefaction analysis confirmed the sequencing depth after rarefying being sufficient to cover the diversity in each sample for both methods.Our results thus suggest that the true richness was higher in the nematode-sample DNA, which likely reflects the much larger amount of soil used for isolating nematodes from soil (up to 200 g) than for extracting the DNA from soil (10 g).Accordingly, Wiesel et al. (2015) also reported different community compositions depending on the amount of soil used for extraction.For methods extracting intact nematodes, >100 g soil is typically used (Wiesel et al., 2015), whereas the amount of soil used for DNA extraction directly from soil is typically 0.25-10 g (Peham et al., 2017;George et al., 2019).Our results suggest that 10 g soil is not enough to cover the local diversity such that more taxa are recovered using larger amounts.This highlights the need to develop direct soil DNA extraction methods that permit using larger amounts of soil.For instance, existing protocols for soil amounts as large as 100 g (Yeates et al., 1998;Taberlet et al., 2012;Gorny et al., 2018) could be combined with a commercial purification kit to remove PCR-inhibiting soil substances that can hamper amplification of soil DNA (Peham et al., 2017).In addition, the relationship between the recovery of additional species and the amount of soils should be studied more systematically, to determine the minimal amount sufficient to represent nematode diversity and to adapt DNA extraction protocols accordingly.

Effect of DNA source on community structure and ecological relationships
We found significant differences in nematode community structure among sites with both extraction methods (soil-extracted DNA vs.nematode-sample DNA), while the effect of biome was significant in the soil-extracted DNA only.Interestingly, differences in community structure across the whole dataset were much smaller for the nematodesample DNA than for the soil-extracted DNA.One explanation is that nematode isolation from soil leads to converging community structures in dissimilar samples.This could be a result of size and/or immotile nematode exclusion effect, as stated above (Verschoor and de Goede, 2000;EPPO, 2013), and indicates that across samples similar taxa are favoured by the extraction.Alternatively, lower amounts of soil (10 g of soil) for soil-extracted DNA may have resulted in an incomplete representation of the nematode community in each sample compared to the classical isolation of intact nematodes (from around 200 g of soil) and thus led to differences between samples.
With both extraction methods (DNA sources), total soil C concentration appeared as an important driver of nematode community structure, which is in agreement with findings in previous studies (Stone   et al., 2016b;Bongiorno et al., 2019;van den Hoogen et al., 2019).Using direct extraction of DNA from the soil, in accordance with previous studies (Nielsen et al., 2014;van den Hoogen et al., 2019), we further found soil texture (sand, silt and clay percentages) and bulk density to associate with nematode community structure.As nematodes are too small to move soil particles, their range of motion depends on pore spaces and aggregate structure, which in turn depend on the soil texture (Sechi et al., 2018).Soil structure likely selects nematode groups based on their size and feeding preferences, as the structure of soil influences access to food sources (Mikola and Sulkava, 2001;Neher, 2010).In addition, we found soil moisture to be strongly associated to community structure.Nematodes are aquatic organisms living in water films surrounding soil particles, and moisture also influences the availability of food sources, such as bacteria, fungi and plant roots (Neher, 2010).Our results thus confirm previous findings that soil moisture and precipitation can be important drivers for nematode communities (Nielsen et al., 2014;Sylvain et al., 2014;Chen et al., 2015;Wilschut et al., 2019) and that response of soil moisture to temperature might explain the temperature effects across latitudinal gradients (Song et al., 2017).

Selection of nematode species by the DNA source
In our comprehensive dataset, most of those nematodes that differed in read abundance between the two DNA extraction methods (DNA sources) were more abundant in the nematode-sample DNA.In the class Chromadorea, those species that differed between the extraction methods were predominantly small, with a weight <1 μg.Exceptions with a greater body weight were Rhabditis brassicae, three species belonging to the genus Anaplectus and Pristionchus lheritieri, all of which were underrepresented in the soil-extracted DNA, and Phasmarhabditis sp., which was overrepresented in the soil-extracted DNA.In the class Enoplea, no overrepresented taxa were found in soil-extracted DNA.In contrast, this class had many underrepresented taxa of larger body weight, the heaviest being Anatonchus tridentatus.Overall, we did not find an unequivocal relationship between the size of species (using weight as a proxy) and difference in abundance between the two extraction methods, i.e., species of similar size were both over-and underrepresented.If the isolation of intact nematodes prior to DNA extraction did select against taxa of certain sizes, one possible reason for not finding a clear difference in their relative abundances between the two extraction methods might be that both very large and very small nematodes can be excluded during the different steps of the Oostenbrink method (Verschoor and de Goede, 2000;EPPO, 2013).Another possible explanation is that the bias between the two extraction methods is related to the motility of nematodes rather than their size: in the final step of the Oostenbrink dish method, nematodes migrate actively through a cotton-wool filter, such that eggs, nematodes in non-motile resting stages, and relatively immotile or otherwise inactive nematodes are excluded (Verschoor and de Goede, 2000;EPPO, 2013;Cesarz et al., 2019).Conversely, soil-extracted DNA targets any nematode DNA, independent of life stage and activity status.
We found a strong relationship between the rank abundance and the differential abundance of taxa between the two methods, i.e. rare taxa were more underrepresented in the soil-extracted DNA.This could be a result of lower amounts of soils and a large fraction of sequences of nonnematode origin in the soil-extracted DNA approach leading to a lower sensitivity to detect rare taxa.Consequently, rare taxa seem to be underrepresented compared to the nematode-sample DNA.These findings additionally emphasize the need to develop methods to extract DNA from larger amounts of soil and to design nematode-specific primers.
In agreement with findings from Griffiths et al. (2018), we found a lower abundance of taxa affiliated with the families Aporcelaimidae, Cephalobidae, Paratylenchidae, Plectidae and Tylenchidae and a higher abundance of taxa affiliated with the family Monhysteridae in the soil-extracted DNA.However, opposite patterns were observed for taxa within the families Aphelenchoididae, Diphtherophoridae and Microlaimidae.One reason for the contrasting results could be the taxonomic resolution, i.e., family level in Griffiths et al. (2018) and species level in our study.For instance, within the family Rhabditidae, we found both over-and underrepresented species in soil-extracted DNA.Differences could also be attributed to different soil types (Verschoor and de Goede, 2000;Cesarz et al., 2015), the sequencing depth and platform used, the bioinformatic processing of the sequences, and the database used for taxonomic classification.

Effects of the DNA source on nematode feeding types
Nematode feeding types are important indicators of soil food-web dynamics and have been widely studied across different soil habitats (Nielsen et al., 2014;Stone et al., 2016b;Song et al., 2017;Treonis et al., 2018;van den Hoogen et al., 2019;Vazquez et al., 2019;Wilschut et al., 2019;Li et al., 2020).In this study, we found a highly significant difference in the abundance of all feeding types between the two extraction methods (DNA sources).In agreement with findings from previous studies, bacterivores were the most abundant group with both extraction methods, followed by fungivores, whereas omnivores and predators

Table 2
Effects of site, biome and extraction method on nematode ASV richness and Shannon diversity index.The effects of site and biome were assessed in separate subsets for the two DNA extraction approaches.To test the site effect, a one-way ANOVA was used, and to test the effects of biome and the DNA extraction method, linear mixed models with site as a random factor were used.Samples with <100 sequences were removed: when the analyses were conducted for data subsets of the two DNA extraction methods, samples with <100 sequences were removed for the respective method only, and when both DNA extraction methods were included, samples with <100 sequences were removed for both methods.For ANOVA results, subscripted indices indicate numerator and total degrees of freedom; for results of the linear mixed models, subscripted indices indicate numerator and denominator degrees of freedom.Significant values are in bold letters.N sites = 29 (4 alpine, 3 arctic, 5 atlantic, 4 boreal, 9 continental, 4 mediterranean), n = 4. were the least abundant (Bongiorno et al., 2019).The abundance of feeding types was highly variable among sites within the same biome, and also across biomes.In the soil-extracted DNA, herbivores (ecto-and endoparasites) differed significantly among biomes, whereas in the nematode-sample DNA, this was the case for ectoparasitic herbivores only.Van den Hoogen et al. (2019) reported that the absolute abundance of all feeding types followed the total abundance of nematodes with a peak at high latitudes, indicating that the relative abundances of different feeding types did not exhibit a latitudinal trend in their study.Conversely, we found an increase in the proportion of bacterivores with increasing latitude in the soil-extracted DNA.In accordance with Song et al. (2017), we found soil moisture to be significantly associated with the abundance of bacterivores.In addition, we also found a significant association between soil and litter C concentration and bacterivore abundance as well as between litter N concentration and fungivores.These soil properties have previously been identified as important drivers of bacterial and fungal communities (Delgado-Baquerizo et al., 2016;George et al., 2019;Smith et al., 2021), thus having potential to affect the abundance of bacterivores and fungivores.The effect of soil properties on herbivores, which was revealed with both extraction methods, is likely to be mediated by the plants these nematodes feed on as plant characteristics have earlier been identified as primary drivers of herbivorous nematodes (van den Hoogen et al., 2019).Collectively, our findings suggest that feeding type abundances and their associations with environmental variables are influenced by the extraction method (DNA source), highlighting the need to take into account the potential methodological biases when making ecological interpretations of DNA metabarcoding data.

General implications
Across all analyses, we observed differences between the DNA extraction methods across the large spatial extent of our study.A previous methodological comparison that aimed at validating the use of molecular approaches for the ecological assessment of nematodes demonstrated a good agreement between morphological, qPCR-based and high throughput sequencing-based analyses of nematode communities (Geisen et al., 2018).As a result, the authors deemed molecular characterization of nematode communities a valid approach and highlighted the increased diversity and taxonomic resolution compared with the morphological approach.Notably, all the approaches compared were based on the Oostenbrink extraction.Overall, the molecular characterization of nematode communities is increasingly used, but the vast majority of studies have first isolated the nematodes from soil and only very few studies have included direct DNA extraction from the soil (environmental DNA).Griffiths et al. (2018) found differences in nematode communities depending on the extraction method applied at a small scale.Here, we were able to demonstrate that the choice of the DNA source material (nematode-sample DNA vs.soil-extracted DNA) can lead to significantly different ecological interpretations at a large spatial scale including sites with a wide range of edaphic and climatic conditions.This suggests that previous studies using isolation of intact nematodes from soil might have led to different conclusions if conducted using direct DNA extraction from the soil.One possible source of differences between approaches compared here, is the preferential extraction of some taxa in the Oostenbrink method depending on size and activity.(Verschoor and de Goede, 2000;EPPO, 2013).Such a bias could be significantly reduced using, a combination of Oostenbrink extraction method followed by a sugar centrifugal flotation including active, inactive and dead nematodes (McSorley and Frederick, 2004; Table 3 Effects of site, biome and extraction method on nematode community structure were tested using permutational multivariate analysis of variance (PERMANOVA) with 9999 permutations.The effect of site and biome were analysed in separate subsets for each extraction method, while the effect of extraction method was analysed using the full dataset.In the subsets, when testing the effect of biome, samples were not permuted across sites to account for increased similarity among replicates within sites.In the full dataset, the effect of extraction method was tested together with the effect of site, but only the effect of extraction is shown.Data was permuted in pairs of samples for nematode-sample DNA and soil-extracted DNA to account for the dependence between DNA extracted with the two different methods in the same plot.The function permutest with 9999 permutations was used to assess homogeneity of variance.Significant values are in bold letters.N sites = 29 (4 alpine, 3 arctic, 5 atlantic, 4 boreal, 9 continental, 4 mediterranean), n = 4.

Table 4
Effects of site and biome on the abundance of nematode feeding types for nematode-sample DNA and soil-extracted DNA, as well as the effect of the extraction method.The site effect was estimated using a generalized linear model with a negative binomial distribution, based on read abundances normalized at the amplicon sequence variant (ASV) level using DESeq2.The effects of biome and extraction method were estimated using a generalized linear mixed model with site as a random factor.Subscripted indices indicate numerator and total degrees of freedom.N sites = 29 (4 alpine, 3 arctic, 5 atlantic, 4 boreal, 9 continental, 4 mediterranean), n = 4. Fig. 6.Abundances of nematode feeding types as a function of latitude for the two extraction approaches (nematode-sample DNA vs.soil-extracted DNA).Amplicon sequence variants (ASVs) were classified into feeding types based on their taxonomic affiliation at the genus level.Abundances were normalized at the ASV level using DESeq2.Lines represent regressions using a generalized linear mixed model with significant regressions in solid.Mean ± SD (n = 4) for each site is shown.Cesarz et al., 2019;Gattoni et al., 2022).Soil-extracted DNA appeared to be biased by a lack of sensitivity to detect rare taxa resulting from low amounts of soil used for extraction and a large fraction of non-nematode sequences.Thus, this approach would significantly benefit from extraction protocols for large amounts of soils (>100 g) and from the need of designing nematode-specific primers.Geisen et al. (2018) pointed out that sequencing methods are not capable of retrieving full quantitative information on the abundance of taxa or functional groups, and for this reason, van den Hoogen et al. (2019) used morphological data in their global modelling study.Using taxon-specific qPCR primers, as used in some recent studies (Vervoort et al., 2012;Quist et al., 2016;Geisen et al., 2018) may contribute more quantitative information on particular groups of interest.Also, it is important to remember that the Oostenbrink extraction approach relies on active migration of individuals and thus targets only living nematodes, whereas direct DNA extraction from the soil includes both active and inactive/dead nematodes.Research questions targeting both life stages could be tackled by targeting ribosomal RNA that is rapidly degraded in dead cells, an approach that has been applied in numerous amplicon and metagenome-based microbial studies (Carini et al., 2020;Donhauser et al., 2020).
In conclusion, nematode communities assessed with two different molecular approaches (soil-extracted DNA and nematode-sample DNA) were comparable although with some differences in nematode community composition at a large spatial scale.Direct extraction of DNA from the soil can be a high-throughput screening approach to target nematode diversity and community structure, which can complement laborious classical methods.In particular, for monitoring programmes where DNA is directly extracted from soils to study bacterial and fungal communities using soil-extracted DNA (Gschwend et al. 2021a(Gschwend et al. , 2021b) ) will allow to analyse nematode community structure at a similar throughput and to assess cross-kingdom trophic structures.A key problem with direct soil extraction is that a large amount of sequences are of non-nematode origin, pointing to a need for developing nematode-specific primers (Kawanobe et al., 2021) and that low amounts of soil used for extraction reduced sensitivity to detect rare taxa.If this obstacle can be overcome, metabarcoding of nematode communities directly extracted from soil can be a powerful, efficient approach that could be implemented in large-scale monitoring programs using nematode-based indicators of soil health.

Fig. 1 .
Fig. 1.The location of the twenty-nine Europe-wide sites involved in the study.

Fig. 3 .
Fig. 3. Alpha-diversity.(a) Observed ASV richness and (b) Shannon diversity for nematode-sample DNA and soil-extracted DNA samples rarefied together.If a sample had <100 sequences for one of the extraction methods, it was removed for both methods.Mean ± SD per biome, N sites = 29 (4 alpine, 3 arctic, 5 atlantic, 4 boreal, 9 continental, 4 mediterranean), n = 4 per site.(c) Rarefaction curves for nematode-sample DNA (left) and for soil-extracted DNA (right).

Fig. 4 .
Fig. 4. Relationships between nematode community structure and environmental variables.(a) Nematode-sample DNA, (b) Soil-extracted DNA.PCoAs were conducted based on Bray-Curtis distances of relative read abundances.Vectors represent a regression of environmental parameters against the PCoA axes.Significant parameters are represented by solid lines.(c) Comparison of community structure across biomes for nematode-sample DNA and soil-extracted DNA.n = 4 per site.

Fig. 5 .
Fig. 5. (a) Log 2 -fold changes in species abundances between extraction methods (soil-extracted DNA vs.nematode-sample DNA), controlling for the effect of site.Significant log 2 -fold changes (P adjusted < 0.05) are shown.Normalization in DESeq2 was conducted at the ASV level.Error bars represent the estimated standard error for the log 2 -fold changes from the model (DESeq2).Positive log 2 -fold changes indicate that taxa were overrepresented in soil extracted DNA, while negative log 2 -fold changes indicate that taxa were overrepresented in nematode-sample DNA.The size of the black points represents the mean normalized read abundance of the respective species.The size of the grey points represents the weight of the species.C-p/P-p represents colonizer-persister classes for free-living and plant-associated nematodes.C-p 1/P-p 1 = highest growth rate; C-p 5/P-p 5 = lowest growth rate.(b) Relationship between log 2 -fold changes and weight of the species.The trendline was obtained by loess smoothing.N sites = 29 (4 alpine, 3 arctic, 5 atlantic, 4 boreal, 9 continental, 4 mediterranean), n = 4 per site.

Table 1
Geographical, climatic and soil physico-chemical properties of the forest sites.MAT = mean annual temperature, MAP = mean annual precipitation, BD = bulk density, OM = organic matter, SM = soil moisture.NA = not analysed; n = 4 per site.