Nematode species were raised on cultured agar plates (1.7%) dosed with cholesterol and E. coli and synchronised regarding life stages. In total, 18 species were used and organisms were individually removed with needles from plates and rinsed in a water drop in order to remove excess agar parts, bacteria and fungi. Organisms were then transferred into an Eppendorf tube containing lysis buffer (Machery & Nagel, Nucleo Spin Tissue XS Kit), according to predefined mock-community ratios (Table 1). Each tube contained 198–200 nematodes from the 18 species, resulting in 9 different mock-communities with 3 replicates for each mock-community. Mock-communities represented ratios that are found in nature (Peters and Traunspurger 2005; Kazemi-Dinan et al. 2014), as well as evenly distributed ratios in order to examine sequencing success. Subsequently, 8 µl of proteinase-K was added and samples were lysed at 56 °C within a rocking water bath for 8 hours. Subsequently, DNA was extracted, according to the NucleoSpin Tissue Kit XS protocol (Macherey & Nagel, Hilden). Prior to amplification, DNA concentration was checked with PicoGreen. Samples were amplified with three different primer pairs that had been commonly used for metabarcoding of microscopic animals previously. The 18S rRNA gene was amplified with F04/R22 (5′-GCTTGTCTCAAAGATTAAGCC-3′/5′-GCCTGCTGCCTTCCTTGGA-3′ (Fonseca et al. 2010), amplifying a ~350 bp fragment of the V1-V2 region and 3NDf/C_1132f (5’GGCAAGTCTGGTGCCAG 3’/5’TCCGTCAATTYCTTTAAGT 3’), amplifying a ~530 bp fragment of the V4 region (Geisen et al. 2018). A ~520 bp fragment of the D3-D5 region of the 28S rRNA gene was amplified with the primer pair 1274/706 (5’–GACCCGTCTTGAAACACGGA-3’/5’-GCCAGTTCTGCTTACC-3‘) introduced by Markmann and Tautz (2005), which has proven to work extremely well for freshwater nematodes (Ristau et al. 2013; Schenk et al. 2016). Primers pairs are subsequently simplified to 18S_V1 (F04/R22), 18S_V4 (3NDf/C_1132f) and 28S_D3 (1274/706). PCR conditions were 30 cycles with 60 s pre-denaturation at 96 °C, followed by 15 s at 96 °C, 30 s at 58 °C and 90 s at 70 °C. In a second PCR with the same conditions, Illumina indices were attached in 10 cycles. Metabarcoding was carried out at LGC Genomics (Berlin) on an Illumina MiSeq (V3 2 × 300 bp). About 20 ng of each sample were pooled for sequencing. Samples were delivered demultiplexed.

Individual species were amplified, Sanger-sequenced and used as reference in the bioinformatic pipeline (see below). PCR conditions followed Schenk et al. (2016). Sequences are deposited at NCBI GenBank under the Accession numbers: MK543220–MK543233, MK541669–MK541684, MK541653–MK541668. The closely related species Acrobeloides tricornus and Acrobeloides cf. nanus, Plectus aquatilis and Plectus cf. acuminatus and Plectus velox, as well as Acrostichus nudicapitatus and Acrostichus sp. could not be distinguished, based on their sequences for the 18S rRNA gene regions. For the 28S rRNA gene region, Plectus aquatilis and Plectus cf. acuminatus, as well as Acrostichus nudicapitatus and Acrostichus sp. could not be distinguished. However, as natural communities contain closely related species as well, this was also maintained in our study. For two species, Panagrolaimus thienemanni and Aphelenchoides parietinus, no Sanger reference sequence could be generated for any of the markers (Suppl. material 1: Table S1). For the latter, NCBI sequences were available for the 28S rDNA gene region (MF325173.1, MF325174.1)

Mean biomass (wet weight) of every species was calculated following (Andrássy 1956) for several individuals (n = 10–15 adults) that were heat-fixated and measured (Table 1). Furthermore, the proportion of generated reads for 18S_V4 and 28S_D3 was plotted against the relative biomass and the relative abundance together with linear regression in MatLab (MATLAB User’s Guide 1998). The correlation for absolute biomass and absolute abundance against the read numbers are given in Suppl. material 2: Fig. S1.

Except for the taxonomic classification and primer removal, the MiSeq standard operation procedure using mothur (Kozich et al. 2013; Schloss et al. 2009) was used. Demultiplexed reads were combined to reach a longer paired-end read for higher phylogenetic resolution by using the make.contigs function of mothur with default settings, thereby correcting sequencing errors in the overlapping region. Primer sequences were then removed from the combined reads using cutadapt with a default error rate of 0.1 (Martin 2011). Reads with ambiguous bases, homopolymers larger than 10 bases and of unexpected short or long length (allowed range: 333–367 for 18S_V1, 514–597 for 18S_V4 and 471–516 for 28S_D3) were removed from the dataset. An overview about distribution of read numbers in bioinformatics is given in Suppl. material 1: Table S2. The remaining sequences were aligned with the SILVA reference (release 132) alignment (Martin 2011) to determine the spanned 18S or 28S rRNA gene region within the alignment in order to remove sequences not spanning this region and to remove overhangs outside of this region as well. The sequences were then clustered with approximately 99% identity which corresponds to a maximum difference of 4 (18S_V1) or 5 (18S_V4 or 28S_D3) bases. An alternative clustering-free approach (using DADA2) by Callahan et al. (2016) was also tested, but as extremely similar results were obtained, this is not discussed further (Suppl. material 2: Fig. S2). Chimeras were then removed using UCHIME (Edgar et al. 2011) with a de novo constructed reference database. All these filter steps reduced the dataset from 448,687 to 404,187 (for 18S_V1), from 315,232 to 241,822 (18S_V4) and from 453,691 to 264,010 (28S_D3). The median of read length was 532-bp for 28S_D3, 583-bp for 18S_V4 and 405-bp for 18S_V1. OTUs were taxonomically classified based on top BLAST hits in the NCBI nt database using an identity cutoff of 97%. Contradicting equal scores were not encountered and OTUs represented by less than 10 reads were discarded. Besides reference sequences obtained from GenBank, Sanger-sequenced reference sequences of the mock-community species were added to support taxonomic classification of the taxa used in our mock-communities (via BLAST with 97% identity cutoff as well). These reference sequences were used instead of the NCBI nt based classification, when they had a higher similarity than the hits in the NCBI nt database or when the hits in the NCBI database were poorly annotated and contained the keywords “uncultured”, “environmental” or “unidentified”. Raw reads were deposited under PRJNA513984.

For each primer pair, read numbers for each species in each replicate were averaged to decrease variation and communities were compared to relative abundance and relative biomass, based on the number of inoculated specimens. For statistical analyses of communities, Bray Curtis similarity was applied with non-transformed data. Non-metric multidimensional scaling (NMDS) plots were created with MathLab (MATLAB User’s Guide 1998), while Permutational multivariate analysis of variance (PERMANOVA) was conducted with PRIMER-E v6 (Clarke and Gorley 2006). PERMANOVA was computed with 9,999 permutations and pairwise comparisons amongst relative abundance, biomass and the read proportions for the three primer pairs. The samples were used as the nested factor in the treatments (e.g. abundance, biomass or primer used). A list of all pairwise comparisons is given in Suppl. material 1: Table S3. The effect of amplicon length and mismatches on sequencing success was analysed with Analysis of Variance (ANOVA) in R (R Core Team 2013) in order to check if there was an significant effect by the two factors on the sequencing success. Sequencing success was defined as expected to encountered sequencing ratios based on biomass proportions.

The three primer pairs had a different amplification success (Table 2) with the 28S_D3 marker resulting in the highest recovery rate of species present in the mock-communities (n = 14), followed by 18S_V4 (n = 12), while the 18S_V1 returned the lowest OTU number (n = 10). Different amplification success could partly be attributed to mismatches in the primer regions, such as for Diploscapter coronatus, which could not be amplified within the metabarcoding approach for any of the markers, although it could be amplified with Sanger-sequencing (Suppl. material 1: Tables S1, S4). ANOVA showed that variations in length of the amplified region of different nematode taxa significantly affected amplification success of 18S_V4 (p < 0.001), the primer pair that differed most in length between nematode taxa (shortest amplicon of 492 bp in 531 bp). For the other two markers, amplicon length (28S_D3: p = 0.41, 18S_V4: p = 0.21) and mismatches (28S_D3: p = 0.51, 18S_V4: p = 0.6) did not significantly affect amplification success (28S_D3: 512 bp – 525 bp, 18S_V1: 333 bp – 364 bp, Suppl. material 1: Table S5). We focused our analyses on all sequences assigned to the nematode taxa that were added, as contamination with non-target sequences was low (0.5–2.1%, Table 2).

The community structure differed between the primer pairs, with 28S_D3 resulting in the most similar community depiction to the initial inoculated communities, while depiction for the 18S_V1 primer pair substantially differed. Nematode communities amplified with the 28S_D3 primer pair grouped together with nematode mock communities presented as biomass data in the NMDS plot (Fig. 1a), which was also supported by pairwise comparisons between 28S_D3 and biomass showing similar patterns (t-value = 0.8193, p = 0.5293). In contrast, nematode abundances differed from the community composition as depicted by the primer pair 28S_D3 (Fig. 1a) and supported statistically (t-value = 2.2339, p = 0.0076). The communities for the 18S_V4 primer pair were different from mock-community data presented as biomass (t-value = 1.9353, p = 0.0279) and even more to abundance (t-value = 3.6898, p < 0.001, Fig. 1b). Communities revealed with the primer pair 18S_V1 were most different to biomass and abundance (t-value = 5.9545, p < 0.001 and 6.7753, p < 0.001, Fig. 1c) and were therefore not evaluated further. The community structure of the inoculated communities differed when represented as biomass or abundance (t-value = 2.281, p = 0.0068). Primer pairs 28S_D3 and 18S_V4 revealed nematode communities in a similar way (t-value = 1.278, p = 0.1638). In contrast, all communities amplified with the primer pair 18S_V1 were strikingly different, being characterised by a clear dominance of few species, for example, Rhomborhabditis regina and Plectus aquatilis (Suppl. material 2: Fig. S3).

Figure 1.

NMDS plots for a) the 28S_D3 region, b) the 18S_V4 region and c) the 19S_V1 region depicting the 9 samples. Each sample is shown with three replicates, together with the relative proportion of biomass (black) and the relative abundance of the original community (cyan).

