All four field sites were located around the towns of Plant City and Wimauma in Hillsborough County, Florida. The region around Plant City has a long history of strawberry production going back to the late 1800s and, to this day, accounts for most of Florida’s commercial winter strawberry production. Three of the fields were commercial strawberry farms known to produce berries for at least 25 years (Fields C1, C2 and O1) and the fourth field belonged to the University of Florida’s Gulf Coast Research and Education Center where strawberries were first planted in 2016 (Field R1). The field selection was based on a previous nematode survey conducted between 2017 and 2021 to identify the extent of recent infestation with Meloidogyne hapla, a nematode likely introduced via out-of-state strawberry transplants (Desaeger 2019; Oliveira 2022). Two of the commercial fields (Fields C1 and C2), as well as the UF research field (Field R1), were conventionally managed, while the other field was organically managed (Field O1). Details of management practices regarding nematicides and cover crops varied amongst fields and are given in Suppl. material 1.

Soil samples for this study were collected during June 2021, about 2 months after terminating the strawberry crop. In fields R1, C2 and O1 a cover crop, sunny hemp (Crotalaria juncea), was growing at the time of sampling, while field C1 had a pepper crop growing instead (Suppl. material 1). Within each field, a M. hapla nematode hot spot was previously identified, based on above-ground symptoms and root galling. Each nematode hot spot was divided into two sites and each site was marked with one GPS point. Within 3.6 m from the GPS point, two replicate soil samples were randomly collected with a gas-powered auger at four separate soil depths: 0 – 25 cm, 25 – 50 cm, 50 – 75 cm and 75 – 100 cm (one core per depth). The auger was equipped with reinforced regular 2–3/4” (6.9 cm) quick connector soil augers (AMS, Inc. American Falls, ID), collecting approximately 1 litre of soil per sample. Soil samples were placed into plastic bags and stored at 4 °C for a maximum of 48 h before nematode extraction. At the UF Field R1, we also measured the temperature at four different depths using a Hobo ProV2 onset data logger (Suppl. material 2) (Onset Computer Corporation, Bourne, MA, USA) confirming that sampling during the June month was in line with the hottest months of the year (average 29.41 °C), comparable only to July. This could not be done in the commercial fields due to constant machinery activities. At the time of sampling, upper soil layers were, on average, lower in temperature than the deeper soil layers (0 – 25 cm: 28.90 ± 0.55 °C, 25 – 50 cm: 29.27 ± 0.75 °C, 50 – 75 cm: 29.58 ± 1.04 °C, 75 – 100 cm: 29.86 ± 1.94 °C). An extra subsample of 200 cm³ from each of the soil samples was sent for soil analysis at the Urban Soil and Water Quality Lab of the University of Florida, Wimauma, FL. Results demonstrated that soil from the different depths collected at the farms were primarily Myakka sand with soil texture ranging from 95 –98% sand, 2 – 3.9% clay and 0 – 1.3% silt. Refer to Suppl. material 3 for further details on the soil type and particle size of the samples.

Nematodes were extracted from a 200 cm³ subsample of soil using the centrifugal flotation technique (Jenkins 1964). After collecting over a 25 µm sieve, nematodes were transferred in water into plastic scintillation vials and stored at 4 °C for a maximum of two weeks prior to counting on a ZEISS Axio Vert.A1 inverted microscope (Carl Zeiss, Gottingen, Germany). A calibrated counting dish was used to count and examine 50% of a nematode suspension and, afterwards, multiplied by two to adjust to the final nematode counts (Suppl. material 4). Plant-parasitic nematodes were identified to genus and selected individuals (see below) to the species level by using the pictorial key proposed by Mai and Mullin (1996). Morphometrics of representative specimens of Nanidorus minor were taken as described in Oliveira et al. (2023a). Polymerase Chain Reaction – Restriction Fragment Length Polymorphism (PCR – RFLP) was used to confirm the identify of representative specimens of Meloidogyne hapla as described in Oliveira et al. (2023b). As endemic to Florida, all specimens resembling Belonolaimus longicaudatus were assumed to be that species. As per standard practice of morphological evaluation of nematode communities in agricultural settings, all free-living nematodes were identified only at the trophic group level following Yeates et al. (1993). Nematodes were carefully returned to the plastic vials and kept at 4 °C before DNA extraction to ensure nematode viability and DNA integrity (Christoforou et al. 2017).

To remove excessive soil and organic matter debris, counted nematodes were passed through stacked 425 µm over 25 µm sieves, then centrifuged, reduced to 0.5 ml and transferred to Qiagen® PowerSoil Kit (Qiagen Inc., Toronto, Canada) bead-beating tubes for DNA extraction following manufacture’s protocols (Mauger et al. 2020). Extracted DNA concentrations were quantified using a Qubit Fluorometer (ThermoFisher Scientific Inc.).

