A COI DNA barcoding survey of Pratylenchus species in the Great Plains Region of North America

Soil and root samples analyzed in this study were obtained from USDA Cooperative Agricultural Pest Survey Program (CAPS), Wheat and Corn Disease Surveys in Kansas and Nebraska, and field samples submitted to the University of Nebraska-Lincoln Disease Diagnostic Clinic, representing different crops and geographic regions primarily within the Northern Great Plains of North America. Five statewide surveys associated with the CAPS program were conducted by Departments of Agriculture in Kansas, Montana, Nebraska, North Dakota, and Wyoming. A minimal number of samples were acquired from South Dakota and no samples were collected from Iowa. Sample collection sites (county locations), nematode identification (NID) numbers, and host information for each analyzed specimen are presented in Supplementary Table 1 and Figure 1.

Figure 1:

Nematodes were extracted from 100 cm³ field soil using a modified flotation-sieving and sugar centrifugation method (Jenkins, 1964) and from the host-root material using root incubation (Russell, 1987; Todd and Oakley, 1996). A majority of the samples recovered from Kansas was obtained from root extracts.

Nematodes extracted from soil and roots were first evaluated under a stereo dissecting microscope and select specimens belonging to the genus Pratylenchus were handpicked for light microscopy examination and DNA extraction. Following immobilization of live specimens by heating, individual nematode specimens were mounted on temporary glass slides, measured, and photographed with a Leica DMLB light microscope with Differential Interference Contrast and a Leica DC300 video camera. Each nematode was assigned a unique Nematode Identification number (NID). The NID number links images, measurements and DNA of each individual nematode specimen. Usually five individual specimens from each sample were mounted, measured, and photographed in this manner. Images were archived in the database system of the Nematology Laboratory at the University of Nebraska-Lincoln. For PCR amplification, image-vouchered specimens were removed from temporary slides and smashed in an 18 µl drop of sterile deionized distilled water using a sterile, transparent micropipette tip. Smashed specimens were transferred to PCR reaction microfuge tubes and stored at −20°C until PCR amplification.

The mitochondrial COI gene was amplified by PCR using primer sets of COI-F7b-Prat (F7bP) (5′-GGDTGRACWTTHTAYCCNCC-3′) developed by the UNL nematology laboratory, and COI-JB5 (5′-AGCACCTAAACTTAAAACATAATGAAAATG-3′) (Derycke et al., 2005), resulting in a fragment length of 727 to 739 bp for genetic analysis after trimming the primers from the amplified product. On occasion, forward primer JB3 (5′-TTTTTTGGGCATCCTGAGGTTTAT-3′) (Derycke et al., 2005) was used in combination with the JB5 primer. The D2-D3 domains of 28S rDNA were amplified using the primer set D2A (5′-ACAAGTACCGTGAGGGAAAGTTG-3′) and D3B (5′-TCGGAAGGAACCAGCTACTA-3′) (De Ley et al. 1999). PCR was conducted in a total volume of 30 µl reaction mix consisting of 1.2 µl molecular biology grade water (Sigma-Aldrich), 2.4 µl of each primer, and 15 µl of 2× JumpStart REDTaq ReadyMix Reaction Mix (Sigma-Aldrich), and 9 µl DNA template. Amplification conditions were as follows: initial denaturation (modified hot start) at 94°C for 5 min followed by 45 cycles of denaturation at 94°C for 30 sec, annealing at 50°C for 30 sec, extension at 72°C for 90 sec, and a final extension at 72°C for 5 min. Annealing temperature was at 48°C for the amplification of D2–D3 domains. Amplification success was evaluated in 1% agarose gels using 0.5X TBE and ethidium bromide.

High-quality PCR products in a 0.7% TAE agarose gel were extracted by x-tracta Tool (USA Scientific), and cleaned using Gel/PCR DNA Fragments Extraction Kit (IBI Scientific). Cleaned PCR products were shipped to the University of California-Davis DNA Sequencing Facility for sequencing in both directions.

DNA sequences were edited using CodonCode Aligner version 8.0.1 (www.codoncode.com) and used in a BLAST search. Sequences were submitted to GenBank, with accession numbers for COI and 28S sequences provided in Supplementary Table 1. Multiple sequence alignment was conducted using MUSCLE (Edgar, 2004) with a gap opening penalty −400 and a gap extension penalty −200, in MEGA version 7 (Kumar et al., 2016). The best DNA Model tool in MEGA 7 was used to determine a best-fit substitution model: general time reversible (GTR) with a gamma distribution (G) and proportion of invariable sites (I). The substitution model (GTR + G + I) was used for maximum-likelihood analysis with 200 bootstrap replications using software MEGA 7 and for Bayesian inference analysis (BI) using the software MrBayes 3.1.2 (Huelsenbeck and Ronquist, 2001). The COI haplotype data set was reduced by removal of redundant sequences using software Jalview.2.10 (Waterhouse et al., 2009).