Relative biomass positively correlated with relative read numbers for the primer pair 28S_D3 (R² = 0.90539, p < 0.001) and 18S_V4 (R² = 0.81396, p < 0.001, Fig. 2 a–b). At the level of individual mock-communities, especially samples with dominant species were often positively correlated with relative biomass proportions for the 28S_D3 and 18S_V4 primer pair, for example, Sample 4 (28S_D3: R² = 0.9995, p < 0.001; 18S_V4:R² = 0.9656, p < 0.001). The same was found for Sample 7 (28S_D3: R² = 0.9885, p < 0.001 18S_V4: R² = 0.923, p < 0.001) or Sample 9 (28S_D3: R² = 0.9162, p < 0.001; 18S_V4: R² = 0.3069, p = 0.02597). In general, species with a high contribution to the total biomass resulted in the highest read proportions, while species with a low biomass were clearly less amplified. The biomass of juvenile and adult specimens differed by a factor of 5.00 (Suppl. material 1: Table S6 and Suppl. material 2: Fig. S4), while the amplification of different life stages of R. regina also resulted in varying read proportions, with the adult nematodes generating significantly more reads than the juvenile nematodes (Wilcoxon-Mann-Whitney test, W = 81, p < 0.001). For the 28S_D3 primer pair, 17.14 times more reads were generated, while the 18S_V4 primer pair resulted in 4.65 times more reads and 5.81 times

