Intraspecific variation in phenotypic and phylogenetic features among Pratylenchus penetrans isolates from Wisconsin, USA

Abstract Pratylenchus penetrans is a common and important agricultural pest in Wisconsin, a USA state with a diverse agriculture. We compared populations from around the state to each other and to data published for populations around the world to gain insight on the variability of features important for identification of this cosmopolitan species. Thirteen isolates from samples collected in soybean fields in ten Wisconsin counties were established in monoxenic cultures. Analysis of morphological features revealed the least variable feature for all isolates collectively was vulva percentage. Features less variable within than among isolates were body width, lip region height, and stylet length. Some isolates showed only the smooth tail tip phenotype and others had a mix of smooth and annulated tail phenotypes. A suite of features provided sufficient pattern to group isolates into four clusters according to hierarchical agglomerative clustering and canonical discriminative analyses, but not with enough distinction to be useful for classification. Haplotype analysis based on the COI mitochondrial gene of the 13 cultured isolates, 39 Wisconsin field populations, and published sequences representing five additional USA states and six countries revealed 21 haplotypes, 15 of which occurred in Wisconsin. Ten haplotypes represented in Wisconsin were shared with populations from Europe, South America, Africa, or Asia. Five haplotypes were unique to Wisconsin, six were unique to The Netherlands, and one was unique to Japan suggesting that even more COI diversity will be revealed when more COI sequences for P. penetrans become available. The maximum pairwise sequence variation was 6% and the SNPs did not alter amino acids, indicating cryptic biodiversity within the species worldwide. The cosmopolitan to localized scale of distribution of COI haplotypes could be due to frequent and ongoing dispersal events, facilitated by life history traits and the broad host range of P. penetrans. Regions of diverse agriculture, like Wisconsin, show promise for studying this important pest and our study confirms the utility of the COI mtDNA gene for studying variation within a species.

common pest of fruit, vegetable, grain, and forage crops.
Identification of P. penetrans is supported by multiple species keys for the genus Pratylenchus (Loof, 1960;Corbett, 1970;Handoo and Golden, 1989;Ryss, 2002;Castillo and Vovlas, 2007). Consensus morphological characters are three lip annuli, round spermatheca, presence of males, smooth tail tip, and four incisures in the lateral field. Morphological measurements and indices of morphological features are also important, particularly the proportional position of the vulva (V%) and relative length of the esophagus and tail.
Morphological features and measurements of P. penetrans have been recorded many times since its original description as Tylenchus penetrans by Cobb (1917). We found records of four populations from the USA and 32 populations from other countries (Sher and Allen, 1953;Taylor and Jenkins, 1957;Loof, 1960;Román and Hirschmann, 1969;Tarte and Mai, 1976a;Vovlas and Troccoli, 1990;Townshend, 1991;Troccoli et al., 1992;Mekete et al., 2011;Mokrini et al., 2016;Janssen et al., 2017;Ozbayrak, 2019;Rusinque et al., 2020). The range in values for quantitative features barely overlaps among some populations. Intraspecific variation in tail morphology was noted by Taylor and Jenkins (1957), Román and Hirschmann (1969), Tarte and Mai (1976a), Townshend (1991) and Rusinque et al. (2020), including the observation that tail features as well as body size could be modified by environmental factors such as plant host and culture media (Tarte and Mai, 1976a;Townshend, 1991).
There is consensus that nematode identification today should be supported by molecular analyses so sequence data are rapidly accumulating for multiple gene fragments including the cytochrome c oxidase subunit1 (COI) mitochondrial gene (Janssen et al., 2017), rDNA internal transcribed spacers (De Luca et al., 2011), D2/D3 expansion region of 28S rDNA, and 18S rDNA (Subbotin et al., 2008). Janssen et al. (2017) used a combination of morphological and sequence data of COI to successfully distinguish species among the Penetrans group of Pratylenchus spp; P. penetrans, P. fallax, P. convallariae, P. dunensis, P. pinguicaudatus, P. brachyurus, P. arlington, and P. oleae.
Intraspecific genetic variation among P. penetrans populations was addressed recently by two studies (Mokrini et al., 2016;Janssen et al., 2017). Mokrini et al. (2016) found relatively little variation among 13 Moroccan isolates of P. penetrans for the D2/ D3 expansion region of 28S rDNA. Janssen et al. (2017) reported 16 haplotypes based on COI for P. penetrans isolates collected around the world.
Their analysis revealed 14 of the 16 haplotypes were represented in The Netherlands. Three of the 16 haplotypes were shared by populations from The Netherlands, Africa, and Asia. North America was represented by only one population of P. penetrans and it was the only member of one haplotype. The study by Janssen et al. (2017) also revealed cryptic biodiversity within P. penetrans, as morphologically indistinguishable populations from The Netherlands were assigned to different haplotypes.
We recently surveyed Wisconsin for species of Pratylenchus with males to determine the prevalence of P. penetrans in crop rotations including soybean. That project (Saikai, 2019) resulted in first reports of P. fallax  and P. alleni  in Wisconsin, but also showed the vast majority of populations with Pratylenchus males to be P. penetrans based on sequence data from 28S rDNA, 18S rDNA, and a species-specific primer. Our objectives in this study were (i) to assess the extent of intraspecific variability in P. penetrans in Wisconsin by comparing 13 populations of P. penetrans for morphological variation and 52 populations for sequence variation in the COI mitochondrial gene, and (ii) to determine the extent of diversity in COI for P. penetrans by combining our data with published COI sequences from populations around the world in phylogenetic and haplotype analyses.

