The various Helicotylenchus spp. populations were recovered from rhizospheric soil samples collected from 9 golf courses. Detailed information regarding the sample localities, host plant, and laboratory codes is summarized in Table 1. The nematodes were extracted from sod samples using a combination of modified Cobb’s sieving and Baermann funnel methods (Whitehead and Hemming, 1965). With the aid of a Nikon SMZ 1000 stereomicroscope (Nikon, Tokyo, Japan), nematode specimens belonging to Helicotylenchus spp. were picked out from the collected nematode suspension and subsequently characterized to species level based on inferences from morphometrics and DNA sequence data.

A defined number of specimens from each population were heat-killed, fixed, and mounted to pure glycerin (Seinhorst, 1959). Morphometric data were taken using a Zeiss imager Z2 microscope (Carl Zeiss, Jena, Germany) fitted with Axio-vision software Material Science Software for Research and Engineering (Carl Zeiss). Morphological species identification was done following the identification keys and species descriptions compiled by Boag and Jairajpuri (1985), Fortuner (1984), Marais (2001), Mwamula et al. (2020b), and Subbotin et al. (2015).

Ribosomal and mitochondrial DNA was extracted from single female and male specimens using the DNA extraction kit WizPure, as described by Iwahori et al. (2000). The 18S-rRNA and COI gene of the previously published Helicotylenchus populations (Mwamula et al., 2020b) were amplified from their stored DNA. The nearly full-length 18S-rRNA gene was amplified as 2 partially overlapping fragments using 2 primer sets: 988F (5′-CTCAAAGATTAAGCCATGC-3′) and 1912R (5′-TTTACGGTCAGAACTAGGG-3′), 1813F (5′-CTGCGTGAGAGGTGAAAT-3′), and 2646R (5′-GCTACCTTGTTACGACTTTT-3′) (Holterman et al., 2006); D2A (5′-ACAAGTACCGTGAGGGAAAGTTG-3′) and D3B (5′-TCGGAAGGAACCAGCTACTA-3′) (Nunn, 1992) amplified the D2-D3 expansion segment of 28S-rRNA; TW81 (5′-GTTTCCGTAGGTGAACCTGC-3′) and AB28 (5′-ATATGCTTAAGTTCAGCGGGT-3′) (Curran et al., 1994), and S-ITS1 (5′-TTGATTACGTCCCTGCCCTTT-3′) and Vrain2R (5′-TTTCACTCGCCGTTACTAAGGGAATC-3′) (Vrain et al., 1992) amplified the partial ITS-rRNA region. The partial COI gene was amplified using a combination of several primer sets: COIF5/COIR9; JB3/JB5 and COIHF/COIPR as detailed in Table 2. Polymerase chain reaction (PCR) was performed with a PCR cycler (T100, Bio-Rad, Hercules, CA, USA), the PCR program with specific primers being set as follows: initial denaturation at 95°C for 5 min, 35 cycles at 95°C for 30 s, followed by an annealing step at 58°C for 30 s (S-ITS1/Vrain2R), 54°C for 40 s (COIHF/COIPR), 56°C for 40 s (COIF5/COIR9); 72°C for 60 s (S-ITS1/Vrain2R), 72°C for 80 s (COIHF/COIPR and COIF5/COIR9), and finally one cycle at 72°C for 5 min (COIHF/COIPR) and 72°C for 10 min (S-ITS1/Vrain2R and COIF5/COIR9). The cycle for D2A/D3B and JB3/JB5 primer sets was as described by Mwamula et al. (2023a) while the thermal profiles for TW81/AB28, 988F/1912R and 1813F/2646R primer sets were set as detailed by Mwamula et al. (2023b). The PCR products were purified using a PCR purification kit (Qiagen, Hilden, Germany) and quantified using a quickdrop spectrophotometer (Molecular Devices, San Jose, CA, USA). The purified products were subsequently used for direct sequencing in both directions using the same primers as specified above. DNA sequencing was performed at Macrogen. The edited DNA sequences were submitted to the GenBank database under the accession nos.: OR843993-OR844017 (18S-rRNA); OR844019-OR844030 (28S-rRNA); OR852617-OR852632 (ITS-rRNA); OR844045-OR844080 and OR855927-OR855934 (COI gene).

