Morphological Characterization and Assessment of Genetic Variability of Tylenchulus semipenetrans Populations from Southern Iran

Nearly 60 citrus nematode-infested orchards from different citrus growing areas with different soil properties in Fars province, Iran were randomly sampled from August to October 2018 and 2019. The citrus orchards in Fars province are mainly located in six regions known collectively as the citrus belt (Fig. 1). The GPS coordinates for each sample were recorded using a GPS device (Garmin ETrex 32x) and are shown in Table 1. The majority of host trees for citrus nematode populations in the present study were orange trees (Citrus sinensis L.), sour orange (C. aurantium L.), sweet lemon (C. limetta R.), mandarin orange (C. reticulata L.) and bitter orange (C. bigaradia L.). Citrus trees (Citrus spp.) in Fars province are usually grafted onto sour orange (Ebadi et al., 2014). Indian bael (Aegle marmelos) is also one of the rootstocks for citrus trees in the province. In addition, some citrus trees, usually sour orange, are grown on their roots. The samples were taken from the top layer of soil, up to 30 cm below the tree canopies. To minimize the influence of the dried topsoil, the surface soil (1–3 cm) was discarded to ensure that no roots of herbaceous plants were included.

Figure 1:

Geographical location of the sampled areas. The majority of the samples were collected in the citrus-growing regions (known as the citrus belt) of Fars province, southern Iran (shown in green: Kazerun, Shiraz, Ghir, Jahorm, Fasa and Darab). Two samples were also collected from Mazandaran province in northern Iran (shown with red border).

Soil samples (approximately 1 kg each) were placed in plastic containers, labeled, and taken to the nematology laboratory at the School of Agriculture, Shiraz University. The collected samples were thoroughly mixed by hand and prepared by sieving (1.25 μm openings). Then the males and second-stage juveniles (J2) of the citrus nematode in 200 cm³ subsamples were extracted for 48 hours using the tray method (Whitehead & Hemming, 1965). Because of the long period between nematode extraction and molecular analysis (DNA extraction), the collected soil samples were inoculated onto sour orange rootstocks under greenhouse conditions. DESS solution (0.25 M disodium EDTA at pH 8.0, 20% dimethyl sulfoxide [DMSO], and saturated NaCl) was used to preserve the individual nematodes until molecular analysis was performed at the Nematology Laboratory of the Czech University of Life Sciences Prague (Yoder et al., 2006; Perry et al., 2020).

For light microscopy, the extracted J2 of T. semipenetrans from the rhizosphere of citrus trees were killed and fixed in hot formaldehyde-acetic acid (4:1), and processed to anhydrous glycerol using Seinhorst’s method (1962). The morphological examination was carried out with an Olympus BX41TF microscope equipped with a camera at up to 1,000x magnification. The nematode population in the infected root was visualized and examined by staining the roots in fuchsin acid and subsequent clarification in acidified glycerol (Eisenback & Hunt, 2021). The morphometric indices of the males and J2s were also measured. The citrus nematode species was verified using Tylenchulus identification keys (Inserra et al., 1988; Tanha Maafi et al., 2012). To find out whether the populations have different morphometrics, a Principal Component Analysis (PCA) of the morphometric characters of males and J2 females of the nematode was performed in R version 3.5.1 (Fig. 2).

Figure 2:

Principal Component Analysis (PCA) performed on populations of Tylenchulus semipenetrans collected from citrus orchards in Fars province, focusing on the morphometric characteristics of the second-stage juveniles (A) and males (B).

To extract genomic DNA, a single J2 of T. semipenetrans was hand-picked from the DESS solution (Yoder et al., 2006) using a pricking needle under a stereomicroscope and then washed thoroughly in distilled water. The nematode was placed in a 1.5-mL Eppendorf centrifuge tube containing 20 μL proteinase K (600 μg/ml) and 180 μL tissue lysis buffer (ATL) (two replicates for each population) and incubated overnight at 56°C. The extraction process was then continued according to the manufacturer’s instructions for the DNeasy Blood and Tissue Kit (QIAGEN, Hilden, Germany). The extracted DNA was either used directly or stored at −20° C in the DNA database of the Nematology Laboratory of the Czech University of Life Sciences, Prague.