Nematode isolates
Nematodes were recovered from soil samples collected from 11 farms in nine Wisconsin counties for a commodity-sponsored testing program. The samples were randomly selected from a pool of stored samples that were known to contain Pratylenchus males. Nematodes were extracted from root fragments sieved from the soil and incubated in water for 48 hr. One female from the incubated fraction was placed on an excised root culture of 'IO Chief' sweet corn (Zea mays L.) grown in Gamborg's B5 media, with the assumption that insemination occurred in the field. Cultures were maintained in the dark at 24°C and transferred to new Gamborg's B5 media every five months to increase nematodes. Each population was verified as P. penetrans by morphology, the species-specific primer PPEN (Al-Banna et al., 2004), and matching sequences of the D2/D3 expansion region of 28S rDNA to accessions reported in GenBank. Two additional isolates of P. penetrans in culture for more than 25 years (Portage1 and Portage2), originally extracted from potato, were included in the study for a total of 13 isolates, each named for the county of origin (10 counties total).

Morphological characterization of isolates
Nematodes from at least three Petri dishes were harvested and pooled for morphological characterization, given the reported sensitivity of some features to environmental conditions (Tarte and Mai, 1976a). Fresh specimens were mounted following the protocol of Eisenback (2012), visualized using differential interference contrast optics, and measured and photographed using NIS-Elements AR Software on a Nikon Eclipse Ti-E inverted microscope (Nikon Instruments, Melville, NY). Tail terminus morphology, tail shape, and spermatheca shape were recorded. Data were collected from a minimum of 25 females of each isolate, except for Chippewa2 (n = 18 females), for 26 morphometric and morphological variables. The same morphometrics, except for those that are specific to females, were measured for males of each isolate but only spicule length and gubernaculum are listed in this study. The following indices were derived: V% (anterior end to vulva length/total body length), and a (total body length/maximum body width), b (total body length/esophageal length), c (total body length/tail length), b′ (total body length/distance from anterior end to the end of esophageal glands), and c′ (tail length/body width at anus). Tail shape and tail tip appearance were assigned according to the categories used by Frederick and Tarjan (1989). The coefficient of variation (CV) was computed for each morphological variable by isolate and for pooled data representing all isolates.