Molecular species delimitation was assessed using Automatic Barcoding Gap Discovery (ABGD) (Puillandre et al., 2012), the Generalized Mixed Yule Coalescent (GMYC) method (Pons et al., 2006), and by statistical parsimony networks (TCS) (Clement et al., 2002), all applied to a non-redundant, 143 specimen COI data set.

ABGD analyses were performed at the online webserver (wwwabi.snv.jussieu.fr/public/abgd/abgdweb.html) using the Kimura (K80) distance model with default parameters.

GMYC is a model-based maximum likelihood approach that uses an ultrametric tree to delimit species and determine diversification and coalescence events based on branching patterns. An ultrametric tree was constructed using Beast 2.5.2 (Bouckaert et al., 2014) with the GTR model with Gamma distribution (G) and four gamma categories for nucleotide evolution. The tree was estimated with a strict clock model or a lognormal relaxed clock model assigned to a yule or a coalescence branching model. All other parameters were set as defaults and the tree was not dated for GMYC analysis. Markov Chain Monte Carlo (MCMC) chain was run at 30 to 250 M generations and sampled every 3,000th to 5,000th generation based on model selection. Tracer V1.7 (Rambaut et al., 2018) was used to visualize convergence and evaluate effective sample size (ESS > 200) of traces. A maximum clade credibility tree was produced by TreeAnnotator V2.5.2 (Rambaut et al., 2018) with common ancestor node heights, after discarding 25% of the trees as burn-in. The tree was visualized using FigTree V.1.4 (tree.bio.ed.ac.uk), and Beast analyses were run on the XSEDE server of the CIPRES Science Gateway (Miller et al., 2011). The GMYC analysis was conducted using the single threshold option (Fujisawa and Barraclough, 2013) using the Species’ limits by threshold statistics (SPLITS) R Package Version 1.0–19  (http://R-Forge.Rproject.org/projects/splits/) and APE package (Paradis et al., 2004) available for R v3.5.2 (www.R-project.org). After analysis all outputs (estimated threshold time and the list of ML clusters and entities) were exported from R.

TCS is a clustering method that calculates a distance matrix for pairwise comparison of haplotypes to recognize species boundaries, while calculating mutational differences at an assigned cutoff probability. Parsimony criterion applies before the mutational differences reach the cutoff percentage (Templeton et al., 1992). In this study, we evaluated two connection limits (90 and 95%) for species delimitation.

A haplotype network analysis was used to visualize the relationships among COI haplotypes of Pratylenchus species based on their geographic and host information. Haplotype networks were calculated using TCS plug-in (Templeton et al., 1992; Clement et al., 2002) in the software PopART 1.7 (Leigh and Bryant, 2015) using a 95% connection limit.

Molecular divergence time was estimated with a molecular clock analysis using the non-redundant COI data set in BEAST v2.5.2 (Bayesian evolutionary analysis by sampling trees) (Bouckaert et al., 2014). Because no fossil or geological calibration points are available for dating the Pratylenchus phylogeny, we used two mutation/clock rates from the literature. The first was a widely used substitutions/site/my/lineage clock rate of 0.0115 corresponding to a common invertebrate COI mitochondrial pairwise sequence divergence rate of 2.3% per site per million years (Brower, 1994). The second was the mitochondrial substitution genome rate of 7.2 × 10⁻⁸ per site per generation experimentally calculated for the nematode Caenorhabditis briggsae (Dougherty & Nigon, 1949) Dougherty, 1953 (Howe et al., 2010). For Pratylenchus, an assumption of two generations per year was applied to the analyses. The BEAST analysis was carried out under a fixed strict clock model with a yule speciation model as tree prior, all other parameters were set as defaults. Markov Chain Monte Carlo (MCMC) chain was run at 30 M generations and sampled every 3,000th generation. Tracer V1.7 (Rambaut et al., 2018) was used to visualize convergence and evaluate effective sample size (ESS > 200) of traces. A maximum clade credibility tree was produced by TreeAnnotator V2.5.2 (Rambaut et al., 2018) with common ancestor node heights, after discarding 25% of the trees as burn-in. The MCMC tree and node ages were visualized using FigTree V.1.4 (tree.bio.ed.ac.uk). Beast analysis was conducted on the XSEDE server of the CIPRES Science Gateway (Miller et al., 2011).

A total of 860 soil samples were assayed during the growing seasons of 2017 and 2018. These samples represented statewide surveys of seven major agronomic crops (wheat, corn, dry beans, barley, alfalfa, sugar beets, and potatoes) in Colorado, Kansas, Montana, Nebraska, North Dakota, and Wyoming, as well as four additional crops (soybean, cotton, apple, and vineyard) opportunistically sampled across the Great Plains Region. In total, Pratylenchus species were recovered from approximately 71% of all samples. Prevalence of Pratylenchus spp. were 89.30% for corn fields (86.40% in Nebraska, 63.60% in Montana, 94.90% in Kansas, and 87.50% in Colorado), 53.15% of wheat fields (100% in Nebraska, 20% in Montana, 37% in North Dakota, 80% in Wyoming, and 80.30% in Kansas), and 41.25% of potato fields (56.50% in Nebraska, 23.80% in North Dakota, and 93.75% in Wyoming).

For DNA analysis, nematodes from fields in Nebraska, Kansas, Wyoming, North Dakota, Montana, and Colorado were chosen to maximize geographic and host coverage across the region (Table 1). Samples from outside the Great Plains were added for phylogenetic context including topotype specimens of P. alleni isolated from a soybean field five miles north of Eldorado City, in Saline County in Illinois. A total of 915 specimens of Pratylenchus were sequenced for a 727 to 739-bp fragment of the COI gene. The sequence length was 730-bp for most of the specimens, except P. neglectus (727-bp), and Pratylenchus sp. 9 and Pratylenchus sp. 10 (739-bp). The initial set of 915 Pratylenchus COI haplotypes was reduced to a 143 unique-sequence data set by the removal of redundant sequences. Maximum likelihood, neighbor-joining, and Bayesian phylogenetic trees each identified 19 distinct COI haplotype groups that were well-supported by bootstrap values (BS: 96-100% based on tree construction method), intra- vs interspecific genetic distance and posterior probabilities (PP:100 for all groups) (Fig. 2). These haplotype groups were tentatively labeled as P. neglectus, P. scribneri, P. thornei, P. hexincisus, P. alleni, P. penetrans, P. zeae Graham, 1951, P. crenatus Loof, 1960, P. vulnus Allen & Jensen, 1951, Pratylenchus sp. 1 to sp. 10, including one unnamed singleton, NID 6402 from Kansas. Species identity was based on a combination of factors including distinctive morphological features, agreement with GenBank vouchers, congruence between independent markers, and in the case of P. alleni, topotype specimens. Four haplotype groups in Figure 2 consisted exclusively of specimens collected outside the Great Plains. These were P. crenatus from Ireland, Poland, and Puerto Rico, P. vulnus from California, and species 9 and 10, both exclusively represented by specimens from Arkansas. Among the 19 haplotype groups, males were recorded in seven haplotype groups including P. penetrans, P. alleni, P. crenatus and haplotype groups sp. 1, sp. 2, sp. 9, and sp. 10. In total, 14 haplotype groups and one singleton (NID 6402) were associated with agroecosystems of the Great Plains Region.

Figure 2:

DNA sequences of the D2-D3 region of the 28S rRNA were evaluated in a phylogenetic tree using 24 specimens previously sequenced for the COI gene together with 70 additional sequences retrieved from GenBank (Supplementary Figure 1). Of the 24 sequences, 10 were placed within well-supported clades representing named species of Pratylenchus. These include P. alleni (NIDs 10848, 10850, 3717), P. crenatus (NID 8539), P. thornei (NIDs 7566, 7567), P. neglectus (NID 10756), and P. penetrans (NID 7091, 6260, 6261). Sequences of the D2–D3 region for the Pratylenchus species sp. 1, sp. 2, sp. 5, sp. 8, and sp. 9 were not closely associated with known Pratylenchus species in the NJ, ML, or Bayesian phylogenetic trees.

A BLAST search on GenBank with Pratylenchus COI sequences generated in this study helped to determine species identity through near identical matches (99-100%) for seven of the COI haplotype groups (Table 2). These well-supported haplotype groups corresponded to P. neglectus, P. penetrans, P. thornei, P. alleni. P. zeae, P. crenatus, and P. vulnus. Pratylenchus scribneri, the second most abundant Pratylenchus species in this survey, provided conflicting results in GenBank BLAST searches. For specimen NID 7839, the percentage identity of the top 9 matches ranged from 98.5% (P. scribneri KY424092) to 99.5% (P. scribneri MH016378) and bracket two close matches to P. hexicusus of 98.7% (KY828320 and KY828322). Similarly, for the D2-D3 domains of the 28S rRNA, the top 25 matches for P. scribneri (percentage identity range from 99.7% to 98.1%), encompassed a range that included three P. hexincisus entries of 99.7%. Specimens hypothesized to be P. hexincisus had no close COI match identified as P. hexincisus (Table 2; Ozbayrak et al. in prep). For the D2-D3 domains, those specimens exhibited a relatively close match to accessions representing P. hexincisus (KY828291–98.9% similarity) and P. scribneri (KX842632–100% similarity). In total, 10 haplotype groups with no clear taxonomic affinities based on DNA sequences were solely identified as Pratylenchus ‘sp.’ These specimens had low percentage identity scores for the COI gene and moderate or ambiguous scores for the D2–D3 domains.

Intragroup and intergroup genetic distance matrices are presented in Table 3. Estimated mean genetic raw distance (p-distances) for all haplotype pairs was 0.20. Pairwise genetic distances within most haplotype groups ranged from 0 to 0.048. Haplotypes groups Pratylenchus sp. 2 and Pratylenchus sp. 9 exhibited the highest intraspecific diversity with distance values of 0.077 and 0.048, respectively. Interspecific variability among haplotype groups was exceptionally high, ranging from a low of 0.079 between haplotypes Pratylenchus sp. 5 and Pratylenchus sp. 6 to 0.409 between haplotypes P. zeae and Pratylenchus sp. 2.

Molecular species delimitation using ABGD suggested 19 to 61 haplotype groups in the recursive partitioning based on differing intraspecific divergence priors (Fig. 3). The last three partitions settled on 19 groups with P-values of 0.0359, 0.0599, and 0.100, respectively. ABGD derived groupings at these partitions were congruent with the haplotype groups on phylogenetic trees.

Figure 3:

Both GYMC and TCS identified multiple subdivisions within the 19 groups derived by ABGD and the highly supported clades in the phylogenetic trees. The single threshold GMYC analysis revealed 22 ML clusters (CI 22–23) and 39 ML entities (CI 39–42) when using a strict molecular clock and yule tree prior (Fig. 3). The likelihood for the null model was 972.166 whereas the maximum likelihood of the GMYC model was 998.4352. The likelihood ratio (52.5385) test rejected the null hypothesis for the models tested, and thus the assumption that all sequences belonged to the same species (LR test: P < 0.001). TCS analysis resulted in 36 and 32 haplotype networks at connection limits of 95 and 90%, respectively, for the complete data set (Fig. 3).

The estimation of divergence time of the major haplotype groups was investigated using BEAST using two different molecular clock rates to calibrate dates of lineage divergence (Fig. 4). The nematode molecular clock rate derived from Caenorhabditis briggsae, using an assumption of two generations per year for Pratylenchus, resulted in younger node ages than the insect substitution rate. Divergence time estimation based on the nematode clock rate revealed an early split of Pratylenchus clades occurred at an estimated 6.28 Mya (CI: 5.12–7.59 Mya) in the late Miocene epoch. Pratylenchus neglectus may have diverged from P. penetrans as early as 4.73 Mya (CI: 3.56–5.61 Mya), and the split of P. thornei from P. zeae at an estimated 5.39 Mya (CI: 3.38–5.47 Mya).

Figure 4:

Later Pratylenchus divergences occurred in the Pliocene epoch at approximately 2.7 Mya, prior to the onset of the ice ages. The most recent common ancestor (MRCA) of P. hexincisus and P. scribneri emerged in the mid-Pleistocene approximately 1.46 Mya (CI: 1.24–1.69 Mya). The MRCA of P. scribneri, a relatively young clade, emerged in the late Pleistocene epoch at 0.65 Mya (0.48–0.73 Mya). Three lineages of sexually reproducing species, Pratylenchus sp. 9, Pratylenchus sp. 10, and P. alleni, all diversified during early to mid-stages of the ice ages between 2 million and 500,000 years ago (Fig. 4).

The most frequently sampled haplotype group in the Great Plains was Pratylenchus neglectus, comprising 53% of all specimens. P. neglectus was identified in 96, 90, and 83% of all wheat, potato, and dry bean fields, respectively. On the other hand, this haplotype group was encountered in corn fields less frequently, collected in 42% of the corn fields across the Great Plains (Fig. 5). Four distinct haplotype networks labeled A, B, C, and D were recognized for P. neglectus at both the 95 and 90% connection limit. These four P. neglectus subgroups were also detected in GMYC delimitation methods. Network A consisted of eight unique haplotypes (or subgroups), occurring in 10 states of North America and Canada (Fig. 6A). The most abundant haplotypes were neg1, neg2, and neg3 within network A. Haplotype neg1 specimens were found in six states and were most often associated with corn, but also associated with dry bean, wheat, alfalfa, potatoes, and cereal rye. Haplotype neg2 and neg3 specimens were found in 10 and 9 states, respectively, as well as Canada, and were most frequently recovered from wheat, but also associated with other crops (Fig. 6B). A total of six haplotypes from networks B, C, and D were located in western and northwestern states of Colorado, Wyoming, Idaho, Montana, and North Dakota with only a few specimens from Nebraska and Kansas included in these networks (Fig. 6A). Specimens in networks B, C, and D were associated with barley, potatoes, alfalfa, dry beans, corn, wheat, and sugar beets (Fig. 6B).

Figure 5:

Figure 6:

The second most abundant haplotype group in the Great Plains, comprising 30% of all specimens, was P. scribneri. This haplotype group appeared as a single TCS network comprised of 17 closely related haplotypes (or subgroups) (Fig. 7). These haplotypes were recovered from 104 fields (45 counties) and were associated most often with corn fields in 4 states; 78 in Nebraska, 14 in Kansas, 2 in South Dakota, and 1 in Montana. Also, they were recorded from four potato fields in Nebraska, one wheat field in both Nebraska and Texas, and a sugar beet field in Colorado (Fig. 5).

Figure 7:

Pratylenchus thornei was comprised of a single haplotype distributed across six different states (Fig. 8A). P. thornei was primarily associated with wheat in Kansas but was also collected from corn and alfalfa in Montana, corn in Oklahoma, a single plot in a cover crop (cereal rye) experiment in Nebraska, in a vineyard in California, and a sugar beet field rotated with barley in Colorado.

Figure 8:

Pratylenchus hexincisus specimens were split into two separate groups in TCS networks and GMYC methods (Fig. 8B), but ABGD only recognized a second network at P-values between 0.008 and 0.022 (Fig. 8). Pratylenchus hexincisus was primarily associated with corn fields in eastern Kansas, Nebraska, and South Dakota, as well as from dry beans in Wyoming and wheat in North Dakota.

Four haplotype groups with male specimens, P. penetrans, Pratylenchus sp. 2, sp. 9, and sp. 10, displayed diverse TCS networks and multiple haplotype subgroups in GMYC methods. The final partition in ABGD, however, recognized these four sexual groups as four distinct entities. For P. penetrans, ABGD revealed two distinct groups at P-values of 0.0359, 0.0599, and 0.100. Six haplotypes were detected, suggesting three distinct networks (networks not shown) in TCS analyses at both 90 and 95% cutoff values, and three entities were revealed in GMYC analysis. Pratylenchus penetrans was not common in the agronomic crops sampled in the Great Plains region. It was recovered from one south central cornfield in Nebraska, an apple orchard in eastern Nebraska, and two corn fields in Montana. Species less commonly found in this Great Plains survey included P. zeae and P. alleni, both collected from single corn fields in Nebraska.

Eleven Pratylenchus sp. 2 haplotypes comprised six networks at 95 and 90% connection limits, respectively, and GMYC methods exhibited nine entities. Pratylenchus sp. 9 and sp. 10 had four networks at 95% cutoff and two at 90% cutoff, respectively. GMYC analysis revealed four entities for both groups. Species delimitation results for other groups exhibited one network per haplotype per entity due to the existence of few specimens in these groups.

Nearly all haplotypes representing unknown Pratylenchus species were associated with corn, with the exception of one wheat, one soybean, and one cotton field. Barcoding also revealed that 44 of the 439 Pratylenchus-infested fields had a mixed population of at least two Pratylenchus species. Mixtures were recorded in 19.5, 10.6, and 27.8% of the corn, wheat, and dry bean fields, respectively. Five corn fields in Kansas and 23 corn fields in Nebraska had a mixture of different haplotype groups. The most common combination of species was P. neglectus and P. scribneri, which were recovered together from approximately 55% of all mixed fields (Fig. 9).

Figure 9: