Multi-locus phylogenetic analyses uncover species boundaries and reveal the occurrence of two new entomopathogenic nematode species, Heterorhabditis ruandica n. sp. and Heterorhabditis zacatecana n. sp.

Abstract Species of the nematode genus Heterorhabditis are important biological control agents against agricultural pests. The taxonomy of this group is still unclear as it currently relies on phylogenetic reconstructions based on a few genetic markers with little resolutive power, specially of closely related species. To fill this knowledge gap, we sequenced several phylogenetically relevant genetic loci and used them to reconstruct phylogenetic trees, to calculate sequence similarity scores, and to determine signatures of species- and population-specific genetic polymorphism. In addition, we revisited the current literature related to the description, synonymisation, and declaration as species inquirendae of Heterorhabditis species to compile taxonomically relevant morphological and morphometric characters, characterized new nematode isolates at the morphological and morphometrical level, and conducted self-crossing and cross-hybridization experiments. The results of this study show that the sequences of the mitochondrial cytochrome C oxidase subunit I (COI) gene provide better phylogenetic resolutive power than the sequences of nuclear rRNA genes and that this gene marker can phylogenetically resolve closely related species and even populations of the same species with high precision. Using this gene marker, we found two new species, Heterorhabditis ruandica n. sp. and Heterorhabditis zacatecana n. sp. A detailed characterization of these species at the morphological and morphometric levels and nematode reproduction assays revealed that the threshold for species delimitation in this genus, using COI sequences, is 97% to 98%. Our study illustrates the importance of rigorous morphological and morphometric characterization and multi-locus sequencing for the description of new species within the genus Heterorhabditis, serves to clarify the phylogenetic relationships of this important group of biological control agents, and can inform future species descriptions to advance our efforts towards developing more tools for sustainable and environmentally friendly agriculture.