Phenotype analyses of isolates
Intraspecific variation in the morphometrics among the P. penetrans isolates were evaluated with Hierarchical Agglomerative clustering (HAC) and Canonical Discriminant analysis (CDA) for the following female features: L, V%, ratios of a, b, c, c′, maximum body width, esophageal length, length of anterior end to excretory pore, lip region height, stylet length, DGO, PUS, body width at anus, tail length, vulva to anus distance, lateral incisures, and labial incisures. Based on high correlation coefficients with other features, the b′ ratio, length of head end to posterior gland and body width at vulva were excluded from phenotype analyses. The shape of the tail and annulation of the tail terminus were included as binary variables: subhemispherical tail = 0, not subhemispherical tail = 1; smooth tip = 0, crenate tip = 1.
For HAC, Euclidean distance coefficients were computed on normalized mean values to evaluate the (dis)similarity between each pair of isolates using the 'dist' function with the Euclidean method option in the 'STATS' R package (R Core Team, 2020). The normalized means were computed by ′ = − − X X X X X min max m in , so that all the morphometrics ranged between 0 and 1. HAC was performed using the 'hclust' function of the same R package with Ward's method. The optimum number of clusters was determined using the Elbow and Silhouette methods.
Canonical discriminant analysis (CDA) of the clusters suggested by HAC was performed using the CANDISC procedure in SAS (SAS Institute, Cary, NC) to determine which variables discriminated clusters. Univariate statistics (ANOVA) were also conducted for each morphometric variable to test equality of means among the 13 isolates.

COI sequencing
Nematodes from the 13 cultured isolates and 39 field populations collected from 39 farms in 27 Wisconsin counties were characterized molecularly for mitochondrial cytochrome c oxidase subunit 1 (COI) (Fig. 1). The 52 populations were named for the county of origin and multiple populations from the same country were distinguished by numbers and individual specimens by letters (Table 3). Nematode DNA was extracted using DNeasy Blood and Tissue Kit (Qiagen, Valencia, CA). The DNA extracts were stored at −20°C if PCR was not conducted immediately. Five microliters of DNA were suspended in a 25 μ l reaction volume composed of 2.5 μ l ddH 2 O, 2.5 μ l of forward and reverse primers, and 12.5 μ l of GoTaq Green Master Mix x2 (Promega, Madison, WI). The primer used for COI was forward JB3 (5′-TTTTTTGGGCATCCTGAGGTTTAT-3′) and reverse JB4.5 (5′-TAAAGAAAGAACATAATGAAAATG-3′) (Derycke et al., 2010). PCR cycling conditions were: 94°C for 5 min, 5 cycles of denaturation at 94°C for 30 sec, annealing at 54°C for 30 sec and temperature decreasing with 1°C for each cycle, and extension at 72°C for 30 sec followed by 35 cycles of 94°C for 30 sec, 50°C for 30 sec, and 72°C for 30 sec, and a final extension of 10 min at 72°C. Five micro liters of the PCR products were resolved by electrophoresis on 1.0% agarose gel mixed with 2.5 μ l of SYBR Safe DNA Stain (Invitrogen, Waltham, WA) at 70 V for 1 hr. The presence of DNA in the PCR products was confirmed by Gel Logic 200 Imaging System (Kodak alaris, Rochester, NY). The crude PCR products were sent to Functional Biosciences (Madison, WI) for purification and Sanger-sequencing.
The DNA sequences were edited with Finch TV (Geospiza Inc, Seattle, WA), and aligned with MUSCLE in MEGA 7 (Kumar et al., 2016) using default parameters with the gap opening penalty set at 15 and the gap extension penalty at 6.66. Alignments were further adjusted by inspection. A complemented sequence was constructed by overlapping sequence reads with forward and reverse primers for each specimen. Basic Local Alignment Search Tool (BLAST) from the National Center for Biotechnology Information (NCBI) was used to compare the complemented sequences to the GenBank sequence database. Sequences were deposited to GenBank and assigned accession numbers (Table 3).

Haplotype network and phylogenetic analyses
Redundant COI sequences among multiple individuals sequenced from the same population were dropped, leaving 72 sequences for the haplotype analysis, 13 from cultured isolates, and 59 from field populations. In total, 48 published sequences from populations in other USA states (Ozbayrak et al., 2019) and the populations used in the Janssen et al.'s (2017) study were also included for a total of 120 sequences used in the analyses. The haplotype network was calculated and visualized using TCS network in the software PopART Version 1.7. Geographical origin of the sequences was used as trait information in the analysis.
For phylogenetic analysis, maximum-likelihood (ML) trees were constructed based on the aligned sequences of COI using the default program settings in MEGA 7 with 1000 bootstrap replications. The best DNA model of the nucleotide substitution tool option in MEGA 7 selected the Hasegawa-Kishino-Yano model with a gamma distribution and proportion of invariable sites (HKY + G + I) for COI. The COI sequences of all the Wisconsin populations and CA85, T716, T723, T726, T713, T132 as well as the other species in the Penetrans group and P. crenatus were included in the analysis. Meloidogyne hapla (JX683719) and M. camelliae (KM887147) were included as outgroup taxa per the study by Janssen et al. (2017).

