Integrative taxonomy and mitogenome characterization of the root-knot nematode Meloidogyne silvestris

Root-knot nematodes (RKNs) of the genus Meloidogyne are the most important plant-parasitic nematodes in agriculture, posing a major threat to global food security by affecting a wide range of economically important crops¹. To date, more than 105 validated Meloidogyne species have been described. Among these, M. incognita, M. javanica, M. arenaria, and M. hapla are considered the most widespread and damaging². These species parasitize a broad range of monocot and dicot hosts, including both herbaceous and woody plants³.

In Europe, several Meloidogyne species; including M. ardenensis, M. arenaria, M. baetica, M. hispanica, M. javanica, M. lusitanica, and M. hapla have been documented parasitizing woody plants across boreal, temperate, and mediterranean regions⁴. Some of these, such as M. baetica, appear largely restricted to woody hosts within their native ranges⁵. A more recently described member of this group is Meloidogyne silvestris, a species showing similar preference for woody plant hosts.

M. silvestris was first described from Ilex aquifolium host plant in Spain⁶. Our study documents its first occurrence beyond the type locality, with populations identified in two Dutch municipalities (Hilversum and Bennekom). According to the European and Mediterranean Plant Protection Organization (EPPO) quick scan report, these populations exhibited an extended host range, parasitizing Ulmus hollandica, U. glabra, and alongside the original I. aquifolium host plants⁷.

The systematics and phylogeny of Meloidogyne have undergone major revision over the last two decades, driven largely by the integration of molecular data. Early phylogenetic analyses based on the internal transcribed spacer (ITS) region revealed extensive cryptic diversity and incongruence with traditional morphology-based classifications⁸. More recently, comprehensive multigene and phylogenomic approaches have provided a robust framework for understanding evolutionary relationships within the genus, revealing multiple deeply divergent lineages and frequent cases of morphological convergence^9,10.

While morphological identification remains central to nematode taxonomy, it demands significant expertise and is frequently complicated by phenotypic plasticity, overlapping morphological and morphometric traits among closely related species^11,12,13. Isozyme methods such as malate dehydrogenase (Mdh; EC 1.1.1.37) and esterase (Est; EC 3.1.1.1.), while valuable for specific diagnostic applications, are limited to particular developmental stages typically young adult females and may yield inconsistent results across laboratories due to variation in protein extraction, sample storage and electrophoretic conditions¹⁴.

Molecular markers have greatly enhanced the precision and reliability of nematode identification. Nuclear ribosomal regions, including 18S rDNA, 28S rDNA, and the internal transcribed spacer (ITS), are widely used for species barcoding and community-level analyses^15,16. However, these markers also have limitations: the 18S gene is often too conserved to distinguish closely related species, while ITS and 28S regions may exhibit high intraspecific variability^17,18.

Mitochondrial DNA (mtDNA) markers, particularly protein-coding genes (PCGs) such as cox1, cox2, and cytb, have become increasingly important in diagnostics and phylogenetic studies^19,20. Advances in high-throughput sequencing have further expanded their utility in nematode diagnostics, phylogenomics, and comparative genomics^21,22,23. However, mitochondrial markers are not without caveats, as maternal inheritance patterns and potential introgression events can obscure phylogenetic interpretations, necessitating careful consideration in taxonomic assessments^24,25.

Given the respective strengths and limitations of each identification method, an integrative approach combining morphological, biochemical, and molecular data is increasingly recognized as essential for robust nematode diagnostics and taxonomy. Here, we applied such framework with three objectives: (1) report for the first time the occurrence of M. silvestris outside its type locality; (2) to expand the species phenotypic and molecular characterization through morphometrics, isozymes and multi-locus DNA analyses, and (3) present its first near-complete mitogenome for characterization and evolutionary insights.

Field sampling and nematode extraction

Field sampling was conducted in October 2022 at two locations in the Netherlands: Hilversum (52°23’53.9″N, 5°18’42.3″E) and Bennekom (51°59’38.1″N, 5°40’34.1″E). The Hilversum site was targeted following the EPPO quick scan that reported M. silvestris in the area, whereas the Bennekom site was selected based on the presence of host plants showing clear root-knot nematode feeding symptoms. At the Hilversum site, root samples were collected from Ilex aquifolium (common holly), while Euonymus fortunei (wintercreeper) roots were sampled in Bennekom.