Relative species abundances were positively correlating with read numbers for the primer pairs 28S_D3 (R² = 0.53694, p < 0.001) and 18S_V4 (R² = 0.33077, p < 0.001), but less strong than for relative biomass-read number correlations (Fig. 2c–d). Correlations of absolute biomass and absolute abundance against absolute read counts showed a similar pattern, but with overall weaker correlation strengths (Suppl. material 2: Fig. S1).

Figure 2.

Correlations of relative biomass and relative abundance against relative read proportions. a–b Correlation of 239 relative biomass proportions in μg (x-axis) against the relative number of generated reads (y-axis). Given is the slope of the 240 correlation and the according R2-value for the 28S_D3 and the 18S_V4 marker. c–d Correlation of relative abundance 241 proportions in μg (x-axis) against the relative number of generated reads (y-axis). Given is the slope of the correlation and 242 the according R2-value for the 28S_D3 and the 18S_V4 marker. 24

A total of 78% of species diversity present in the samples was recovered for the 28S_D3 primer pair, while, for the other primers, resolution was lower (18S_V1:55%; 18S_V4: 67%). A higher diversity was recovered considering indistinguishable taxa (28S_D3 and 18S_V4: 89%; 18S_V1: 78%). Those that could not be distinguished shared the same marker gene sequence and, as such, cannot be distinguished, based on sequencing (Table 2). The rRNA genes were recently suggested to be incapable of identifying many nematode taxa to species and genus level resolution was consequently suggested as a reliable lower-resolution alternative (Creer et al. 2016; Sahraean et al. 2017). It should, therefore, be considered in metabarcoding studies of unknown species composition, when the full species diversity is aimed at identification in a sample. However, species that are closely related and, as such, share the same genetic sequence commonly can share similar morphological, physiological and functional traits (Potapov et al. 2019). For example, nematode species have proven to be equally sensitive to pollutants within the same genus (Höss et al. 2017). Therefore, the relative-read abundance links at lower taxonomic resolution, which likely links to functional information, still holds. From a total of 18 species used, one species could not be recovered with any of the primer-pairs and other species such as Aphelenchoides parietinus were under-represented in relative sequence reads. This can be mainly explained by mismatches present in primer sites. Accordingly, a 10 fold drop in amplification efficiency was reported to be caused by only one nucleotide mismatch in the priming region (Piñol et al. 2015). Due to amplicon length variations, Panagrellus redivivus (531 bp) might not be amplified with 18S_V4, as the primer length is at the maximum of the capacity our Illumina approach can sequence (2× 300 bp). This is supported by the fact that P. redivivus was found in unmerged forward and reverse reads (Suppl. material 1: Table S7). Furthermore, we demonstrate that longer amplicon sizes negatively influenced amplification success for the 18S_V4 marker (Suppl. material 1: Table S4).

The community depiction for this study varied due to primer-induced differences. The 18S_V1 marker was, in our case, not suitable for distinguishing communities due to a low species recovery and an extreme over-amplification of mainly two species (Rhomborhabditis regina and Plectus aquatilis) and is, therefore, not discussed further. In turn and in line with most accurate diversity representations, the 28S_D3 primer pair could reliably distinguish between communities, followed by the 18S_V4 marker. The difference for the V4-marker is likely introduced by higher length-variation in this marker (Suppl. material 1: Table S4). Especially, communities with dominant taxa were depicted accurately and uncertainties were mostly detected from less abundant taxa for both primer pairs. Many metabarcoding studies reported a problem in detecting rare taxa, which can be masked by the more abundant taxa during PCR or sequencing (Evans et al. 2016) or even filtered out during bioinformatics analyses due to too low read abundances (Elbrecht et al. 2018).

Relative biomass and read proportions for almost all mock-community species at the species level were not different for the 28S_D3 primer pair, suggesting that relative biomass can reliably be depicted by relative read numbers obtained by our nematode metabarcoding approach. In support, species that were missing or amplified at low proportions for the 28S_D3 and 18S_V4 markers in the dataset were those that also had a low overall biomass (Panagrolaimus thienemanni and Diploscapter coronatus, Table 1). These results imply that upscaling with more taxa and more complex communities might be possible. This study covered many of the extremely ubiquitous Rhabditiade, as well as other members of the nematode taxa as presented by van Megan et al. (2007), but more detailed studies using a wider phylogenetic range of nematode taxa are needed to confirm these findings across a wider species-range within Nematoda. If relative sequencing data actually represents relative biomass data, the information provided by metabarcoding might be implemented in studies investigating food-web structures (Sechi et al. 2018). This relationship, although by far less strong, has been observed for nematodes recently (Wilschut et al. 2019), but was limited to a 18S rDNA gene region. The here observed biomass-read number link was further confirmed by differences in read numbers between juvenile and adult specimen of the same species, such as shown before in fishes (Maruyama et al. 2014). As many ecosystem functions rely on biomass estimations, such as production, our results reveal that metabarcoding can be used to assess biomass distribution of previously hardly studied organisms, such as minute organisms living in soils and sediments. While we here show this potential for nematodes, we are convinced that this approach can be applied and provide biomass estimates also for other groups of organisms. Thorough calibration efforts are now needed to reliably link diversity and (relative) biomass of taxa in a sample.