Morphology and morphometrics of isolates
All 13 Wisconsin isolates possessed the diagnostic characters of the species; three lip annuli, round spermatheca filled with sperm, four incisures in the lateral field, and a majority of specimens had a smooth tail tip (Tables 1 and 2). Stylet basal knobs were usually round but sometimes cupped anteriorly. The CV for the pooled data of all isolates was greatest for PUS (24.8%) followed by lip region height (17.4%) and c′ (17.3%). The morphometrics with the smallest variability were V% (1.8%), stylet length (6.0%), and length of head end to the end of the esophageal gland (9.0%). The maximum value for the CV for individual isolates was lower than the pooled CV for body width, lip region height, and stylet length.
Mean morphometric values for the isolates were all within the range reported for P. penetrans by Loof (1960) except three isolates were lower than Loof's range (5.3-7.9) for b, one was higher than Loof's range (15-24) for c, and one was lower than Loof's range (15-17) for stylet length (Tables 1 and 2). All isolates had four lateral incisures but individuals in seven isolates had striae in the middle band in the vicinity of the vulva. Nine isolates had a subhemispherical tail with a smooth tip. Individuals that deviated from this typical tail were detected in four isolates; the frequency of truncate tails ranged from 10 to 30% for three isolates and the frequency of annulated tail tips was ca 10% for three isolates. Males were common in all the isolates.