Plant identification was carried out by a senior botanist at the Netherlands National Plant Protection Organization. The sampled species are listed as Least Concern on the International Union for Conservation of Nature (IUCN) Red List and are not included in the Convention on International Trade in Endangered Species of Wild Fauna and Flora (CITES); therefore, no special permits or licenses were required for sampling.

Fine root segments (< 1 cm length) were processed using the mist chamber extraction method²⁶ to isolate second-stage juveniles (J2s) and males of M. silvestris. Adult females were carefully extracted from root galls using fine dissection tools under a Leica stereomicroscope at 40× magnification. For diagnostic purposes, only J2s and females were used for primary taxonomic identification, as males of Meloidogyne species typically provide limited taxonomic information. Nevertheless, to complement our morphometric analyses, adult males were also examined.

Preparation for morphology, morphometry and biochemical analysis

Freshly extracted J2s were individually selected using a nematode picking needle under a Leica stereomicroscope and temporarily mounted for examination under a Leica DM2500 differential interference contrast (DIC) microscope (Leica Microsystems GmbH, Wetzlar, Germany) at 50×, 100×, 200×, 400×, 1000× magnification. High-resolution images were captured using a Leica DFC450 camera coupled with IMAGIC-IMS Professional Imaging software v21H1 (Bildverarbeitung AG, Glattbrugg, Zurich, Switzerland).

From each population, 20 J2 individuals were analyzed, with 14 primary morphological measurements and 5 derived morphometric indices recorded per specimen (19 variables in total). Additionally, 12 adult males were analyzed, with 15 primary morphological measurements and 5 derived morphometric indices recorded per individual (20 variables in total; Table 1). Morphologically characterized specimens were preserved at -20 °C for subsequent molecular analysis. For perineal pattern preparation, the anterior region (head and neck) of young females was removed, and the posterior end was cleared of residual body tissues to isolate the perineal pattern. Temporary slides were then prepared for observation.

Biochemical analysis

Isozyme analysis targeting Mdh; EC 1.1.1.37 and Est; EC 3.1.1.1, was conducted on all populations. The methodology for the Hilversum population followed Karssen et al., 1995²⁷. Briefly, ten young females were isolated and maintained in 0.9% NaCl. They were then desalted in distilled water on an ice bath for no longer than 5 min and individually macerated in 0.6 µL extraction buffer using a metal rod. The homogenized protein lysate was separated via polyacrylamide gel on a PhastGel system (Cytiva, Uppsala, Sweden).

For the Bennekom population, the same procedure was followed using ten young females, except that Novex™ WedgeWell™ Tris-Glycine Mini Gels (Invitrogen™, Thermo Fisher Scientific, Carlsbad, CA, USA) were used instead of the PhastGel system. In both cases, the gels were stained for Mdh and Est activity for 30 and 60 min, respectively. After staining, gels were rinsed with distilled water and fixed in a solution of 10% glycerol, 10% acetic acid, and distilled water. M. javanica was included as a reference at the center of the gels.

DNA extraction and Illumina sequencing

ClearDetection kit (Wageningen, The Netherlands) was used for DNA extraction of nematodes from both populations. Briefly, samples were lysed with proteinase K and Dithiothreitol (DTT) in 50 µL of 2x extraction buffer, incubated at 65 °C for 30 min with shaking at 800 rpm, followed by incubation at 95 °C for 5 min. The condensate was briefly centrifuged for (30 s at 1500 × g) and stored at -20 °C.

About 30 µL of extracted genomic DNA (gDNA) from each sample was sent to GenomeScan (Leiden, The Netherlands) for 150 bp paired-end sequencing using Illumina technology. Prior to library preparation, the gDNA was fragmented with a Bioruptor Pico (Diagenode, Belgium). Library preparation was conducted using the NEBNext Ultra II DNA Library Prep Kit (New England Biolabs, Ipswich, England), following the manufacturer’s instructions.

DNA quality and concentration before and after library preparation were assessed with 2100 Bioanalyzer, developed by Agilent, California, USA. Commercially available Cat (Felis silvestris catus L) was included as Positive process control (PPC) to detect possible contamination during the library prep.