Abstract
Species of the nematode genus Heterorhabditis are important biological control agents against agricultural pests. The taxonomy of this group is still unclear as it currently relies on phylogenetic reconstructions based on a few genetic markers with little resolutive power, specially of closely related species. To fill this knowledge gap, we sequenced several phylogenetically relevant genetic loci and used them to reconstruct phylogenetic trees, to calculate sequence similarity scores, and to determine signatures of species-and population-specific genetic polymorphism. In addition, we revisited the current literature related to the description, synonymisation, and declaration as species inquirendae of Heterorhabditis species to compile taxonomically relevant morphological and morphometric characters, characterized new nematode isolates at the morphological and morphometrical level, and conducted selfcrossing and cross-hybridization experiments. The results of this study show that the sequences of the mitochondrial cytochrome C oxidase subunit I (COI) gene provide better phylogenetic resolutive power than the sequences of nuclear rRNA genes and that this gene marker can phylogenetically resolve closely related species and even populations of the same species with high precision. Using this gene marker, we found two new species, Heterorhabditis ruandica n. sp. and Heterorhabditis zacatecana n. sp. A detailed characterization of these species at the morphological and morphometric levels and nematode reproduction assays revealed that the threshold for species delimitation in this genus, using COI sequences, is 97% to 98%. Our study illustrates the importance of rigorous morphological and morphometric characterization and multi-locus sequencing for the description of new species within the genus Heterorhabditis, serves to clarify the phylogenetic relationships of this important group of biological control agents, and can inform future species descriptions to advance our efforts towards developing more tools for sustainable and environmentally friendly agriculture.
Nematodes of the genus Heterorhabditis Poinar, 1976 are soil-dwelling organisms that parasitize and kill certain small arthropods, mainly insects (Kaya and Gaugler, 1993). Their lifestyle is particularly interesting as they establish an obligated, mutualistic symbiosis with entomopathogenic bacteria of the genus Photorhabdus (Clarke, 2020;Machado et al., 2018). Nematodes colonize their prey, and upon sensing unknown chemical cues, they release their symbiotic bacterial partners inside the bodies of the infected organisms (Ciche et al., 2008;Dillman et al., 2012). The bacteria establish, multiply and produce an arsenal of immunosuppressors, lytic enzymes, and toxins that kill the infected organism and predigest its tissues, which serve as food for the bacteria and the nematodes (Shankhu et al., 2020;Tobias et al., 2016;Vlisidou et al., 2019). The nematodes grow, reproduce, and, upon resource depletion, reestablish symbiosis with Photorhabdus bacteria, and abandon the consumed cadavers in search for new prey (Somvanshi et al., 2012). Given this peculiar lifestyle, this deadly symbiotic pair is commonly used as a biocontrol agent in agricultural settings (Kajuga et al., 2018;Paddock et al., 2021;Toepfer and Zellner, 2017;Zhang et al., 2019). In addition, given the enormous biosynthetic capacity of Photorhabdus bacteria, they are of great medical, agricultural, and biotechnological importance (Blackburn et al., 1998;Bode, 2009;Hill et al., 2020;Joyce and Clarke, 2003;Lacey and Georgis, 2012;Machado et al., 2018;Machado et al., 2020;Tobias et al., 2018).
The number of described species of the genus Heterorhabditis is steadily growing, mainly boosted by recent advances in genomics. Up to now, the genus includes between 16 and 21 valid species, several synonymized species and some species inquirendae (Boemare et al., 1993;Hunt and Nguyen, 2016;Maneesakorn et al., 2011;Sudhaus, 2011;Tóth and Lakatos, 2008). Given the discrepancy in the number of recognized valid species and the increasing number of synonymized species, a throughout revision of the current literature related to the description, synonymisation, and declaration as species inquirendae of Heterorhabditis species may help to determine the actual number of valid species in this genus. As some species were described prior to the discovery of modern molecular techniques, and therefore the sequences of phylogenetically relevant gene markers are not available, morphological characters play an important role in this context (Andaló et al., 2006;Edgington et al., 2011;Hunt and Nguyen, 2016;Li et al., 2012;Malan et al., 2008;Malan et al., 2014;Nguyen et al., 2004Nguyen et al., , 2006Nguyen et al., , 2008Pereira, 1937;Phan et al., 2003;Poinar and Veremchuk, 1970;Poinar, 1971Poinar, , 1976Poinar et al., 1987Poinar et al., , 1992Poinar, 1990;Stock et al., 2002).
Ribosomal RNA (rRNA) gene sequences such as ITS sequences and the sequences of the D2-D3 expansion segments of the 28S rRNA are traditionally used for identification purposes and for novel taxonomic status descriptions of the species of the genus Heterorhabditis (Adams et al., 1998;Campos-Herrera et al., 2011;Li et al., 2012;Malan et al., 2008;Nguyen et al., 2008;Rana et al., 2020;Spiridonov and Subbotin, 2016). As a recently evolved group, marginal variations in the rRNA gene sequences are expected in this genus, which limits the use of these genetic markers for taxonomic purposes, especially of closely related species (Blaxter et al., 1998;Blouin, 2002;Haag et al., 2018). In addition, the use of sequences containing several ambiguous nucleotides, potentially arisen from sequencing errors and/or poor quality-control, leads to erroneous taxonomic affiliations, as it is exemplified by the relatively high number of synonym species in the genus Heterorhabditis (Dhakal et al., 2020;Hunt and Nguyen, 2016). The use of mitochondrial DNA such as COI sequences, the gold standard taxonomic marker for species delimitation in the Kingdom Animalia, may help to overcome the taxonomic limitations of rRNA gene sequences. However, this taxonomic marker has been used only sporadically for identification purposes, barely used for taxonomy, and never used to describe new Heterorhabditis species (Chaubey et al., 2016;Hebert et al., 2003;Joyce et al., 1994a;Kuwata et al., 2007). As a consequence, the availability of COI sequences for this genus re mained very limited for several years, limiting our understanding of the phylogenetic relationships of this genus (Chaubey et al., 2016;Dhakal et al., 2020;Kuwata et al., 2007).
To improve our understanding on the phylogenetic relationships of Heterorhabditis nematodes, to determine the most suitable genetic markers for the rapid and reliable identification of the species of this genus, specially of closely related species, and to determine species boundaries in this genus, we generated nucleotide sequences of several phylogenetically relevant gene markers and used them to reconstruct phylogenetic trees, to calculate sequence similarity scores, and to determine signatures of species-and population-specific genetic polymorphism. To improve our understanding on the taxonomic relationships of Heterorhabditis nematodes, we revisited the current literature related to the description, synonymisation, and declaration as species inquirendae of Heterorhabditis species to compile taxonomically relevant morphological and morphometric characters, characterized new nematode isolates at the morphological and morphometrical level, and conducted self-crossing and cross-hybridization experiments. Our study illustrates the importance of multi-locus sequencing for the characterization of new species within the genus Heterorhabditis, serves to clarify the phylogenetic relationships of these important biological control agents, and can inform future species descriptions to advance our efforts towards developing more tools for sustainable and environmentally friendly agriculture.

Nematode origin
Heterorhabditis nematodes used in this study were collected by us during different nematode collection campaigns carried out in Rwanda, Mexico, and India, or were collected by different collaborators at different locations around the world (Table S1) (Bai et al., 2013;Bhat et al., 2021b;Carrera, 2015;Fallet et al., 2020;Mukuka et al., 2010;Rana et al., 2020;Yan et al., 2016).

Self-crossing and cross-hybridization experiments
Self-crossing and cross-hybridization experiments were carried out on lipid agar plates as described by Dix et al. (1992). Heterorhabditis ruandica n. sp. Rw14_N-C4a and H. zacatecana n. sp. MEX-39 were self-crossed, hybridized with each other and with H. bacteriophora CH21 (Rana et al., 2020). For this, one second-generation male and one second-generation virgin female were placed on lipid agar plates (35 mm diam.) and incubated at 27°C. Ten independent plates per crossing type were set. Progeny production was observed daily for a period of five consecutive days. Experiments were repeated three times under the same conditions.

Nematode molecular characterization and phylogenetic relationships
Genomic DNA from about 10 to 20 thousand nematodes was extracted using the genomic DNA isolation kit following manufacturer's instructions (Norgen Biotek Corp., Thorold, Ontario, Canada). The following genes/genomic regions were amplified by polymerase chain reaction (PCR): the D2-D3 expansion segments of the 28S rRNA, the internal transcribed spacer (ITS) region of the rRNA, the cytochrome c oxidase I (COI), the thin filament (F-actin)-associated protein (unc-87), and the calmodulin 1 (cmd-1). Primers used were selected based on previous publications (Dhakal et al., 2020;Joyce et al., 1994b;Regeai et al., 2009;Subbotin et al., 2006) (Table S2). PCR reactions consisted of 1 µL of genomic DNA, 12.5 µL of EmeraldAmp GT PCR Master Mix (Takara Bio, Shiga, Japan), 0.5 µL of both forward and reverse primers at 10 mM and 10.5 µL of dH 2 O. The PCR reaction was performed using a thermocycler (Mastercycler nexus gradient, Eppendorf, Germany) with the following settings: (i) for ITS and D2-D3, 1 cycle of 1 min at 98°C followed by 35 cycles of 10 sec at 98°C, 30 sec at 50°C, 1 min 30 sec at 72°C, and by a single final elongation step at 72°C for 10 min; (ii) for cmd-1 and unc-87, 1 cycle of 1 min at 98°C followed by 40 cycles of 10 sec at 98°C, 30 sec at 50°C, 30 sec at 72°C, and by a single final elongation step at 72°C for 10 min; (iii) for COI, 1 cycle of 1 min at 98°C followed by 40 cycles of 10 sec at 98°C, 30 sec at 40°C, 30 sec at 72°C, and by a single final elongation step at 72°C for 10 min. PCR was followed by electrophoresis (45 min, 100 V) of 5 μ l of PCR products in a 1% TBA (Tris-boric acid-EDTA) buffered agarose gel stained with SYBR Safe DNA Gel Stain (Invitrogen, Carlsbad, California, USA). PCR products were purified using the FastGene Gel/ PCR extraction kit (Nippon Genetics Co., Japan) and sequenced using reverse and forward primers by Sanger sequencing (Microsynth AG, Balgach, Switzerland). Obtained sequences were manually curated and trimmed and deposited in the NCBI under the accession numbers given in Table S3   To complete this data set and to obtain genomic sequences of nematodes that belong to all the validly described species of the genus Heterorhabditis, we searched the database of the National Center for Biotechnology Information (NCBI) by the Basic Local Alignment Search Tool (BLAST) using the accession numbers of the sequences obtained previously (Dhakal et al., 2020) (Table S3). Resulting sequences were used to reconstruct phylogenetic relationships by the Maximum Likelihood method based on the following nucleotide substitution models: Hasegawa-Kishino-Yano (HKY + I) (cmd-1), Tamura-Nei (TN93 + G + I) (COI), Kimura 2-parameter (K2 + G) (D2-D3 and ITS), and Tamura 3-parameter (T92) (unc-87). To select the best substitution model, best-fit nucleotide substitution model analyses were carried out in MEGA 7 (Hasegawa et al., 1985;Kimura, 1980;Kumar et al., 2016;Nei and Kumar, 2000). Sequences were aligned with MUSCLE (v3.8.31) (Edgar, 2004). The trees with the highest log likelihood are shown. The percentage of trees in which the associated taxa clustered together is shown next to the branches. Initial tree(s) 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. In some cases, a discrete Gamma distribution (+G) was used to model evolutionary rate differences among sites and the rate variation model allowed for some sites to be evolutionarily (+I). The trees are drawn to scale, with branch lengths measured in the number of substitutions per site. Graphical representation and edition of the phylogenetic tree were performed with Interactive Tree of Life (v3.5.1) (Chevenet et al., 2006;Letunic and Bork, 2016).

Symbiotic relationships
The Photorhabdus entomopathogenic bacteria associated with H. ruandica n. sp. Rw14_N-C4a and H. zacatecana n. sp. MEX-39 nematodes were isolated as described by Machado et al. (2019), (2021b). Briefly, Galleria mellonella larvae (Lepidoptera: Pyralidae) were exposed to 100 nematode infective juveniles. Three to four days later, insect cadavers were surface-sterilized and cut open with a blade. Bacteria-digested internal organs were spread onto LB agar plates and incubated at 28°C for 24 to 48 h. Photorhabdus-like colonies were sub-cultured until monocultures were obtained. A single primary form colony was then selected and used for further experiments. Bacteria primary forms were determined by examining colony characteristics and by examining pigments uptake on NBTA plates (LB agar plates supplemented with 25 mg l -1 bromothymol blue and 4 mg l -1 triphenyl-2,3,5-tetrazolium chloride). The strains were further sub-cultured and maintained on LB agar plates at 28°C. To establish their taxonomic identities, we reconstructed phylogenetic relationships based on whole genome sequences of the isolated bacteria and all the different species/ subspecies of the genus Photorhabdus (Machado et al., 2021a, b). To obtain genomic sequences, genomic DNA was extracted and purified using the GenElute Bacterial Genomic DNA Kit (Sigma-Aldrich, Switzerland) following manufacturer's instructions. The resulting DNA was used for library preparation using the TruSeq DNA PCR-Free LT Library Prep (FC-121-3003) kit. Indexed libraries were then pooled at equimolar concentrations and sequenced (2 × 150 bp) on an Illumina HiSeq 3000 instrument. Genomes were assembled using the Bactopia pipeline (Petit and Read, 2020). Briefly, the raw Illumina reads were quality trimmed using Trimmomatic 0.39 (Bolger et al., 2014). The resulting reads were assembled with SPAdes 3.14.1 (k-mer sizes of 31, 51, 71, 91, and 111 bp) (Bankevich et al., 2012). Scaffolds with a mean read-depth smaller than 20% of the median read-depth of the longer scaffolds (≥5,000 bp) as well as scaffolds that were shorter than 200 bp were removed. The final assemblies were polished using Pilon 1.22 (Walker et al., 2014). Genome sequences were deposited in the National Centre for Biotechnology Information. Accession numbers are listed in Table S4. Phylogenetic relationships were reconstructed based on the assembled genomes and the genome sequences of all validly published species of the genus (Machado et al., 2021a, b). For this, core genome alignments were created using Roary 3.6.2 (Page et al., 2015). Using this alignment, a maximum likelihood tree was constructed using Fasttree 2.1.10 based on the Jukes-Cantor + CAT nucleotide evolution model (Price et al., 2009(Price et al., , 2010. These analyses were carried out in Galaxy (Afgan et al., 2018). Whole genome sequence similarities were calculated by the digital DNA-DNA hybridization (dDDH) method using the recommended formula 2 of the genome-togenome distance calculator (GGDC) web service of the Deutsche Sammlung von Mikroorganismen und Zellkulturen GmbH (DSMZ) (Auch et al., 2010a(Auch et al., , 2010bMeier-Kolthoff et al., 2013.

Results and discussion
Heterorhabditis ruandica n. sp.

Males
Body 0.65 to 0.86 mm long, C-shaped after fixation. Cuticle almost smooth, with transversal striae poorly developed. Lateral field not visible. Lip region with six lips developed but not fused bearing six acute labial papillae at oral margin and four rounded cephalic papillae at the base of lips. Oral opening almost rounded with thick margins. Amphidial apertures pore-like, ovoid and located posterior to lateral labial papillae. Stoma rhabditoid type, 1.2 to 2.3 times the lip region width, with short cheilostom with poorly refringent rounded cheilorhabdia, short gymnostom with refringent barlike rhabdia, and long stegostom surrounded by the pharyngeal collar and bearing bar-like pro-mesorhabdia and small poorly refringent meta-telorhabdia. Pharynx poorly developed with robust corpus without differentiated metacorpus, short and slightly narrow isthmus and pyriform bulb with poorly visible valvular apparatus. Nerve ring encircling the isthmus at 58 to 75% of neck length, just anterior to basal bulb. Excretory pore located at basal bulb level, located at 61 to 97% of neck length. Cardia poorly developed, surrounded by intestinal tissue. Intestine without differentiations. Cardiac anterior end with thin walls. Genital system monorchic, laterally reflexed. Spicules well-developed, separate, with small angular manubrium, calamus poorly developed, and robust lamina with acute tip, scarcely prominent dorsal hump and poorly developed ventral velum. Gubernaculum with manubrium straight and slightly ventrally curved corpus, 40 to 50% of spicule length. Tail conoid with acute tip, ventrally curved posteriorly, flanked by the bursa. Bursa peloderan, with nine pairs of genital papillae (1 + 2/3 + 3), one of them probably the phasmid: three pairs pre-cloacal (GP1-GP3) and six pairs post-cloacal being three pairs at mid tail length (GP4-GP6) and three pairs (GP7-GP9) terminal; GP1 and GP2 more spaced, GP2 and GP3 closely spaced (Figs. 1-4).

Hermaphroditic females
Body 2.91 to 4.12 mm long, arcuate with general morphology similar to male, having labial papillae very acute and prominent. Nerve ring encircling the isthmus at 56 to 78% of neck length. Excretory pore located at or posterior to basal bulb, located at 67 to 103% of neck length. Genital reproductive system didelphic-amphidelphic with ovaries well developed, reflexed, oviducts and uteri not well visible, vagina very short and vulva small having transverse slit opening. Rectum slender, 0.8 to 1.3 times longer than the anal body diameter. Anus with prominent lips. Tail conoid with acute tip lacking mucro, having cellular part simple at its junction with the hyaline part. Phasmids inconspicuous (Figs. 1-4).
Heterorhabditis ruandica n. sp. IJs can be distinguished from the IJs of H. bacteriophora by the distance between the anterior end and the excretory pore (67-90 vs. 87-110 µm) and the distance between the anterior end and the nerve ring (52-64 vs. 72-93 µm), and by the tail length (49-65 vs. 83-112 µm). The males of Heterorhabditis ruandica n. sp. can be distinguished from the males of H. bacteriophora by the distance from the anterior end to the excretory pore (61-109 vs. 114-130 µm) and by the lower D% value (61-97 vs. 117). Hermaphroditic and amphimictic females also show various morphometric differences (Tables 3-6).

Type host and locality
The type hosts are unknown as the nematodes of this genus can be hosted by different insect species and were isolated from soil samples by the Galleria baiting technique (Bedding and Akhurst, 1975;White, 1927). Nematode strains H. ruandica n. sp. Rw18_M-Hr1a and Rw18_M-Hr1b were collected in the district of Karongi, Western province of the Republic of Rwanda (Decimal degrees coordinates: -2.131500, 29.325467) in a moist habitat along a river bench covered with sweet potato plants. Heterorhabditis ruandica n. sp. Rw14_N-C4a nematodes were collected in a ploughed cropland on terraces in a hilly area near Kanyirandori village, Tare sector, Nyamagabe district, Southern province of the Republic of Rwanda (Decimal degrees coordinates: -2.500000, 29.483333).

