High genetic diversity and geographic subdivision of three lance nematode species (Hoplolaimus spp.) in the United States

Lance nematodes (Hoplolaimus spp.) feed on the roots of a wide range of plants, some of which are agronomic crops. Morphometric values of amphimictic lance nematode species overlap considerably, and useful morphological characters for their discrimination require high magnification and significant diagnostic time. Given their morphological similarity, these Hoplolaimus species provide an interesting model to investigate hidden diversity in crop agroecosystems. In this scenario, H. galeatus may have been over-reported and the related species that are morphologically similar could be more widespread in the United States that has been recognized thus far. The main objectives of this study were to delimit Hoplolaimus galeatus and morphologically similar species using morphology, phylogeny, and a barcoding approach, and to estimate the genetic diversity and population structure of the species found. Molecular analyses were performed using sequences of the cytochrome c oxidase subunit 1 (Cox1) and the internal transcribed spacer (ITS1) on 23 populations. Four morphospecies were identified: H. galeatus, H. magnistylus, H. concaudajuvencus, and H. stephanus, along with a currently undescribed species. Pronounced genetic structure correlated with geographic origin was found for all species, except for H. galeatus. Hoplolaimus galeatus also exhibited low genetic diversity and the shortest genetic distances among populations. In contrast, H. stephanus, the species with the fewest reports from agricultural soils, was the most common and diverse species found. Results of this project may lead to better delimitation of lance nematode species in the United States by contributing to the understanding the diversity within this group.

Of the twenty-nine species described in the genus (Handoo and Golden 1992), H. galeatus is the most commonly reported in the United States (Lewis and Fassuliotis 1982). Morphologically, Hoplolaimus galeatus belongs to a group of species that Fortuner (1991) called "ancestral", characterized by three esophageal gland nuclei, four incisures in the lateral field, excretory pore posterior to the hemizonid, and the presence of abundant males. Siddiqi (2000) classifies this group as the subgenus Hoplolaimus. This "ancestral" subgenus includes H. magnistylus, H. concaudajuvencus, and H. stephanus, among a few other species. Hoplolaimus magnistylus and H. concaudajuvencus differ from H. galeatus by the possession of a longer stylet (Table 1, Fig. 1). Hoplolaimus concaudajuvencus has more definitely tulip-shaped stylet knobs and a second-stage juvenile with conically pointed tail. Hoplolaimus stephanus can be distinguished from H. galeatus by the 24-28 longitudinal striations on the basal annule of the lip region compared to 32-36 in H. galeatus, shorter spicules, less areolation of the lateral field, and shorter body length (Sher 1963) (Table 1, Fig. 1). The "derived" species, which constitute the subgenera Basirolaimus and Ethiolaimus and include species such as H. columbus and H. seinhorsti, possess six esophageal gland nuclei and fewer than four incisures in the lateral field, the excretory pore is anterior to the hemizonid, and the mode of reproduction is mostly parthenogenetic, with males rare or absent (Fortuner 1991).
Morphometric values of the "ancestral" species overlap considerably (Table 1, Fig. 1) (Sher 1963;Vovlas et al. 1991;Handoo and Golden 1992), and useful morphological characters for their discrimination require high magnification and significant diagnostic time. Given their morphological similarity, these Hoplolaimus species within the "ancestral" clade provide an interesting model to investigate hidden diversity in crop agroecosystems, considering the wide host range and distribution of H. galeatus in the United States (Fortuner 1991). In this scenario, H. galeatus may have been over-reported and the related species within the subgenus Hoplolaimus (Siddiqi 2000) that are morphologically similar could be more widespread in the United States that has been recognized thus far. The main objectives of this study were (1) to delimit Hoplolaimus galeatus and morphologically similar species using morphology, phylogeny, and a barcoding approach and (2) to estimate the genetic diversity and population structure of the species found. With this study, we expect to understand the diversity within this group and contribute to the elucidation of the delimitation of species of lance nematodes in the United States.

Materials and Methods
Nematode sampling, identification, and DNA isolation Nematode populations were obtained from soil samples collected in 2011-2013 in agricultural fields, golf courses, and lawns of different regions in the United States (Table 2). Nematodes were extracted from soil by sugar centrifugal flotation (Jenkins 1964). Lance nematode specimens from each soil sample were morphologically identified using the key by Handoo and Golden (1992) and the original descriptions of the species (Sher 1963;Golden and Minton 1970;Robbins 1982). DNA from individual nematodes was extracted using the Sigma Extract-N-Amp kit (XNAT2) (Sigma, St. Louis, MO) as reported by Ma et al. (2011), and DNA was stored at À20°C until use.