Whole-genome sequencing libraries were prepared using the Kapa BioSystems HyperPlus Kit (KR1145 –v.3.16), optimised for low-input DNA, multiplexed and enriched for mitochondrial DNA with capture/hybridisation probes specifically designed to target nematode mitochondrial sequences (Sevigny et al. 2021). DNA sequencing was completed on a NovaSeq 6000 instrument to produce 250-bp paired-end reads. Raw sequencing reads were demultiplexed using the Illumina bcl2fastq conversion Software v.1.8.4. FastQC v.0.11.5 was used to remove low-quality reads and sequencing adapters (Andrews 2010; Gendron et al. 2023).

To recover the maximum number of nematode contigs, they were assembled using three different protocols: 1) de novo assembly followed by COX1 gene annotation, 2) read mapping to the reference database of COX1 gene followed by de novo assembly and 3) read mapping to the reference database of all mitochondrial genes followed by de novo assembly and gene annotation. Regardless of the assembly protocol, all sequences were first assembled using the metaSPAdes pipeline (Nurk et al. 2017; Prjibelski et al. 2020). Assembled contigs were tentatively identified as mitochondrial in origin using BLAST against the NCBI nucleotide database at ≥ 75% ID match and filtered of non-mitochondrial sequences. Mitochondrial contigs, generated under the 1^st assembly protocol, were scanned for open reading frames using the ‘orfipy’ python script tool using the invertebrate mitochondrial codon table as the reference (https://github.com/urmi-21/orfipy; Singh and Wurtele 2021). Open reading frames were then BLASTed against the nematode mitochondrial database (Nema-mtDB; see Gendron et al. (2023) for database composition and overview) and the best BLAST hits (based on the highest bit scores and lowest e-values) were used for initial taxonomy and gene identification. Under the 2^nd assembly protocol, sequencing reads were first mapped (at 100% match over a short raw read sequence) to the Nema-mtDB COX1 reference sequences using the Burrows-Wheeler Aligner (BWA; Li and Durbin (2010)) and then assembled through de novo assembly as described above. Finally, the same process applied to the 3^rd assembly protocol, except that reads were first mapped to all genes in the Nema-mtDB mitochondrial reference sequences.

Due to the capacity of nematode communities for a high degree of gene flow and migration potential, the likelihood of population level diversity is high (Derycke et al. 2013). In order to differentiate between true and erroneous assemblies, we examined all recovered haplotypes amongst identified taxa. First, all filtered contigs were tagged with sample origin, accession ID and % ID at the species level, based on matches to the Nema-mtDB (Gendron et al. 2023). As evaluation of haplotypes, based on % similarity, does not effectively identify errors, contigs and their best matched reference sequences were aligned using the ‘mafft’ aligner under the settings ‘—auto’ and ‘–adjustdirectionaccurately’. Alignments were used to develop Maximum Likelihood trees under the general time-reversible model in ‘FastTree’ with bootstrap support values based on 1000 replications (Price et al. 2010). Potential haplotypes and errors were initially examined using trees at the suborder level (i.e. Tylenchina, Rhabditina and Spirurina) and promising haplotype clades further followed using specific trees at the genus level (Meloidogyne, Bursaphelenchus, Diploscapter, Acrobeloides, Mesocriconema, Litoditis and Panagrellus). Contigs with suspect validity (i.e. poor alignment and placement on trees) were considered as errors and removed from further analyses. While this approach currently provides for the most conservative diversity estimates, as the mtMG database expands, future work should confirm the validity of this haplotype analysis.

The assembled and filtered contigs from each of the three contig assembly protocols were further filtered for nematode ID assignment using the highest scoring hit passing previously established cutoffs (> 100 bp length, > 91% ID similarity to the Nema-mtDB reference) (Gendron et al. 2023). A further cut-off of > 31% query coverage was set for contigs to improve retention of only reliable matches. Identified contigs passing BLAST cutoffs were equated with positive taxon detections. Samples generating multiple contigs assigning to the same taxon with identical matches were characterised as single hits; hence, all downstream analyses were based on the presence/absence data. Detected taxa were used to construct a presence/absence community table and classified to trophic groups following Yeates et al. (1993). The annotated COX1 approach far outclassed the two read-mapping approaches in recovery of nematode diversity and so we used it for further statistical analysis and proper comparison to the morphologically derived communities. However, communities based on the two read mapping protocols (read mapping to COX1 and read-mapping to all mitochondrial genes) have been included in Figures to emphasise the power of annotation in identifying taxa.

Assembled, identified and validated contigs from the annotated COX1 approach communities were analysed for richness and compositional variation and compared to the standard practice of morphologically identified communities. Annotated COX1 richness was calculated for the entire nematode community at the trophic group level as well as at the species level. In addition, to compare the mtMG method to morphology for PPN, richness at the genus/species level for this trophic group only was also calculated. All calculations were done using base R v.4.0.3 (Suppl. materials 4, 20). For the morphology-based communities, richness was calculated at the trophic group level for the whole community and genus/species level for plant parasites. Due to uncertainty of whether multiple contigs with the same classification reflected relative abundance or an artefact of mitochondrial enrichment and assembly, compositional differences were evaluated, based on presence/absence by calculating Jaccard dissimilarity distances using the R function ‘vegdist’ in the R package ‘vegan’ v.2.5-7. We assembled dissimilarity matrices for the entire community at the trophic level and plant parasitic community at the species level across both methods. In addition, we also developed a matrix for the entire annotated COX1 nematode community at the species level.

In order to evaluate potential relationships between and amongst both the morphological and mtMG community data, we tested them with a pairwise Pearson correlation using the ‘cor’ function in base R 4.3.2. We specifically examined morphology-based counts and richness, recovered contig numbers, number of identified mitochondrial contigs, DNA extraction concentrations and annotated COX1 richness (Suppl. material 4; Stanish et al. (2013)). To determine if significant Pearson correlations were potentially indicative of direct linear relationships, variables (particularly indicators of biomass, such as concentrations of DNA, numbers of assembled contigs and richness) were further investigated using post-correlation linear models. Nematode community richness for both the morphological and mtMG data was tested for method, field and sampling depth as the explanatory factors using paired T-tests, from the base R function ‘t.test’ and Generalised Linear Models using the ‘glm’ function in the R package ‘stats’. Post-hoc GLMs were then run on richness for both methods separately in order to better identify trends present within each method. Community compositional diversity using Jaccard distances was tested with PERMANOVA models using the ‘adonis’ function in the R package ‘vegan’ v.2.5-7 with method, field and sampling depth as the explanatory factors. Post-hoc PERMANOVAs were used to identify specific fields and depths that significantly explained community composition blocked by either Field or Depth. To compare the effects of field and sampling depth on richness and composition between morphology and mtMG-based data, we ran analyses for the whole community at the trophic level as well as species level (or genus where species ID was not possible, i.e. Tylenchorynchus spp.) but only for plant parasites. We also tested the effect of field and sampling depth on richness and composition of the entire community at the species level, but only using the mtMG data from the annotated COX1 approach.

Overall, there was very little evidence of relationships between morphology-based counts and metrics of DNA biomass (e.g. number of total contigs and DNA concentrations) and when there was, correlations were generally weak (e.g. bacterivore counts and total number of contigs) (Suppl. material 5). Significant correlations were observed between morphological counts for plant parasites and the BLAST measures for query coverage and contig match lengths, but these relationships were also relatively weak. The number of recovered contigs (both total and mitochondrial) significantly correlated with DNA extraction concentrations (Suppl. material 6). Out of the richness metrics, only trophic level richness (annotated COX1 and morphology) showed a significant correlation with DNA extraction concentrations. Post-correlations GLMs confirmed these relationships with relatively low R² values (Suppl. material 6, R² 0.069–0.134).

Out of 301 contigs that passed initial classification, as expected, most belonged to Rhabditida. Out of these contigs, 35 were removed from the detection table based on poor alignments to their classified reference sequences indicative of misidentification. The removed contigs predominately belonged to Tylenchina (29 contigs). Examination of the Tylenchina contigs against their reference sequences generally showed strong clustering by genus, except for Acrobeloides and Panagrellus redivivus contigs (along with one of its reference sequences) which clustered together. Panagrellus redivivus demonstrated well-clustered clades, but they also showed some of the largest branch distances between the contigs and their reference sequences of any of the detected taxa (Suppl. material 7). These clades were retained, in contrast to many P. redivivus contigs outside this cluster that scattered throughout Tylenchina which were removed (12 contigs with average 92% similarity to the reference sequence). One clade of M. javanica (9 contigs, average 99% ID similarity) demonstrated poor alignment with both its reference sequences and the other Meloidogyne species and, thus, was removed. Lastly, based on misplacement on the tree or low % ID similarity, the following Tylenchina contigs were removed: Steinernema longicaudum (1 contig, 91% similarity), Steinernema sp. (LT963444, 1 contig, 92%), Meloidogyne hapla (3 contigs, average 98%), Meloidogyne luci (1 contig, 99%), Meloidogyne incognita (1 contig, 92%) and Bursaphelenchus mucronatus (3 contigs, average 93%).

The remaining filtered contigs (6 contigs) belonged to Rhabditina. These taxa included: Pristionchus pacificus (1 contig, 92%), Litoditis aff. (2 contigs, average 92%) and Allodiplogaster seani (3 contigs, average 93%). Two unique clades of Litoditis aff. clustered separately from each other as well as their reference sequences. One clade (19 contigs, average 94% ID) aligned with the Plectus outgroup and the other (15 contigs, average 94% ID) with the reference sequences of Rhabditis blumi (Suppl. material 8). These clades were retained largely because the lack of clear placement within Rhabditina reflects most likely the absence of their expected reference sequences (and their high sequence similarity 92–95%). Another feature of note was a clade of long branch lengths within one of the Diploscapter coronatus clades. This clade (4 contigs, average 94% ID) was also retained because of its clear clustering with other D. coronatus contigs belonging to the expected species (92–96%).

To reduce potentially erroneous estimates of population level diversity, we examined our sequences for the presence of haplotypes. Based on well supported nodes amongst clades (> 0.8 bootstrapped node value), the contig trees highlighted a potential for haplotypes amongst the genera of Meloidogyne, Bursaphelenchus, Diploscapter, Acrobeloides, Mesocriconema, Litoditis and Panagrellus. All potential haplotypes were detected in Tylenchina (Suppl. material 7) and Rhabditina (Suppl. material 8) suborders, but no haplotypes were detected within recovered Spirurina taxa (Suppl. material 9). At the suborder level, Acrobeloides clustered well as part of clade with Panagrellus redivivus (Suppl. material 7). However, at the genus level, Acrobeloides contigs did not show consistent clustering of species into distinct clades and, instead, the five recovered species showed high similarity (average 95%) and were intermixed with each other (Suppl. material 10). Bursaphelenchus contigs (Suppl. material 11) fell well into a monophyletic clade as predicted by assigned taxonomy to B. arthuri with at least two distinct haplotypes (average 93%) distinguished by high node supports. Diploscapter was only represented by D. coronatus (average 94% ID) and showed evidence of two highly-related clades (1.000 node support) when treed with sequences at both the suborder and genus level (Suppl. material 12). Litoditis (Suppl. material 13) was represented by contigs clustering into two distinct clades. However, one clade was much more distant from the reference sequences (both with an average 94% similarity) than the other and the distance between contig clades and the references were most pronounced when examined within the context of the rest of Rhabditina (Suppl. material 8). There is a good chance that these clades have been assigned the wrong taxonomy, but these still likely belong to Nematoda and do not produce better alignments to known taxa. Meloidogyne showed clear definitions of clades at both the genus and species level (Suppl. material 14). The genus was represented by M. hapla (avg. 99% ID), M. javanica (avg. 99% ID), M. luci (avg. 99% ID) and M. incognita (avg. 99% ID). Two potential haplotypes of M. hapla, M. javanica and M. incognita were detected, based on clade node supports. Mesocriconema was represented by two species that showed clear phylogenetic separation (Suppl. material 15). M. ornatum (avg. 100% ID) showed evidence of two potential haplotypes, of which only one was represented by a single contig and a reference sequence (91% ID). At the suborder level, Panagrellus redivivus contigs (92 – 95% ID) and their reference sequences clustered within the Acrobeloides clade (Suppl. material 7). At the genus level, contigs clustered into at least two well-supported clades (Suppl. material 16).

With the addition of tree-based contig filtering, the three mtMG approaches recovered different numbers of taxa with read-mapping methods underestimating diversity and the annotated COX1 approach producing the closest replication of the morphological data (Fig. 1). Based on this annotated COX1 approach, the detection of bacterial-feeding (BF) and plant parasitic (PP) taxa was the most consistent with morphology, but only at the trophic level (Fig. 1). Based on morphology, the presence of BF nematodes was observed in all samples and depths, but their identity remained unknown. In addition to the presence of bacterivores, the annotated COX1 approach tentatively identified 14 BF species in nine genera (Fig. 1), including the frequent detection of Diploscapter, Halicephalobus, Litoditis and Panagrellus and less frequently the detection of Acrobeloides, Panogrolaimus, Pellioditis, Rhabditidoides and Allodiplogaster. Although detected at all depths, BFs were found mostly in the upper soil layers or at the deepest (75 – 100 cm) depths, especially in Fields C1 and C2 where they were also the most diverse.

Figure 1.

Nematode Presence/Absence Detection Method Comparison Heatmap. Binary heatmap of detected taxa and trophic groups, based on morphology, cox1 read-mapping, all mitochondrial genes read-mapping and annotated COX1 (black = detected, white = not detected). Trophic groups detected in fields included Bacterial Feeders (BF; blue), Plant Parasites (PP; Green), Fungal Feeders (FF; brown), Animal parasites (AP; red), Entomopathogenic Nematodes (EPN; yellow), Predators (PR; black) and Omnivores (OM; grey). Detections are split by field of origin, site sampled, site replicates and depth. Annotated COX1 detections were the most similar to the morphology data of any mtMG method and so only the detected species for the Annotated COX1 data are presented.

Despite the lower detection consistency across all replicates and sites than the morphology approach, the annotated COX1 approach generally agreed with morphology in detection of Meloidogyne species across all fields. In addition, it also provided detailed insights into specific species, confirming that M. hapla was present in Fields C1, C2 and O1, but also detecting the co-occurrence of other Meloidogyne species including M. javanica, M. incognita and M. luci in Field O1. However, the less frequently detected plant-parasitic taxa showed more discrepancies between the mtMG and morphology approaches. Either the detection of these taxa was limited to only microscopy, but not detected with the mtMG approach (Nanidorus sp., Tylenchorynchus sp. and Belonolaimus longicaudatus) (Fig. 1, Suppl. material 17) or they were detected with the mtMG approach, but not microscopy (Mesocriconema sp. and Xiphinema brevicollum). Interestingly, Mesocriconema species were detected almost exclusively in Field R1 (with a single detection in Field C2) and Xiphinema in a single sample in Field C2.

Based on morphology, fungal-feeding nematodes (FF) were observed in almost every sample and depth, albeit of unknown identity. The annotated COX1 approach primarily detected FF taxa in Field C2 with sparse distributions in other fields. The most commonly detected FF taxa in Field C2 consisted of Bursaphelenchus species in contrast to Aphelenchus sp. detected only in two samples in Field R1.

The remaining trophic groups including animal parasitic (AP), entomopathogenic (EPN), omnivory (OM) and predatory (PR) nematodes showed no overlap between morphology and mtMG approaches. The largest difference between morphology- and mtMG-based datasets was the detection of AP with no taxa detected morphologically, but 10 molecularly. Although the two read-mapping approaches detected APs in very few samples, when detected, the results agreed across all assembly approaches with the highest number of taxa from the highest number of samples being detected in Field C2. Although the presence of APs was detected at all depths, 0 –25 cm generated the most positive hits, especially in Field C2. Almost all detected taxa represented terrestrial mammal parasites, albeit the % ID matches to their reference sequences in the Nema-mtDB were relatively low (average 93% ID). EPN taxa (Steinernema diaprepesi), similar to AP taxa, were only detected with the mtMG approach, but only in Field R1. In contrast, both PRs and OMs were only observed with morphology.

Richness of trophic groups, based on morphological counts, significantly varied by both field and sampling depth, with the highest richness present in the upper soil depths and with the highest average across depths in Field O1 (Field: P = 0.021, Depth: P = 0.033; Fig. 2A, Suppl. material 18). Plant parasite richness also significantly varied by field and depth (Field: P = 0.034, Depth: P = 0.018; Fig. 2B, Suppl. material 18), with the highest average plant parasite richness in Field O1. The annotated COX1 mtMG data also showed significant variation amongst fields and depths for trophic richness (Field: not significant, Depth: P = 0.008), plant parasitic richness (Field: not significant, Depth: not significant) and total community richness (Field: P = P < 0.001, Depth: not significant), but the exact patterns were significantly different from the morphology-based data; Figs 2A, B, 5, Suppl. material 18). In terms of depths, morphology-based trophic richness was significantly explained (P = 0.033) across all depths, while the annotated COX1-based trophic richness was significantly explained (P = 0.008) by only a single (75 – 100 cm) depth.

