Molecular Identiﬁcation and Phylogenetic Diversity of Native Entomopathogenic Nematodes, and Their Bacterial Endosymbionts, Isolated from Banana and Plantain Crops in Western Colombia

: With the increasing negative impacts on worldwide food production caused by pests, the recovery of native entomopathogenic nematodes (EPNs) is relevant, since they are adapted to local environments, entomofauna, and signiﬁcant virulence. Therefore, the present study was designed to recover and understand the phylogenetic diversity of EPNs and their associated bacterial endosymbionts, from banana and plantain crops, as alternatives for the control of weevil species. An extensive sampling of western Colombia covered 325 ha, yielding the recovery of three EPNs’ isolates (0.49% of the samples). The molecular characterization included four mitochondrial and nuclear loci, which, after merging with the sequences of 48 species, conﬁrmed the presence of Steinernema carpocapsae , the ﬁrst report of S. costaricense in South America, and monophyly in most of the Steinernema clades. The tree topologies were consistent for the nuclear loci but not for mitochondrial, probably due to the high nucleotide substitution rate, deﬁcit in the number of species available for these loci, and incomplete lineage sorting. The endosymbiotic bacteria associated with S. carpocapsae were identiﬁed as Xenorhabdus nematophila . However, the S. costaricense bacterial symbiont presented a genetic similarity to X. koppenhoeferi and X. khoisanae , which are still uncertain in their classiﬁcation. The identiﬁcation of S. costaricense in South America indicates the wide range distribution of this species in the Americas and its ability to persist in different soil types. For the ﬁrst time, EPN isolation and phylogenetic characterization are directed to plantain and banana crops. Leveraging EPNs’ diversity promises novel applications for crop protection, while the genetic resources from the bacterial endosymbionts may provide metabolites with a wide spectrum of uses, either for agricultural or medicinal purposes.