Sequence alignment
Contigs were assembled in Sequencher 5.1 (Genes code corp., Ann Arbor, MI). All sequences were checked and edited manually, and chromatograms were inspected to confirm base calling and to identify recombination sites. Consensus DNA sequences were then aligned using Clu-stalW (Thompson et al. 1997) including three out-group taxa: H. columbus sequences obtained from this study and sequences of Rotylenchus robustus (JX015440) and Rotylenchus paravitis (JX015415) from GenBank. The original alignment for ITS consisted of 1050 bp, but several indels were detected. Therefore, divergent and ambiguously aligned positions were removed and conserved blocks selected using the software Gblocks v0.91b (Castresana 2000) with default values. The resulting dataset comprised 550 bp of the ITS1 portion of the gene. For the mitochondrial region, the alignment consisted of 347 bp. The   new generated haplotypes for both genes were deposited in GenBank (Table 1). For the COI marker, the aligned sequences were well defined and chromatograms had no double peaks, ambiguous positions or indels. Nonetheless, we tested for the occurrence of stop codons that could denote the presence of nuclear copies of mitochondrial-derived genes (numts) or COI pseudogenes (Zhang and Hewitt 1996;Song et al. 2008;Moulton et al. 2010). Numts are copies of mitochondrial genes moved to the nuclear genome that become nonfunctional and noncoding. Consequently, these numts can confuse phylogenetic analyses (Song et al. 2008;Moulton et al. 2010;Baeza and Fuentes 2013). To check for the presence of numts, we followed Song et al. (2008) and did a basic local alignment search (blast) of all COI sequences in NCBI (National Center for Biotechnology Information) against the database nucleotide collection (nr/nt) and optimized for highly similar sequences to include only haplotypes that showed E-values ≥ 1.0e-45 and similarity ≥ 90% with plant-parasitic nematodes. All retrieved sequences were of plant-parasitic nematodes most commonly of the genera Rotylenchus and Scutellonema followed by Heterodera, Punctodera, and Meloidogyne. After this, the COI haplotypes were translated using the invertebrate mitochondrial code in Mega v.5 (Tamura et al. 2011) to verify the protein coding frameshifts and nonsense codons for each of the six putative reading frames in DNAsp (Librado and Rozas 2009).

