Comparative metabarcoding study performed on the mock soil nematode communities, composed of DESS-preserved and freshly extracted nematodes

Nematode samples were collected mostly in years 2019 and 2020, from various environments and localities in Poland and laboratory cultures, which are specified in Table 1. The nematodes were extracted from the soil samples by decantation and sieving method followed by the centrifugal flotation method (van Bezooijen, 2006). In the case of fresh leaf and tomato root samples, the material was cut into smaller pieces and left on moist nematode filters in an extraction sieve until sufficient nematodes migrated into the water.

After the nematodes were extracted from soil, roots, or leaves and placed in water suspensions, they were heat-killed in tap water. Subsequently, they were observed and sorted using a stereoscopic microscope. Sorted and labelled nematodes were either fixed in a triethanolamine formalin (TAF) water solution or in a DESS mixture (Yoder et al., 2006), or placed directly in composed (mock) community test samples used for metabarcoding analyses. Temporary (DESS) and permanent (TAF) slides were made for further detailed morphological identification according to the method of Seinhorst (1959). After the photographic documentation was prepared using Leica DFC 500 camera on a Leica DM 5000B microscope (Leica Microsystems, Wetzlar, Germany) and after identification of the family, genus, or species levels, the nematodes from the temporary slides were rinsed of the preservative solution and placed in the freezer (the washing procedure as in Rybarczyk-Mydłowska et al., 2022). They constituted the base of reserve material for additional molecular analyses including the acquisition of new DNA fragments for reference sequence databases as well as for supplementary molecular identification. The permanent slide specimens are stored at the Laboratory of Nematology at the Museum and Institute of Zoology PAS.

After the preliminary identification of selected nematode taxa, two nematode sets (Set 1 and Set 2), composed of DESS-treated (at least for a year) and DESS-untreated (fresh) individuals, were prepared (Fig. 1, Table 1). These sets differed in taxa composition and number of individuals and were composed by transferring each nematode into 2 ml Eppendorf tube containing 150 µl of Milli-Q water or 180 µl of ATL buffer (Qiagen, Hilden, Germany), depending on the DNA extraction method used. Six replicates of each set, DESS-treated and untreated communities, were prepared: three were suspended in Milli-Q water (a–c) and three in the ATL buffer (d–f). Prior to the DNA isolation procedures, these DESS-treated mock communities (Set 1) were washed with phosphate buffered saline (PBS) solution and Milli-Q water as in the study by Rybarczyk-Mydłowska et al. (2022). The second set (Set 2) of the mock communities consisted of nematodes picked from water suspensions, directly after they were extracted from soil (roots or leaves), heat-treated, and sorted under a stereo microscope.

Figure 1

Simplified scheme of (mock) community sample preparation used for metabarcoding analyses. (a) Extracted heat-killed nematodes. (b) DESS preservative treatment or lack of treatment. (c) Transfer of the individual representatives of the segregated and identified nematode taxa into the test tubes. (d) Replicates containing milliQ water (marked in blue) or ATL buffer (yellow). Sets 1 and 2 differed in taxa composition (Table 1). Each replicate within the same set comprised of the same taxa and their amounts. (e) DNA isolation either with “CD” or “Q” method. As a result, 12 DNA samples were prepared and further analysed: 1aCD-1cCD, 1dQ-1fQ, 2aCD-2cCD, and 2dQ-2fQ.

Three replicates of each nematode set were used for DNA isolation using one of the two tested isolation kits (“CD” – ClearDetections, Wageningen, The Netherlands or “Q” – Qiagen, Hilden, Germany; Fig. 1). As a result, 12 samples were prepared (representing 2 sets of soil nematode taxa, isolated with 2 methods in 3 replicates each): 1aCD-1cCD, 1dQ-1fQ, 2aCD-2cCD, and 2dQ-2fQ, which were subsequently stored in a freezer (−20^oC) until the total DNA isolation process. Five nematode taxa (Acrobeles sp., Anaplectus sp., Aphelenchoides ritzemabosi, Meloidogyne hapla, and Rotylenchus uniformis) having the same number of individuals were included in both of the mock community sets (see the first five nematode taxa in Table 1) in order to better assess the influence of DESS on the tested metabarcoding approach. Four additional nematode taxa, with different abundances, were also included in both sets (see the nematode taxa from 6 to 9 in Table 1). As we aimed for the biggest nematode diversity possible, and the repeated sampling after 1 year allowed us to obtain the same taxa only partially, the remaining nematode taxa varied in both their identities and abundances across the two sets. However, these set up was sufficient to test the efficiency of the method itself and to evaluate if a specific nematode taxon, represented by a specific number of individuals, will be retrieved. The total number of individuals for Set 1 and Set 2 were similar (273 and 304, respectively). Finally, to compare the efficiency of the two chosen DNA isolation methods combined with the metabarcoding approach, the nematode sets were constituted in a way that included representatives of all trophic groups, varying in numbers, sizes, and cuticle structures. Our aim was not to test all possible DNA extraction kits, as it would further complicate the complex set up of the experiment. The selection of the two kits was subjective and chosen as one of them (“Q”) was successfully used in previous research (Ahmed et al., 2019, Waeyenberge et al., 2019) and the other one is specifically designed by the manufacturer for works with nematode suspensions (“CD”).

Isolation of DNA from the two sets of the artificially composed samples destined for metabarcoding studies (Table 1 and Fig. 1) was done using either the “Nematode DNA extraction & purification kit for nematode suspensions and multiple cysts” (Clear^®Detections, “CD”) or the “DNeasy^® Blood & Tissue Kit” (Qiagen, “Q”), following manufacturers’ instructions, with minor modifications to the original protocols. Previous results of our studies (data not published) showed that the initial direct resuspension of nematodes in the buffer ATL, instead of water, increases the efficiency of tissue lysis for the “Q” method. For this reason, the samples destined for the DNA isolation with this kit (1dQ-1fQ and 2dQ-2fQ) were transferred to ATL buffer prior to the extraction. In the “CD” method, instead of the recommended 200 μl, the total nematode lysate (300 μl) was transferred to the “pre-purification tube.” Additionally, due to such increased sample volume, two (instead of one) extraction columns were used per sample so that the final volume of the obtained sample was 200 μl. Similarly, the final eluate volume for the Q method was changed from 400 to 200 μl by adding half of the recommended amounts of the AE buffer in the last two elution steps. During the enzymatic lysis phase of both methods (200 rpm; 56°C or 65°C depending on the method), the incubation of the samples in the lysis buffers was prolonged to 24 h; however, it was interrupted several times for short periods of observations to examine the degree of nematode tissue degradation. These observations were conducted under a stereomicroscope, through plastic wall of a sample tube, without opening the lid. After the isolation process, DNA concentration was measured using the NanoDrop^TM instrument.

Single nematode individuals selected for additional molecular studies, for the purpose of supplementing the reference sequence database or confirming the morphological identifications, were transferred to separate 0.2 ml PCR tubes containing 25 μl of sterile water. An equal volume of the lysis buffer (Holterman et al., 2006) was added to each tube. The lysis conditions were 65°C for 3 h followed by a 5 min incubation at 100°C.

Two reference sequence databases were created: 18S rDNA and 28S rDNA. They included GenBank available nematode sequences as well as sequences acquired in this study (Accessed: 14 September 2022). The search term used for the 28S rDNA was: “((“Nematoda”[Organism] OR nematode[All Fields]) AND 28S[All Fields]) NOT 5.8S[All Fields] NOT uncultured[All Fields] NOT (“unidentified”[Organism] OR unidentified[All Fields]) OR ((“Nematoda”[Organism] OR nematode[All Fields]) AND lsu[All Fields]) NOT 5.8S[All Fields] NOT uncultured[All Fields] NOT (“unidentified”[Organism] OR unidentified[All Fields]) NOT chromosome[All Fields] NOT wgs[All Fields] NOT refseq[All Fields] AND (“300”[SLEN]: “3500”[SLEN]).” A similar search term was used for the 18S rDNA, with 28S and LSU (large subunit) replaced with 18S and SSU (small subunit), respectively. Sequences of the selected single nematodes and GenBank-downloaded sequences together formed the databases comprising of 38,955 18 rDNA and 12,733 28S rDNA sequences in total.

In order to improve the reference databases, rDNA was amplified from the four investigated Hoplolaimidae taxa (H. pseudorobustus, Helicotylenchus sp., Rotylenchus goodeyi and R. uniformis), Geocenamus longus, Mylonchulus brachyuris, and Tylocephalus auriculatus (Table 1). The detailed DNA amplification and sequencing conditions are described in previous studies focusing on characterizations of G. longus and R. goodeyi (Rybarczyk-Mydłowska et al., 2022; Singh et al., 2022). The DNA sequences obtained for the rest of the mentioned species were submitted under the GenBank accession numbers: OQ377005-OQ377012 (18S rDNA) and OQ377013-OQ377022 (28S rDNA).

From each of the isolated DNA sample replicates (1aCD-1cCD, 1dQ-1fQ, 2aCD-2cCD, 2dQ-2fQ), two amplicons were generated using primer combinations targeting 18S and 28S rDNA genomic regions (Table 2). The D2A/D2 primer combination was chosen in this study (and for the first time for metabarcoding applications) as it targets approximately 100 bp longer region in comparison to the frequently used D3FA/D3BR combination (Ahmed et al., 2019, Waeyenberge et al., 2019, Zapałowska et al., 2023). Each primer was ligated at its 5′ end either with the 5′-TCGTCGGCAGCGTCAGATGTGTATAAGAGACAG-3′ (forward) or the 5′-GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAG-3′ (reverse) Illumina adapter sequence. Each reaction mix contained 12.5 μl of 2× KAPA HiFi HotStart ReadyMix, 2.5 µl of each primer (5 µM), 2 μl of DNA, and 5.5 µl of Milli-Q water. The cycling conditions were 95°C for 3 min, followed by a specific number of cycles at 95°C for 30 s, 62°C for 30 s, 72°C for 30 s, and a final step of 72°C for 5 min. The specific number of cycles was determined via an optimization qPCR, which identified 30 cycles as the optimal for samples that were isolated with the “CD” kit and 25 cycles for samples that were isolated with the Q method. The optimization qPCR reaction contained 10 µL of Sso Advanced Universal SYBR Green Supermix (Bio-Rad), 1.5 µl of each primer (5 µM), and 2 μl of DNA. The qPCR conditions were 98°C for 3 min, (95°C for 15 s, 62°C for 30 s, and 72°C for 30 s) ×40 and 60–95°C for the melting curve.

After the amplicon PCR was performed, the acquired products were sent for NGS sequencing to Genomed company (Warsaw, Poland). The PCR products containing the adapter sequences were purified and indexed with the Nextera XT indexes. Reaction conditions were based on manufacturer’s recommendations. Sequencing was conducted on the Illumina MiSeq 2× 300 sequencer, using a paired-end approach. Automatic preliminary data analysis was carried out on the MiSeq using the MiSeq Reporter v2.6 software. The analysis consisted of two stages: automatic sample demultiplexing and generating FASTQ files containing raw reads. The obtained sequencing data are available from NCBI SRA database under the Bioproject PRJNA940653.

Bioinformatic analysis was carried out using USEARCH v11.0.667 (Edgar, 2010). The paired raw reads were first merged using the command -fastq_mergepairs. Pairs with more than 10 base mismatches within their region of overlap were discarded. The resulting contigs were filtered using the command -fastq_filter, to remove ones with “expected error” of more than 1 or contig lengths shorter than 250 bp. Before clustering, the filtered reads were reduced to unique sequences (dereplicated) using the command -fastx_uniques. The dereplicated sequences each have their number of copies appended to their description lines in the dereplicated FASTA file. Clustering was performed using the -unoise3 command, which is based on the UNOISE algorithm (Edgar, 2016). This method clusters reads into amplicon sequence variants-like groups referred to as Zero-radius Operational Taxonomic Units (ZOTUs), instead of the traditional OTUs. ZOTU-based clustering results in the recovery of a higher number of correct biological sequences than the OTU-based approach (Edgar, 2016). Taxonomy assignments of the ZOTUs were performed using blast matches against the two custom reference databases described above, one for the 18S rDNA and the other for the 28S rDNA (suppl. data S1).