The obtained and edited sequences (18S-rRNA, 28S-rRNA, ITS-rRNA and COI gene) were aligned using ClustalX (Thompson et al., 1997) along with the sequence data sets of Helicotylenchus species published in GenBank (Divsalar et al., 2020; Ehtesham et al., 2021; Etongwe et al., 2020; Mohammadi Zameleh et al., 2020; Mwamula et al., 2020b; Palomares-Rius et al., 2018; Riascos-Ortiz et al., 2020; Rybarczyk-Mydłowska et al., 2019; Shokoohi et al., 2018, 2021; Subbotin et al., 2011, 2015; Xia et al., 2022). Only a representative number of 18S-rRNA sequences were included in the data set for analysis. The data sets for 28S-rRNA, ITS-rRNA, and COI gene comprised essentially all the comparable sequences of Helicotylenchus species available in GenBank database at the time of the current analysis. However, the too divergent sequences possibly belonging to different genera or deemed to contain errors were omitted from phylogenetic reconstruction. The partial sequences of Rotylenchus urmiaensis (KP718969) and Scutellonema bradys (AJ966504) were selected as the outgroup taxa for the 18S-rRNA gene; Hoplolaimus galeatus (EU626787) and Hoplolaimus seinhorsti (DQ328752) for 28S-rRNA; Scutellonema bradys (JX472067) and Rotylenchus pumilus (JX015436) for ITS-rRNA; and Bursaphelenchus mucronatus (AB067765) and Bursaphelenchus sinensis (AB232163) for COI gene. The generated alignments were analyzed with Bayesian inference (BI) using MrBayes 3.2.7 (Ronquist et al., 2012) under the GTR + I + G model. BI analysis for each gene was initiated with a random starting tree and run with four chains for 1 × 10⁶ generations. The Markov chains were sampled at intervals of 100 generations. After discarding burn-in samples, a consensus tree was generated with a 50% majority rule. Trees were subsequently visualized and edited using FigTree v1.4.4 software. Posterior probabilities (PP) exceeding 50% are given on appropriate clades. Interspecific and intraspecific sequence variations were examined using PAUP* v4.0a169 (Swofford, 2003).

Identification of the newly recovered populations was made using an integrative approach considering morphological characteristics and the results of sequence analyses and phylogenetic inferences. A total of 97 sequences (25 nearly full-length 18S-rRNA, 12 partial 28S-rRNA, 16 partial ITS-rRNA, and 44 partial COI gene sequences) were newly obtained from the various populations of the four studied species in the present study. All the newly recovered populations were identified morphologically as representatives of the same species published by Mwamula et al. (2020b). The morphological and morphometric data of the current populations of H. microlobus agree well with those of Mwamula et al. (2020b), Siddiqi (1972, 2000), Subbotin et al. (2015), and Tzortzakakis et al. (2018). Helicotylenchus dihystera population YJ2 and H. caudatus Sultan 1985 population PJ3 are morphologically and morphometrically similar to the topotypes detailed by Sher (1966) and Sultan (1985), respectively. Helicotylenchus asiaticus is herein recorded for the first time outside its type locality. The morphometric data for the newly recovered populations are detailed in Table 3. The partial COI gene sequences of the studied species (H. microlobus, H. dihystera, and H. asiaticus) are herein supplied for the first time. Similarly, the partial 18S-rRNA gene sequences of H. caudatus and H. asiaticus are also reported for the first time.

The amplification of the nearly full-length 18S-rRNA gene from all the studied species and their respective populations yielded fragments of approximately 1,650-1,700 bp. The partial 18S-rRNA gene sequence alignment was 1,685 bp in length and contained 103 sequences of Helicotylenchus and 2 sequences of the outgroup taxa (Scutellonema bradys and Rotylenchus urmiaensis). Phylogenetic relationships between Helicotylenchus species, as inferred from Bayesian phylogenetic analysis of the dataset are shown in Fig. 1. Seven moderately supported major clades were distinguished within the analyzed taxa. Sequence variation within all taxa including outgroups reached 114 bp (6.8%) and for Helicotylenchus, it reached 91 bp (5.4%). Minimal to moderate sequence variation (0-45 bp [0.0-2.7%]) was evident among taxa belonging to Clade I (PP = 84%). The clade comprised eleven putative taxa including H. caudatus, H. indicus, H. multicinctus, H. digitiformis, H. paraplatyurus, H. pseudorobustus, H. certus, H. microlobus, H. crenacauda, H. dihystera, H. rotundicauda and two unidentified populations (Helicotylenchus sp. 1 JH-2014 and Helicotylenchus sp. 5 JH-2014).

All four newly obtained sequences of H. caudatus were identical with no recorded intraspecific variation and differed from the closely clustered H. indicus sequences by 23-26 bp (1.4-1.6%). Intraspecific sequence variation for the newly obtained H. microlobus sequences varied by 0-2 bp and these differed from H. microlobus populations from China; OP225352 and ON117514 by 0-9 bp (0.0-0.6%) and from the closely clustered taxa (H. digitiformis, H. crenacauda, H. pseudorobustus, H. certus, Helicotylenchus sp. 5 JH-2014, Helicotylenchus sp. 1 JH-2014, including H. dihystera accessions: KJ934127 and JX069950) by 0-26 bp (0.0-1.6%). An intraspecific sequence variation of 0-2 bp was recorded among the obtained H. dihystera sequences. The new sequences differed from other H. dihystera populations by 0-24 bp (0.0-1.5%) and from the closely clustered H. rotundicauda (KM603516 and OP225358) by 2-8 bp. Helicotylenchus asiaticus sequences formed a moderately supported independent clade (PP = 87%), with no recorded intraspecific sequence variation.

The amplicons of the D2-D3 expansion segment of 28S-rRNA gene were approximately 750 bp in length. The partial 28S-rRNA gene sequence alignment constituted all the comparable sequences of Helicotylenchus published in GenBank at the time of this analysis. A few incomparable sequences were excluded from the analysis (see remarks on omitted sequences). The alignment contained 471 sequences of Helicotylenchus including 2 sequences of the outgroup taxa (Hoplolaimus galeatus and Hoplolaimus seinhorsti). Twelve highly and moderately supported clades were distinguished within the taxa. Phylogenetic relationships between Helicotylenchus species, as inferred from Bayesian phylogenetic analysis of the dataset are shown in Fig. 2A-E. The general topology of the resultant Bayesian tree agrees with the phylogenetic data presented by Subbotin et al. (2015). Intraspecific sequence variation for the three newly obtained sequences of H. microlobus was 0-1 bp, and the intraspecific sequence variation for all H. microlobus sequences varied from 0-10 bp (0.0-1.8%); for H. pseudorobustus − 0-11 bp (0.0-2.0%); and for H. paxilli − 0-5 bp (0.0-0.8%). Interspecific sequence diversity within H. microlobus, H. paxilli, and H. pseudorobustus sequences reached 29 bp (6.0%) (15-29 bp [2.6-6.0%]). Intraspecific sequence variation for H. digonicus “Type A” varied from 0-15 bp (0.0-2.6%); and for H. broadbalkiensis − 1-6 bp (0.1-1.0%). Sequence variation within all comparable H. dihystera sequences (except ON117414) varied from 0-15 bp (0.0-2.6%). High sequence divergences of 11-21 bp (1.9-3.6%) were recorded in pairwise comparison of ON117414 with all other H. dihystera accessions. Intraspecific sequence variation for H. paraplatyurus including Helicotylenchus sp. X-1 (KM506841-KM506843) varied by 0-3 bp (0.0-0.5%) and for H. californicus − 0-1 bp. Sequences of H. caudatus from Korea (including the newly obtained sequences) varied by 1-6 bp (0.1-0.8%) but differed from the shorter sequences of H. caudatus from Ethiopia (OP647413-OP647416) by 12-16 (2.7-3.7%). Intraspecific sequence variation within H. asiaticus was 0-2 bp and these differed from the closely clustered unidentified species from Korea (KY484829-KY484832) by 5-9 bp (0.6-1.2%). Sequence diversity within H. multicinctus varied by 0-13 bp (0.0-2.5%); for H. labiodiscinus − 3-9 bp (0.5-1.5%); for H. oleae − 0-8 bp (0.0-1.1%); for H. cuspicaudatus − 3-12 bp (0.3-1.8%); for H. microcephalus including Helicotylenchus sp. VII − 0-6 bp (0.0-0.8%); for H. digonicus “Type B” and H. minzi − 0-13 bp (0.0-2.2%); and for H. vulgaris Type A and H. canadensis − 0-6 bp (0.0-1.0%).