Type material
Rw14_N-C4a nematodes are the type material for Heterorhabditis ruandica n. sp. Holotype male, and 15 paratype hermaphrodites, males and amphimictic females and 15 third stage juveniles were deposited in the National Nematode Collection of India, IARI, New Delhi, India. Additional specimens were deposited at the nematode collection of the Department of Animal Biology, Plant Biology and Ecology of the University of Jaén, Spain, under the following slide numbers: Rwa001-01 to -12 (25 hermaphrodite females and 6 juveniles), Rwa002-01 to -05 (8 amphimictic females and 9 males), and Rwa003-01 to -02 (8 juveniles). Nematode cultures are maintained in the Institute of Biology, University of Neuchatel, Switzerland and in the Rwanda Agriculture and Animal Resource Development Board, Rubona, Rwanda.

Etymology
The specific name refers to the country, the Republic of Rwanda (Africa), where the type material, Heterorhabditis ruandica n. sp. Rw14_N-C4a nematodes, used to phenotypically characterize the species, were collected.
Heterorhabditis zacatecana n. sp. Males Body 0.81 to 0.91 mm long, J-shaped after heat killing and body arcuate posteriorly. Cuticle almost smooth, with transversal striae poorly developed. Lateral field not visible. Lip region with six lips poorly developed bearing six acute labial papillae at oral margin and four rounded cephalic papillae at base of lips. Oral opening almost rounded with thick margin. Amphidial apertures pore-like, ovoid and located posterior to lateral labial papillae. Stoma rhabditoid type, 0.9 to 1.6 times the lip region width, with short cheilostom with poorly refringent rounded cheilorhabdia, short gymnostom with refringent bar-like rhabdia, and long stegostom surrounded by the pharyngeal collar and bearing bar-like pro-mesorhabdia and small poorly refringent meta-telorhabdia. Pharynx poorly developed with robust corpus without differentiated metacorpus, short and slightly narrow isthmus and robust pyriform bulb with poorly visible valvular apparatus. Nerve ring encircling the isthmus at 61% to 96% of neck length, just anterior to basal bulb. Excretory pore located at or posterior to the basal bulb, located at 78% to 134% of neck length. Cardia poorly developed, surrounding by intestinal tissue. Intestine without differentiations. Genital system monorchic, laterally reflexed. Spicules well developed, separate, with more or less rounded manubrium, calamus poorly developed, and thinner and slender lamina with acute tip, scarcely prominent dorsal hump and poorly developed ventral velum. Gubernaculum with manubrium slightly ventral curved and straight corpus, 40% to 60% of spicule length. Tail conoid with acute tip, ventrally curved posteriorly, flanked by the bursa. Bursa peloderan, with nine pairs of genital papillae (1 + 2/3 + 3), one of them probably the phasmid: three pairs pre-cloacal (GP1-GP3) and six pairs post-cloacal being three pairs at mid tail length (GP4-GP6) and three pairs (GP7-GP9) terminal; GP1 and GP2 more spaced, GP2 and GP3 closely spaced (Figs. 5-8).

Hermaphroditic females
Body 4.41 to 6.18 mm long, arcuate with general morphology similar to male, having labial papillae more acute and prominent. Genital reproductive system didelphic-amphidelphic with ovaries well developed, reflexed, oviducts and uteri not well visible, vagina very short and vulva small having transverse slit opening. Rectum slender, about 1.5 times longer than the anal body diameter. Anus with prominent lips. Tail conoid with acute tip lacking mucro, having cellular part simple at its junction with the hyaline part. Phasmids inconspicuous (Figs. 5-8).

Amphimictic females
Body similar to, but usually smaller than hermaphrodites, 1.95 to 2.80 mm long. Rectum very long, about twice longer than the anal body diameter. Anus with posterior lip more prominent. Tail conoid with acute tip lacking mucro, having cellular part simple at its junction with the hyaline part (Figs. 5-8).