Sequence assembly, annotation and barcoding

Raw Illumina reads were processed using CLC Genomics Workbench v23.0.4 (Qiagen, USA) with an automated pipeline²⁸. Contigs were quality-filtered and visualized using Krona²⁹.

For ribosomal DNA (rDNA), only 18S and 28S genes were retained for downstream analysis. Fragmented rDNA sequences were de novo assembled in Geneious Prime v2020.0.5 (Biomatters Ltd., Auckland, New Zealand) (http://www.geneious.com/) to generate complete contigs when necessary.

Putative mtDNA contigs were identified via BLASTn against a custom database and annotated using MITOS and MITOS v2^30,31. PCGs were manually curated for start/stop codons. Near complete mitogenome was plotted using OGDraw³².and nucleotide composition was analysed with strand asymmetry calculated as:

rDNA and mtDNA barcodes were generated computationally (in silico) in Geneious Prime using primer sets (Table 2), and alignments and similarity estimates were computed with MAFFT v7.490³³. Secondary structures of predicted tRNAs were visualized using R2DT framework³⁴.

Phylogenetic and mitogenomic analysis

All phylogenetic trees were generated using Maximum Likelihood (ML) in MEGA v11.0.8³⁵ with 1000 bootstrap replicates³⁶. Nuclear rDNA phylogenies were inferred from 18S and D2–D3 28S rDNA alignments (2,105 bp and 1,106 bp, respectively). The alignment length reflects sequence length variation among taxa and gaps introduced during alignment. The best-fit substitution models, selected based on the Bayesian Information Criterion (BIC), were GTR + G+I for both the 18 S and 28S rDNA datasets. The generated trees from 18S rDNA dataset (48 Meloidogyne sequences) and 28S rDNA dataset (44 Meloidogyne sequences) were both rooted with Pratylenchus penetrans and Pratylenchus vulnus. GenBank accession numbers of gene sequences from Meloidogyne species used in the phylogenetic analysis are provided in Table S1.

For mitogenomic reconstruction, amino acid sequences of 12 mitochondrial PCGs from 33 nematodes (Chromadorea and Enoplea) were aligned using MAFFT in MEGA and concatenated. Poorly conserved or uninformative regions were removed with Gblocks³⁷, resulting in a final alignment of 1,284 aa. The best-fit substitution model, LG + G+I + F, was selected based on the BIC. The phylogeny was rooted using two arthropod outgroups: Limulus polyphemus and Lithobius forficatus. Mitochondrial synteny was analyzed across nematode taxa (including M. silvestris) to identify gene order rearrangements. GenBank accession numbers of mitochondrial gene sequences used in these analyses are provided in Table S2.

Morphology, morphometry and biochemical analysis

Females and J2s of M. silvestris from both populations exhibited morphological and morphometric uniformity to the original description, with slight variation in J2s body length, maximum body diameter, “a” ratio and stylet length. Among all measured traits, body length showed the highest variability (Table 1). Perineal patterns were similar between the two examined Dutch populations and consistent with the original description (Fig. S1d).

Biochemical profiles were also similar with the original species description⁶. M. silvestris from both populations showed a single A1 esterase (Est) band and an N1c malate dehydrogenase (Mdh) phenotype. Reference species M. javanica exhibited J3 Est and N1 Mdh phenotypes³⁸(Fig. 1).

DNA extraction for Illumina and de novo assembly

Measured DNA concentrations ranged in each sample ranged from 1.52 to 2.06 ng/µL. Illumina sequencing yielded a total output of 2Gb raw data and approximately 5.2 to 60.4 million reads per sample, resulting in 5,372 to 130,950 assembled contigs. The majority of assembled contigs did not produce BLASTn hits or aligned primarily to bacterial and fungal sequences. However, BLASTn searches against a locally installed NCBI nr/nt database produced 5–23 ribosomal DNA (rDNA) contigs and 2–6 mitochondrial DNA (mtDNA) contigs per sample.

Ribosomal and mitochondrial DNA Sequences

From both Dutch populations, a total of eleven annotated rDNA contigs were obtained. For the 18S rDNA region, both partial (< 1000 bp) and complete (1730 bp) barcode sequences were retrieved. 18S rDNA contigs from the two Dutch populations showed 99.9% similarity to each other and ~ 99.36% similarity to the Spanish population. For 28S rDNA (768 bp), the D2–D3 region showed 99.7% similarity within Dutch populations and ~ 99.35% similarity to the Spanish population. Barcode sequences of 450 bp (cox1) and 432 bp (cox2) were obtained, with representative alignments showing 100% identity between the two populations (Fig. S2 c-d).

Phylogenetic and mitogenomic analysis

Phylogenetic analysis of 18S rDNA placed M. silvestris sequences from both sampling sites in a clade with conspecifics from the type locality, albeit with low bootstrap support (44%). This clade also included M. ardenensis, M. dunensis, M. dutysi, M. graminis, M. hapla, M. maritima, M. marylandi, M. microtyla, M. partityla, M. spartelensis and M. spartinae (Fig. 2). For 28S rDNA, sequences of M. silvestris from both sampling sites and the type locality formed a highly supported monophyletic clade (100% bootstrap). This clade also showed a sister-group relationship with M. dunensis and M. spartelensis though with moderate support (76% bootstrap) (Fig. 3).

The amino acid ML mitogenomic tree grouped M. silvestris with M. hapla, while tropical (M. incognita, M. javanica, M. arenaria, M. enterolobii) and temperate (M. oryzae, M. graminicola, M. chitwoodi) RKNs formed two distinct, well-supported subclades (bootstrap support > 90%). Notably, the tree also showed that P. vulnus clustered more closely with Meloidogyne spp. than with H. glycines (moderate bootstrap support; Fig. 4, left). Comparative synteny analysis showed that the PCG order in M. silvestris was conserved across other Meloidogyne species and shared more similarity (four genes at the same position) with P. vulnus than with H. glycines (one gene at the same position; Fig. 4, right).

Mitochondrial genome of M. silvestris

The nearly complete mitochondrial genome of M. silvestris was assembled as a circular DNA molecule of 15,048 bp (Fig. 5). Its nucleotide composition was: A = 27.3%, T = 49.3%, C = 6.3%, G = 17.1% resulting in a strong A + T bias of 76.6%. The AT skew and GC skew were − 0.2866 and 0.4623, respectively.

Twelve protein-coding genes (PCGs) were annotated, including: cox1, cox2, cox3 (cytochrome c oxidase subunits), cob (cytochrome b) nad1-6, nad4L (NADH dehydrogenase subunits) and atp6 (ATP synthase subunit) comprising 10,053 bp residues. The longest genes were cox1 (1525 bp) and nad5 (1533 bp); the shortest were nad4L (176 bp) and nad3 (324 bp) (Table 3).

Two gene overlaps were observed: 31 bp between nad5 and atp6, 18 bp between nad3 and cob. The longest intergenic spacer (382 bp) was located between cox2 and 16S rRNA gene. Eighteen tRNA genes were identified ranging from 53 bp to 59 bp, all encoded in the same orientation as PCGs (Fig. 5). The predicted secondary structures of these tRNAs are provided in Fig. S3. Four tRNAs (trnA, trnH, trnQ, and trnR) were not recovered. The 16S rRNA gene (rrnL, 468 bp) was located between cox2 and nad3, while the 12S rRNA gene (rrnS, 403 bp) was found between cox1 and nad1, both showing similar nucleotide composition. A long non-coding region (2118 bp) was observed between trnD and trnL1.

Our combined morphological, biochemical, and molecular analyses provide a strong support that geographically separated populations of M. silvestris are conspecific, consistent with the original type population from Spain⁶. This is further supported by the consistently high sequence similarity observed across all markers.

Females and J2s from both Bennekom and Hilversum populations were largely morphologically uniform to the original description, with only slight variations observed in J2s body length, maximum body diameter, “a” ratio and stylet length. These levels of intraspecific variation are comparable to those reported in previous studies^39,40. Isoenzyme analyses further supported species identity: both populations exhibited an A1 esterase (Est) band and an N1c malate dehydrogenase (Mdh) phenotype⁴¹, clearly distinct from reference species M. javanica.

ML phylogenetic trees based on 18S and 28S rDNA sequences complemented the phenotypic findings and were congruent with previously published maximum parsimony (MP) analyses^6,42. As expected for conserved markers, 18S rDNA-based topologies exhibited limited species-level resolution, reflected in low bootstrap support for M. silvestris. By contrast, the D2–D3 region of 28S rDNA, which contains hypervariable domains, provided stronger phylogenetic signal, yielding a monophyletic M. silvestris clade with M. dunensis and M. spartelensis as a sister lineage. This pattern aligns with previous studies showing that 18S evolves slowly and is suitable for broad phylogenetic placement, whereas 28S evolves more rapidly and is highly informative for resolving closely related species^43,44.

Phylogenetic reconstructions from concatenated PCGs, together with synteny comparisons, provided robust support for the species-level distinctiveness of M. silvestris and its broader placement within the genus. Consistent with patterns inferred from nad5 markers⁴⁵, our topologies recovered distinct tropical and temperate clades of root-knot nematodes, with M. silvestris occupying an intermediate position between these lineages.

Comparative synteny analysis showed that the PCG order in M. silvestris is largely conserved across Meloidogyne species and is more similar to Pratylenchus vulnus than to Heterodera glycines, consistent with previous studies indicating a closer evolutionary relationship between Pratylenchus and RKNs^46,47. This suggests that mitochondrial gene order carries phylogenetic signals that complement sequence-based analyses and may help resolve relationships among plant-parasitic nematodes⁴⁸. Within the limited dataset analysed in this study, Chromadorean nematodes appear to exhibit a more conserved mitochondrial gene order than the Enoplean taxa included. Nevertheless, the small number of Enoplea representatives prevents broader generalization across the clade.

The nearly complete mitochondrial genome of M. silvestris is a 15,048 bp circular molecule, encoding 12 protein-coding genes (cox1-3, cob, nad1- 6, nad4L, atp6), 18 tRNAs, and two rRNAs, consistent with other Meloidogyne species, which lack the atp8 gene but otherwise exhibit a conserved PCG repertoire^49,50,51. Genome size variation among root-knot nematodes is often driven by differences in the length and composition of non-coding regions⁵². Although mitogenome sizes of ~ 15 kb falls within the broader range reported for plant-parasitic nematodes, published mitochondrial genomes of Meloidogyne species are typically larger (> 17 kb). The relatively small size of the 15,048 bp mitogenome of M. silvestris therefore represents a deviation from the common size range for the genus and may reflect reduction in non-coding regions. The non-coding region is probably involved in regulatory functions analogous to the control regions described in M. incognita and M. chitwoodi⁴⁹. Such regions often contain tandem repeats and secondary structures involved in replication, and their variation may contribute to intraspecific diversity⁵³.

The strong A + T bias (76.6%) and pronounced AT (− 0.29) and GC (0.46) skews mirror patterns described for other chromadorean nematodes⁴⁸. Such biases are thought to reflect asymmetric mutational processes during replication and transcription⁵⁴. Consequently, amino acid alignments and models that account for compositional heterogeneity should be employed in phylogenomic analyses^55,56.

Two gene overlaps (nad5-atp6 and nad3-cob) and a reduced tRNA complement (missing trnA, trnH, trnQ, and trnR) illustrate mitogenome compaction and may result from stop/start codon sharing⁵⁷. To the best of our knowledge, no studies have reported a complete absence of nematode tRNAs; however, reduced or aberrant tRNAs are relatively common and may complicate genome annotation⁵⁸. Functional compensation by nuclear-encoded tRNAs has been suggested in cases of true loss⁵⁹. Manual inspection of predicted structures or experimental validation could help resolve whether these tRNAs are absent or simply undetected by computational tools.

Across multiple lines of evidence, Meloidogyne silvestris populations showed striking uniformity, with concordant morphological, biochemical, and molecular data, including identical cox1 and cox2 barcodes, near-identical rDNA sequences, and a conserved mitogenomic architecture. These results confirm the taxonomic integrity of M. silvestris and demonstrate the power of integrative datasets for species delimitation in root-knot nematodes. The near-complete mitochondrial genome presented here expands genomic resources for the genus and provides a foundation for future studies of population structure, evolutionary dynamics, and the development of robust diagnostic markers.