The amplification of the partial ITS-rRNA gene yielded single amplicons of approximately 1,200 bp. The ITS-rRNA gene sequence alignment contained 185 sequences of Helicotylenchus including 2 sequences of outgroup taxa (Scutellonema bradys and Rotylenchus pumilus). Phylogenetic relationships between Helicotylenchus species, as inferred from Bayesian phylogenetic analysis of the dataset are shown in Fig. 3A-C. Ten highly and moderately supported clades were distinguished within the taxa. ITS-rRNA gene sequence variation for H. microlobus populations varied from 0-48 bp (0.0-4.9%); for H. digonicus − 4-16 bp (0.4-1.6%), differing from the closely clustered H. ciceri (MN194203 and MN194204) by 85-96 bp (8.8-9.8%); for H. broadbalkiensis − 12-14 bp (1.2-1.5%), differing from H. scoticus (MN194205) by 40-50 bp (4.0-5.0%); for H. pseudorobustus sensu stricto populations − 0-5 bp (0.0-0.6%); and for H. paxilli − 1-15 bp (0.1-1.5%). Sequence diversity within all but four (ON123797-ON123799 and MH142618) H. dihystera sequences varied from 0-33 bp (0.0-3.2%). The accessions ON123797-ON123799 and MH142618 contain a significant number of nucleotide deletions, insertions, and substitutions which might have been caused by sequencing errors. Pairwise comparison with other H. dihystera sequences showed sequence variation of 21-59 bp (2.5-5.4%). Despite the disparities, these accessions clustered together with other H. dihystera sequences (Fig. 3A). The top BLASTN hit for the H. caudatus ITS-rRNA gene sequences was H. crenacauda (DQ309586) with a percentage identity of 85% and 100% query coverage. The newly obtained H. caudatus sequences varied by 0-2 bp and differed from the closely clustered H. abunaamai (MW488036, MW488037), H. crenacauda (DQ309586), H. caudatus from Ethiopia (OP650246) and H. indicus (ON128730, ON128731), by 29-59 (7.2-9.2%), 92-158 bp (11.4-13.8%), 93-130 (14.9-15.6%), and 114-158 (16.6-19.5%), respectively. Intraspecific sequence variation within H. multicinctus varied from 1-37 bp (0.08-3.0%); for H. oleae − 5-24 bp (0.9-3.0%); for H. cuspicaudatus − 1-29 bp (0.07-2.3%); for H. microcephalus − 0-26 bp (0.0-2.3%); and for H. asiaticus − 0-8 bp (0.0-0.7%), differing from closely clustered Helicotylenchus spp. (KY484812-KY484814, KY484818-KY484822, KY512778, KY512779, and OR159472) by 42-56 bp (4.3-5.2%).

Amplification of the partial CO1 gene using the primer sets: COIF5/COIR9, COIHF/COIPR, and JB3/JB5 yielded single amplicons of ca 700, 400, and 350 bp in length, respectively. By using the primer set COIHF/COIPR, we were able to easily obtain amplicons for all four species. However, due to limited DNA template for H. caudatus and H. dihystera, amplicon sequencing was only done for two species (H. microlobus and H. asiaticus) and their respective studied populations. JB3/JB5 successfully amplified the target region for two (H. microlobus and H. dihystera) of the four species while COIF5/COIR9 was successful for only H. microlobus (Table 4). We also included the recently designed primer set M2F (5′-ATTGGiGSTTTTGGTAATT-3′) and RH1R (5′-CCAACAATGAATATATGATG-3′) (Rybarczyk-Mydłowska et al., 2019). However, amplification of the target region was not successful despite M2F/RH1R being the most successful primer set in the molecular studies performed by Rybarczyk-Mydłowska et al. (2019). The partial COI gene sequence alignment used in phylogenetic reconstruction was limited to 398 bp in length and contained 63 sequences of Helicotylenchus and two sequences of Bursaphelenchus used as outgroup taxa. Relationships between Helicotylenchus species, as inferred from Bayesian phylogenetic analysis of the dataset are shown in Fig. 4. Pairwise sequence comparison between the contentious H. microlobus and H. pseudorobustus was done using a longer alignment of 691 bp. Helicotylenchus microlobus differed from H. pseudorobustus by 13.4% (84 bp of the 628 aligned nucleotides) to 13.7% (64 bp of the 465 aligned nucleotides). Intraspecific sequence variation within all the newly obtained sequences of H. microlobus varied by 0-3 bp. Helicotylenchus microlobus and H. pseudorobustus clustered together with the new sequences of H. dihystera in a well-supported clade (PP = 100%). Intraspecific sequence variation within H. dihystera varied by only 0-1 bp. The newly obtained sequences of H. asiaticus also varied by 0-1 bp, and clustered together with H. oleae, H. varicaudatus, H. canadensis, and H. digonicus.

