Molecular identification of a Heterorhabditis entomopathogenic nematode isolated from the northernmost region of South Africa

Entomopathogenic nematodes (EPNs) are among the frequently used and commercialized bio-pesticides. However, they are restricted in their infectivity, persistence, storage, and cost of production. One of the methods used to improve this is a continuous search for new isolates with significant behavioral and physiological characteristics. A novel EPN isolate, Heterorhabditis zealandica strain ETL, isolated from South Africa (GPS co-ordinates − 24.849721 and 28.336980) is described and studied against late-instar of Galleria mellonella (L.) and Tenebrio mollitor (L.) larvae. The morphological and molecular studies indicated this isolate as a Heterorhabditis strain. The comparison of sequences of the internal transcribed spacer (ITS) region, the 18S rRNA gene, and the D2-D3 region of the 28S rRNA gene with available sequences of other described species within the genus indicate the isolate as a new Heterorhabditis zealandica strain. The phylogenetic analysis of the sequence data places strain ETL, closest to H. zealandica strain Bartow (GU174009) in the Heterorhabditis group. This EPN was lethal to G. mellonella and T. mollitor; as infections occurred within 24–96 h. Fifty percent of the larvae population were killed within 24 h and 100% after 96 h.


Background
Entomopathogenic nematodes (EPNs) are highly pathogenic to numerous insect pests (Malan and Ferreira 2017). Effective biocontrol agents of insects have been mainly isolated from the genera Steinernema and Heterorhabditis, which are mutually associated with bacteria of the genera Xenorhabdus and Photorhabdus, respectively (Malan and Ferreira 2017). Recently, Oscheius EPN genus associated with symbiotic bacteria of Serratia genus was reported (Torres- Barragan et al. 2011). To date, EPNs are used and commercialized as bio-control agents of insect pests in large-scale agricultural and individual home gardens in many countries (Azazy et al. 2018). South Africa (SA) demonstrates a rich EPN fauna with biocontrol potential (Hatting and Malan 2017). Few data on the diversity and activity of EPNs from different provinces of SA are available.
The EPN characterization is particularly useful to (1) match the appropriate EPN species with a targeted pest (or environmental conditions) for an effective biocontrol potential and/or integrated pest management (IPM) programs (Peat et al. 2009), (2) to prevent the importation of species (and/or strains) that can outcompete native species and compromise the local biodiversity, and (3) to select species with beneficial properties for commercial bio-prospecting purposes, such as the use of patents or other legal issues (Peat et al. 2009). Therefore, identifying the EPN species accurately is critical.
The initial morphological characterization of the EPN taxonomy is difficult because of the relatively conservative morphological traits of nematodes (Peat et al. 2009). Molecular systematics, compounded by the improvement of both polymerase chain reaction (PCR) and DNA sequencing techniques, have proven adequate and played a key role in advancing the EPN taxonomy, which contributed to the expansion of knowledge on EPN biodiversity, geographic distributions, host ranges, ecology, behavior, and co-evolution (Peat et al. 2009). Molecular characters are considered as the most suitable approach for species identification and systematics of nematodes, especially for taxonomic ambiguities, identification of members of a species complex, and the differentiation of morphologically similar species (Peat et al. 2009). Thus, they have been widely used (Peat et al. 2009;Campos-Herrera 2015).
Sequence data from nuclear and mitochondrial genetic loci, and whole genomes (more recently) are used in EPNs studies, as well as nuclear ribosomal DNA (rDNA) genes, partial messenger RNA (mRNA) copies of protein-coding genes (expressed sequence tags (ESTs)) are used to infer the phylogenetic relationship (Peat et al. 2009). Among these, nuclear genes have been used extensively; each has numerous copies and evolves at different rates. Their variable and conserved regions render it possible to differentiate taxonomic levels and delimitate the nematode taxa (Stock 2009). They include the external non-transcribed spacer (NTS); the small subunit (SSU) or 18S; the internal transcribed spacer (ITS) region that separates rDNA-coding regions (Adams et al. 2006) composed of ITS-1, 5.8S, and ITS-2; and the large subunit (LSU) or 28S (Stock 2009).
The rRNA gene sequences can vary in length and nucleotide composition, insertion, and deletion events, sometimes involving blocks of multiple nucleotides occur frequently (ITS sequence length differences of > 100 bp can be observed), and can result in rDNA size differences between sequences (Peat et al. 2009). These, unlike protein-coding genes, do not compromise the function of the ribosome (Nguyen et al. 2001), but can render dubious the homological position of characters (useful for delimitating species) during the phylogenetic reconstruction (Peat et al. 2009). Thus, inaccurate tests of homology statements during alignment, compounded by a high number of taxa (which drastically increases the number of possible phylogenetic solutions) can lead to spurious phylogenies (Peat et al. 2009). The aim of this study was to isolate and identify a native EPN species that can be suitable as a biocontrol agent by testing it against Galleria mellonella L. and Tenebrio mollitor L. larvae under laboratory conditions.

Soil sampling and isolation of entomopathogenic nematodes
Eighty-seven soil samples were collected from both agricultural disturbed and natural undisturbed soils from the locality of Bela-Bela in Limpopo province of South Africa according to the method of Kaya and Stock (1997). Accordingly, samples were taken from a depth of up to 20 cm and placed in plastic containers (2000 ml) with lids and transported to the laboratory for processing. In the laboratory, isolation of EPNs from soil samples was performed, following the live bait method-using laboratory-reared late-instar larvae of the greater wax moth, Galleria mellonella L. (Bedding and Akhurst 1975). Briefly, after adjusting the soil humidity to about 10% using sterilized water (to encourage the potential EPNs movements) (Ruan et al. 2018), 10 late-instar larvae were placed in each container, and the container was inverted. The containers were incubated at room temperature (± 25°C) and checked daily for larval infection and death based on color change. Dead larvae were collected and cultured according to the White-trap method (White 1927). Infective juveniles (IJs) emerging from the cadavers of G. mellonella were collected from the White-traps and used for re-infecting fresh larvae to ensure the maintenance of the nematode's population, and some of the IJs were used for the subsequent determination of the pathogenicity against G. mellonella and T. mollitor larvae. The most virulent isolate (isolated from GPS co-ordinates − 24.849721 and 28.336980) was selected for subsequent analyses.

EPN pathogenicity test
Nematode initial pathogenicity was evaluated according to Glazer and Lewis (2000), using late instar of G. mellonella and T. mollitor larvae in a Petri dish with 10-cm diameter lined with 2 layers of moist filter paper (Whatman No. 1). Accordingly, each Petri dish contained 10 g of sterile sandy loam soil, four live larvae, and approximately 250 of IJs, and control dishes had only distilled water instead of IJs. The Petri dishes were incubated at 25°C in the dark and observed daily for mortality of the larvae for 96 h. Dead larvae were collected daily and placed on White-traps for collection of emerging IJs. The experiment was conducted in quintuplets and repeated 4 times.

Morphological characterization of EPN
An initial morphological characterization of EPN IJs (Stock and Hunt 2005) was done using a scanning electron microscope (SEM). Accordingly, IJs and 1st-generation adults (males, hermaphrodite, and/or amphimictic females) were heat-killed at 60°C for 2 min in Ringer's solution. Infective juveniles were rinsed with Ringer's solution 3 times with a 5-min interval between each rinse. These were then fixed in 8% glutaraldehyde overnight (glutaraldehyde 25% EM grade, diluted in Ringers' solution). For post fixing, these were rinsed by sterile distilled water 3 times and dehydrated with 30, 50, 70, 90, 95, and 100% ethanol at 10 min intervals. These were then critically point dried with CO 2 for 2 h and mounted on SEM stubs, coated with gold and palladium, and viewed using TESCAN VEGA3 microscope (TESCAN ORSAY HOLDING, a.s., Brno -Kohoutovice, Czech Republic).

Molecular characterization of EPN Amplification and sequencing
The EPNs were inbred 13 times to eliminate heterozygosity and the total genomic DNA extracted. Protocol from ISOLATE II Genomic DNA Kit (Bioline USA Inc., 305 Constitution Drive, Taunton, MA, 02780, USA) was followed to extract and purify the total genomic DNA. To identify the EPN, 18S rRNA gene, internal transcribed spacer (ITS) region, and the D2-D3 region on the 28S rRNA gene were used (Nguyen and Hunt 2007). The ITS region was amplified using universal primers TW81 and AB28 (Joyce et al. 1994), 18S rRNA gene using primers (18S-26S) by Subbotin et al. (2006) and primers (D2A-D3B) described by Vovlas et al. (2006) were used for D2-D3 region on the 28S rRNA gene. The PCR solution had a final volume of 25 μl, including 12.5 μl of 2× PCR Master Mix, 4 μl genomic DNA, 3 μl of each primer with 10 μM concentration, and 2.5 μl of nuclease-free water. The PCR reactions were conducted using a Bio-Rad T100™ Thermal Cycler, (Bio-Rad Laboratories, Inc., Hercules, California 94547, USA). This consisted of an initial denaturation at 94°C for 4 min, followed by 35 cycles of 94°C for 1 min, 50-55°C for 30 s, and 72°C for 2 min; and a final extension at 72°C for 10 min. The PCR product was electrophoresed at 80 V for about 30 min on agarose gel (1% v/v: agarose/1X TBE [Tris base, boric acid and EDTA]) buffer stained with ethidium bromide (EtBr). PCR products were purified and concentrated using the QIAquick PCR Purification kit before sending for sequencing at a commercial service provider, Inqaba Biotechnical Industries (Pty) Ltd., Pretoria, South Africa.

Phylogenetic analysis of the EPNs
BioEdit v 7.2 (Hall 2005) software was used for editing and creating the consensus sequences. The consensus DNA sequences for each marker were searched against the NCBI (National Center for Biotechnology Information) nucleotide collection database, using BLAST (Basic Local Alignment Search Tool) for closely related species with high similarity percentages as well as low e values. All phylogenetic analyses post-BLAST were performed in MEGA 7 (Kumar et al. 2016). The selected species from the BLAST search and present isolates in the study were aligned, using MUSCLE (MUltiple Sequence Comparison by Log-Expectation) (Edgar 2004). Maximum likelihood trees of the homologous sequences obtained were inferred with 1000 replications, using the Jukes-Cantor model (Jukes and Cantor 1969). Bootstrap analysis with 1000 replicates was used to determine branch support (Nguyen and Hunt 2007).

Statistical analysis
The EPN pathogenicity method was statistically validated by a one-way analysis of variance (ANOVA) using IBM SPSS Statistics 22 (SPSS/IBM, Chicago, IL, USA). Mean values of the number of dead G. mellonella and T. mollitor larvae recorded daily post-EPN application were compared. Significant differences at p ≤ 0.05 level of probability were reported.

Isolation and characterization of EPNs Isolation and morphological characterization
EPNs were isolated in four samples (5%). Of the 4 samples, only one sample displayed complete virulence against the larvae allowing subsequent IJ production for further study. Although this represents a low recovery, this study highlights the need for more EPN surveys and studies in the northern parts of SA. From the scanning electron microscope, the isolated IJ had similar characters as those described by Stock and Hunt (2005) within the family Heterorhabditidae. Accordingly, the cuticle had longitudinal ridges throughout most of its body length and a tessellate pattern in the anterior-most region, lateral field with 2 ridges, head with prominent cuticular dorsal tooth, excretory pore located posterior to the basal bulb, and tail short, conoid and tapering to a small spike-like tip (Additional file 1: Figure S1).

Molecular characterization
The sequences of the ITS, 18S rRNA gene, and the D2-D3 region of the 28S rRNA gene were analyzed. The BLAST analysis of the ITS consensus gene (819 bp) attributed 99.88 and 99% of the similarities and query coverage with both H. zealandica strain Bartow (GU174009) and H. zealandica isolate WS20 (KY055373), and 99.51 and 99% of similarities and query coverage with both H. zealandica strain NZH3 (EF530041) and H. zealandica (AY32148) ( Table 1) The name Heterorhabditis zealandica strain ETL was assigned to the isolate, and the consensus sequences were submitted to NCBI-GenBank database and accession numbers MH443371 (for ITS), MH443381 (for D2-D3 region on 28S rRNA gene), and MN341010 (18S rRNA gene) were assigned.
The evolutionary history was inferred by using the maximum likelihood method, using Jukes-Cantor model (Jukes and Cantor 1969). Trees with the highest loglikelihoods (− 310.50, − 409.01, and − 1505.68 for the ITS, 28S, and 18S gene, respectively) are shown in Fig. 1. The percentage of trees in which the associated taxa clustered together is shown next to the branches. Initial trees for the heuristic search were obtained automatically by applying Neighbor-Join and BioNJ algorithms to a matrix of pairwise distances estimated using the maximum composite likelihood (MCL) approach and then selecting the topology with superior log likelihood value. The trees were drawn to scale, with branch lengths measured in the number of substitutions per site. The analysis involved 18, 15, and 18 nucleotide sequences for the ITS, 28S, and 18S taxonomic study, respectively. Codon positions included were 1st + 2nd + 3rd + noncoding. All positions containing gaps and missing data were eliminated. There were 152, 131, and 412 positions in the final dataset of the ITS, 28S, and 18S taxonomic studies, respectively.
The 18S rRNA sequence was detected to be suboptimal at resolving the taxonomic conflicts at species level, as it had fewer taxon, which is reported to evolve at a slower rate (Stock 2009). Reportedly, this is because of its slow evolutionary rate (Stock 2009) and has been used to differentiate the monophyletic origins of nematode groups (Peat et al. 2009). The 5.8S rRNA region of the ITS is short and highly conserved than the ITS-1 and ITS-2 regions (Stock 2009), but evolve more rapidly than the 18S and 28S genes. Thus, they are ideal (i) for the EPN taxonomic studies at species (population) levels and (ii) for population genetic studies (Stock 2009). The ITS-1 region is reportedly sufficient at differentiating species and assessing their evolutionary relationships, particularly among the Heterorhabditis spp. (Stock 2009;Peat et al. 2009). The 28S rRNA varied more rapidly than the 18S rRNA and had fewer positional ambiguities during alignment than ITS (Stock 2009). However, the 28S rRNA has been attested more informative and suitable for the assessment of phylogenetic relationships, delimitation of terminal taxa, and for diagnostic purposes among Steinernema spp. (Stock and Hunt 2005;Stock 2009).

Entomopathogenic nematode pathogenicity test
The EPN pathogenicity test was repeated (4×), and a number of dead G. mellonella larvae (mean ± RSD) 1.72 ± 0.5, 2.75 ± 0.5, 3.5 ± 0.5, and 4 ± 0 (p 0.006) were collected daily after EPN application, respectively. The data are presented in Table 2. No statistical difference was found with the G. mellonella quadruplet EPN pathogenicity experiments. The scatter-gram and line of best fit were obtained from the time vs. mortality of larvae for each experiment (Fig. 2a). A coefficient of linearity (r) of 0.95 was obtained for all replicated experiments. As for G. mellonella, the number of dead T. mollitor larvae (mean ± RSD) 1 ± 0, 2 ± 0, 2.75 ± 0.5, and 4 ± 0 (p > 0.001) were recorded at the 1st, 2nd, 3rd, and 4th day post-EPN application, respectively (Table 2). Again, a non-statistical difference was found. The scatter-gram and line of best fit (Fig. 2b) showed similar r values (> 0.96) for all replicated experiments. Heterorhabditis zealandica strain ETL was capable of infecting G. mellonella and T. mollitor within 24-96 h (Fig. 2a, b), which was within the criterion (5 days) for entomopathogenicity suggested by Dillman et al. (2012b). About 50% of the larvae population was killed within 24 and 48 h for G. mellonella and T. mollitor, respectively. After 96 h, mortality reached 100% and the nematodes emerged from the cadaver 2-3 days after death. The color change of dead larvae observation (Fig. 2c) disclosed a biocontrol potential of H. zealandica strain (See figure on previous page.) Fig. 1 Maximum likelihood phylogenetic tree of Heterorhabditis zealandica strain ETL and its closest neighbors based on a the ITS region, Caenorhabditis elegans strain RBS-CE KX572972 was included as an out-group; b D2-D3 region of the 28S rRNA gene, C. remanei strain EM464 AY602174 was included as an out-group; and c the 18S rRNA gene, C. elegans NR_000055 was included as an out-group. The numbers on nodes indicate bootstrap values after 1000 replicates expressed in percentages ETL against insects. This could be due to synergetic effects of cell clumping, protease activities, crystal protein production, and/or anthraquinone production (pigmentation) by EPN symbiotic bacteria within insect cadavers (Eckstein and Heermann 2019). This is an interesting observation and it revealed a biocontrol potential of H. zealandica strain ETL against insects from the Lepidoptera and Coleoptera orders and possible use against other insects. One approach to enhance the EPN efficiency in fields has been isolating EPN species from different areas that could display a significant behavioral and physiological adaptation (Lu et al. 2016). The recovery of native H. zealandica strain ETL from the province of Limpopo, an unexplored region in South Africa, further demonstrates the richness of South African EPN fauna with biocontrol potential (Hatting and Malan 2017), suggesting further research and extensive sampling.
The species H. zealandica strain ETL was successfully characterized based on morphological characters of the IJ and based on its molecular sequences of three (ITS, 18S, and 28S) genetic markers. The usefulness of molecular data to distinguish these species was documented (Campos-Herrera 2015). The comparison of the ITS, 18S, and 28S rDNA sequences with available sequences of other described species indicated the isolate as an EPN isolate among the Heterorhabditis group. The phylogenetic analysis data to infer the relationships of the isolated strain with other Heterorhabditis found that it was in a polyphyletic clade with 100% support (Fig.  1c) and was placed closest to H. zealandica strain Bartow (GU174009) in the Heterorhabditis group (Fig. 1b) supported with 99.88 and 99% of similarities and query coverage of a 819-bp sequence. This was demonstrated by the 18S rRNA sequence (Fig. 1c) and ITS rRNA sequence (Fig. 1a), respectively, while the 28S rRNA sequence ( Fig. 1b) provided more information on its proximity to strain Bartow (GU174009). Moreover, the 18S rRNA sequence (Fig. 1c) revealed the proximity of strain ETL to other strains from Thailand, and other H. zealandica isolates recovered from Australia, Lithuania, New Zealand, Russia, and USA) confirm the widespread occurrence of these strains; although they have been rarely recovered in SA (Malan et al. 2011). Future surveys are needed to better assess the distribution of these strains in South Africa. Other H. zealandica strains isolated from neighboring provinces within South Africa (De Waal et al. 2018) such as Mpumalanga, KwaZulu-Natal, Eastern Cape, and Western Cape provinces demonstrate a behavioral and physiological adaptation that could be beneficial for their bio-control potential. Strain ETL was lethal to G. mellonella and T. mollitor; infections occurred within 24-96 h. Fifty percent of the larvae population were killed within 24 h and 100% after 96 h.

Conclusion
This is among the first documented record of species of EPNs isolated from the Limpopo Province of South Africa. It was identified as Heterorhabditis zealandica strain ETL and is a potential biocontrol agent against insects from the Lepidoptera and Coleoptera orders, thus indicating a potential use against other agricultural industry key insect pests.
Additional file 1: Figure S1. Morphological features of the isolate showing a tessellate pattern in anterior most region (a). Lateral field with two ridges (b). Prominent cuticular dorsal tooth present (c, d, and e).