Figure 2.

Trophic and Plant Parasite Community Richness. A comparison between mtMG Annotated COX1-based data (A; grey) and morphology-based data (M; light blue) for A. Trophic community richness and B. Plant Parasite community richness using box and whisker plots. Richness metrics presented are split by field of origin and depth sampled with sites and replicates for each field-depth pooled together.

Overall, the mtMG data described simpler and more similar communities than the morphology data. Trophic community composition, independent of method, reflected the richness results, with both field and sampling depth significantly explaining variation in Jaccard distances (Field: R2 = 0.04 and P = 0.004, Depth: R2 = 0.04 and P = 0.002, Field*Depth Interaction: R2 = 0.06 and P = 0.037; Fig. 3A, B, Suppl. materials 18, 19), but, in general, there was high variation amongst both fields and depths. Significant differences between morphology and mtMG approaches were also observed for community composition (R2 = 0.35 and P = 0.001; Suppl. material 18) at the trophic group level. Post-hoc analysis showed that significance of field was driven by fields C1, O1 and C2 (P = 0.018, 0.003 and 0.004, respectively) for the morphology trophic data, while by fields C2 (P = 0.003) and O1 (P = 0.003) for the annotated COX1 trophic data.

Figure 3.

Trophic Community Composition. PCoA plots of Jaccard-based community compositional dissimilarity for trophic-level communities, based on A. Morphology data and B. Annotated COX1 data. Dissimilarity explained by axes are presented in the axes labels and colours correspond to field (R1 = black, C1 = red, O1 = green, C2 = blue) with increasing depth represented by lighter shades of each field colour.