All statistical analyses were conducted using R and RStudio software (R Core Team 2023; Posit team 2024, v. 2024.4.2.764). For the mock community diversity analyses, mainly ggplot2 (Wickham, 2016) and phyloseq (McMurdie and Holmes, 2013) R packages were used. The 18S and 28S rDNA ZOTUs tables, the taxonomy (including information on trophic preferences) tables and the samples table (supp. data S1) were transformed into phyloseq objects. The number of reads from each of the two sets and each of the two primer combination sequencings were normalized using the median sequencing depth ( T = median ∑ i = 1 n x i ; x i ′ = round T· x i ∑ j = 1 n x j ({T}=\hspace{0.25em}{\rm{median}}\hspace{0.25em}\left(\hspace{0.25em}{\sum }_{i=1}^{n}{x}_{i}\right){\rm{;}}\hspace{0.25em}{x}_{i}^{^{\prime} }={\rm{round}}\left({T\cdot }\frac{{x}_{i}}{\mathop{\sum }\limits_{j=1}^{n}{x}_{j}}\right) ) before the calculation of the abundances of the retrieved taxa, creating heatmaps of the most abundant sequences, the determination of the ZOTU richness, and the ordination plots. The latter were created using the non-metric multidimensional scaling and Bray-Curtis dissimilarity index.

Finally, tables of the summarized reads, recovered from each taxon and all taxa in total, comparisons of sample means, and the relative abundance table were created (suppl. data S2). To evaluate the possible influence of factors such as DESS treatment, different primer combination use, and the method of the DNA extraction on the obtained number of reads, the PERMANOVA tests were performed, using the lmPerm (Wheeler and Torchiano, 2025) package in R. Based on whether the reads were derived from a sample that was DESS-treated or untreated, an amplification of the 18S or 28S rDNA and the “CD” or the “Q” DNA extraction, eight possible treatment combinations were established, with three replicates each. Prior to running the PERMANOVA tests, the reads were transformed by dividing the number of reads designated to the taxon by the specific number of nematodes placed in samples (as specified in Table 1 and the suppl. data S2). The untransformed data were used only for the total read PERMANOVA analyses. The default 5,000 permutations were used for all tests. The test for indication of significant differences between the eight treatments, as well the combinations of four and two treatments, was performed for the taxa that were present in both the sets (see the first 9 taxa in Table 1) and for the total amount of reads obtained for all taxa together. This test aimed to check for joint and separate influence of the three applied factors (DESS preservative treatment, primers, and isolation) on the number of the obtained reads (results are given in Table 3A). Five additional tests were performed for each taxon separately that was present either in the Set 1 (DESS-treated) or Set 2 (fresh) and for the total amount of reads obtained for all taxa together (results are given in Table 3B and 3C, respectively). Significant differences were tested considering either all four treatments, that could correspond to the effects of joint influence of primers and isolation, or two of them in all possible combinations that would correspond with the influence of isolation within the 18S rDNA or the 28S rDNA and the influence of primers where the “CD” or the “Q” method was applied. In some cases, M. hapla, A. ritzemabosi, and M. xenoplax were excluded from the analyses due to the absence of reads for these taxa.

Bioinformatics procedures and the R code are available from the corresponding author upon reasonable request.

Twenty-four nematode taxa were chosen for this study (Table 1). Morphological observations, and in some cases, the molecular analyses as well, resulted in the identification of 13 taxa to the species level, 8 to the genus, and the other 3 to the family level.

After the mock samples were composed (Fig. 1, Table 1), the DNA extraction using two different DNA isolation methods (“CD” and “Q”) was performed. Initial observations of tissue breakdown during the lysis step showed that the most difficult individuals to disintegrate were the very large nematodes such as the representatives of Dorylaimida as well as the nematodes with thick cuticular structures (Hoplolaimidae and Criconematidae). Based on the microscopic observations, a slightly stronger digestive effect was observed in the case of the “CD” isolation kit. However, the observations showed that none of the kits guaranteed 100% lysis efficiency of the prepared mixtures of about 300 nematodes after 24 h of incubation in the two tested lysis buffers. The DNA isolations resulted in substantially higher DNA concentrations for the “CD” method (suppl. data S3). Additionally, the results of the performed qPCRs showed that the amplification of rDNA fragments from the samples isolated with the “CD” kit was delayed by about 5–8 cycles in comparison to the samples isolated with the “Q” method, where the Ct values were relatively lower (data not shown).

Approximately 80,000–100,000 reads were obtained for the 18S rDNA and 40,000–80,000 for the 28S rDNA samples (suppl. data S2; suppl. Fig. S1). The obtained relative abundance for both the 18S and 28S rDNA regions did not exactly reflect the number of individuals put in the samples. The obtained reads for the particular taxa in the tested replicates constituted different fractions of the total amounts (Fig. 2). All the taxa were represented by reads from multiple ZOTUs (Fig. 2; suppl. data S1). Regardless of whether the nematodes were treated with DESS prior to DNA isolation or not, reads were recovered from almost all taxa tested (Fig. 2; suppl. data S1 and S2). Nevertheless, instead of the assumed Qudsianematidae (other than Discolaimus), the current representatives of the Dorylaimidae family (Labronema – 28S rDNA, earlier also in Qudsianematidae; Mesodorylaimus – 18S rDNA; Prodorylaimus – 18S rDNA) were detected (Figs 2 and 3, suppl. data S1). Therefore, they were treated as misidentification, and these reads were counted in as the taxon initially defined as Qudsianematidae for the read summary table (suppl. data S2). Moreover, the reads were not recovered from every genomic region to the same degree. For the five nematode taxa represented in both sets in similar amounts (Table 1), the recovery was usually higher in case of 18S rDNA from fresh nematode samples (Fig. 2; suppl. data S2). The lowest amounts of reads or no detection at all was obtained in case of the 28S rDNA from M. hapla or A. ritzemabosi in the DESS-treated samples (suppl. data S2). M. xenoplax was not detected from any of the samples where only four individuals of this species were placed. In fresh samples where 20 M. xenoplax were included, this species was detected in most of the samples (suppl. data S2). The strongest read recovery in comparison to other taxa was obtained from large-sized nematodes of the order Dorylaimida (e.g. identified as Aporcelaimellus, Labronema, Longidorus, Mesodorylaimus, and Sectonema), even with only few of them present in the sample, especially in the DESS-treated samples (Fig. 3; suppl. data S2). However, the number of reads obtained for some taxa showed an association with the method of DNA isolation. For example, the recovery from thick-sheathed nematodes such as Mesocriconema or Rotylenchus was stronger in the samples isolated with the “CD” kit (suppl. data S2) than in the ones isolated with the Q method. The “CD” method recovered more than expected reads of Longidorus, Mylonchulus, or Paratylenchus, while Acrobeles, Alaimidae, Anaplectus, Aporcelaimidae, Coomansus, Discolaimus, Geocenamus, or Paratrichodorus were recovered better with the “Q” method. At the same time, the mean overall number of reads for the samples isolated with the “CD” and “Q” methods remained close; however, in most cases, slightly higher for the “Q” method: ∼97,000 vs ∼95,000 (18S rDNA, DESS-treated), ∼57,000 vs 70,000 (28S rDNA, DESS-treated), 85,000 vs 90,000 (18S rDNA, fresh), and 53,000 vs 58,000 (28S rDNA, fresh). False positive reads were recorded for four taxa: Aporcelaimidae, Discolaimus major, Paratrichodorus pachydermus, and Tylenchorhynchus sp. (suppl. data S2).

Figure 2

Taxa abundance and composition in two mock nematode community sets. (a–c) Set 1, DESS-pretreated nematodes. (d–f) Set 2, DESS-untreated, freshly extracted nematodes. (a) Exact numbers of nematode taxa included in the samples 1aCD-1cCD and 1dQ-1fQ. (d) Exact numbers of nematode taxa, included into the samples 2aCD-2cCD and 2dQ-2fQ. Taxa initially occurring in the two nematode community sets in the same amounts (Table 1) are marked with an asterisk. Taxa occurring in the two nematode community sets in the different amounts are marked with two asterisks. (b and e) Read numbers and retrieved genera after 18S rDNA sequencing. (c and f) Read numbers and retrieved genera after 28S rDNA sequencing. The obtained false-positive reads were not included into the visualization and the corresponding data are available in the suppl. data S1. “CD” and “Q” samples were isolated with the Clear Detection or Qiagen DNA extraction kits, respectively.

Figure 3

Heatmaps of the most abundant sequences obtained after metabarcoding sequencing of sample replicates (a–f), representing two mock nematode communities (Sets 1 and 2), isolated with two DNA extraction kits (“CD” and “Q”). (a) Most abundant nematode species sequences after 18S rDNA sequencing. (b) Most abundant nematode species sequences after 28S rDNA sequencing.

The total 18S and 28S rDNA sequence (ZOTUs) richness among the 12 analysed samples, representing the two analysed sets of species (and preservation), was estimated with Chao1 and Shannon indices (Fig. 4). The latter balances richness with evenness. The obtained alpha diversity measures were similar for the samples representing the two sets. Mean Chao1 values for 18S rDNA were: 189.33 (Set 1 isolated with “CD”), 293.48 (Set 1 isolated with “Q”), 238.13 (Set 2 isolated with “CD”), and 299.93 (Set 2 isolated with “Q”); and for 28S rDNA: 85.57, 142.69, 67.67, and 200.25, respectively. Correspondingly, mean Shannon values for 18S rDNA were: 3.53, 4.15, 4.04, and 4.37; and for 28S rDNA: 2.15, 2.78, 2.23, and 2.46. However, the sequence diversity (ZOTUs richness) was generally lower in case of the 28S rDNA data and was strongly influenced by the extraction method, in both the 18S and the 28S rDNA. The samples isolated with the Q method exhibited higher diversity and more uniform results among the three replicates. This was also confirmed with the ordination analysis of the particular samples (suppl. Fig. S1).

Figure 4

ZOTUs richness estimated with Chao1 and Shannon indices obtained after 18S (a) and 28S rDNA (b) metabarcoding sequencing of sample replicates (a–f), representing two mock nematode communities (Sets 1 and 2), isolated with two DNA extraction kits (“CD” and “Q”).

The reads obtained for each taxon separately and for all taxa together (the total amount) were used for PERMANOVA tests (suppl. data S2; Table 3). Combinations of up to eight treatments, with three replicates each, were tested to examine joint or separate effects of: DESS preservative treatment, primers used, and the isolation method applied. Consequently, the obtained statistically significant results of the PERMANOVA tests indicate significant differences between the tested treatments, suggesting possible influences of the tested factors on the retrieved amounts of reads.

The initial PERMANOVA tests performed for each of the nine taxa, that were included in both of the prepared sets, indicated significant differences among the eight and four treatment groups (Table 3A; suppl. data S2). This suggests the joint influence of preservation and isolation method, preservation and primers, and all three of these factors on the read abundance. However, for the main effect of the DESS preservation, some of the tests did not indicate significant differences for all taxa except the representatives of Hoplolaimdae (Helicotylenchus sp. and R. uniformis).

PERMANOVA tests performed on the reads obtained from the taxa derived from the DESS-treated Set 1, indicated a significant difference in at least one among the five treatment groups for all the taxa except for T. arenicola (Table 3B). Significant results that might indicate the influence of the isolation method were observed for all taxa except: Aporcelaimidae, M. hapla, T. arenicola, T. auriculatus, and Tylenchorhynchus, in case of the amplification with the 18S rDNA primers, and Acrobeles sp., C. parvus, Helicotylenchus sp., Longidorus sp., P. projectus, Qudsianematidae, and R. uniformis, in case of the amplification with the 28S rDNA primers. The results indicating the type of primer combination used might have had an influence on the reads obtained for A. ritzemabosi, Helicotylenchus sp., Longidorus sp., M. hapla, Plectus sp., R. uniformis, and T. auriculatus when DNA was extracted using the “CD” method, and for all taxa except: Aporcelaimidae, Acrobeles, C. parvus, Longidorus sp., and T. arenicola, when “Q” extraction method was used.

PERMANOVA tests performed on the reads obtained from the taxa of Set 2 mock community samples indicated a significant difference in at least one among the five treatment groups for all of the taxa. For the 18S marker treatments, no differences were found between the isolation method region only as far as M. hapla, Longidorus sp., and Qudsianematidae are concerned. Likewise, for the 28S marker, no differences were found between the isolation methods for A. ritzemabosi, Helicotylenchus, H. pseudorobustus, M. hapla, and Qudsianematidae. When “CD” method was used for the extraction, the marker used had no influence on Alaimidae, Anaplectus sp., G. longus, Longidorus sp., Pratylenchus sp., and Qudsianematidae. For mock community samples extracted using the “Q” method, Acrobeles, D. major, G. longus, Longidorus, M. xenoplax, and P. pachydermus were the taxa for which the marker region used for the amplification had no effect at all.