Type host and locality
The type hosts are unknown as the nematodes of this genus can be hosted by different insect species and were isolated from soil samples by the Galleria baiting technique (Bedding and Akhurst, 1975;White, 1927). Type material MEX-39 nematodes are the type material for Heterorhabditis zacatecana n. sp. Holotype male, 15 paratype and 15 third stage juveniles were deposited in the National Nematode Collection of India, IARI, New Delhi. Additional specimens were deposited in the nematode collection of the Department of Animal Biology, Plant Biology and Ecology of the University of Jaén, Spain, under the following slide numbers: Mex001-01 to -03 (6 hermaphrodite females), Mex002-01 to -04 (8 amphimictic females and 3 males), and Mex003-01 to -04 (14 juveniles). Nematode cultures are maintained in the Institute of Biology, University of Neuchatel, Switzerland.

Etymology
The specific name refers to the Mexican state, Zacatecas, where the type material, Heterorhabditis zacatecana n. sp. MEX-39 nematodes, used to phenotypically characterize the species were collected.

Nematode molecular characterization and phylogenetic relationships
Phylogenetic reconstructions based on nuclear and mitochondrial genes (ITS, D2-D3, COI, umc-87, and cmd-1), either individually or concatenated, confirm that the nematodes of the genus Heterorhabditis are grouped into three major clades: the "Megidis-group", the "Indica-group" and the "Bacteriophora-group", which is consistent with previous studies (Dhakal et al., 2020) (Fig. 9, Fig. S1). The clade of the "Bacteriophoragroup" is, in turn, separated into five subclades. Three of them are composed of already described species: H. beicherriana, H. georgiana, and H. bacteriophora, and two of them are composed of two new, undescribed species, which we named here H. zacatecana n. sp., and H. ruandica n. sp. (Fig. 9, Fig. S1). Clearer phylogenetic separations within the species of the clade of the "Bacteriophora-group" were observed when phylogenies were reconstructed based on COI, ITS, or on concatenated sequences of COI, ITS, and D2-D3 (Fig. 9, Fig. S1). Closer inspection at the ITS, D2-D3 and COI sequences reveals unambiguous genetic differences between the nematodes of the "Bacteriophora-group" (Fig. 10). Sequence similarity scores and nucleotide difference counts show a closer relationship between H. bacteriophora, H. ruandica n. sp., and H. zacatecana n. sp. nematodes (Fig. 11 and Figs. S2-S6). Heterorhabditis ruandica n. sp. and H. bacteriophora share 99.1% and differ in 6 nucleotide positions in the ITS sequences flanked by primers TW81 and AB28, share 99.8% and differ in 1 nucleotide position in the D2-D3 sequences flanked by primers D2A and D3B, and share 94.1 to 94.7% and differ in 18 to 19 nucleotide positions in the COI sequences flanked by primers HCF and HCR ( Fig. 11 and Figs. S2-S6). Heterorhabditis zacatecana n. sp. and H. bacteriophora share 99.4% and differ in 4 nucleotide positions in the ITS sequences flanked by primers TW81 and AB28, share 99.8% and differ in 1 nucleotide position in the D2-D3 sequences flanked by primers D2A and D3B, and share 94.1 to 94.4% and differ in 19 to 20 nucleotide positions in the COI sequences flanked by primers HCF and HCR ( Fig. 11 and Figs. S2-S6). Heterorhabditis ruandica n. sp. and H. zacatecana share 99.7% and differ in 2 nucleotide positions in the ITS sequences flanked by primers TW81 and AB28, share 100% and differ in no nucleotide position in the D2-D3 sequences flanked by primers D2A and D3B, and share 97.6% to 98.2% and differ in 6-8 nucleotide positions in the COI sequences flanked by primers HCF and HCR ( Fig. 11 and Figs. S2-S6). Noteworthy, we observed almost no intraspecific variations within the nematodes of the "Bacteriophora-group" at different genetic loci (Figs. 10,11,. However, the sequences of the COI gene show very interesting signatures of population-specific polymorphism (Figs. 10D-F, 11). Specifically, Heterorhabditis ruandica n. sp. Rw18_M-Hr1a and Rw18_M-Hr1b nematodes that were collected in the same western Rwandan region differ from the Heterorhabditis ruandica n. sp. Rw14_N-C4a nematodes collected in a southern Rwandan region in a transitional nucleotide change (g.1212A > G) (Fig. 10D). Moreover, H. zacatecana n. sp. MEX-39 and MEX-40 nematodes collected in north-central Mexico and H. zacatecana n. sp. MEX-41 nematodes collected in central Mexico differ in three transitional nucleotide changes (g.1257T > C, g.1324T > C, and g.1464A > G) (Fig. 10D-F). Hence, due to its highly conserved species-specific polymorphism, and the consistent population-specific polymorphic patterns, the COI gene emerges as an important phylogenetic marker also for the genus Heterorhabditis, in a similar manner as it is for many other taxonomic groups (Hebert et al., 2003;Pentinsaari et al., 2016).

Interspecific genetic variability within the H. bacteriophora clade
In a recent study, Dhakal et al. (2020) Table S3.
(Figs. S7 and S8). Phylogenetic recon structions show a clear phylogenetic separation between all these haplotypes (Fig. S8). Hence, some of the haplotypes described by Dhakal et al. (2020) represent new species, closely related to H. bacteriophora, and some others likely represent new species, which highlights the power of statistical parsimony network analyses to uncover undescribed species of the genus Heterorhabditis, and supporting previous hypothesis regarding the taxonomic status of these nematode isolates Dhakal et al., 2020;Fallet et al., 2020).