The plant parasite composition, independent of method, was significantly explained by field (R2 = 0.20 and P = 0.001); however, depth was only significant for morphology (Field: R2 = 0.12 and P = 0.001, Depth: R2 = 0.11 and P = 0.031; Fig. 4A, B, Suppl. materials 18, 19) and not for mtMG (Field: R2 = 0.50 and p = 0.001, Depth: R2 = 0.05 and p = 0.860; Fig. 4, Suppl. materials 18, 19). The significance of field was driven by all fields for both the morphology and mtMG approaches (P = 0.001 – 0.005), but both approaches showed that Field R1 accounted for the majority of explained variation (mtMG R2 = 0.29, morphology R2 = 0.15) and showed little variation in composition compared to the other fields. In terms of depths, morphology-based PP richness was only explained by the single 25–50 cm depth (R2 = 0.05 and P = 0.006), with no effect of depth detected for the mtMG data.

Figure 4.

Plant Parasite Community Composition. PCoA plots of Jaccard-based Plant Parasite community compositional dissimilarity, based on A. Morphology data and B. Annotated COX1 data. Dissimilarity explained by axes are presented in the axes labels and colours correspond to field (R1 = black, C1 = red, O1 = green, C2 = blue) with increasing depth represented by lighter shades of each field colour.