The pairwise sequence comparisons and phylogenetic analyses revealed several possible misidentifications or mislabeling of several Helicotylenchus species in GenBank database. Suggested changes are indicated with double asterisk signs on the phylogenetic tree Figs. 2A-E and 3A-C. A comparison of the D2-D3 sequence data accessioned ON005156 and identified as H. dihystera with H. microlobus sequences revealed that it belongs to H. microlobus. Pairwise sequence comparison revealed intraspecific sequence variation of 0-9 bp (0.0-1.4%). The accession also significantly differed from all H. dihystera sequences by 17-30 bp (2.2-4.2%). Accessions MF996708, MF996709, MG925220, MK358143, and OK256081-OK256084 (all identified as H. pseudorobustus) clustered with H. microlobus. Sequence variation varied from 0-18 bp (0.0-3.1%); with the highest variation of 13-18 bp being recorded in sequence comparison with the short sequences accessioned MF996709 and OK256082. The accession OP297979 and OP297982 identified as H. microlobus and Helicotylenchus sp., respectively, probably belong to H. paxilli (differed from H. paxilli by 1-3 bp and 4-5 bp, respectively); OK256086-OK256088, and OK256092 identified as H. pseudorobustus probably belong to H. digonicus “type A” (differed from H. digonicus “type A” by 0-8 bp [0.0-1.4%]); MN783728-MN783732 (Helicotylenchus sp.) belong to H. pseudorobustus (differed from H. pseudorobustus by 0-10 bp [0.0-2.0%]); ON014758 (identified as H. microlobus) belong to H. dihystera (showed intraspecific variation with H. dihystera of 0-9 bp, and differed from H. microlobus by 18-33 bp [3.1-5.7%]); accession OK275497, identified as H. dihystera is herein designated as Helicotylenchus sp. III-4 based on the grouping of Subbotin et al. (2015). It differed from H. dihystera sequences by 22-34 bp (2.9-4.5%); OP297983 (identified as Helicotylenchus sp.), MN922339, and OP297977 (identified as H. dihystera) probably belong to Helicotylenchus sp. III-1 (intraspecific variation of 2-13 bp [0.3-2.2%] with Helicotylenchus sp. III-1 but differed from H. dihystera by 22-36 bp [2.9-4.8%]); KM506841-KM506843 (designated as Helicotylenchus sp. X-1 by Subbotin et al., 2015) belong to H. paraplatyurus (intraspecific variation of 1-3 bp); KM506847 (Helicotylenchus sp. V-2) belongs to H. caudatus (intraspecific variation of 0-4 bp); HM014301 (Helicotylenchus sp. V-II) belongs to the same taxa as H. microcephalus (intraspecific variation of 1-5 bp); and DQ328755 (Helicotylenchus sp. IX-3) belongs to the same taxa as H. varicaudatus “Type B” (MN783723-MN783725) (intraspecific variation of 3-5 bp).

Based on the ITS-rRNA gene sequence analysis, accession FJ427209 identified as H. dihystera from Malaysia probably represents unidentified species. It significantly differed from all other H. dihystera sequences by 53-101 bp (7.6-12.8%) and formed a separate independent clade from other H. dihystera accessions (Fig. 3A). Similarly, accession GQ906354 identified by Zhao et al. (2010) as H. dihystera is also too divergent from all Helicotylenchus spp. (up to over 200 bp, with 80% coverage) and formed an independent clade (Fig. 3C). Accession GQ906356 also identified by Zhao et al. (2010) as H. crenacauda most likely belongs to H. microlobus. The sequence evidently contains multiple nucleotide insertion and deletions (indels), including nucleotide substitutions. However, the accession clustered together with H. microlobus sequences and differed by 13-60 bp (1.4-6.7%). The alterations might be a result of sequencing errors, misassembles or presence of cloning vectors.

