Reliability and Utility of Standard Gene Sequence Barcodes for the Identification and Differentiation of Cyst Nematodes of the Genus Heterodera

Analyses here are based on sequence data obtained from the NCBI public database, GenBank, up to 10 August 2021. We downloaded all available partial or complete sequences of the 18S rDNA, ITS, 28S rDNA, and coxI gene regions for all named species of Heterodera. We note that ITS is situated between the 18S and 28S genes and comprises the ITS1, 5.8S, and ITS2 genes; we considered a sequence as “ITS” if it included a partial fragment or complete sequence of one of the latter three genes. We searched for all valid species individually and also searched using junior synonyms. We created a database of these sequences where each line of data records the GenBank accession number, species, gene region, geographic collection information, and sequence author(s) or publication reference. In those cases where some of this information was not included with the sequence record on GenBank, we sought it from the publication listed and/or through web searches of scholarly literature aiming to identify if the GenBank accession number in question had been published or referenced. Where collection locality could not be determined and where possible, collection locality was inferred based on the institutional addresses of the author(s) listed in the respective sequence record. Base statistics were calculated in Microsoft Excel, and a map depicting the geographic spread of available sequences was generated using ggplot2 (Wickham, 2009) in R (https://www.R-project.org) and edited in Adobe Illustrator CS6.

Our analyses indicate that several of the most commonly employed molecular markers for the characterization of species of Heterodera lack the necessary resolution for distinguishing between some species, including a number of plant pests of significant global biosecurity concern. Among those for which molecular identification issues were detected, 13 species (H. avenae, H. carotae Jones, 1950a, H. ciceri Vovlas, Greco & Di Vito, 1985, H. cruciferae Franklin, 1945, H. daverti Wouts & Sturhan, 1978, H. elachista Ohshima, 1974, H. filipjevi (Madzhidov, 1981) Stelter, 1984, H. glycines Ichinohe, 1952, H. goettingiana Liebscher, 1892, H. hordecalis Andersson, 1975, H. oryzae Luc & Brizuela, 1961, H. schachtii and H. trifolii) are listed as regulated pests in at least one country (Singh et al., 2013). This has numerous implications for routine molecular diagnoses of species of Heterodera, such as the potential for confusing a major pest with a minor or non-pest species. For example, H. avenae is a major pest of cereals in temperate regions and causes significant annual yield losses throughout its range (e.g., Meagher, 1982; Smiley et al., 2005; Nicol and Rivoal, 2008), but our analyses showed that it could not be reliably distinguished from H. arenaria Cooper, 1955 or H. pratensis using the most commonly employed molecular marker, ITS. Both latter species parasitize non-crop grasses (Subbotin et al., 2018); thus, there is potential for confusing H. avenae with one of these non-pest species. Such cases of mistaken identity could result in misdiagnosed infestations and novel incursions. In addition, mixed infestations of some closely related species which parasitize similar crops, such as members of the Schachtii group, may go undetected.

Although a combination of morphological data and sequences from multiple gene regions is best for identification of species of Heterodera, our results suggest that where resources or expertise are limited, the coxI region should be favored over 18S, ITS, and 28S for basic diagnostic purposes. It is significant that a number of economically important species of Heterodera cannot be delineated using ITS; as to date this marker has been used more extensively than any other for both identification and phylogenetic Purposes (e.g., Subbotin et al., 2001, 2017; Tanha Maafi et al., 2003). Furthermore, more species of Heterodera have been characterized with ITS than any other marker, and for many species, this is the only gene region for which sequences are available. Thus, when developing molecular diagnostic tests for plant-parasitic nematodes like Heterodera spp., the trade-off between delineation power and species coverage needs careful consideration. There is currently a great amount of interest in employing high-throughput sequencing methods, such as metabarcoding and eDNA, for detection and monitoring of pest species (e.g., Abdelfattah et al., 2018; Valentin et al., 2018; Ruppert et al., 2019; Hardulak et al., 2020; Young et al., 2021). Many metabarcoding studies have utilized short fragments of nuclear genes such as 18S for species detection and identifications (e.g., Macheriotou et al., 2019; Ruppert et al., 2019; Giebner et al., 2020), but for Heterodera, fragments of the nuclear genes tested here lack the power to distinguish between all species. Those using such tools will need to be aware of the characteristics and shortcomings of the marker(s) employed and may need to follow-up species detections with other methods to confirm identifications.

The coxI gene seems robust for species delineation in the genus Heterodera; however, we did observe two instances in which this marker showed potential shortcomings. The first issue was detected in our coxI phylogenetic analyses of the Avenae group, in which sequences of H. arenaria formed a clade within the larger H. avenae species cluster. The sequences of H. arenaria differ by just two bp from those of H. avenae, well below the intraspecific variation observed for the latter species. However, these two bp appear to be unique changes, suggesting that H. arenaria can at least be distinguished from H. avenae from a barcoding, if not a phylogenetic, perspective. This is consistent with the findings of Subbotin et al. (2018) where H. arenaria and H. avenae were shown to be distinct in a haplotype network but formed a polyphyly in a tree derived from Bayesian analysis. Subbotin et al. (2018) speculated that H. arenaria represents a species recently diverged from a European population of H. avenae and remarked that, from a phytosanitary perspective, it is best to retain the species status of H. arenaria because it parasitizes coastal grasses of no economic importance, rather than cereals as in H. avenae.

The second instance in which coxI showed a possible shortcoming relates to the case of distinguishing H. carotae from H. cruciferae. These two species have previously been reported as indistinguishable using PCR-ITS-RFLP profiles and ITS gene sequences (Subbotin et al., 2000, 2001). In a study of H. carotae in Mexico, Escobar-Avila et al. (2018) produced novel coxI sequences of that species but also several coxI sequences for H. cruciferae, which were the only coxI sequences available for the latter species at the time of writing. Escobar-Avila et al. (2018) found that H. carotae and H. cruciferae could not be distinguished using coxI sequences and concluded that to differentiate these species an integrated approach including morphology and a test of host range was necessary. Although it is entirely possible that H. carotae and H. cruciferae are indistinguishable using coxI, considering the former species appears restricted to carrots as hosts (Jones, 1950; Winslow, 1954; Mugniery and Bossis, 1988) and the latter to brassicas and a few species of the Lamiaceae (Winslow, 1954; Baldwin and Mundo-Ocampo, 1991), it is somewhat surprising that these species would exhibit so little divergence in such a rapidly evolving gene. Escobar-Avila et al. (2018) provided three sequences of H. cruciferae, two from specimens from the USA and one from Russia, but did not provide a morphological account of these specimens, presumably because they were all consumed in molecular analyses. The Russian sequence of H. cruciferae differs from those from the USA by 15 bp, whereas the USA sequences of H. cruciferae differ from sequences of H. carotae from multiple countries by as little as 1 bp. Thus, the current intraspecific variation for H. cruciferae is greater than the difference between some sequences of H. carotae and H. cruciferae. Because these two species are fairly similar morphologically and could easily co-occur in mixed or rotated vegetable fields, it is possible that some of the putative specimens of H. cruciferae used by Escobar-Avila et al. (2018) were misidentified. At present, there are simply too few coxI sequences of H. cruciferae to be certain of the utility of this marker, or lack thereof, for the H. carotae–H. cruciferae species pair. Heterodera carotae is an important pest of carrots throughout its range (Greco et al., 1994; Esquibet et al., 2020), and although H. cruciferae is largely considered a minor pest, there have been a few reports of significant damage to crops infested by this species (Lear et al., 1965; Sykes and Winfield, 1966). It is critical then that a suitable molecular marker be identified which can reliably distinguish between these two species. Using a collection of microsatellite markers, Esquibet et al. (2020) demonstrated a high level of genetic divergence between H. carotae and H. cruciferae and remarked that microsatellites could be used to develop a diagnostic test for these species. If further study confirms that the coxI gene cannot reliably distinguish between H. carotae and H. cruciferae, other molecular methods such as microsatellites (e.g., Gautier et al., 2019; Esquibet et al., 2020) can fill this diagnostic gap.

Although the number of obviously erroneous sequences detected was low, for several species of Heterodera, we observed high levels of intraspecific variation within one or more of the molecular markers evaluated, primarily ITS. Unsurprisingly, this was most apparent in those species/marker combinations with a large pool of sequences, such as H. avenae, H. glycines, and H. latipons Franklin, 1969, each of which had over 100 ITS sequences available. In these species, we observed intraspecific maxima of 80–122 bp (12%–13%) in the ITS gene. Additionally, some species had large intraspecific variation despite lower sample size, such as H. goettingiana with 15 ITS sequences and an intraspecific maximum of 81 bp (11%) and H. hordecalis with 37 ITS sequences and an intraspecific maximum of 68 bp (9%). None of these intraspecific maxima are unbelievable when considering that these datasets include sequences of individuals sourced from many geographically disparate populations (e.g., Subbotin et al., 2003, 2018). However, these large levels of intraspecific variation do clash with the patterns observed for most other members of the genus and notably in several other species with medium to large sequence pools, such as H. schachtii and H. filipjevi, both of which exhibit intraspecific maxima of less than 5% in ITS. We think it is likely that for several species the overall level of intraspecific variation observed for some molecular markers is inflated. In our assessment of sequence reliability, in many instances, we could not be certain if the intraspecific variation observed was due to real population-level genetic variability or artifacts of poor-quality sequencing results and/or sequence editing errors. We flagged sequences that diverged greatly from their congeners as problematic; however, in larger sequence pools, relatively small levels of variation between individual sequences become amplified, resulting in the very large levels of overall intraspecific variation detected in pairwise comparisons. This is an important limitation of our analyses and of the overall pool of publicly available sequences.

A related issue is species for which only a few sequences are available, but one or more of those sequences diverge significantly from the others. Again, it can be difficult to determine if such divergence represents natural variation or an artifact such as misidentification or sequence editing errors. Furthermore, the large number of sequences submitted to GenBank that are not associated with a published manuscript leaves many issues related to sequence identity ambiguous as there is no account of the morphology of specimens utilized. A good example of both of the above relates to the 28S sequences available for H. cruciferae. Sasanelli et al. (2013) performed a thorough study of H. cruciferae in Italian cabbage, identified specimens using morphology, and characterized them with ITS and 28S rDNA gene sequences. However, BLAST results of the 28S sequence of H. cruciferae from that study suggests that it is probably representative of H. schachtii, another cyst nematode widespread on brassicas (Subbotin et al., 2010b). There is one other 28S sequence of H. cruciferae available, but it is not associated with a published manuscript (KP114546; Jabbari et al., unpublished); this sequence is very close to H. carotae but not identical. This suggests that this “unpublished” sequence is accurate, but without a morphological account or other sequences with which to compare, it is still uncertain. This results in a situation where, although two 28S sequences are available for H. cruciferae, we still cannot be confident that either are truly representative of that species. Thirty species of Heterodera lack molecular data and thus cannot be diagnosed using barcoding methodologies. Of the 57 species of Heterodera for which molecular data are available, eight are represented by just a single sequence and over half of the remainder have just one sequence for one or more of the markers evaluated here. Where only a single sequence of a particular marker is available, it is best treated with caution when used for diagnostic purposes as additional sequences, ideally from independent studies, are needed for confidence that such sequences are truly representative of the species they are purported to be of.

There is little doubt that erroneous sequences we observed on GenBank were uploaded in good faith. However, to avoid errors, it is important that authors carefully compare their novel sequences with those available in public databases via BLAST or other phylogenetic methods before upload. Original sequencing results should be scrutinized for base-call quality, and poor-quality sequences should be discarded. Where possible, multiple sequence replicates of each gene should be produced and compared prior to upload to ensure base calls are consistent. It is also advisable to sequence multiple gene regions from individual isolates to ensure that molecular based identification is consistent across markers. We also encourage authors who have uploaded sequences which later prove to be mislabeled or problematic to correct them.

High-throughput sequencing techniques utilizing new and more accurate markers are already superseding other molecular-based diagnostic methods for cyst nematodes in many laboratories (e.g., Gautier et al., 2019; Esquibet et al., 2020). Despite this, we foresee the four markers evaluated here remaining in use for identification of plant parasitic nematodes for many years to come. With that in mind, of the four markers evaluated, coxI shows the greatest utility for identification and delineation of species of Heterodera. However, the cox1 gene might not be best for phylogenetic interpretation (Subbotin et al., 2015, 2018); so, going forward, a combination of coxI and ITS, plus 28S where possible, seems ideal, especially as coxI data are presently available for only around a third of species of Heterodera. There is a real possibility of hybridization between species of Heterodera (Potter and Fox, 1965); thus, a combination of a mitochondrial and nuclear marker is recommended for areas where the range of multiple species overlap. Molecular diagnoses should be based not only on multiple molecular markers used in concert but also on a combination of sequence matching and tree-based phylogenetic methods. Lastly, it is important to keep in mind that no one technique is likely to be a panacea for the identification of species of Heterodera. New species are being described at a fairly steady rate (e.g., Li et al., 2020; Phougeishangbam et al., 2020; Jiang et al., 2022), and sequence data continue to accumulate and provide further insight into the population genetics, host associations, and phylogeography of Heterodera (Subbotin et al., 2018; Esquibet et al., 2020; Oro and Tabakovic, 2020). It is telling that so few sequences have been generated from countries in South America and Sub-Saharan Africa as this seems an obvious result of lack of resources, rather than lack of Heterodera species richness. Thus, there is still a great need for traditionally trained taxonomists and diagnosticians that can employ a range of techniques to identify these problematic nematodes.