Figure 5.

Annotated COX1 Community Richness. Species richness box and whisker plots for each field and depth based on pooled replicate site samples (R1 = black, C1 = red, O1 = green, C2 = blue).

When examining the species level for entire community composition, based on annotated COX1 mtMG data, only field was a significant factor in explaining the variation (R2 = 0.16 and p = 0.001; Fig. 6, Suppl. materials 18, 19). This significance was driven by all fields (Field R1: R2 = 0.06 and P = 0.001, Field C1: R2 = 0.03 and P = 0.023, Field O1: R2 = 0.07 and P = 0.001, Field C2: R2 = 0.07 and P = 0.001). Generally, Fields R1 and C1 were characterised by high variation that overlapped much more than the more tightly clustered variation of Fields C2 and O1. Furthermore, while not significant, the composition of communities across different depths generally followed the pattern of the field.

Figure 6.

Annotated COX1 Community Composition. PCoA plot of Jaccard-based community compositional dissimilarity for the mtMG Annotated COX1-based data. Dissimilarity explained by axes are presented in the axes labels and colours correspond to field (R1 = black, C1 = red, O1 = green, C2 = blue) with increasing depth represented by lighter shades of each field colour.

We observed significant differences in the detected taxa amongst all different methods. First, all the targeted mtMG approaches detected bacterial feeders (BF), plant parasites (PP), fungal feeders (FF), animal parasites (AP) and entomopathogenic parasitic nematodes (EPN). However, mtMG assemblies, based on read-mapping, produced ~ 40 fewer nematode taxa than the de novo assembled annotated COX1 approach. The greater richness recovered with the annotated COX1 approach was particularly evident for FF and BF taxa. In addition, although the annotated COX1 protocol detected higher total PP richness, it did so with less consistency. Specifically, while the annotated COX1 protocol detected multiple Meloidogyne spp., the read mapping approach retrieved multiple M. hapla haplotypes.

