Genome Characterization and Development of Real-Time PCR Assays for Ditylenchus dipsaci and D. weischeri

Nematodes were identified morphologically. The following nematode isolates, as mixed-stage populations, were previously collected and propagated in culture (summarized in Table 1): three nematode isolates were identified as D. dipsaci (C-100, E-105, G-137), two nematode isolates were identified as D. weischeri (O-100, S-100), two nematodes isolates were identified as D. destructor from (O-101, J-100), one nematode isolate was identified as D. africanus (D-afri-South-Afr), one nematode isolate was identified as a Ditylenchus sp. (Dity-sp-Jord-ONT), and one nematode isolate was identified as Litylenchus crenatae (Lity-cren-ONT-01). Mixed-stage nematodes were collected using the Baermann funnel technique (Baermann, 1917).

DNA was extracted using QuickGene-810 instrument with QuickGene DNA Tissue Kit S (Kurabo, Osaka, Japan), following the manufacturer’s protocol with a modification: prior to overnight incubation in lysis buffer, tissue lysis buffer MDT was added to the nematodes in sterile 1.5-ml tubes, and the nematodes were crushed using a micropestle for 1–2 min. The quality and quantity of DNA were determined using agarose gel electrophoresis, Qubit^® 2.0 Fluorometer (Invitrogen, Carlsbad, CA), Nanodrop^® ND-1000 UV/ Vis spectrophotometer (NanoDrop Technologies, Wilmington, DE), and Agilent 4200 TapeStation (Agilent Technologies, Santa Clara, CA).

Additionally, DNA was extracted from individual nematodes based on the proteinase K method described in Kumari and Subbotin (2012). One, five, or ten nematode specimens were put in 15 μl of DNA extraction buffer, containing 0.25 mg proteinase K (BP-1700-500, Thermo Fisher Scientific, Waltham, MA) per 1 ml of 1χ PCR buffer (BP6112, Thermo Fisher Scientific). The nematodes were crushed under a dissecting microscope. The suspension with nematodes was then transferred into a sterile 1.5 ml centrifuge tube containing 15 μl of DNA extraction buffer. Tubes were incubated at 65°C for 1.5 hr, followed by an incubation at 95°C for 15 min to inactivate proteinase K. All nematode isolates were barcoded using 28S and ITS to confirm species identity (data not shown).

Nanopore and Illumina sequencing was performed on four of the isolates: D. dipsaci C-100, D. dipsaci G-137, D. weischeri O-100, and D. weischeri S-100. For nanopore sequencing, the DNA was sheared to 15–30 kb using Covaris g-TUBE (Covaris, Woburn, MA). The SQK-LSK109 kit was used to prepare the library from genomic DNA following the manufacturer’s instructions (ONT, Oxford, UK). Nanopore sequencing was performed on the Oxford Nanopore MinION R9.4 flow cell for 48 hr, producing “1D” reads. Basecalling was performed during the run using the high-accuracy calling option in Guppy v. 3.0.3. Paired-end Illumina sequencing from genomic DNA of all four isolates was performed on an Illumina NextSeq 550 Instrument at the Molecular Technologies Laboratory (MTL) at the Ottawa Research and Development Centre, Agriculture and Agri-Food Canada. Additionally, paired-end Illumina sequencing (2χ 101 bp) from total RNA of D. dipsaci grown on plant material was performed on an Illumina HiSeq 2500 instrument at the National Research Council Canada in Saskatoon, Saskatchewan, Canada (Yu et al., 2014).