Phenotype analysis of isolates
HAC separated the 13 P. penetrans isolates into four clusters based on the Euclidean distance in canonical space ( Fig. 2): (I) Chippewa1 and Marquette1, (II) Chippewa2, Marathon2, Iowa and Portage1, (III) Buffalo and Grant, and (IV) Calumet, Sheboygan, Portage2, Marathon1, and Wood. There was no association between clusters and geographical origins of the isolates. Isolates in HAC clusters I and II had smooth and round tail tips, whereas both isolates in the HAC clusters III and some isolates in the HAC cluster IV (i.e. Calumet and Sheboygan) had a mix of both characters.
The clusters determined by HAC were confirmed by the CDA analysis (Wilks' λ < 0.001) (Fig. 3). The means of all the morphometrics, except for esophageal length (P = 0.06), were not equal among the 13 isolates (P < 0.01). The first two canonical variates for the CDA models explained 88% of the variability in morphometrics for clusters based on the HAC groups. The variables with the greatest influence on the discriminant functions (vector loading > 0.5) were lip region height (0.69), maximum body width (−0.61), stylet (−0.55), body width at anus (−0.53), and a (0.50) for the first canonical variate; maximum body width (0.64), total body length (0.61), length from anterior end to excretory pore (0.52) for the second canonical variate; and vulva percentage (0.71) and tail length (0.54) for the third canonical variate. There was reasonable separation of clusters I, II, and III, but the distinction of cluster IV from the other clusters was not resolved (Fig. 3).

Molecular characterization
The 402 bp of nucleotide sequence was read without gaps by the COI primer. The pairwise sequence identity between the 13 Wisconsin isolates ranged from 94.2% (Portage1 and Marathon2) to 100% similarity (Table 4). There were 30 single nucleotide polymorphisms (SNPs) on COI sequences among the 13 P. penetrans isolates. The most common mutations were between C and T. Similar mutations were also observed at the same sites among the COI sequences of our field populations and those obtained from GenBank. The largest nucleotide divergence in COI between a Wisconsin population and a population elsewhere was 92.7% for Sheboygan-a and the T716 population (Stramproy) from the Netherlands (Janssen et al., 2017), but the species identification of the Dutch population was not confirmed by a gene in addition to COI. The largest nucleotide divergence in COI with nematodes confirmed to be P. penetrans by 28S rDNA was between Marathon2-a and the V1B population (Meijel) from the Netherlands (Janssen et al., 2017), at 94.5% similarity. There was no variation in amino acid sequences of P. penetrans populations as a result of SNPs on COI.

Haplotype network analysis
A total of 120 sequences, 72 from Wisconsin, revealed 21 haplotypes for P. penetrans based on COI (Fig. 4).  The Wisconsin populations were assigned to 15 of the 21 haplotypes: three shared with other populations from the USA, seven shared with populations from other countries, and five that were unique to Wisconsin. The most common haplotype in Wisconsin was Pp1, which included 32 Wisconsin populations as well as populations from Nebraska, USA, Colombia, The Netherlands, and Japan (Fig. 4). The second most frequent haplotype in Wisconsin, Pp10, included 19 Wisconsin populations and one population from Idaho, USA. In 17 cases, multiple specimens from the same population were assigned to different haplotype groups (Table 3). No relationship was observed between haplotype and HAC assignment for the 13 Wisconsin isolates.

Phylogenetic analysis
The ML tree based on COI validated that the 72 Wisconsin populations were conspecific with P. penetrans (Fig. 5). The bootstrap value supported the assignment of P. penetrans populations to a clade distinguished from the other species in the Penetrans group and P. crenatus. Within the P. penetrans clade, T716, the Dutch population from Janssen et al. (2017), was assigned to a sub-clade that was distant from all others, but it should be noted that species confirmation by 28S rDNA or other genes is not available for that population. Haplotypes Pp1 to Pp4 were distinguished from others and grouped together in one clade, as were Pp6 to Pp8. The evolutionary relationship of other haplotypes was not clearly resolved in the analysis.

Discussion
Our study demonstrated considerable intraspecific variation in morphological and molecular characters for P. penetrans populations from Wisconsin that was not associated with geography. Haplotype analysis showed differences among individuals collected from the same field population, as was the case for P. penetrans populations in the Netherlands (Janssen et al., 2017). Cultured isolates originating from samples collected in the same county (Chippewa) also differed for haplotype and showed morphological variation detected by HAC analysis. Intraspecific variation at local scales was juxtaposed to similarity in haplotypes that spanned continents. Wisconsin populations shared the same haplotypes as populations from North America, South America, Europe, Asia and Africa. Populations of P. penetrans from the Netherlands were sympatric for COI with populations from Europe, South America, Asia, and Africa (Janssen et al., 2017), so this is the first report of P. penetrans COI haplotypes shared by Europe and North America. Six haplotypes unique to the Netherlands (Janssen et al., 2017)   with our analysis and we detected five haplotypes unique to Wisconsin, so it is likely that even more intraspecific variability will be revealed as more populations of P. penetrans are sequenced. Three morphological features used to distinguish P. penetrans as a species, three lip annuli, a round spermatheca filled with sperm (i.e. presence of males), and four lines in the lateral field, were robust for our Figure 4: COI (mitochondrial cytochrome c oxidase subunit 1) haplotype network of Pratylenchus penetrans using TCS network. Circle size and hash marks between haplotypes represent numbers of specimens included in each haplotype and mutations, respectively. Geographical location of the population is designated by color. Median vectors are shown as black circles.
isolates. Most keys for identifying Pratylenchus spp. assign a smooth tail tip to P. penetrans (Loof, 1960;Corbett, 1970;Handoo and Golden, 1989;Castillo and Vovlas, 2007), but we and others (Taylor and Jenkins, 1957;Tarte and Mai, 1976a;Townshend, 1991;Rusinque et al., 2020) noted some individuals with annulated tail tips. Tarte and Mai (1976a) demonstrated that the frequency of annulated tail tips was influenced by the plant host, but they also found the only case of exclusively smooth tail tips in an isolate started with a single smooth-tailed female (Tarte and Mai, 1976b). They speculated that tail tip phenotype is at least partly under genetic control. Nine of our isolates had a 100% frequency of the smooth tail tip phenotype, but four did not and all were reared in vitro in the same stable environment, supporting the possibility that tail tip phenotype is determined genetically.
Data for 49 isolates and populations of P. penetrans (Sher and Allen, 1953;Taylor and Jenkins, 1957;Loof, 1960;Román and Hirschmann, 1969;Tarte and Mai, 1976a;Vovlas and Troccoli, 1990;Townshend, 1991;Troccoli et al., 1992;Ozbayrak, 2019;Ryss, 2002;Mekete et al., 2011;Mokrini et al., 2016;Janssen et al., 2017;Rusinque et al., 2020), including our study, provided opportunity to review morphometrics of the species presented in seminal publications by Sher and Allen (1953) and Loof (1960). Loof's range for V% (75-84) was a closer fit to the data (n = 39) than that of . The minimum mean value for V% was 75% shared by two isolates from Morocco (Mokrini et al., 2016) and the maximum was 82% for a population from Nebraska (Ozbayrak, 2019). In total, 49% of the isolates and populations had a mean value for V% of less than 80%, suggesting that the range of 80 to 85% used in a recent tabular key (Castillo and Vovlas, 2007) is not representative of many P. penetrans, including three isolates from Wisconsin. The smallest value for mean a index (n = 46) was 21.1 in a population from Algeria (Troccoli et al., 1992) and the maximum was 33.1 in two isolates from Morocco (Mokrini et al., 2016), which was above the maximum set by Loof (31.9) and Sher and Allen (30). The largest mean value of the b index (n = 43) was within the maximum set by Loof (7.9), but well above the maximum of Sher and Allen (6.5) and the minimum value was 4.5 for a population from the Netherlands (Janssen et al., 2017) which was below the minimum of both Loof (5.3) and Sher and Allen (5.7), as were Haplotypes are labeled and also designated by the vertical red bar. Information about the haplotypes appears in Table 3. three Wisconsin isolates. Three isolates (n = 46), one each from Wisconsin, Morocco (Mokrini et al., 2016) and the Netherlands (Janssen et al., 2017), had a mean value for the c index that exceeded the upper limit published by Loof (23.8) and Sher and Allen (21). Seven Wisconsin isolates and a population from Canada (Townshend, 1991) had a mean stylet length below 15 μ m and one population from Portugal had a mean stylet length of 18.5 μ m (Rusinque et al., 2020), falling outside the maximum published by Loof (17 μ m) and the minimum published by both Loof (15 μ m) and Sher and Allen (17 μ m). Most of the populations outside of the commonly recognized ranges were identified using molecular as well as morphological criteria, so they are likely to be valid P. penetrans which suggests there is even more intraspecific variability within the species than previously considered.
We detected enough pattern in morphological variation to assign our 13 Wisconsin isolates of P. penetrans to four clusters. Two features, stylet length and height of the cephalic lip region, showed less variation within than among isolates and were identified as important by CDA analysis. Contrary to our results, these features were similar among 13 Moroccan P. penetrans isolates reared in vitro (Mokrini et al., 2016). Based on the relatively weak separation of some of our clusters and the possibility that our results are idiosyncratic to Wisconsin, we do not consider our clusters to represent morphotypes.
Our haplotype analysis revealed up to 5.8% divergence in COI DNA among populations confirmed to be P. penetrans by multiple genes. As was the case in the Netherlands (Janssen et al., 2017), the nucleotide divergence did not change amino acid sequences, confirming that they were still conspecific with P. penetrans. Our ML tree showed the same four groups presented by Janssen et al. (2017) with Wisconsin isolates assigned to three (i.e., Janssen's A, B, and C) groups. The most common haplotype in Wisconsin, Pp1, was the second most common haplotype of the Dutch populations (Pe11) and included populations from South America. The second most common haplotype in Wisconsin, Pp10, was comprised exclusively of P. penetrans populations from the northern USA. In total, 11 haplotypes were unique to a single country. A gradient ranging from cosmopolitan to localized distribution suggests that dispersal events for P. penetrans are frequent and ongoing.
The widespread incidence and high degree of intraspecific variability of P. penetrans coupled with the state's edaphic and agricultural diversity makes Wisconsin an interesting place to study this cosmopolitan plant pest. The success of P. penetrans is most likely due to its life history traits and agricultural activities. The species, highly polyphagous and capable of long-term survival in an anhydrobiotic state (Townshend, 1984), was likely introduced to Wisconsin in vegetatively propagated crops such as potato or in grain contaminated with soil peds. It would be interesting to compare Wisconsin and German isolates, as 37% of the state's citizens in 1890 originated from Germany and most were farmers (Wyman, 1968). There are many approaches and characters useful for studying intraspecific variability for nematodes. Our results echo previous studies (Janssen et al., 2017;Ozbayrak et al., 2019) in demonstrating the utility of the COI mtDNA gene for distinguishing populations of P. penetrans as well as their relationships.