Symbiotic relationships
Up to now, the bacterial genus Photorhabdus Boemare, Akhurst and Mourtant 1993 contains 27 taxa, including species and subspecies . Phylogenetic relationship reconstructions based on whole genome sequences show that the bacterial symbionts isolated from H. zacatecana n. sp. MEX-39 and H. ruandica n. sp. Rw14_N-C4a nematodes, named here as MEX-39 and RW14-46, respectively, show high similarity with two of the already described Photorhabdus species: Photo- Figure 12: Phylogenetic reconstruction based on core genome sequences of Photorhabdus bacterial strains. Numbers at the nodes represent SH-like branch supports. Bar represents average nucleotide substitutions per sequence position. Accession numbers of the genome sequences used for the reconstruction are shown in Table S4.
On the synonymization and declaration of species inquirendae of some species We revised the original publications of all syno nymized species and based on their morphology and molecular data (when available), we reinforce the synonymized status of most of them (Khan et al., 1976;Wouts, 1979;Stock, 1993;Gardner et al., 1994;Liu, 1994;Stock et al., 1996;Plichta et al., 2009;Maneesakorn et al., 2015;Hunt and Nguyen, 2016;Shahina et al., 2017;Dhakal et al., 2020). However, the original description of H. bacteriophora provided by Poinar (1976) shows males with very anterior GP1 while in its synonymized species H. heliothidis (Khan, Brooks & Hirschmann, 1976) Poinar, Thomas & Hess, 1977(=Chromonema heliothidis Khan et al., 1976) the GP1 appears more posterior (Khan et al., 1976;Poinar, 1976). Hence, it is likely that both species are not conspecific. Therefore, we declare H. heliothidis (Khan, Brooks & Hirschmann, 1976) Poinar, Thomas & Hess, 1977 as species inquirenda. Heterorhabditis hoptha and H. poinari were poorly described (Turco, 1970;Kakulia and Mikaia, 1997). Original descriptions lack differentiated description of all diagnostic characters of adult and larval stages. According to this, both species should remain in the list of species inquirendae. Heterorhabditis egyptii and H. hambletoni were described showing all diagnostic characters of adults and larvae stages. According to this, both species are considered valid herein (Pereira, 1937;Abd-Elgawad and Ameen, 2005). The lack of molecular data, however, impairs their inclusion in future phylogenetic studies. Nevertheless, new species description should contrast morphological characters with these species. An updated dichotomous key to identify the species of the genus Heterorhabditis is provided (Fig. 13, Tables 3-6).

On the species of the genus Heterorhabditis
Considering the results of this study and the analyses of all the literature that describes new species of the genus Heterorhabditis, the updated list of the species of the genus, including their status, is as follows. N. Zealand Poinar (1990) as      Poinar (1990) as      Poinar et al. (1992)