The contrasts amongst the mtMG protocols emphasises the influence of the differences of reference databases on the results. While the annotated COX1 approach produced the most comprehensive results, the read-mapping approach highlighted the presence of a diverse range of M. hapla sequences that matched those in the reference database, demonstrating the ability to assemble many of rare haplotypes that otherwise would have been missed (required 100% matching) with less targeted approaches. In contrast, since our reference database is heavily weighted towards PP and AP taxa, the read-mapping approach missed many or all of the under-represented taxa in other trophic groups (i.e. BF, FF, OM and PRs) due to a lack of close reference sequence matches. Although the annotated COX1 approach was much more capable of capturing this diversity, likely due to the flexibility of the ORF search, it was more prone to missing rare taxa during read consolidation within the assembly process (Bankevich et al. 2012; Lapidus and Korobeynikov 2021). This suggests that, to increase the recovery of rare and underrepresented taxa, it may be necessary to include a variety of mtMG data analysis approaches. Furthermore, without a significant expansion of reference sequences, the mtMG approach for nematode diversity surveys may be limited (Gendron et al. 2023), particularly with regards to identification of variant haplotypes (e.g. resistant to nematicides) (unpublished data from Braun et al. in the conference presentation Plant pathogenic nematodes behave badly: Uncovering the first case of field-relevant nematicide resistance at the Symposium of the European Society of Nematologists, Cordoba, Spain 2024, April 15–19).

In addition to differences amongst the three mtMG protocols, we also observed significant departures from the traditional morphology-based method, albeit similarities at the trophic group level were present. Differences were largely driven by AP, OM and PR taxa. While the morphology-based approach failed to detect any animal parasitic species in any sample, the mtMG-based approach failed to detect any omnivorous or predatory taxa. In contrast, BF, PP and FF trophic groups were similarly detected by both methods, albeit with the species ID accessible only with the mtMG approach. In terms of richness and identity of plant parasitic species, morphology and mtMG-based approaches showed significant differences. The mtMG method not only detected different Meloidogyne species, but also multiple haplotypes of M. hapla that were not differentiated with morphology, thus demonstrating a promising benefit in the capture of population level diversity in nematode communities. Unfortunately, except for Meloidogyne, none of the other PPNs detected by morphology (Tylenchorhynchus, Nanidorus and Belonolaimus) were recovered by the mtMG method. The differences also applied to PPNs detected with the mtMG approach, but not morphology (Mesocriconema and Xiphinema). There are several factors that could account for these differences. For instance, extracted DNA quality may affect the contig assembly process. However, nematode DNA has been shown stable for long periods of time, even after nematode death, as long as it is not exposed to natural soil conditions (Kunin et al. 2008; Christoforou et al. 2017). A more significant factor is likely the lack of mitochondrial reference sequences for these taxa, preventing detection of many PP, OM and PR taxa by the mtMG approach (Gendron et al. 2023).