Phylogenetic analysis
For phylogenetic analysis, the most appropriate evolutionary model was selected for each gene dataset using the Akaike information criterion (AIC) in the software Modeltest v3.7 (Posada and Crandall 1998). For both genes, the best-fit model was GTR with invgamma-shaped rate variation ( Phylogenetic relationships were constructed for each gene separately using maximum-likelihood (ML) and Bayesian inference (BI). Maximum-likelihood analysis was performed in Treefinder (Gangolf et al. 2004) using the default parameters. Branch support was based on 1000 bootstrap pseudoreplicates (Felsenstein 1985), and clades were considered as well/strongly supported when bootstrap was >70%. Bayesian inference was implemented in the software MrBayes 3.1.2 (Ronquist and Huelsenbeck 2003). The analysis was conducted for 6 million generations, and trees were sampled for every 100th generation from the Markov Monte Carlo chain (MCMC) analysis. A burn-in period was set to discard the first 1250 trees with nodal support defined as posterior probabilities, and clades were considered strongly supported when values were > 0.95 (Alfaro et al. 2003). Additionally, for the COI sequences, a neighbor-joining (NJ) tree was constructed in MEGA v.5.0 (Tamura et al. 2011) using the Kimura two-parameter (K2P) model under the default settings. From this dataset, a pairwise distance matrix among haplotypes was generated to calculate intra-and interspecific genetic divergence among haplotypes, as it is typically performed for DNA barcoding studies (Hebert et al. 2004).

Genetic diversity
Genetic diversity analyses were based on mitochondrial DNA data because COI gene trees provided better resolu-tion than ITS1. Mitochondrial haplotypes networks were constructed for each species using median-joining (MJ) networks in Network 4.5.1.0 (http://www.fluxus-engineering.com/sharenet.htm). Analysis of molecular variance (AMOVA) was conducted in Arlequin v. 3.11 (Excoffier et al. 2005) to infer genetic structure within species. Values of Fst were also calculated in Arlequin v. 3.11 to estimate genetic differentiation among clades and tested for significance by permuting haplotypes between species/ populations (10,000 replicates).

Phylogenetic relationships
Phylogenetic reconstruction using the mtDNA sequences with ML and BI yielded trees with highly similar topologies and revealed five major clades (Fig. 2) Within clades, phylogenetic analyses using the three methods (ML, BI, and NJ) for the COI region showed that Hoplolaimus individuals were structured according to locality (Figs. 2, 3). In the H. stephanus clade, for example, COI sequences were subdivided into two major lineages: one composed mainly of different populations from NC and a sister clade mostly comprised of individuals from     NE. The other H. stephanus subclade included three geographically clear lineages: individuals from IA that formed a single lineage closely related to the KS population and a third distinct lineage corresponding to populations from OH. Similar topology was revealed within the H. magnistylus clade, with a single clear subclade formed by two populations from IL and another lineage mostly composed of TN populations. Specimens of H. concaudajuvencus were only collected on grasses from Texas. The H. galeatus clade was the only one that did not show geographic structure, with two major subclades with broadly overlapping geographic ranges, which included populations from FL, SC, and AL.
For the nuclear marker, the ITS1 trees using ML and BI methods identified similar clades as for the COI sequences, with the difference that H. concaudajuvencus individuals were included in the H. magnistylus clade (Fig. 3), and trees estimated from nuclear data gave less resolution, with little or no structure within clades. The only clade that showed fragmentation with ITS1 was H. stephanus, in which IA, KS, and NE populations formed one subclade, and OH and NC populations appeared as a closely related group (Fig. 3). In general, the BI tree for ITS1 was supported by strong posterior probabilities (0.9-1), whereas the ML tree had bootstrap values from 35 to 100%. Another difference was that sequences from H. sp. 1 appeared as a sister clade of H. galeatus although with low bootstrap support (37%) (Fig. 3).
Genetic distances (K2P) for COI between Hoplolaimus groups (clades) ranged from 11.65% to 23.21% with a mean of 16.64%. The maximum interspecific distances (22.78% and 23.21%) were found between H. galeatus and H. columbus populations, and the minimum values (11.65% and 11.99%) were found between H. stephanus and H. magnistylus populations from NC and IL, and OH and IL. Within groups (clades), genetic divergence varied from 0 to 12.83% with an average of 6.26%. Therefore, a clear overlap between intra-and interspecific genetic distances was detected. The deepest levels of intragroup divergence that caused this overlap mainly occurred between populations of the H. stephanus clade: OH and NE (10.93-12.83%), NC and IA (10.55-12.83%), NC and KS (11.31-12.44%), and NE and KS (11.29-12.44%). Additional high levels of intraspecific genetic divergence (10.91-12.04%) occurred in the H. magnistylus clade between some populations from IL and TN. In general, the H. galeatus clade showed lower levels of intraspecific divergence ranging from 0 to 7.01%.

COI mtDNA population structure
Haplotype networks for COI were congruent with the phylogenetic reconstructions, showing segregation of the mitochondrial haplotypes according to locality (Fig. 4). Sixteen haplotypes were identified in the H. stephanus clade, distributed in three main geographic groups: one mostly composed of populations from NC and NE, which share two haplotypes among them; the second with unique haplotypes from OH populations; and the third comprised of two haplotypes from IA and one from KS populations. In the H. magnistylus network, the 12 haplotypes identified were divided into two groups: one composed of haplotypes from populations collected in IL, and the second of haplotypes from TN populations. For H. concaudajuvencus and H. galeatus, fewer haplotypes were identified (4 and 7, respectively). For H. concaudajuvencus, all were unique haplotypes from TX populations. Hoplolaimus galeatus was the only species that did not show correlation with geographic regions, having a single haplotype shared by several individuals from SC, AL, and FL, with the remaining haplotypes being from populations collected in SC.
Frequencies and distribution patterns of COI haplotypes showed that H. stephanus was the most widespread and with the most diverse lineage, distributed from NE, KS, IA, and OH to NC (Fig. 4). Hoplolaimus galeatus was the least diverse species with a predominant haplotype present in populations from AL, one from FL and one from SC. Haplotypes for the remaining species were predominantly location specific, conforming haplogroups concordant with sampling locations. The highest nuclear diversity was observed in H. stephanus (p= 0.068) and H. magnistylus (p=0.052). Haplotype diversity (Hd) was similar for H. stephanus and H. magnistylus (0.908 and 0.909, respectively) ( Table 2 and Fig. 4).
Results from the 112 haplotypes analyzed with AMOVA suggest strong genetic differentiation in all species (clades) as explained by the high F ST value (0.96%, P < 0.0001) ( Table 3). AMOVA also revealed that 51.11% of the mtDNA variation was attributed to differences among the Hoplolaimus clades, followed by the variation among populations within each clade (43.07%) and very low variation within populations (5.82%). However, higher variation (76.79%) was recovered when AMOVA was performed across geographic regions, suggesting strong genetic differentiation among sampling locations. With regard to each morphospecies, the strongest genetic subdivision was found within the H. stephanus populations

2939
(85.11%), whereas populations of H. galeatus and H. magnistylus were less differentiated. The strong genetic differentiation in H. stephanus and the high intraspecific variability caused overlap between intra-and interspecific K2P genetic distances. Therefore, AMOVA was performed for H. stephanus with different combinations of populations to determine which grouping better explained the variability among groups. Results from the AMOVA showed the highest variation (83.10%) when populations of H. stephanus were partitioned as NC, NE, OH, and KS+IA. Based on this combination of populations, new K2P genetic distances were estimated considering these populations as different genetic entities (interspecific variability), and a barcoding gap was detected for this species.

Discussion
Molecular studies on Hoplolaimus species have resolved some phylogenetic relationships within the genus using the internal transcribed spacer (ITS1) (Bae et al. 2008), the 28S ribosomal DNA (Bae et al. 2009), and the actin gene (Ma et al. 2011). Mitochondrial markers had not been used before for molecular analysis of this genus. The protein coding cytochrome oxidase 1 (COI) gene has been proposed for species delimitation in several taxa by comparing pairwise genetic divergence among individuals of the same species versus individuals of different species (Hebert et al. 2003). Although the reliability of genetic distances for species delimitation has been subject of criticism (De Ley et al. 2005;DeSalle et al. 2005), COI barcoding can be a useful tool that complements morphological and molecular identification methods in an integrated taxonomy approach (Ferri et al. 2009).
In the COI haplotype distribution map, Hoplolaimus species and populations appear allopatrically distributed (geographically isolated). Genetic divergence between allopatric populations depends on the timing and level of gene flow between them (Mayr 1963;Nosil 2008). Geographic barriers seem to be the main source of variation for nematodes and soil-dwelling arthropods due to the low dispersal capacity of nematodes in soils, causing high levels of population fragmentation, which eventually lead to diversification of the species (Picard et al. 2004). Low gene flow among populations may have influenced the deep structure detected in the COI sequences. However, several authors have detected high gene flow over large distances in plant-parasitic nematodes such as H. schachtii, G. pallida, and Bursaphelenchus mucronatus (Picard et al. 2004;Plantard and Porte 2004;Pereira et al. 2013). For C. elegans, a free-living soil nematode, long-range dispersal was also suggested (Koch et al. 2000). For H. schachtii, G. pallida, and B. mucronatus, the high gene flow was attributed to transport of soil by farm machinery, sewage farms around sugar factories, by water through irrigation, flood or drainage, and /or by wind. For H. galeatus, we hypothesize that the application of nematicides could be responsible for the reduction in genetic diversity within populations, considering that this species was identified mainly on golf courses that require high maintenance and nematicide applications.
Hoplolaimus stephanus showed deep intraspecific variation that might be an indication of cryptic speciation. Strong population structure and genetic differentiation interpreted as cryptic speciation has been identified for marine and free-living nematodes using the same set of mitochondrial primers used in this study (Derycke et al. 2005;Ristau et al., 2013). The mitochondrial data evidenced the presence of three main groups in the H. stephanus haplotype network (NC + NE, OH, KS + IA). These groups also showed large K2P genetic distances that overlapped with the interspecific genetic distances, and AMOVA was better explained when the populations were partitioned in the same groups. However, the This study is the first to use molecular phylogenetic analysis based on mitochondrial data for species within the genus Hoplolaimus. The COI gene resolved species in highly supported clades and showed deep variation within species, mainly correlated with geographic origin of the populations, indicating that this portion of the COI gene is adequate to define species and to study the genetic structure. Furthermore, we observed no heteroplasmy (presence of more than one mitochondrial genome in an individual) as evidenced by clean chromatograms with the absence of double peaks. Heteroplasmy has been reported in plant-parasitic nematodes like M. chitwoodi (Humphreys-Pereira and Elling 2013). However, Kiewnick et al. (2014) studied the feasibility of mitochondrial markers COI and COII on several Meloidogyne populations and did not detect heteroplasmy.
With ITS sequences, high genetic structure was only detected for H. stephanus. Previous studies on the phylogenetic relationships of Hoplolaimus species (Bae et al. 2008(Bae et al. , 2009 found that the intron ITS1 showed higher variability among species within this genus. Our data showed that ITS1 provided poor resolution within clades, showing low variability for the majority of lineages. Other studies on the genetic diversity of plant-parasitic nematodes such as reniform nematode also detected a lack of genetic differentiation in the ITS region (Agudelo et al. 2005) and recently, Fu et al. (2011) for species discrimination within the Heterodera genus. This confirms the need to use different markers to infer phylogenetic relationships and for genetic population structure analysis. Incomplete lineage sorting could be a possible explanation for the high level of mitochondrial differentiation detected compared to the nuclear sequences (Gebiola et al. 2012), which is an expected consequence of populations undergoing speciation. The end result within a phylogeny is a subset of characteristics that have discordance within the species tree (Derycke et al. 2005).
It is possible that the inconsistency between nuclear and mitochondrial differentiation is attributable to a higher rate of molecular evolution for the latter. Because haploidy and uniparental inheritance reduce the effective number of mitochondrial genes to about one quarter of that of nuclear (Piganeau and Eyre-Walker, 2009), divergence can be accomplished earlier than for nuclear genes after a recent speciation event. Alternatively, the failure to detect nuclear DNA variation could be an artifact of the use of an inappropriate marker, although ITS rRNA markers have proven useful for species diagnosis in many nematode taxa (De Ley et al. 2005). The primers used in this study were also tested by Toumi et al. (2013) for discriminating several Heterodera and Punctodera species with successful results, which suggest that this portion of the COI gene can be a suitable marker for plant-parasitic nematode species. A nuclear marker that could be considered in future analysis is elongation factor-1a (EF-1a), a conserved nuclear coding gene that can be used to investigate recent divergences due to the presence of rapidly evolving introns (Kawakita et al. 2003). An ideal marker to confirm the lack of congruence between nuclear and mitochondrial differentiation is microsatellites, which are more sensitive in detecting recent or ongoing speciation events (Michel et al. 2010).
Morphological and molecular data using nuclear and mitochondrial genes revealed the presence of four recognized Hoplolaimus morphospecies in the soil samples collected: H. galeatus, H. magnistylus, H. concaudajuvencus, and H. stephanus. A fifth species (H. sp. 1), currently undescribed, was collected on maple. This is the first report for H. concaudajuvencus in TX on bentgrass and for H. magnistylus on corn in IL. Ours is also the first report for H. stephanus in OH and KSon bentgrass, and in IO and NE on corn. Our group had recently published a first report of H. magnistylus in TN on soybean, corn, and cotton (Donald et al. 2013).
To our surprise, H. galeatus, the species most frequently reported in crops in the United States, was only identified in SC, FL, and AL on turfgrasses. This species also exhibited low genetic diversity (low haplotype and nucleotide diversity) and showed the shortest genetic distances among its populations. In contrast, H. stephanus, the species with the fewest reports from agricultural soils, was the most common and diverse species found in this study. Sher (1963) described H. stephanus from specimens collected in Nicols, SC, and in the description included mention of an additional population from New Jersey. Nearly three decades later, Vovlas et al. (1991) published morphological observations from individuals collected in Raleigh, NC, from an unspecified host. Two decades later, Ma et al. (2011) found three populations of H. stephanus, one in Pennsylvania and two in SC from which they published nuclear sequences and designed species-specific primers for molecular diagnosis. The population from Pennsylvania (Ma et al. 2011) constitutes the only report of this species on a grass (Poa pratensis) host. These three publications total the extent of the research available for H. stephanus. Considering it appears to be the most widely distributed lance nematode species in agricultural soils in the United States, further work is needed to study its biology, ecology, pathogenicity, and economic thresholds.