Some of the notable excluded sequences included the following accessions. For the 28S-rRNA gene, accessions MH087063, MH087064, MH084946, MH087062, and MF996772 (all identified as H. multicinctus from India) were omitted due to their significant divergence from all H. multicinctus sequences (differences of up to 200 bp) and all other Helicotylenchus spp. These accessions probably represent inaccurately edited sequence data. Accession MT321730 (H. dihystera) was omitted due to the existence of unusual multiple degenerate DNA letter codes. For the ITS-rRNA gene, MW404014 belongs to Merlinius sp. but was mistakenly labeled as H. digonicus (see Abdulsalam et al., 2021); GQ906351-GQ906353 identified by Zhao et al. (2010) as H. digonicus were also omitted. These accessions are highly divergent from all other Helicotylenchus spp. and are closer to Merlinius spp; accession MZ706973 identified as H. multicinctus from India is significantly divergent from H. multicinctus sequences, with biased substitution patterns and indels despite the 98% query coverage; ON117615 and ON117616 were omitted due to limited query coverage of below 30% when compared with other Helicotylenchus spp.; MW812365 and MN262451 identified by Abdulsalam et al. (2021) and Lamula (2020), respectively, as H. dihystera probably contain sequencing errors in form of multiple nucleotide insertion and deletions (indels), including multiple nucleotide substitutions and hence, when included in the analysis data, they generated a poor alignment that bore a substantial effect on phylogenetic reconstruction (data not shown).

The genus Helicotylenchus comprises several valid but cryptic species. The diagnostic characters of these closely related species are marred by phenotypic plasticity with no well-defined boundaries (Subbotin et al., 2011, 2015). Linking molecular markers to morphological data is ideal as it augments and provides a better understanding of the morphologically unresolved species within the genus. However, in most cases, molecular taxonomists tend to extract DNA from random unmeasured specimens taken from larger target population(s) putatively identified based on morphometrics of a few specimens to represent the species of interest. But it is common for two or more Helicotylenchus species to co-exist as mixed populations. For instance, in the current study, H. caudatus and H. asiaticus populations (PJ2) were recovered from the same subsample in a ratio of 1:6. Random selection of specimens from such populations may result in amplifying DNA from unidentified species and erroneously linking such DNA barcodes to morphological data of rather a morphometrically identified species. We envisage this to be one of the main causes of incorrect linkages of morphospecies and molecular barcodes that frequently appear in the public databases.

Case studies of sequence mislabeling in other groups such as Pratylenchus spp. have already been detailed, including assigning sequence data of bacterial feeding Acrobeloides to a quarantine species Pratylenchus goodeyi (Janssen et al., 2017). Acrobeloides and Pratylenchus commonly occur as mixed populations in nematode communities (Lü et al., 2020). Unfortunately, apart from the obvious differences in stoma and pharyngeal corpus, some species of Acrobeloides such as A. nanus generally present body postures that are similar to Pratylenchus, including morphometric characters such as general body length, position of vulva, and tail shape, among others. It is therefore possible to mistake Acrobeloides individuals for Pratylenchus in a mixed population, especially when observed at a low magnification under reflected illumination of a stereo microscope that tends to obscure the details of the internal organs. We therefore suggest that morphometric data of specimens to be used in DNA barcoding should be taken prior to DNA extraction to avoid mis-linking barcodes of unidentified populations to putative species. In addition, as already noted by Janssen et al. (2017), several other factors are evidently implicated in contributing to the unsuitable nature of several sequences in public databases. For instance, several sequences contain PCR amplification or sequencing errors, cloning vectors, chimera sequences, unexcised primer sequences, and bioinformatics analysis errors, among others. Therefore, to ensure correct integrative nematode diagnosis, it is crucial to mitigate these potential pitfalls during the various steps involved in DNA barcoding to avoid the consequential cascade of erroneous interpretations in diagnosis.