These differences highlight the limitations and advantages of each method. On one hand, the targeted mtMG approach can capture a wide array of diversity, including APs likely present in samples as eggs and morphologically ambiguous dauers, while also reducing time and cost. However, the assembly process faces challenges in assembling and identifying contigs from rare taxa due to low read depth (Lapidus and Korobeynikov 2021) and limited reference database coverage (Gendron et al. 2023). Although read-mapping can help remedy contigs assembly difficulties, it relies on a database populated with reference sequences matching the reads (Xu et al. 2022). In contrast, the morphology approach excels at detecting rarer taxa, but falls short in identifying species and population level dynamics. Nevertheless, our study demonstrates that both methods capture factors predicting nematode communities and their combined use is complementary.

We can observe the contrast of strengths and weakness for each method, as well as the benefits of combining them, when analysing the diversity metrics. For instance, we found that nematode communities were primarily structured by the field, although the importance of this relationship varied across fields. By comparing the results of both methods, we were able to better identify which fields and specific variables drove these relationships. First, the influence of the field on community richness (i.e. trophic, PP and total richness) differed between the mtMG and morphological approaches. While the mtMG methods did not identify the role of the field in trophic or PP richness, they allowed us to pinpoint the field as the driver of the total nematode community richness. In comparison, morphology-based trophic and PP richness were mostly explained by the field, but, due to its limitations, we could not test the effect of field on the total community. Overall, our combined methods highlighted Field C2 and Field O1 as the drivers of differences in richness across our samples.

Both fields demonstrated some of the highest richness of the four fields, with Field C2 showing higher total richness (based on mtMG), but lowest PP richness (similar to Field C1; based on mtMG and morphology) and Field O1 some of the highest PP richness indicating that the richness of these two fields is driven by different trophic groups and taxa. However, Field O1 also showed the greatest difference between mtMG and morphology approaches. Specifically, PP richness, based on the targeted mtMG approach, was much more variable and trophic richness was significantly lower. As discussed above, the strengths and weakness of our methods explained why we observed these differences. The PP richness, based on morphology, resulted from the detection of M. hapla and B. longicaudatus, while based on mtMG from the detection of multiple Meloidogyne species and M. hapla haplotypes. However, the advantage of morphology in detecting the rarer species in these fields (i.e. Belonolaimus longicaudatus, Nanidorus minor and Tylenchorhynchus spp.) (123 combined counts out of 7239 total counts; Suppl. material 4; Bello et al. (2020)) would quickly diminish once their reference sequences become a part of the mtMG database (Gendron et al. 2023). In contrast, the mtMG approach excelled at the detection of different Meloidogyne species (i.e. M. luci and M. javanica) and haplotypes amongst its closest morphologically indistinguishable relatives.

Although the two methods exhibited some similarities at the trophic group level, they diverged at the species level. Much of the variation in richness between the two methods can be attributed to the detection (or lack thereof) of FFs, APs, EPNs, OMs and PRs. Field C2 and O1 already differed at the trophic level. For example, trophic richness in Field O1 showed the most pronounced differences between the two methods. Morphologically, Field O1 differed from C2 due to higher detection of OM and FF taxa. This increased detection in Field O1 likely resulted from organic management practices, leading to higher organic matter and the absence of broadcast pesticides compared to the conventionally managed fields (Sánchez-Moreno et al. 2008; Neher 2010). In contrast, targeted mtMG analysis revealed no OM taxa and few FF taxa, but multiple AP taxa were recovered. Trophic richness in Field C2 showed the smallest differences between methods, primarily driven by substantial overlap in the detection of FF and PP taxa, although OM, PR and AP taxa still exhibited disparities. The lack of molecular detection of OM and PR likely stems from their lower abundance relative to other trophic groups, resulting in low sequencing coverage and reduced likelihood of assembly. Moreover, the scarcity of reference sequences for OM and PR contributes to this limitation. While the low recovery of OM and PR taxa could potentially result from our use of double sieves to filter organic debris prior to DNA processing, the successful recovery of other large-bodied taxa (e.g. Xiphinema sp.) makes this scenario unlikely. Notably, the mtMG approach excelled in detecting AP taxa (i.e. Ancylostoma caninum, Heligmosomoides polygyrus, Ortleppascaris sinensis, Parastrongyloides trichosuri, Strongyloides venezuelensis, Trichostrogylus axei and Trichuris skrjabini). These taxa, often present in soils as eggs or unidentifiable juveniles (Qazi et al. 2020), are likely overlooked by morphology due to a lack of distinguishing characteristics. Furthermore, AP taxa have been extensively studied and their sequences dominate nematode databases, often overshadowing other trophic groups (Gendron et al. 2023). Although not statistically significant in our richness models, similar considerations apply to EPN taxa (i.e. Steinernema diaprepesi) which were detected by the mtMG method, but missed by morphology.