Raw nanopore reads were filtered using filtlong v0.2.0 (https://github.com/rrwick/Filtlong) where reads shorter than 5 kb (D. weischeri O-100 and D. dipsaci G-137) or 10 kb (D. weischeri S-100 and D. dipsaci C-100) were discarded. The remaining reads were assembled using canu v1.8 (Koren et al., 2017) with default settings with an estimated genome size of 200 Mb (genomeSize = 200 m). Assembled contigs were subjected to Purge Haplotigs (Roach et al., 2018), with the following parameters for read depth cut-off step: -l 1 -m 17 -h 50, to improve the haploid representation by identifying and reassigning allelic contigs and to reduce duplication arising from both haplotypes of a given part of the genome being assembled as separate primary contigs.

Sequences were verified for contaminants. Assembled contigs were searched against a local whole-genome database of plants, fungi, oomycetes, and bacteria (downloaded from trusted sources such as NCBI RefSeq, JGI Mycocosm, and Ensembl) with BLASTn v.2.2.31+ (Camacho et al., 2009) using a threshold E-value of 1 χ 10^-100. Any reported hits were manually inspected. Contigs with multiple and long matches to plants, fungi, oomycetes, or bacteria were removed from subsequent analyses.

Illumina reads were checked with fastqc v0.11.8 (ht tps://www.bioinformatics.babraham.ac.uk/ projects/fastqc/). These reads were trimmed with bbduk.sh script (ref=adapters qtrim=rl trimq=20 minlength=36 ktrim=r forcetrimleft=15 tossjunk=t) as part of the BBMap package v38.22 (https:// sourceforge.net/projects/bbmap/). The trimmed reads generated from genomic DNA were mapped to the canu assembly with BWA v0.7.17 (Li and Durbin, 2009), and errors were corrected with Pilon v1.23 (Walker et al., 2014). Genome assembly statistics were calculated with QUAST v5.0.2 (Gurevich et al., 2013).

The trimmed RNA reads were assembled with Trinity v2.8.5 (Haas et al., 2013) with default settings. The genome assembly of D. dipsaci E-105 (Mimee et al., 2019) was downloaded from NCBI. It was annotated along with the four isolates sequenced in this study. Genome annotation was performed with funannotate v1.5.2 (https:// zenodo.org/record/2576527#.YR2AcC2caJY).Sequences shorter than 1,000 bp were removed followed by repeat masking with RepeatMasker using “ditylenchus” species. Gene prediction and annotation steps were run with the nematoda_odb10 Benchmarked Universal Single Copy Orthologs (BUSCO) dataset (–busco_db nematoda_odb10); Caenorhabditis as the seed species (–busco_ seed_species caenorhabditis); 50 minimum training models (–min_training_models 50); GeneMarkES prediction set to non-fungal organism (–organism other); and the Trinity assembled transcriptome from D. dipsaci was used as transcript evidence. To evaluate the completeness of the genome annotation procedure, BUSCO analyses, using the nematoda_odb10 database, were conducted with BUSCO v5.3.2 (Manni et al., 2021) running in the protein mode.

All genome assembly and genome annotation statistics are summarized in Table 2. Raw sequence data and assemblies were uploaded to NCBI under BioProject PRJNA847237.

Following assembly and annotation of Ditylenchus genomes, orthologous group analysis was performed with OrthoFinder v2.3.12 (Emms and Kelly, 2019), using default settings with input protein sequences from D. dipsaci E-105, D. dipsaci G-137, D. dipsaci C-100, D. weischeri O-100, and D. weischeri S-100, to find species-specific genes. Using BLAST + v2.11.0 (Camacho et al., 2009), these genes were compared against the genome assemblies of the non-target species using BLASTn and an E-value threshold of 1 χ 10^-5. Genes that still have some similarities will return hits and those that are completely unique will have no hits. This process further narrowed the list of species-specific genes for D. dipsaci and D. weischeri. Next, to find sites that were conserved between the different isolates of the same species, MUSCLE v3.8 (Edgar, 2004) was used to align sequences of the genes in this list of narrowed down candidate genes. Genes with perfect or near-perfect pairwise alignment between the different isolates of the same species were selected for further consideration.

These genes were further filtered based on low similarity to other sequences within the same species. Candidates that resulted in BLAST hits with >85% sequence similarity and length of >150 bp were discarded. Using Geneious R10 (Biomatters Ltd., Auckland, New Zealand), the remaining candidates were mapped to the genome assemblies of the same species to verify that the candidates were located on a contig with several other predicted genes and that they were not located at the periphery of the contig. Finally, after all these considerations, two candidate genes for D. dipsaci and two for D. weischeri were chosen for primer and probe design. The protein products of the final candidate genes were scanned for functional domains with InterProScan on the web (https://www.ebi.ac.uk/interpro/, accessed June 21, 2022).

Primers and probes were designed as described in Nguyen et al. (2019) using Geneious R10 (Biomatters Ltd., Auckland, New Zealand). Stringent parameters were used: amplicon size between 100 bp and 150 bp, primers and probes inside an exon, primers with a melting temperature (T_m) of approximately 60°C and roughly 23 bp long, and probes with a T_m of approximately 70°C and roughly 27 bp long. Geneious DNA Fold was used to check the products for secondary structure. Out of several equivalent primer-probe sets, the ones with least secondary structure and resulting in shorter products were chosen. Primer and probe sequences are summarized in Table 3.

The optimal annealing temperature for real-time PCR primers was tested using end-point PCR with temperature gradient of 56°C to 63°C. All reactions were prepared using the Titanium Taq PCR Kit (Takara Bio USA, Mountain View, CA), with 1 μl of approximately 1 ng/μl of DNA, 2.5 μl of 10χ PCR buffer, 1.0 μl of 2.0 mM dNTPs, 0.4 μl 10 μM primers, and 0.2 μl of Titanium Taq DNA polymerase.

Real-time PCR assays were optimized for primer concentrations. Primer concentrations were tested at 100 nM, 200 nM, 400 nM, 600 nM, and 800 nM final concentrations in the reaction mix. The optimized primer concentrations are summarized in Table 3. PrimeTime^® qPCR probes with 5’ 6-FAM^TM fluorophore and ZEN–3’ Iowa Black^® FQ double quencher were used (Integrated DNA Technologies, Coralville, IA). Real-time PCR reactions of 10 μl contained 5 μl 2χ LightCycler^® 480 Probes Master (Roche Diagnostics, Laval, Canada), primers at the optimized final concentrations of 100–800 nM, probes at 100 nM, and 1 μl of template DNA diluted to approximately 1 ng/μl.

Real-time PCR conditions were as follows: polymerase activation at 95°C for 10 min, followed by 40 cycles of denaturation (10 sec at 95°C), primer annealing (30 sec at 60°C), and extension (10 sec at 72°C), with a final cooling at 40°C for 10 sec. All real-time PCR assays were run with the LightCycler^® 480 Instrument (Roche Diagnostics, Laval, Canada). C_q values were determined with the LightCycler^® 480 Software using the second derivative maximum method. Sanger sequencing of qPCR products confirmed specific amplification of the intended targets (data not shown).

Standard curves were generated to determine amplification efficiencies and reproducibility of the real-time PCR assays (Figs. 1,2). DNA from different isolates of the target species was pooled in equimolar ratios as follows: (1) D. dipsaci C-100, D. dipsaci E-105, and D. dipsaci G-137 and (2) D. weischeri O-100 and D. weischeri S-100. The pooled DNA samples were serially diluted 3-fold, starting at approximately 1 ng/μl. Three separate real-time PCR runs were carried out for each of the primer sets to assess reproducibility.

Figure 1

Figure 2

To test primer specificity, real-time PCR assays were performed using DNA from the following nematode isolates: D. dipsaci C-100, D. dipsaci E-105, D. dipsaci G-137, D. weischeri O-100, D. weischeri S-100, D. destructor J-100, D. destructor O-101, D. africanus D-afri-South-Afr, Ditylenchus sp. Dity-sp-Jord-ONT, and Litylenchus crenatae Lity-cren-ONT-01 (Table 4). Ditylenchus dipsaci, D. weischeri, and D. destructor DNA extracted using QuickGene was used at a concentration of approximately 1 ng/μl as determined by Qubit fluorometer. DNA prepared using the proteinase K method was used, from 5 nematodes each of D. africanus and Ditylenchus sp. Dity-sp-Jord-ONT, and 1 or 10 nematodes of Litylenchus crenatae Litycren-ONT-01. The specificity of the two primer sets was further verified using a mixture of DNA. Pooled DNA of each of the three species was further combined in equal amounts as follows: (1) D. dipsaci and D. weischeri, (2) D. dipsaci and D. destructor, (3) D. destructor and D. weischeri, and (4) D. dipsaci, D. weischeri, and D. destructor.

Assay sensitivity was tested using DNA extracts from 10, 5, and 1 nematode of D. dipsaci E-105 and D. weischeri S-100 (Table 5).

After basecalling, nanopore sequencing yielded approximately 8.5 Gb (1,321,464 reads) for D. dipsaci C-100, 6.8 Gb (1,147,869 reads) for D. dipsaci G-137, 3.7 Gb (813,509 reads) for D. weischeri O-100, and 7.8 Gb (1,329,588 reads) for D. weischeri S-100. Illumina sequencing of genomic DNA yielded 18.4 Gb for D. dipsaci C-100, 15.6 Gb for D. dipsaci G-137, 16.3 Gb for D. weischeri O-100, and 17.5 Gb for D. weischeri S-100 in both directions of the paired-end sequencing. The reads assembled into 2,256 contigs (228.2 Mb genome size) for D. dipsaci C-100, 5,530 contigs (239.5 Mb genome size) for D. dipsaci G-137, 1,715 contigs (177.0 Mb genome size) for D. weischeri O-100, and 2,012 contigs (196.3 Mb genome size) for D. weischeri S-100. The assemblies were not scaffolded. Mimee et al. (2019) obtained a scaffolded genome size of 227.2 Mb for D. dipsaci E-105. It appears that D. dipsaci has a larger genome size compared with D. weischeri. The sequencing coverage, determined by taking the total number of bases sequenced using Nanopore and Illumina and dividing by the assembly size, resulted in coverages ranging from 93χ to 129χ. The GC content of each assembly was comparable (~37.5%) to the D. dipsaci E-105 assembly from Mimee et al. (2019). A total of 21,403 to 27,365 putative genes were predicted. Using the protein sequences from the genome annotations, ~50%–65% complete BUSCOs were found, which is in agreement with the results from Mimee et al. (2019). Genome statistics are summarized in Table 2.

Orthologous group analysis found 1,539 single-copy genes completely unique to D. dipsaci and 2,306 unique to D. weischeri. Using BLAST, species-specific genes were compared against the genome assemblies of the non-target species to further weed out similar sequences. This process narrowed the list of species-specific genes down to 140 candidate genes for D. dipsaci and 327 candidate genes for D. weischeri. Finally, there were only 8 candidate genes for D. dipsaci and 108 candidate genes for D. weischeri with perfect or near-perfect pairwise alignment between the different isolates of the same species. After checking for the possibility of amplifying multiple targets within species, two genes for D. dipsaci (OG0014442 and OG0014782) and two genes for D. weischeri (OG0018807 and OG0020570) were chosen for primer and probe design.

Using InterProScan, no functional domains were detected in OG0014442 and OG0020570, whereas OG0014782 is a protein kinase and OG0018807 contains signal peptide at the N terminus and a non-cytoplasmic domain that spans through the centre of the protein sequence. For all primer sets, 60°C was chosen as the optimal annealing temperature, at which only single bands of the expected size were obtained, and there was no amplification in the negative controls. All primer sets produced single bands of expected size.

The lowest concentration that could be reliably detected with all the primer-probe sets was at 81-fold dilution of the DNA template, corresponding to approximately 12 pg of DNA. Average amplification efficiencies ranged between 95.9% and 100.2% (Figs. 1,2). For all the primer-probe sets, no amplification was observed in the negative control reactions. Both D. weischeri primer-probe sets amplified only in the target species and not in the non-target species tested (Table 4). Both D. dipsaci primer-probe sets amplified in the target species but also produced signal in the non-target D. destructor isolates at C_q values >35 (Table 4).

The specificity of the two primer sets was further verified using a mixture of DNA. The D. weischeri primer-probe sets amplified successfully in the mixtures containing target nematode species and not in those containing non-target species. The D. dipsaci assays additionally produced a signal with high C_q values (37-38) in mixtures containing D. destructor (Table 4). As expected, C_q values for all the primer-probe sets were higher in mixed than pure Ditylenchus DNA samples, due to the two or three-fold dilution effect of mixing.

Additionally, assay sensitivity was tested using DNA extracts from 10, 5, and 1 nematode of D. dipsaci E-105 and D. weischeri S-100 (Table 5). All the assays showed amplification from the different numbers of nematodes, with C_q values increasing with decreasing number of nematodes, as expected. Amplification from single nematodes resulted in C_q values of approximately 32 to 34 cycles, which is outside of linear dynamic range of the assays as determined from standard curves. Five nematodes resulted in C_q values of approximately 29 to 31 cycles, which is right at the edge of the lower detection limit of the assays.