For the analysis of the genetic diversity of T. semipenetrans, the DNA Sanger sequencing technique was performed. For this purpose, ITS and the D2-D3 segments of the 28S rDNA and COI mtDNA genes of 55 populations of the citrus nematode isolated from different districts were amplified (partial/complete), purified, and sequenced. For amplification of ITS rDNA (including ITS1, 5.8S and ITS2 regions), a universal primer pair of forward 18S and reverse 21S, as described by Marek et al. (2010), and a species-specific primer pair of forward Ts2-IF (TTCGAGAAACTTGGGGATTGGC) and reverse Ts2-IR (CAGGGACCTATGATCAAGTGCT) [presented in our study] were used. Since the amplification of the ITS rDNA of T. semipenetrans with the 18S and 21S primer pair was inefficient (Supplementary Fig. 1. A), this primer set was used as a basis for the development of a new primer pair (Ts2-IF: TTCGAGAAACTTGGGGATTGGC and Ts2-IR: CAGGGACCTATGATCAAGTGCT) with a smaller amplification product [770 bps] (Supplementary Fig. 1. B). The newly designed primer pair was developed based on the alignment of the sequences in the Vector NTI software (ThermoFisher Scientific) and the NCBI Primer-BLAST Tool (National Library of Medicine, Bethesda, MD). In addition, the designed primer pair fulfilled almost all criteria on the PCR Primer Stats website to be considered a suitable primer set (e.g., high GC content, low self-compatibility rate). To test the new primer set, the ITS rDNA sequences of common plant parasitic nematode species occurring in citrus orchards were obtained from NCBI and aligned in Mega 7 (https://www.kent.ac.uk/software/mega-7). The sequences were then searched for annotated primers. The forward primer (Ts2-IF) matched completely with two isolates of the citrus nematode, but at least five mismatches were detected for other taxa (including Hemicycliophora sp., Mesocriconema sp., Paratylenchus sp., Gracilacus sp., Tylenchorhynchus sp., Hoplolaimus sp., and Meloidogyne sp.) (Supplementary Fig. 2A). The result of the reverse primer search (Ts2-IF) in the ITS sequences of the other species showed that there were at least nine mismatches between the reverse primer sequence and the target sequences. Therefore, the ITS region of the above-mentioned taxa was not amplified with this primer set (Supplementary Fig. 2B). The sequences of the primer sets used in this study are listed in Table 2.

The total volumes of all PCR reactions were 20 μL, and contained 1 μL of DNA template; 0.4 μL forward and reverse primers mix 50 pmol μL⁻¹ (Table 2); 10 μL Phusion HSII High Fidelity PCR Master Mix (17 μL ddH₂O; 2.5 μL 10× buffer; 1.5 mM MgCl₂); 1 U Taq DNA polymerase (200 μM each dNTP); and sterile ddH₂O, added to a final volume of 20 μL. A negative control (without a DNA template) was also included in all reactions. Amplification was performed using a C1000 Touch Thermal Cycler (Biorad, Hercules, CA, USA). The thermocycling profile for the ITS rDNA gene consisted of an initial hot-start denaturation at 98°C for 30 s; 35 cycles of denaturation at 98°C for 10 s; 30 s of annealing [at 58°C with the 18S/21S primer pair and 56.5°C, using Ts2-IF/Ts2-IR primer pair]; and 72°C extensions for 35 s, followed by 72°C for 8 min to complete the process. The forward primer COIF5 and the reverse primer COIR9 (Powers et al., 2014) were used for amplification of the COI mtDNA gene (Table 2). PCR conditions for the COI gene consisted of 5 min at 94°C, followed by 40 cycles of 30 s at 94°C, 30 s at 48°C, and 90 s at 72°C, with a final extension at 72°C for 5 min. D2–D3 expansion of 28S rDNA was amplified using D2A and D3B (Subbotin et al., 2006) as forward and reverse primers, respectively (Table 2); the PCR program for this fragment included an initial treatment at 98°C for 30 s, 35 cycles at 98°C for 10 s, 58.4°C for 30 s, and a 72°C elongation step for 35 s, followed by an extension at 72°C for 8 min.

PCR products were run on a 1% TAE-buffered agarose gel (stained with ethidium bromide), visualized, and photographed under a UV transilluminator (80 V, 50 min). The high-volume PCR products were amplified (in 40 μL), run on a 1.5% TAE-buffered agarose gel, cut from the gel, and then purified using the GeneJET Gel Extraction kit (Fisher Scientific Lithuania, Vilnius) according to the manufacturer’s instructions. The proper fragments were sent to Eurofins Genomics (Ebersberg bei München, Germany) for sequencing in both directions with the corresponding primers. All newly obtained sequences were submitted to the GenBank database under the accession numbers listed in Table 1.

Before the alignment, the new DNA sequences obtained in this study were manually trimmed (if necessary) using Mega 7.1.0, and contiguous strands (contigs) were obtained using the online CAP3 sequence assembly program. A simple local alignment search tool in NCBI was used to check the species identity of the DNA sequences. Three datasets were created (one for each of the three sequenced genes) that contained the sequences of all collected populations, as well as closely related taxa. Taxa were selected based on the literature (Holterman et al., 2006; Subbotin et al., 2005, 2006; Bert et al., 2008; Rashidifard et al., 2015b). The available sequences and outgroups were aligned using MUSCLE (Edgar, 2004) implemented in Geneious Prime 2021.2.2 (www.geneious.com). According to jModeTest 2.1.10, General Time Reversible with a gamma distribution (GTR + G) was the best-fitting nucleotide substitution model for the ITS dataset and GTR was the best model for the D2-D3 and COI datasets (Darriba et al., 2012).

Bayesian analysis was performed using MrBayes version 3.1.2 (Ronquist & Huelsenbeck, 2003), which was used in Geneious Prime. For each of the datasets, the chain runs for 3 ×10⁶, after discarding 25% as burn-in samples. The Markov chain Monte Carlo (MCMC) method was used to estimate the posterior probabilities (PB) of the phylogenetic trees (Larget & Simon, 1999) using the 50% majority rule. A principal component analysis was performed in R version 3.5.1. to assess the dissimilarity between the sequences and the population measurements.

For haplotype analysis, the obtained consensus sequences were aligned with sequences from GenBank (if present) and screened for the presence of single nucleotide variations (SNV) and/or single nucleotide polymorphism (SNPs) using DnaSP Version 5.10.01 (Rozas et al., 2010). Further, nucleotide and haplotype diversity and other characteristics of the groups were evaluated. Moreover, the TCS haplotype network analysis and haplotype genealogy graphs, constructed using COI mtDNA genes and ITS sequences, were generated utilizing Hapsolutely version 0.2.2 (Clement et al., 2002; Vences et al., 2024).

The comparative morphometric data of males and J2 females of 30 populations of T. semipenetrans are shown in Tables 4 and 5. The characteristics were generally consistent with those reported for a population of T. semipenetrans from Florida (Inserra et al., 1988) and the previously reported population from Fars province (Rashidifard et al., 2015b). However, the populations in the present study differ from the Florida population and the other population from Fars province by the following morphometric indices means, respectively: slightly shorter male stylet (8.41 vs. 9.3 and 9 μm); shorter J2 stylet (11.4 vs. 12.3 and 13 μm); shorter male body length (347 vs. 362 and 368 μm); and shorter male tail (36 vs. 39.9 and 41.4 μm). Moreover, the body length of J2 was shorter than in the Florida population (332 vs. 363 μm) and slightly longer than in the Fars population (332 vs. 313 μm). Morphological identification of populations from Fars province in the present study was also verified by species-specific PCR using ribosomal DNA from J2. Moreover, the results of morphometric analysis of the characters based on Principal Component Analysis (PCA) showed that some populations have slightly different morphometry – for example, males and J2 females of populations 772 (Darab, orange trees) and 720 (Jahrom, sweet lemon) were located marginally apart from other populations. Also, males of 789 and 746 isolates and J2 of 777 isolates stood far from the other isolates (Fig. 2).

In the present study, a total of 134 new sequences of T. semipenetrans were obtained, including 48 sequences each of ITS rDNA, D2-D3 of 28S rDNA, and COI mtDNA (Table 1). Based on BLAST search, the sequences of ITS and D2-D3 rDNA matched the sequences of T. semipenetrans available in GenBank, with similarities ranging from 96.3% to 100%. There was no deposited sequence of COI mtDNA for the citrus nematode in the database, so the sequences were aligned with other taxa in the superfamily Criconematoidea – e.g., Hemicriconemoides macrodorus (KM577167), Hemicriconemoides promissus (KM577164), Ogma seymouri (MN711327 & MN711328), and Crossonema menzeli (MN710911 & MN710914), with a similarity level higher than 86%. The results indicate a low variability among the analyzed isolates.

The electrophoretic separation of the amplified ITS rDNA resulted in a single product size of 770 bp for all samples. After aligning the sequences of the ITS dataset, 597 bp were used for phylogenetic analysis. Alignment of the ITS consensus sequences of our T. semipenetrans populations with a reference sequence of T. semipenetrans from GenBank (JN112270) revealed the presence of 17 single nucleotide variations (SNVs), 11 of which were identified as SNPs (Table 3). The results showed that on average one nucleotide variation occurred per 31 base pairs. Additional information about the DnaSP analysis and the location of the mutations are shown in Tables 3 and 6. The most frequent mutation was the transition from cytosine to thymine (with seven transitions located at nucleotides 92, 151, 161, 225, 405, 423 and 523), followed by the transition from guanine to adenine (three transitions situated at nucleotides 37, 46 and 139) (Table 6).

The phylogenetic relationships of T. semipenetrans populations based on the ITS rRNA gene are shown in Figure 3. Isolates of five known Tylenchulus species were well-delineated in the tree and formed five distinct subclades. The main clade, which included all species of Tylenchulus, was closely related to another clade that included two isolates of Trophotylenchulus floridensis (JN112261 and JN112262) (PB: 0.73). All newly obtained sequences of T. semipenetrans grouped with those of the same species from GenBank in a strongly supported subclade (BP: 0.97), and T. musicola was the closest species to T. semipenetrans.

Figure 3:

Bayesian phylogenetic tree of Tylenchulus semipenetrans isolates from citrus orchards in Fars province based on the ITS of rDNA sequences, analyzed under the General Time Reversible with a gamma distribution (GTR + G) model. Numbers at nodes are posterior probability values. Sequences with codes in the parentheses generated in this study. The codes in parentheses indicate the haplotype of the relevant population.

PCR of the D2-D3 fragments of the 28S rDNA yielded amplification products 774–777 bp in length which, was obtained from a single J2 of T. semipenetrans. Comparison of 48 sequences of this gene in the present study with a reference sequence of T. semipenetrans from GenBank (KM598334.1) resulted in an alignment 676 bp in length. The alignment showed 24 segregating sites, or SNVs, of which 12 were SNPs (Table 7). In other words, there was an average of one segregating site per 28 bp of the D2D3 28S sequences (Table 3). The most common SNV was cytosine-to-thymine transition (seven, at positions 55, 84, 157, 260, 263, 355, and 497), followed by thymine-to-cytosine transition (four, at positions 59, 94, 238, and 502) and guanine-to-thymine transversion (four, at positions 34, 119, 340 and 369).

Phylogenetic studies for the citrus nematode isolates based on D2-D3 rDNA sequences are shown in Figure 4. Similar to the ITS tree, all isolates of T. semipenetrans formed a separate, maximally supported subclade (BP: 1.00) from the other known species of the genus. Tylenchulus musicola was the closest species to T. semipenetrans based on the D2-D3 sequences. The main clade, which included all species of Tylenchulus, was closely related to another clade comprising Trophotylenchulus floridensis isolates (JN112253 and JN112254).

Figure 4:

Bayesian phylogenetic tree of Tylenchulus semipenetrans isolates from citrus orchards in Fars province based on D2-D3 28S rDNA partial sequences, analyzed by the General Time Reversible (GTR) model. The numbers shown at the nodes are posterior probability values. Sequences with codes in the parentheses generated in this study. The codes in parentheses indicate the haplotype of the relevant population.

This study provided the first sequences of the cytochrome oxidase subunit 1 (COI) gene for T. semipenetrans from Iran. The PCR reaction for the COI mtDNA gene of the citrus nematode populations yielded a single fragment of 790-bp nucleotides. The 603-bp alignment of 48 COI mtDNA gene sequences of our T. semipenetrans isolates revealed 16 single nucleotide variations, corresponding to an average of one variation per 38 bp. Out of 16 segregating sites (or SNVs), 11 were identified as SNPs (Tables 3 and 8). The most common mutation was the transition of adenine to guanine, with five transitions at nucleotides 406, 481, 515, 535 and 574.

The phylogenetic relationships of T. semipenetrans populations based on the COI mtDNA gene revealed that all isolates of T. semipenetrans formed a strongly supported clade (BP: 1.00); this clade was in a well-supported sister relationship with two unidentified species of Trophotylenchulus (MN711381 and MN711382) (BP: 0.97). No COI sequence was available for other known species of Tylenchulus (Fig. 5).

Figure 5:

Bayesian phylogenetic tree of Tylenchulus semipenetrans isolates from citrus orchards of the Fars Province based on COI mtDNA partial sequences, analyzed under the General Time Reversible (GTR) Model. Numbers shown on nodes are posterior probability values. All T. semipenetrans sequences were produced in this study. The codes in parentheses indicate the haplotype of the relevant population.

Based on the occurrence of SNV among the sequences of T. semipenetrans, haplotypes corresponding to each of the genes are presented in Tables 6–8. Preliminary principal component analysis of the SNV level of the haplotypes based on all genes was performed separately. PCA analysis revealed distinct groupings within the COI mtDNA gene (Fig. 6). In contrast, the ITS and D2D3 genes did not exhibit clear clustering patterns, so data for these genes are not presented. The haplotypes of T. semipenetrans based on the COI mtDNA gene were TsI, TsII, TsIII, TsIV, TsV, TsVI, TsVII, TsVIII, TsIX and TsX. The TsI, TsIII and TsII haplotypes were more frequent than others (Table 8 and Fig. 5). The concordance between the PCA biplot and the corresponding phylogenetic tree of haplotypes was observed (Fig. 6). This congruence supported the clustering of TsII, TsIII, TsIV, and TsV, while other haplotypes formed a distinct group (Fig. 6). The TCS network analysis also indicated that there was little variation among the samples (Supplementary Fig. 3).

Figure 6:

A: principal component analysis (PCA) generalized linear modeling of Tylenchulus semipenetrans haplotypes from citrus orchards of the Fars province based on COI mtDNA, and B: the corresponding phylogenetic tree, analyzed under the Hasegawa Kishino Yano (HKY) model in MEGA 7. Numbers shown on nodes are posterior probability values.