In contrast to the richness metrics, nematode community composition showed more consistent agreement between the two methods. For example, both morphology-based and mtMG-based communities (trophic, PP and total) were significantly influenced by field. Post-hoc analysis confirmed the role of field in structuring these communities, consistently identifying it as a significant factor shaping nematode community compositions. Although not all fields exhibited differences in all models, the overall findings remained consistent across the methods. Specifically, trophic community composition yielded similar statistical outcomes regardless of the method. Taken together, the community composition analysis emphasised that, despite variations in detected taxa, both methods converged on the underlying factors shaping nematode community composition. This underscores the reliability of both approaches in providing results representative of real ecological processes. In addition, using both methods allowed for a more detailed recovery of species. Indeed, both methods reliably identified field as the major explanatory factor, supporting our hypothesis of strong field-based structuring of these communities, likely driven by each field’s unique management practices.

While some soil depth results overlapped between the two approaches, they diverged from the field findings, revealing substantial differences. For example, the variation explained by depth for both mtMG- and morphology-based trophic community composition was nearly equal (~ 11.6%). The high degree of similarity in trophic diversity stemmed from both methods’ ability to detect BF and FF taxa. However, the morphology approach indicated that PP community composition was significantly influenced by depth, whereas the mtMG approach showed no significant explanation by depth. This discrepancy likely arises from the morphology-based recovery of Belonolaimus longicaudatus, Nanidorus minor and Tylenchorhynchus spp., which were predominantly abundant in the top soil layers, while detecting Meloidogyne hapla taxa throughout the soil column. Although the mtMG approach consistently recovered Meloidogyne species across multiple soil depths in a similar pattern to the morphology results, it detected few other PP taxa (and none of the morphologically other identified taxa) resulting in a relatively uniform PP community across the soil column. Notably, only annotated COX1 trophic richness and composition metrics derived from the targeted mtMG approach showed a significant relationship with depth, whereas all morphology-derived richness and composition metrics were significant. Additionally, post-hoc analyses revealed that morphology-based trophic community composition variation was influenced by multiple depths (0 – 75 cm), whereas mtMG showed evidence of being driven primarily by the deeper depths (50–75 cm). Overall, the lower abundance of taxa and lower likelihood of presence in lower soil depths may contribute to more variable results, affecting the consistency of analyses. Nevertheless, considering both approaches, these results highlight depth as a key factor structuring nematode communities, with lower depths playing a particularly important role.

Going into the study, we expected to find morphological evidence of M. hapla at lower depths indicative of active parasitism (Wesemael and Moens 2008), in part due to the high summer temperatures during sampling, but the targeted mtMG approach revealed a more complex phenomenon than expected by morphology. During adverse conditions or absence of a plant host, plant parasitic nematodes display a survival behaviour by migrating into deeper soil layers (Pinkerton et al. 1987). Morphology-based trophic and PP richness and composition were largely driven by variation in the upper soil layers (0 – 50 cm) and heavily influenced by the detection of PP taxa that were missing in the mtMG-based communities. As indicated by morphology, Nanidorus minor tended to be associated with the upper soil layers and, although the same pattern applied to Tylenchorhynchus spp. and B. longicaudatus, their presence was much more sporadic. M. hapla was the only plant parasite most consistently recovered by both methods across all soil depths. However, in addition to M. hapla and in contrast to morphology, the mtMG detected other Meloidogyne species. Taken together, these results support past findings showing the capacity of these PPNs (particularly root-knot species) to migrate through soil layers, likely in response to adverse conditions, such as temperature, tillage, fumigation or crop succession (Pinkerton et al. 1987; Mojtahedi et al. 1991; Eo et al. 2007; Wesemael and Moens 2008) and is likely further facilitated by the sand and clay making up the soil composition of the fields (Prot and Van Gundy 1981). This has potential implications for Meloidogyne management with soil layers potentially acting as reservoir from which post fumigation population re-establishment can occur (Giannakou et al. 2005).

The authors have declared that no competing interests exist.

No ethical statement was reported.

This study was funded through the USDA NIFA AWD07689 awarded to DL Porazinska and TO Powers and the Florida Strawberry Growers Association (FSGA) AWD09030 awarded to Johan Desaeger.

CO and JD conceptualised the experimental design, performed field investigations and sample collection. CO and EG performed the laboratory work. EG performed the formal data analysis, wrote the scripts and produced the figures for visualisation. CO and EG curated the data. EG and DP wrote the first draft of the manuscript and all authors reviewed and edited the manuscript. DP and JD provided funding, resources, supervision and overall project administration.

Eli M.S. Gendron https://orcid.org/0000-0002-2552-4658

Johan Desaeger https://orcid.org/0000-0001-5121-5691

Raw sequences and metadata have been deposited into the NCBI Sequence Read Archive (SRA) under the accession numbers SAMN40620041 - SAMN40620104 in BioProject PRJNA1092106. The pipeline and code used to process and analyse the data, along with community tables and corresponding metadata used for the analyses are available at https://github.com/WormsEtAl/mtMG-Analysis-of-Florida-Strawberry-Field-Nematode-Communities.