Introduction
Worldwide food production is affected by approximately 10,000 insect and 30,000 beetles weevil species [1,2], putting the world food supply at risk and representing an estimated crop loss of 18-20% of annual production, reaching an economic value of USD 470 billion [3]. These agricultural losses are more accurate in developing countries [4]. For instance, current crop losses in Latin America account for >14% [3].
Together, bananas and plantains are one of the most important crops worldwide, classified as essential for food security due to their high nutritional value [5][6][7][8]. These crops are grown in more than 135 countries, and the main producers are located in Central and South America, and the Philippines, totalizing 20 million tons of exports [9]. Colombia is one most important banana producers after Ecuador, the Philippines, Guatemala, and Costa Rica. In 2018, this fruit was the fourth main product in the global export trade (8.4% and USD 914,937,000) [10]. However, banana and plantain crops are severely affected by weevil species, generating losses of 2,600,000 tons/yield cycle [11,12]. Five banana weevil species have been reported so far. Cosmopolites sordidus (Germar) is a species with tropical and subtropical distribution around the world, with its center of origin in Southeast Asia [13,14], from where it spread widely over banana-production areas [13]. Metamasius hemipterus (Linnaeus), Metamasius hebetatus (Gyllenhal), and Metamasius submaculatus (Champion) are all distributed in Central and South America, with the Antilles as the plausible center of origin for the genus [15]. Meanwhile, Polytus mellerborgii (Boheman), an Indo-Malaysian native with pantropical distribution, has recently been spread by international trading in the tropical world [16], reaching the Neotropics within the last decade [12]. C. sordidus was reported as the most damaging species of banana and plantain plantations in Colombia [11,12] and worldwide [17][18][19]. The larval stage bores in the corm, creating galleries that reduce nutrient uptake and weaken the stability of the plant [17].
Conventional methods for pest control include chemical agents, which have been used worldwide for decades, bringing negative consequences to the environment, reducing the microbiome diversity in soils, jeopardizing ground water, increasing the flux of greenhouse gases, and accounting for 9.3 billion tons of CO 2 [4]. Besides the high toxicity of these chemicals, the resistance of weevil species (i.e., M. hemipterus) to a broad spectrum of synthetic insecticides has been reported [13], where endosymbiotic bacterial associations could be playing an important role in detoxification mechanisms [20][21][22], meaning the chemical method for pest control will become obsolete over time. Alternative strategies such as biological control are currently capturing the attention as methods to decrease soil pollution, while favoring costs for banana and plantain production. Several studies have applied this method by using insect enemies [23], fungal species such as Bauveria bassiana [24], and entomopathogenic nematodes to control C. sordidus populations [25]. However, these approaches still remain incipient.
Entomopathogenic nematodes (EPNs) are soil-dwelling, obligate insect-parasitic organisms with a cosmopolitan distribution (except the Arctic and Antarctic). Up to now, 96 Steninernema, one Neosteinernema, and 21 Heterorhabditis species have been described [2]. Their efficiency as biological control agents has targeted agricultural pests [26] in many different crops, with outstanding efficiency against Coleoptera [27], Lepidoptera [28], and Diptera [29]. The obligate parasitic nature of EPNs makes them optimal for integrated pest management (IPM) programs, as they are safe for the environment and other organisms compared to insect pests [1,30]. They are symbiotically associated with enterobacteria of the species Xenorhabdus for Steinernema spp. and Photorhabdus for Heterorhabditis spp., which have received increasing interest due to the antibiotic and antifungal properties of their secondary metabolites [2,[31][32][33][34]. EPNs are classified into two major families: Heterorhabditidae and Steinernematidae. All of their species are considered excellent biological control agents against the insect pests affecting many crops worldwide [30]. Based on their high entomopathogenic efficiency [26,28,[35][36][37], EPNs are commercially available as an alternative method to reduce soil and water pollution in agricultural systems [30].
The infection efficiency of EPNs relies on the enhanced pathogenicity of their bacterial symbiont [31,33]. Initially, the nematode enters into the insect's body via natural openings, with a posterior release of the symbiotic bacteria into the host hemocoel, causing septic shock and killing the insect in 24-48 h [26,38-40]. This enables the optimal conditions for EPNs' reproduction and multiplication before the insect's death. The new generation of EPNs feed on the insect's cadaver and acquire the bacterial symbiont [26,30]. Once the resources are depleted, the infective juveniles (IJ) exit the cadaver in search of a new host to infect [26,[41][42][43].
Some attempts of pest control were carried out by incorporating non-native EPNs in major crops worldwide [44,45]. However, these strains failed to parasitize local entomofauna, probably due to environmental conditions or the presence of pollutants in soils [44]. Therefore, it is paramount to isolate and characterize native EPNs, which may be locally adapted to the regional entomofauna. Local EPNs may be more infectious and competent in parasitism at either the adult or larval stage of development [46]. Particularly interesting are the strains that could tolerate intervened soils. In this sense, those EPN strains could provide potential applicability for in situ pest control.
The identification and evaluation of local EPN strains have partly been studied in Colombia with potential application in plantain and banana crop pest control [47,48], directed at the main banana beetle species found in the country [49]. These previous studies contemplated potential negative impacts on the crop [50] as well as the monitoring of reproductive cycles in order to develop target strategies for biological pest control [51]. Nevertheless, EPN molecular identification in the region remains incipient. Given this research gap, the present study aimed to: • Isolate and characterize native entomopathogenic nematodes that are locally adapted to plantain and banana crops, as options for weevil pest biological control; • Identify the symbiotic bacteria associated with EPNs; • Explore EPNs and their bacterial symbionts' phylogenetic diversity.

Sampling Area, EPN Isolation, and Propagation
Sampling area comprised banana and plantain plantations from 67 farms distributed in 12 municipalities located in the north, center, and south of Valle del Cauca province, in southwestern Colombia, which is a major hub of banana and plantain production in the Americas. Altitudes ranged from 7 to 1734 m asl. A total of 603 soil samples were obtained (nine from each farm), covering 325.7 hectares. Using a soil bore corer, the soil samples were taken at a depth of 20 cm and at a distance of 50 cm from the plant corm (rhizome).
The strategies followed for sampling included nine randomly chosen banana or plantain plants per farm. Three soil samples were taken around each plant in a triangle-shaped fashion, mixed into a single bag, and transported to the laboratory. To evaluate the presence of entomopathogenic nematodes (EPNs), 1 kg per soil sample was separated into three subsamples using plastic containers. Ten last instar Galleria mellonella larvae (Lepidoptera: Pyralidae) were employed as hosts and were placed alive into each container, following the method proposed by [52] with minor modifications. The containers were flipped over every 24 h and kept at room temperature (22 • C) for 3-7 days. All larvae cadavers were separated based on the color patterns drawn by the nematodes' infection (brownish for Steinernematids and reddish for Heterorhabditids). Afterward, they were rinsed with distilled water and placed in modified white traps [53,54].
The nematodes that were recovered were stored at 8 • C in a concentration of 2000 EPNs/mL and were reactivated by using a concentration of 200 EPNs/mL for each host, including G. mellonella, Metamasius hemipterus, Cosmopolites sordidus, and Tenebrio molitor larvae, for further analyses. The isolation and laboratory methods were carried out following the conditions described by Stock and Goodrich-Blair [55]. All these processes were conducted at Universidad Nacional de Colombia, Palmira campus.

DNA Extraction and PCR Amplification from EPN
The three obtained Steinernema isolates (coded as UNPR52, UNPR70, and UNPR73) were reared in vivo using last instar G. mellonella larvae. Afterward, a single F1 female was picked for DNA extraction using Quiagen ® DNeasy blood and tissue kit, in accordance with the instructions of the manufacturer. A final volume of 25 µL of DNA solution per isolate was obtained with an average concentration of 20 ng/µL, from which all the downstream molecular analyses were carried out. Four standard invertebrates' barcoding loci (COI and 12S mtDNA and ITS and LSU rDNA) were considered for phylogenetic and DNA pattern analyses; two of them were from the nuclear genome (ribosomal regions ITS and LSU), and two were from mitochondrial origin (COI and 12S). PCR amplifications were performed for each locus, mixing 2 µL of DNA 20 ng/µL, 2.5 µL of PCR Buffer 10× (NH4) 2 SO 4 , 2 µL of MgCl 2 25 mM, 1 µL of dNTP's 20 mM, 0.4 µL of each forward and reverse primer 10 mM, 0.5 µL of BSA 5×, 1 µL of Trehalose 10%, 0.16 µL of Taq DNA polymerase 5 U/µL, and 15 µL of ultra-pure water for a final volume of 25 µL. All sets of primers used are listed in Table 1. The PCR profile was as follows: 95 • C × 5 min, 95 • C × 1 min, 42 • C × 45 s, and 72 • C × 1 min, followed by 32 cycles at 95 • C × 1 min, 40 • C × 45 s, and 72 • C × 1 min, plus a final extension at 72 • C × 10 min. These amplification profiles varied in the annealing temperature: 42 • C for COI, 52 • C for 12S, 55 • C for ITS, and 54 • C for LSU. The PCR products were observed in 1.5% agarose gel and electrophoresed at 120 volts for 45 min in TBE solution 0.5×. All PCR products were purified following the protocol proposed by [56], and 1 µL was electrophoresed to verify the concentration.

DNA Sequencing and Data Analysis
DNA sequencing was performed in both directions (forward and reverse) for all loci. Yet, for ITS and LSU rDNA, internal primers were used to obtain high-quality sequence data. For ITS rDNA, the first primer combination AB28-534 flanked the 5 end of the sequence (~500 bp), whereas the second primer combination 533-TW81 flanked the 3 end (~600 bp). For LSU rDNA, the primer combinations 391-503 (~600 bp) and 502-501 (~670 bp) flanked the 5 and 3 ends of the sequences, respectively (Table 1). Raw sequences were initially assembled following each primer combination separately and were manually edited. PCR primers were anchored and removed from the sequences prior to the whole assembly of 5 and 3 end fragments. For ITS rDNA, a final sequence of~780 bp was obtained after contigs were assembled with an overlapping region of~300 bp, whereas for LSU the sequence comprised~900 bp with an overlapping region of~400 bp. All steps were carried out in Geneious R10 ® [60]. Multiple alignment was performed using Bioedit ® v.7.1.11 [61], including sequences from Steinernema species reported by Hunt and Nguyen [62] for ITS rDNA, LSU rDNA, 12S mtDNA, and COI mtDNA.

EPN Phylogenetic Analysis
Phylogenetic inferences were executed for each locus independently, following the best model of nucleotide substitution as obtained in Mega v.7.0 [63]. Maximum Parsimony (MP) and Maximum Likelihood (ML) analyses were performed using a bootstrap of 1000 replicates, based on General Time Reversible (GTR + G + I) model [64], including invariable sites (I) [65] and rate of variations across sites (G) [66]. A concatenated analysis was also performed, comprising sequences from 48 species reported for both ITS and LSU [62]. These MP and ML analyses comprised 1000 bootstrap replicates based on (GTR + G + I) model. Additionally, Bayesian analyses were carried out in MrBayes v.3.2 [67] under the GTR + G + I model with 10 million generations, sampling every 100 trees.
Caenorhabditis elegans was used as outgroup to root all phylogenies. Codon assignation, reading frames, and amino acid identification for COI mtDNA were performed in Geneious ® R10 [60], considering the genetic code for mitochondrial invertebrates. Each sample was individually translated to verify the correct reading frame. No stop codons were present. S. carpocapsae (AY944007) [68] was used as reference sequence to identify all substitutions at nucleotide and amino acid levels. The quantification of genetic diversity between EPN species relative to the isolates from this study was carried out in Mega 7 [63]. Meanwhile, the number of segregating sites based on the average number of polymorphic nucleotides between sequences and G + C content was computed in DnaSP 5.1 [69]. All nucleotide sequences for EPN isolates were deposited in the GenBank NCBI database ( Table 2). Table 2. NCBI accession numbers for each locus and georeferenced strain voucher isolated in the present study. Species names' assignation was based on molecular approaches. Some loci could not be sequenced for S. costaricense and are indicated by dashes (-).

Bacterial Symbionts Isolation and DNA Extraction
Bacterial recovery was executed four days after nematodes infection (using 200 IJs/ isolates/G. mellonella larvae), using a drop of hemolymph from infected larvae, which was then isolated and plated in NBTA agar incubated at 27 • C, following the method described by [55].
These isolates were coded as UNPB52, UNPB70, and UNPB73. Additional purification process was performed by choosing an individual colony per plate and placing it into a new culture medium, following the conditions described by Stock and Blair [55]. Afterward, a single colony per isolate was picked up and transferred into 15 mL of liquid culture medium (Luria Bertani LB), using oscillating shaking conditions and at 27 • C for 24 h. Bacterial cells were recovered by centrifugation. For the DNA extraction, DNeasy Blood & Tissue kit QUIAGEN ® was used, as in the manufacturer's guide.

Molecular Identification and Phylogenetic Analysis of EPNs' Bacterial Symbionts
Ribosomal 16S rDNA region was employed for molecular characterization using a set of universal primers (Table 1). External primers 27F and 1492R were employed for initial PCR reactions by mixing 2 µL of DNA 20 ng/µL, 2.5 µL of PCR Buffer 10× (NH 4 ) 2 SO 4 , 2 µL of MgCl 2 25 mM, 1 µL of dNTP's 20 mM, 0.4 µL of each forward and reverse primer 10 mM, 0.5 µL of BSA 5×, 1 µL of Trehalose 10%, 0.16 µL of Taq DNA polymerase 5 U/µL, and 15 µL of ultra-pure water, for a final volume of 25 µL. The thermal profile included an initial denaturing step at 94 • C × 4 min, followed by 30 cycles at 94 • C × 1 min, 55 • C × 45 s, and 72 • C × 1.30 min plus a final extension of 72 • C × 7 min. This locus generated a fragment of 1400 bp. However, internal primers were employed for better sequencing results. The primer combination of 27F and 800R flanked the 5 end of the sequence, generating a partial fragment of 760 bp, while the primer combination of 518F and 1492R flanked the 3 end and generated a sequence of 940 bp. Sequence edition and contig assembly was manually and independently carried out in Geneious R10 ® [60] for each isolate. PCR primers (external and internal) were removed prior to fragment assembly. For nucleotide pattern identification and phylogenetic analysis, the conditions described above in Section 2.4. were applied for 16S rDNA region, based on the GTR + G + I model [64] for phylogenetic inferences. All nucleotide sequences from bacterial isolates were deposited in the GenBank NCBI database (Table 3).

Sampling Variation
From the extensive sampling strategy carried out in banana and plantain plantations, covering 12 municipalities of Valle del Cauca province (Figure 1), three soil samples out of 603 (0.49%) were positive for entomopathogenic nematodes (EPNs). These isolates were recovered from the municipalities of Obando (UNPR52), Palmira (UNPR70), and Buenaventura (UNPR73). All samples and isolates were properly coded and stored, following the conditions previously defined at the molecular biology laboratory for the project BPIN 2014000100014. Based on the characteristics of G. mellonella larvae cadavers and the morphology of infective juveniles (IJ) and adults, all isolates were identified as Steinernema spp. Isolates UNPR52 and UNPR73 showed a good performance of pathogenicity on either G. mellonella, M. hemipterus, C. sordidus, or Tenebrio molitor larvae (≥90% mortality), while UNPR70 failed to parasite after a couple of infection cycles and could not be reared under laboratory conditions for any detailed morphological evaluations. All evaluations of infection, multiplication of EPNs, and downstream molecular analyses were performed in the laboratory of molecular biology at the Universidad Nacional de Colombia, Palmira campus.

EPN Nucleotide Diversity
Four loci were sequenced, two ribosomal (LSU and ITS) and two mitochondrial (12S and COI) regions. For isolates UNPR52 and UNPR73, the sequencing covered 2,918 bp, whereas, for the isolate UNPR70, only the LSU region was sequenced, generating 882 bp. Comparisons among the three isolates based on the LSU region showed 115 polymorphic sites. Two nucleotides varied between UNPR52 (accession number MH231232) and UNPR73 (accession number MH231233), whilst most of the nucleotide variation was exclusive to UNPR70 (accession number MH231234). In order to identify the genetic similarity of these isolates to the species reported so far for the Steinernema genus, the analysis included 60 LSU sequences from Hunt and Nguyen [62], generating a total of 1130 sites, of which 190 were polymorphic and 149 were parsimony-informative. The genetic distance indicated no variation between UNPR52 and UNPR73, with respect to the S. carpocapsae isolates reported by Mráček et al. [70] (KJ950293 and KJ950292). On the other hand, UNPR70 showed genetic proximity (0.034) to S. costaricense isolate "Bush-Augusta" (NCBI accession number JX493015), as reported by Ivanova et al. [71], with 14 nucleotide changes. The ITS locus analysis included isolates UNPR52 (accession number MH231235) and UNPR73 (accession number MH231236), without any nucleotide substitution in between. The whole set of sequences reported by Hunt and Nguyen [62] for the ITS region included 85 Steinernema spp. and generated 1293 sites, of which 127 were polymorphic and 78 were parsimony-informative. The genetic distances showed high similarity (0.001) to S. carpocapsae "Breton" isolate (under NCBI accession number AF121049, [72]) and to isolate "NCR" (with accession number KJ950291, [70]), showing a genetic distance of 0.006 (five nucleotide changes). The mitochondrial loci (12S and COI) exhibited contrasting results, with respect to the information obtained for the ribosomal loci. For the 12S mtDNA, 507 site nucleotide positions were considered, of which 210 were polymorphic and 147 were parsimony-informative. The 27 Steinernema spp. sequences reported for this locus by Nadler et al. [68] were included.

EPN Nucleotide Diversity
Four loci were sequenced, two ribosomal (LSU and ITS) and two mitochondrial (12S and COI) regions. For isolates UNPR52 and UNPR73, the sequencing covered 2918 bp, whereas, for the isolate UNPR70, only the LSU region was sequenced, generating 882 bp. Comparisons among the three isolates based on the LSU region showed 115 polymorphic sites. Two nucleotides varied between UNPR52 (accession number MH231232) and UNPR73 (accession number MH231233), whilst most of the nucleotide variation was exclusive to UNPR70 (accession number MH231234). In order to identify the genetic similarity of these isolates to the species reported so far for the Steinernema genus, the analysis included 60 LSU sequences from Hunt and Nguyen [62], generating a total of 1130 sites, of which 190 were polymorphic and 149 were parsimony-informative. The genetic distance indicated no variation between UNPR52 and UNPR73, with respect to the S. carpocapsae isolates reported by Mráček et al. [70] (KJ950293 and KJ950292). On the other hand, UNPR70 showed genetic proximity (0.034) to S. costaricense isolate "Bush-Augusta" (NCBI accession number JX493015), as reported by Ivanova et al. [71], with 14 nucleotide changes. The ITS locus analysis included isolates UNPR52 (accession number MH231235) and UNPR73 (accession number MH231236), without any nucleotide substitution in between. The whole set of sequences reported by Hunt and Nguyen [62] for the ITS region included 85 Steinernema spp. and generated 1293 sites, of which 127 were polymorphic and 78 were parsimonyinformative. The genetic distances showed high similarity (0.001) to S. carpocapsae "Breton" isolate (under NCBI accession number AF121049, [72]) and to isolate "NCR" (with accession number KJ950291, [70]), showing a genetic distance of 0.006 (five nucleotide changes). The mitochondrial loci (12S and COI) exhibited contrasting results, with respect to the information obtained for the ribosomal loci. For the 12S mtDNA, 507 site nucleotide positions were considered, of which 210 were polymorphic and 147 were parsimony-informative. The 27 Steinernema spp. sequences reported for this locus by Nadler et al. [68] were included. UNPR52 (accession number MH231121) and UNPR73 (accession number MH231122) were genetically similar (0.004) to S. carpocapsae isolate "N2" (accession number AY944007, [68]), S. anatoliense isolate "Al-Jubiha" (accession number GU569024, [73]), and S. websteri isolate "Peru" (accession number GU569039, [73]). S. anatoliense and S. websteri were reported as synonyms of S. carpocapsae [62]. For the COI mtDNA analyses, the final dataset considered 31 species. A total of 676 nucleotide positions were included, of which 211 showed polymorphism and 191 were parsimony-informative. The nucleotide composition comprised Thymine (T; 45.6%), Cytosine (C; 14.3%), Adenine (A; 22.9%), and Guanine (G; 17.2%), with a strong A + T bias (68.5%). Nucleotide variation generated a genetic distance of 0.009 between the UNPR52 and UNPR73 isolates, with respect to S. carpocapsae isolate "N2" (accession number AY943981, [68]), 0.004 with S. anatoliense isolate "Al-Jubiha" (accession number GU569062, [73]), and no variation respect to S. websteri isolate "Peru" (GU569074, [73]). This analysis confirms the synonymy of these three species, all of them belonging to the Carpocapsae clade [62].
Following the multispecies clades proposed by Hunt and Nguyen [62], a boxplot for FST analysis indicated no statistical differences between groups for either of the evaluated loci (Figure 2a  A heat map analysis, based on FST values, provided an alternative perspective for clades interaction (Figure 3), where the Affine is distant for both loci and concatenated. Following the superclades classification proposed by Hunt and Nguyen [62] and Nguyen and Hunt [74], we identified that superclade II was consistent for LSU ( Figure 3a) and ITS (Figure 3b) but failed when concatenated (Figure 3c). These results indicate a tendency of similar allele frequencies with possible signatures of positive selection at the LSU locus, conversely to the ITS locus (contrasting selective pressures on both loci). Despite the two markers have enough resolution for species identification, ITS rDNA has better differentiation sensitivity, since the variance in allele frequency is higher, providing additional information about the genetic diversity status at the species level optimal for the identification of new strains and their possible correlation to particular habitats. A heat map analysis, based on F ST values, provided an alternative perspective for clades interaction (Figure 3), where the Affine is distant for both loci and concatenated. Following the superclades classification proposed by Hunt and Nguyen [62] and Nguyen and Hunt [74], we identified that superclade II was consistent for LSU ( Figure 3a) and ITS (Figure 3b) but failed when concatenated (Figure 3c). These results indicate a tendency of similar allele frequencies with possible signatures of positive selection at the LSU locus, conversely to the ITS locus (contrasting selective pressures on both loci). Despite the two markers have enough resolution for species identification, ITS rDNA has better differentiation sensitivity, since the variance in allele frequency is higher, providing additional information about the genetic diversity status at the species level optimal for the identification of new strains and their possible correlation to particular habitats. For amino acid composition analysis (at COI locus), 26 Steinernema spp. sequences were included, obtaining 225 amino acid residues with 15 changes among the species (Table S1). Leucine (Leu) and Serine (Ser) were the most frequent (18.1% and 10.3%, respectively), whereas Glutamine (Gln) and Lysine (Lys) represented the less-frequent amino acid residues (0.8% and 0.9%, respectively). The nucleotide substitutions that occurred among UNPR52, UNPR73, S. carpocapsae (AY944007), S. anatoliense (GU569024), and S. websteri (GU569074) included four transition events (A/G, T/C, C/T, and G/A) and two transversions (T/A and T/G), though all of them were synonymous mutations.
The concatenated analysis (ITS and LSU rDNA regions) included 51 nucleotide sequences, of which 50 belonged to Steinernema spp. in addition of C. elegans as the outgroup for phylogeny rooting. All the sequences were retrieved from the NCBI and reported by Hunt and Nguyen [62]. The final dataset comprised 2022 nucleotides, 361 polymorphic positions, 277 parsimony-informative sites, and a G + C content of 47.3%. There was not genetic distance between UNPR52 and UNPR73 compared to S. carpocapsae (KJ950291-KJ950293), supporting the classification into the latter. Concatenated analyses were unviable due to poor open data for mitochondrial loci (i.e., COI and 12S).

EPN Phylogenetic Analysis
Initially, all loci were independently analyzed to verify the congruence of phylogenetic inferences as well as the clustering of our isolates, with respect to the sequences previously reported. This approach corroborated the accuracy by each locus in a phylogenetic context. For the ribosomal locus LSU rDNA, 64 samples belonging to 58 Steinernema species, including the three isolates obtained here, were analyzed by the Maximum Parsimony (MP), Maximum Likelihood (ML), and Bayesian Inference (BI) methods, using C. elegans as the outgroup (Figure 4). The final dataset comprised 1154 positions, yielding a parsimonious tree of 1935 steps (C.I. 0.45 and R.I. 0.74). Isolates UNPR52 and UNPR73 were grouped into the Carpocapsae clade. Although the posterior probability for this clade was not well-supported (ML 54% and MP 59%), the BI was significant (0.93). Conversely, the isolate UNPR70 was clustered with S. costaricense in the Costaricense clade with a strong bootstrap and posterior supports (99% ML, 100% MP, and BI 1.0), except for S. scarabaei, which was conflicting for all phylogenetic inferences. Following the phylogenetic hypothesis of multispecies clades proposed by Hunt and Nguyen [62], we identified Khoisanae, Feltiae, Cameroonense, Bicornutum, Carpocapsae, and Affine as monophyletic groups for all the three analyses. With the exception of the Carpocapsae clade, all the remaining clades were well supported, with bootstrap values ≥70% for ML and MP, and posterior probability ≥0.93 for BI. For amino acid composition analysis (at COI locus), 26 Steinernema spp. sequences were included, obtaining 225 amino acid residues with 15 changes among the species (Table S1). Leucine (Leu) and Serine (Ser) were the most frequent (18.1% and 10.3%, respectively), whereas Glutamine (Gln) and Lysine (Lys) represented the less-frequent amino acid residues (0.8% and 0.9%, respectively). The nucleotide substitutions that occurred among UNPR52, UNPR73, S. carpocapsae (AY944007), S. anatoliense (GU569024), and S. websteri (GU569074) included four transition events (A/G, T/C, C/T, and G/A) and two transversions (T/A and T/G), though all of them were synonymous mutations.
The concatenated analysis (ITS and LSU rDNA regions) included 51 nucleotide sequences, of which 50 belonged to Steinernema spp. in addition of C. elegans as the outgroup for phylogeny rooting. All the sequences were retrieved from the NCBI and reported by Hunt and Nguyen [62]. The final dataset comprised 2022 nucleotides, 361 polymorphic positions, 277 parsimony-informative sites, and a G + C content of 47.3%. There was not genetic distance between UNPR52 and UNPR73 compared to S. carpocapsae (KJ950291-KJ950293), supporting the classification into the latter. Concatenated analyses were unviable due to poor open data for mitochondrial loci (i.e., COI and 12S).

EPN Phylogenetic Analysis
Initially, all loci were independently analyzed to verify the congruence of phylogenetic inferences as well as the clustering of our isolates, with respect to the sequences previously reported. This approach corroborated the accuracy by each locus in a phylogenetic context. For the ribosomal locus LSU rDNA, 64 samples belonging to 58 Steinernema species, including the three isolates obtained here, were analyzed by the Maximum Parsimony (MP), Maximum Likelihood (ML), and Bayesian Inference (BI) methods, using C. elegans as the outgroup (Figure 4). The final dataset comprised 1154 positions, yielding a parsimonious tree of 1935 steps (C.I. 0.45 and R.I. 0.74). Isolates UNPR52 and UNPR73 were grouped into the Carpocapsae clade. Although the posterior probability for this clade was not well-supported (ML 54% and MP 59%), the BI was significant (0.93). Conversely, the isolate UNPR70 was clustered with S. costaricense in the Costaricense clade with a strong bootstrap and posterior supports (99% ML, 100% MP, and BI 1.0), except for S. scarabaei, which was conflicting for all phylogenetic inferences. Following the phylogenetic hypothesis of multispecies clades proposed by Hunt and Nguyen [62], we identified Khoisanae, Feltiae, Cameroonense, Bicornutum, Carpocapsae, and Affine as monophyletic groups for all the three analyses. With the exception of the Carpocapsae clade, all the remaining clades were well supported, with bootstrap values ≥70% for ML and MP, and posterior probability ≥0.93 for BI.   For the ITS rDNA analyses, a total of 88 nucleotide sequences were considered, of which 85 belonged to the Steinernema species, 2 sequences were from the isolates found in this study (UNPR52 and UNPR73), and C. elegans was the outgroup. The final dataset contained 1691 positions, yielding a parsimonious tree of 7479 steps (C.I 0.28 and R.I 0.60). Highly congruent topologies were identified for all phylogenetic inferences (MP, ML, and BI; Figure 5), under bootstrapping resampling for most of the clades. The isolates UNPR52 and UNPR73 were clustered with S. carpocapsae, with high support for the branch tip (100% and 1.0) and at the Steinernema clade backbone (100% and 0.96). Unlike the LSU rDNA region, the ITS rDNA region revealed moderate (≥75%) to strong (100%) support for the nodes containing the multispecies clades, which were Longicaudatum, Monticolum, Kushidai, Feltiae, Cameroonense, Bicornutum, Carpocapsae, and Affine as the best-supported clades (≥75% for ML and MP and ≥0.68 for BI).
The tree topologies generated by these two loci demonstrated the consistency for the aforementioned clades, except for the Longicaudatum clade that was not supported by LSU rDNA. Additionally, the position of these clades, with respect to one another, was congruent with the topology proposed by Hunt and Nguyen [62], with superclade 2 (Feltiae, Kushidai, and Monticolum) and superclade 3 (Carpocapsae, Bicornutum, and Affine) being well-identified.

Concatenated Analysis for Ribosomal Loci in EPNs
A total of 51 sequences for the concatenated analysis comprised the LSU rDNA and ITS rDNA regions for the same isolates, where 48 belonged to Steinernema spp., two were the isolates UNPR52 and UNPR73, and one was the C. elegans outgroup. From the 2081 positions in the final dataset, a tree length of 5849 steps was obtained (C.I 0.41 and R.I 0.67). Eight clades were covered, five of which were highly congruent in tree topology [(Kushidai (two species), Feltiae (11), Cameroonense (4), Carpocapsae (5), and Affine (4)], with bootstrap values ≥97% and Bayesian posterior probabilities ≥0. 61.
The isolates that clustered into the Carpocapsae clade generated a sub-cluster at the tip of the branches, being strongly supported by the bootstrap values (ML 100%, MP 99%, and BI 1.0; Figure 6). The tree comparison showed conflicting topologies for multispecies clades inferred from both the concatenated and individual datasets, where Glaseri, Khoisanae, and Karii showed unstable topologies. Additionally, the Bicornutum clade was conflicting only in the concatenated dataset, where the species were separately clustered into two clades, with one next to the other (i.e., sister clades). On the contrary, monophyletic clades were identified for either individual locus and/or concatenated loci, from weakly (≥54% LSU Carpocapsae clade) to moderately (≥70% LSU rDNA Feltiae) to strongly (≥99% ITS rDNA and concatenated-Monticolum, Kushidai, Feltiae, Carpocapsae, Cameroonense, and Affine) supported topologies. In general, ITS rDNA showed better resolution for tree topology, bootstrapping scores (≥75%), and posterior probability values (≥0.68) for 10 (Khoisanae, Longicaudatum, Costaricence, Monticolum, Kushidai, Feltiae, Cameroonense, Bicornutum, Carpocapsae, and Affine) of the 12 multispecies clades, as previously reported [62,74].    [62]. Triangles indicate the strains obtained in this study. From left to right on top of the main branches of each cluster: ML, MP, and Bayesian support scores. Caenorhabditis elegans was used as the outgroup.

Concatenated Analysis for Mitochondrial Loci in EPNs
Datasets for mitochondrial concatenated analysis included the 12S mtDNA and COI mtDNA. Twenty sequences included the final dataset, 17 from Steinernema spp. retrieved from Hunt and Nguyen [62], two isolates UNPR53 and UNPR73, and the outgroup. This analysis comprised 1296 positions generating a tree length of 1585 steps (C.I 0.46 and R.I 0.50). MP, ML, and BI supported the relationships (100%, 100%, and 1.0, respectively) between our isolates and S. carpocapsae (Figure 7).  [62]. Triangles indicate the strains obtained in this study. From left to right on top of the main branches of each cluster: ML, MP, and Bayesian support scores. Caenorhabditis elegans was used as the outgroup.

Concatenated Analysis for Mitochondrial Loci in EPNs
Datasets for mitochondrial concatenated analysis included the 12S mtDNA and COI mtDNA. Twenty sequences included the final dataset, 17 from Steinernema spp. retrieved from Hunt and Nguyen [62], two isolates UNPR53 and UNPR73, and the outgroup. This analysis comprised 1296 positions generating a tree length of 1585 steps (C.I 0.46 and R.I 0.50). MP, ML, and BI supported the relationships (100%, 100%, and 1.0, respectively) between our isolates and S. carpocapsae (Figure 7). In general, all trees showed the same topology from weakly (≥51%) to strongly (≥93%) supported bootstrap values at the most recent nodes. For the Carpocapsae clade, four species (S. carpocapsae, S. anatoliense, S. scapterisci, and S. siamkayai) were considered to generate a monophyletic group with high bootstrap values and posterior probabilities support (100% and 1.0). Other monophyletic groups identified were Bicornutum and Feltiae clades, with three species each, from moderately (≥0.85) to strongly (≥92%) supported nodes.

Nucleotide Diversity in EPN's Symbiotic Bacteria
Molecular characterization for the bacterial symbionts considered the whole 16S rDNA region. The isolates were coded the same way as their respective EPN sources: UNPB52, UNPB70, and UNPB73. Initial evaluations comprised 32 sequences belonging to 23 species plus the three isolated bacteria from the present study. A total of 1550 sites were included, of which 118 were polymorphic and 89 were parsimony-informative. Nucleotide composition corresponded to A 24.8%, C 23.1%, G 32,5%, and T 19.6%, with G + C content of 55.6%. Genetic distance reinforced the closeness between UNPB52 (accession number MH249989) and UNPB73 (accession number MH249990) with Xenorhabdus nematophila strain "DSM 3370" (accession number AY278674 [75] at 0.002, and with UNPB73 at 0.003). Additionally, only three nucleotide substitutions were identified between UNPB52 and UNPB73.
For the isolate UNPB70 (accession number MH249991), which belongs to the symbiotic bacteria isolated from S. costaricense, the genetic similarity was toward X. koppenhoeferi strain "USNJ01" (accession number DQ205450, [76]) and X. khoisanae strain "SF362" (accession number JX623978 [77]), having 21 nucleotide substitutions for each comparison with a distance of 0.012. Overall comparisons between Xenorhabdus species showed 147 mutations in specific regions of the 16S rDNA, distributed mainly in three mutational hotspots at the 5′ end of the sequence (from nucleotides 70 to 100 and from nucleotides 460 to 480) and at the 3′ end of the sequence (from nucleotides 1000 to 1050). In general, all trees showed the same topology from weakly (≥51%) to strongly (≥93%) supported bootstrap values at the most recent nodes. For the Carpocapsae clade, four species (S. carpocapsae, S. anatoliense, S. scapterisci, and S. siamkayai) were considered to generate a monophyletic group with high bootstrap values and posterior probabilities support (100% and 1.0). Other monophyletic groups identified were Bicornutum and Feltiae clades, with three species each, from moderately (≥0.85) to strongly (≥92%) supported nodes.

Nucleotide Diversity in EPN's Symbiotic Bacteria
Molecular characterization for the bacterial symbionts considered the whole 16S rDNA region. The isolates were coded the same way as their respective EPN sources: UNPB52, UNPB70, and UNPB73. Initial evaluations comprised 32 sequences belonging to 23 species plus the three isolated bacteria from the present study. A total of 1550 sites were included, of which 118 were polymorphic and 89 were parsimony-informative. Nucleotide composition corresponded to A 24.8%, C 23.1%, G 32.5%, and T 19.6%, with G + C content of 55.6%. Genetic distance reinforced the closeness between UNPB52 (accession number MH249989) and UNPB73 (accession number MH249990) with Xenorhabdus nematophila strain "DSM 3370" (accession number AY278674 [75] at 0.002, and with UNPB73 at 0.003). Additionally, only three nucleotide substitutions were identified between UNPB52 and UNPB73.
For the isolate UNPB70 (accession number MH249991), which belongs to the symbiotic bacteria isolated from S. costaricense, the genetic similarity was toward X. koppenhoeferi strain "USNJ01" (accession number DQ205450, [76]) and X. khoisanae strain "SF362" (accession number JX623978 [77]), having 21 nucleotide substitutions for each comparison with a distance of 0.012. Overall comparisons between Xenorhabdus species showed 147 mutations in specific regions of the 16S rDNA, distributed mainly in three mutational hotspots at the 5 end of the sequence (from nucleotides 70 to 100 and from nucleotides 460 to 480) and at the 3 end of the sequence (from nucleotides 1000 to 1050).

Phylogenetic Analysis of EPN's Symbiont Bacteria
The phylogenetic analyses carried out for EPN's symbiotic bacteria 16S rDNA included 33 nucleotide sequences from the 16S rDNA locus covering 21 Xenorhabdus species, three isolates from the present study, and Photorhabdus luminescens strain DSM 3368 (X82248) [78] as the outgroup. From 1550 sites in the final dataset, 89 were parsimony-informative. A parsimonious tree was yielded with 609 steps (C.I 0.34 and R.I 0.62). Isolates UNPB52 (MH249989) and UNPB73 (MH249990), corresponding to the UNPR52 and UNPR73 EPN isolates, respectively, were grouped in the same cluster with Xenorhabdus nematophila under strong bootstrap support (96%) and Bayesian posterior probability (1.0), whereas the bacterial isolate UNPB70 (MH249991, corresponding to the EPN isolate UNPR70) was clustered with X. koppenhoeferi with very weak bootstrap support (40%) and Bayesian posterior probability (0.63). Tree comparisons showed instability for the most ancient nodes in each phylogenetic inference. Meanwhile, recent nodes showed monophyletic clades (X. nematophila, X. beddingii, X. poinarii, X. bovienii, and X. japonica) strongly supported for both bootstrap resampling (≥96%) and Bayesian inference (≥0.96; Figure 8). The clustering of X. nematophila with the sequences from this study (UNPB52 and UNPB73) in all tree topologies suggests S. carpocapsae as the host nematode, thus clarifying the molecular taxonomic of these isolates.

Phylogenetic Analysis of EPN's Symbiont Bacteria
The phylogenetic analyses carried out for EPN's symbiotic bacteria 16S rDNA included 33 nucleotide sequences from the 16S rDNA locus covering 21 Xenorhabdus species, three isolates from the present study, and Photorhabdus luminescens strain DSM 3368 (X82248) [78] as the outgroup. From 1550 sites in the final dataset, 89 were parsimonyinformative. A parsimonious tree was yielded with 609 steps (C.I 0.34 and R.I 0.62). Isolates UNPB52 (MH249989) and UNPB73 (MH249990), corresponding to the UNPR52 and UNPR73 EPN isolates, respectively, were grouped in the same cluster with Xenorhabdus nematophila under strong bootstrap support (96%) and Bayesian posterior probability (1.0), whereas the bacterial isolate UNPB70 (MH249991, corresponding to the EPN isolate UNPR70) was clustered with X. koppenhoeferi with very weak bootstrap support (40%) and Bayesian posterior probability (0.63). Tree comparisons showed instability for the most ancient nodes in each phylogenetic inference. Meanwhile, recent nodes showed monophyletic clades (X. nematophila, X. beddingii, X. poinarii, X. bovienii, and X. japonica) strongly supported for both bootstrap resampling (≥96%) and Bayesian inference (≥0.96; Figure 8). The clustering of X. nematophila with the sequences from this study (UNPB52 and UNPB73) in all tree topologies suggests S. carpocapsae as the host nematode, thus clarifying the molecular taxonomic of these isolates.

Discussion
The identification of native EPNs that are associated with the plantain and banana rhizosphere is of significant importance for biological control purposes and agricultural land use due to their adaptation to environmental conditions, soil pollution tolerance, and interaction with the native insect community, providing new insights into EPNs' efficiency. In the present study, three native EPNs (two Steinernema carpocapsae and one S. costaricense) could be recovered from soils that were destined for plantain and banana crop plantations. The low rate of EPN recovery (0.49%) could be influenced by several factors such as the high use of agrochemical pollutants in agricultural soils [27,[79][80][81][82], soil sampling strategy, and failure of EPNs to parasite insect hosts [83,84]. Laboratory assays showed differential pathogenicity responses of S. carpocapsae and S. costaricense isolates in parasitizing G. mellonella and weevil species (C. sordidus and M. hemipterus larval stages), with S. carpocapsae being the most pathogenic for all the insect larvae species evaluated, whereas S. costaricense failed to parasite any of the larvae host species.
Several factors have been proposed to explain these differences in pathogenicity, including (1) the EPN body size, with S. carpocapsae being the smallest, which provides efficiency for host penetration [85]; (2) the pathogenicity of symbiotic bacteria types "swarming" and "stimulated" [86], where the former type has a more lethal activity and is present in the dauer juveniles, whilst the second type is a latent immobile stage that has a muchreduced lethal activity in response to environmental adversities [83]; and (3) the immune response of the host, such as the phagocytosis and humoral mechanisms against the symbiotic bacteria [42,83]. It was reported that S. carpocapsae possesses a higher rate of mortality on several species of Hemiptera [83], Coleoptera [27,39,[87][88][89], Diptera [29,90], and Lepidoptera [28], which correlates with our experimental results where S. carpocapsae killed G. mellonella, C. sordidus, and M. hemipterus larvae in 48 h. However, S. costaricense could not parasite after one round of infection. Shapiro-Ilan and Raymond [84] proposed that the main drawback of EPN reduction virulence is due to the need of serial propagation in vivo, which can accentuate competition and reduce the reproductive rate, leading to virulence instability. On the other hand, Adams et al. [91] argued that EPN virulence could be related to strategic dispersal, based on the host's dependence on mobility to reach other potential hosts, which is a strong selection pressure to reduce EPN virulence, while keeping the host alive long enough to provide resources for the life cycle.
In Colombia, the first report on EPN species, from forest and crop soils, was relatively recent [92]. Additionally, López-Núñez et al. [57] reported a new Steinernema species S. colombiense, which, nevertheless, was excluded from the analysis due to inconsistencies in the phylogenetic assignment, demonstrating the potential diversity of these organisms in the country. Several EPN strains were imported to Colombia for biological control strategies that unfortunately failed. Therefore, identifying native EPNs well adapted to specific environmental conditions may convey applications to a broad spectrum of crops in the region. For instance, S. costaricense was initially reported by Uribe-Lorío et al. [93] in forest soils from Costa Rica based on morphological and molecular characterizations (LSU rDNA). Later, Ivanova et al. [71] reported that the same species were isolated from deciduous woodland in Missouri, USA, based on morphology and molecular data (ITS rDNA).
Here, we report for the first time S. costaricense in Colombia, based on the LSU rDNA sequence isolated from agricultural soils, as well as its symbiotic bacteria based on molecular data (16S rDNA gene), showing close genetic similarity to Xenorhabdus koppenhoeferi [76] and X. khoisanae [77], while still being uncertain in its classification. Yet the identification of S. costaricense in South America indicates not only the wide range of this species in the Americas but also the ability to persist in different soil types, an ad hoc hypothesis that deserves further exploration. Also for the first time, EPN isolation is directed at plantain and banana crops as an alternative to identify the natural persistence of these entomopathogens and, eventually, to couple land agricultural practices with native organisms as potential biological control agents directed at insect pests.
Recently, similar approaches were carried out in the isolation of indigenous EPNs for pest control, targeting local adaption for biological control [27,89,94,95]. These agents provide promising results to reduce intense chemical pest management [89,96]. The isolation of environmentally adapted EPNs to cropland and, specifically, from plantain and banana cultures, is a feasible alternative in weevil-species-complex pest control, while mitigating soil, water, and environmental pollution. The strains obtained here and their symbionts showed a remarkable virulence against the two most common weevil pest species (Cosmopolites sordidus and Metamasius hemipterus), killing over 90% of the hosts (to be included in an oncoming report) and demonstrating their utility as biological control agents.

Polymorphism of Gene Barcodes in EPN and Their Bacterial Endosymbionts
Several genetic regions from eukaryotes were analyzed as potential sources of information for species identification [70,72,74,[97][98][99][100][101][102]. These target barcode regions comprise nuclear ribosomal loci (ITS, SSU, and LSU) and mitochondrial genes (COI, 12S, and NAD). Among these, ribosomal loci represent a major tool for investigating and identifying nematodes species [72,74,89,95,99,101,103]. They also represent the main source of information in the databases, providing more reliability for species identification [74,99,101]. Among the ribosomal loci, the coding regions SSU and LSU play an important role in entomopathogenic nematodes' species validation and phylogenetic relationships, due to its conserved nature [74,99,101]. On the contrary, ITS is a rapidly evolving, presumably neutral, noncoding region that provides information about recent evolutionary events, which is useful for intra-species genetic diversity [72,74,99,104].
The patterns of genetic variation between EPN taxa demonstrate a high rate of nucleotide substitution in the ITS region compared to the LSU region across all the Steinernema species analyzed here. Nevertheless, very low genetic variation was observed at the intra-species level for all nuclear loci (ITS and LSU), even for the isolates retrieved by this study, despite the highly intervened soil conditions from which they were recovered. Entomopathogenic nematodes present a higher rate of evolution relative to other metazoans [72,74,105]. However, several factors could constrain certain genomic regions from evolving freely [99]. For rRNA loci, some gene variants have been detected as dominant in nematodes populations, which are often interpreted under purifying selection and concerted evolution [99,106]. Favored variants could in turn act against alternative versions as a mechanism to prevent rare rRNA alleles that could affect their function (e.g., structural modifications leading to folding errors in the secondary structure), maintaining in this way the fitness of the population [99]. This process may be reinforced by the pathogenicity nature of EPNs because virulence against several host species could favor slow evolution. After all, EPNs with a broad spectrum of hosts exhibit higher fitness to spread over populations and distant geographic regions. This scenario would, ultimately, maintain the gene flow over time, keeping the dominance of certain few allelic variants [74]. Given this situation, it is expected that most of the genetic variation may go undetected, limiting the quantification of biodiversity in EPNs. Most of the reported sequences for these taxa were obtained using Sanger sequencing technology, which summarizes the information from the large gene pools of PCR products enriched in abundant genomic rRNA copies, while being limited to specific barcodes [99]. Meanwhile, with the advance of sequencing technologies (high-throughput sequencing), it has been possible to detect small copy-number variants [107] in multiple eukaryotic genomes [106][107][108][109][110][111] and, particularly, in nematodes (i.e., Caenorhabditis elegans, [99]). Therefore, these technologies have offered a more precise approximation of EPN diversity, particularly for native species adapted to pollute environments. More recently, genomic studies on EPN [112][113][114], particularly in S. carpocapsae, found a correlation between non-coding RNA (ncRNA) and its functional role in the pathogenic lifestyle, which is completely absent in free-living nematodes [114]. These findings provide new insights into the molecular mechanisms underlying pathogenicity, as well as into the selective pressures targeting specific genomic regions such as rDNA loci, which could be under control of the gene expression by ncRNA or epigenetic factors.
The search for more precise species-identification genetic barcodes is mandatory because morphological traits are not completely satisfactory for EPN classification (e.g., S. carpocapsae [70]). This misclassification leads to a poor interpretation of the ecology and evolution in these complex taxa. Furthermore, the high copy number of the rDNA loci identified in nematodes [99] could distort the evolutionary hypothesis in phylogenetic studies, leading to incorrect species assignation [70]. This makes it difficult to report new strains with highly divergent rDNA sequence data, to the point that rare individuals may display a type of pre-adaptation to new environments [115] or host-associated habitats [116]. This case was identified in the present study with the Steinernema colombiense isolate SNI-0198 (NCBI accession number EU345420), originally reported by López-Núñez et al. [57] from a coffee plantation in Colombia. However, its correct classification remains uncertain, so it had to be removed from our analysis due to discrepancies in the phylogenies.
Analysis of the 16S rDNA sequences obtained from bacterial symbionts confirmed the species Xenorhabdus nematophila was associated with S. carpocapsae, whilst, for S. costaricense, it was not possible to determine its symbiont at the species level, being the closest taxa to X. koppenhoeferi (mutualistic to S. scarabei) and X. khoisanae (mutualistic to S. khoisanae, S. jeffreyense, and S. saccharii). Steinernema costaricense was first identified by Uribe-Lorío et al. [93], following a detailed morphological description and the sequencing of the LSU rDNA locus. However, its bacterial symbiont was not characterized. Later on, Stock and Blair [35] related X. szentirmaii as the mutualistic bacteria associated with S. costaricense; nevertheless, this association was not supported by other studies [34,74,117]. The isolation and identification of the EPN's bacterial symbionts has received special attention, not only for their entomopathogenic efficiency [27,29,31,33,39,89,[94][95][96][118][119][120][121][122] but also for their potential applicability as antimicrobials [123], even targeting drug-resistant bacteria [34], their richness in metabolites with agricultural interest [32,124], and their abundance of peptides with functionality for cancer treatment [125].
Species identification based on the 16S rDNA locus has been challenging due to intragenomic variable copy numbers, marker heterogeneity, and sequence variation within closely related taxa. Despite certain taxa harbor common 16S rDNA sequences from multiple species [58], genomes with more 16S rDNA copies are expected to carry diverse variants from unrelated origins (i.e., accounting for 2.4% [58] to 10% [126] of the copy number variants). In the present study, we did not analyze the presence of multiple 16S rDNA copy numbers in the EPN bacterial symbionts, since the similarity of this ribosomal locus between different species suggests that this gene may be coupled with multiple symbiotic interactions [74]. Espejo and Plaza [127] proposed the underlying effect of concerted evolution running at the intra-and inter-genomic levels (through horizontal gene transfer) in proteobacteria species (including endosymbionts), leading to rRNA locus homogenization, which in turn could confer advantages for environmental adaptations. Moreover, intergenomic recombination is favored and maintained within clonal populations, overcoming the mutation rate [73]. Multiple EPN species co-infecting the same host has been observed [128,129], providing close proximity in such way that possible horizontal gene transfer could be accomplished with posterior gene conversion, making the 16S rDNA locus still uncertain for species identification and diversity screening [126,127,130,131].
Similarly, microbiome studies carried out in Steinernema species demonstrated that several bacterial symbionts could be related with pathogenicity (not only mediated by Xenorhabdus) and distinctive antimicrobial activity [132], bringing a new perspective to the complexity of symbionts' interactions and the understanding of the pathobiome and its utility for pest control. Anyhow, native EPN isolates can provide a new source for pest control, directed to specific environmental conditions. For instance, the EPN isolates reported here showed high virulence (>90%) against three weevil species (Cosmopolites sordidus, Metamasius hemipterus, and M. hebetatus), which are persistent pests on plantain and banana plantations. EPN recoveries from polluted soils could be under strong selective pressures, favoring the dominance of specific alleles, not only in the EPN but also in the associated bacteria, which is critical for their adaptation and reproductive success.

Phylogenetic Inferences at Ribosomal Loci Enlightens Natural History of EPN
The phylogenetic studies carried out on EPN have proven to be conflicting since the application of several descriptors, morphological [133,134] and molecular [62,[72][73][74]95,104,123], has provided alternative inferences to the understanding of EPN's natural history. Extensive exploration has focused on the search for the optimal loci applicable for species identification and delimitation. Since the studies of Szalanski et al. [135], who identified agreement between morphological and molecular phylogenies, further studies converged in the selection of the ITS and LSU ribosomal loci as the appropriate markers [70,72,95,104,116,136], currently denoted as DNA barcodes for EPN. Mitochondrial loci have tentatively been used as optional barcoding regions, yet significant inconsistencies with morphological and ribosomal hypotheses have limited the factual utilization of these loci. Against these expectations, the present study exhibited agreement of the mtDNA reconstruction with ribosomal phylogenies for Carpocapsae, Bicornutum, and Feltiae. Several drawbacks were identified for the mitochondrial genes, such as the reduced information available in the databases (making mtDNA markers less attractive for EPNs), their haplotypic-inheritance nature, and their high mutation rate. They may perform well at population levels to resolve recent evolutionary events, but they are suboptimal for deep phylogenetic studies, restricting the resolution power to infer a phylogenetic scenario due to marker saturation [137]. Due to this, the whole-genome mitochondrial sequencing of S. carpocapsae demonstrated genetic closeness to Ascaris suum and Caenorhabditis elegans, not supporting the phylogeny based on nuclear rDNA [138]. The present study covered the most important loci selected for EPN identification and phylogeny (nuclear and plastid DNA) to identify new EPN strains from the Neotropics.
Phylogenetic inferences provided stable tree topologies for the ITS and LSU loci at the clade level. The most comprehensive phylogeny for EPNs was provided by Hunt and Nguyen [62], in which all the extant species were included and analyzed through ribosomal loci. Following this hypothesis, our isolates were compared with this whole dataset grouping the Carpocapsae and Costaricense clades as well as the phylogeny obtained by their symbiotic bacteria, confirming the species obtained here (S. carpocapsae and S. costaricense, respectively). Despite the strong agreement and accuracy in species identification and delimitation provided by ITS and LSU, tree topologies were not completely consistent at the clade level [62,70,74]. This discrepancy was also identified in the present study, where Glaseri, Khoisane, and Bicornotum clades were the most unstable clades for both ribosomal loci. DNA sequence variability could be an important factor underlying these disagreements, since the high nucleotide variability found at ITS [74,139] is not of utility in resolving relationships at deep phylogenetic levels. However, it still may be useful in discriminating species within clades [62]. Previous studies reported the instability of tree topologies when analyzed with an ITS ribosomal locus [72]. The discrepancies related to the ITS locus were reported at the intra-individual level, where alternative haplotypes could be concurrently identified within a single genotype, increasing the differentiation of some species in the phylogenetic inferences [139]. In that sense, it is important to consider the sequencing of the ITS and LSU loci in order to achieve a more reliable phylogenetic reconstruction. At the supra clade level, our clusters were concordant with those reported by Hunt and Nguyen [62], in either the ITS, LSU or concatenated datasets, where superclades I (Glaceri, Karii, Longicaudatum, and Khoisanae) and II (Feltiae, Kushidai, and Monticolum) were in complete agreement, while superclade III (Carpocapsae and Bicornutum) genetically diverged from the others.
The limited resolution of the phylogenetic reconstruction for endosymbiotic bacteria through the 16S rDNA locus is already a well know fact [126,140]. In the present study, the endosymbionts of the Kushidai, Feltiae, and Carpocapsae clades were statistically supported (≥98% bootstrap posterior probability). However, most of the bacterial symbionts remain unresolved. This may be caused by several factors such as the host shift due to the dissociation of the endosymbiont and EPN in the insect cadaver, as a consequence of their amphimictic nature, leading to dissociation from the immediate ancestor [73].
Moreover, the spread of bacterial endosymbionts, naturally occurring via infected insects moving to new environments [121,140], or artificially occurring through the trade of commercially available strains for pest biocontrol, enables the spread of conservative traits that provide ecological advantages [73]. This may generate homogenization of the 16S rDNA variants via gene conversion by homologous recombination as a mechanism of concerted evolution [58,126], unleashing the observed complexity for the phylogenetic reconstruction of biotic interactions [141], including the EPNs' bacterial endosymbionts' evolutionary history. Ultimately, these effort to reconstruct the phylogenetic history of the species will prove essential to assist crop protection initiatives, prioritize multi-locality monitoring [142], envision modern plant enhancement strategies [143] for resistance using modern technique [144,145], and improve seed-delivery recommendations [146] together with commercial EPN, among other downstream steps.

•
The identification of indigenous entomopathogenic nematodes in altered environments, such as agricultural soils, sources pest control for local conditions. • The isolation and characterization of Steinernema costaricense is a key finding from our study, taking into account that there are only two reports worldwide, from Central and North America.

•
Here, we report a new S. costaricense strain, which provides clues about the possible wide distribution of this species, not only in the Americas but also in highly intervened environments. • Additionally, the molecular report of its bacterial symbiont was carried out, contributing to the knowledge about the understanding of specific interactions between EPNs and their bacterial endosymbionts.

Perspectives
• This study promotes the still incipient research of EPN in Colombia, where the real diversity of these organisms remains obscured. • An extended sampling would convey more accuracy for comparisons with other studies. Still, EPNs' potential persists for new discoveries, applications for crop protection, and genetic resources from the bacterial endosymbionts that can provide metabolites with a wide spectrum of uses, for agricultural or medicinal purposes. • New approaches should be considered in these organisms, such as next-generation technologies, in order to maximize their utility for pest control management that could be directed to specific phytosanitary issues.