Introduction
The entomopathogenic nematodes (EPNs) Steinernema and Heterorhabditis are classified within the families Steinernematidae and Heterorhabditidae (order Rhabditida), respectively. These nematodes exhibit a unique parasitic strategy for their infective juveniles (IJs) penetrating insect hosts through natural openings, subsequently releasing their symbiotic bacteria – Xenorhabdus (associated with Steinernema) or Photorhabdus (associated with Heterorhabditis) – into the hemocoel (Stuart et al., Reference Stuart, Barbercheck, Grewal, Taylor and Hoy2006). The bacteria rapidly multiply and secrete potent toxins that quickly kill the host insect by inducing septicaemia, after which the nematodes complete their life cycle by feeding on both the bacterial cells and host tissues (San-Blas, Reference San-Blas2013; Stuart et al., Reference Stuart, Barbercheck, Grewal, Taylor and Hoy2006; Stock, Reference Stock2019). Due to their broad host range, rapid life cycle, suitability for mass production, and strong environmental adaptability, these nematodes have become efficient and eco-friendly biocontrol agents in integrated pest management systems (Askary & Abd-Elgawad, Reference Askary and Abd-Elgawad2021; Bhat et al., Reference Bhat, Machado, Abolafia, Ruiz-Cuenca, Askary, Ameen and Dass2023a; Deka et al., Reference Deka, Baruah and Babu2021).
Extensive field surveys have revealed the global distribution of EPNs, particularly in temperate and tropical regions (Bhat et al., Reference Bhat, Chaubey and Askary2020). The documented biodiversity now comprises over 130 valid Steinernema species (Bhat et al., Reference Bhat, Machado, Abolafia, Askary, Půža, Ruiz-Cuenca, Rana, Sayed and Al-Shuraym2023b; Anes et al., Reference Anes, Patil, Babu, Aarthi, Gotyal, Josephrajkumar, Mhatre, Sajan, Gowda and Půža2025; Baniya et al., Reference Baniya, Subkrasae, Ardpairin, Anesko, Vitta and Dillman2024; Pervez et al., Reference Pervez, Eapen and Devasahayam2024) and at least 23 Heterorhabditis species (Machado et al., Reference Machado, Abolafia, Robles, Ruiz-Cuenca, Bhat, Shokoohi, Puza, Zhang, Erb, Robert and Hibbard2025a), highlighting their considerable taxonomic richness. Ongoing worldwide collection efforts continue to isolate novel strains. For instance, a recent survey in Tarim Basin of Xinjiang, China, has discovered over 20 EPN isolates, including a new species, Steinernema tarimense (Zhan et al., Reference Zhan, Tian, Li, Yang, Bao, Zhang, Zhang, Shi, Tomalak, Půža and Guo2025). An integrative taxonomic framework is employed for classifying these nematodes, which incorporates morphological, molecular, phenotypic, and ecological traits and their relationship with symbiotic bacteria.
The advancement of high-throughput sequencing technologies has significantly enhanced our capacity to sequence and analyse both nuclear and mitochondrial genomes (Bernt et al., Reference Bernt, Donath, Jühling, Externbrink, Florentz, Fritzsch, Pütz, Middendorf and Stadler2013; Dierckxsens et al., Reference Dierckxsens, Mardulyn and Smits2017). Mitochondrial genomes (mitogenomes), which are relatively compact (typically 10–20 kb in nematodes), offer distinct advantages for sequencing and phylogenetic reconstruction. As they evolve independently of nuclear genes, mitogenomes provide unique insights into evolutionary history and have become a powerful phylogenetic marker, consistently yielding well-supported phylogenies within Nematoda (Kern et al., Reference Kern, Kim and Park2020; Gendron et al., Reference Gendron, Qing, Sevigny, Li, Liu, Blaxter, Powers, Thomas and Porazinska2024). Moreover, whole-genome sequencing data have demonstrated the ability to resolve phylogenetic relationships with high congruence and robustness, both across deep divergences and at recent taxonomic scales throughout the Nematoda tree of life (Ahmed et al., Reference Ahmed, Roberts, Adediran, Smythe, Kocot and Holovachov2022; Qing et al., Reference Qing, Zhang, Sun, Ahmed, Lo, Bert, Holovachov and Li2025).
Despite their considerable economic importance, both mitochondrial and nuclear genomic resources for Steinernema remain limited. To date, mitogenomes have been reported for only four species: S. carpocapsae (Montiel et al., Reference Montiel, Lucena, Medeiros and Simões2006); S. glaseri, S. kushidai, and S. litorale (Kikuchi et al., Reference Kikuchi, Afrin and Yoshida2016). Similarly, nuclear genomes are available for only eight species: S. feltiae, S. glaseri, S. monticolum, and S. scapterisci (Dillman et al., Reference Dillman, Macchietto, Porter, Rogers, Williams, Antoshechkin and Mortazavi2015); S. carpocapsae (Rougon-Cardoso et al., Reference Rougon-Cardoso, Flores-Ponce, Ramos-Aboites, Martínez-Guerrero, Hao, Cunha and Montiel2016); S. diaprepesi (Baniya et al., Reference Baniya, Huguet-Tapia and DiGennaro2020); S. khuongi (Baniya & DiGennaro, Reference Baniya and DiGennaro2021); and S. hermaphroditum (Schwarz et al., Reference Schwarz, Baniya, Heppert, Schwartz, Tan, Antoshechkin, Sternberg, Goodrich-Blair and Dillman2025).
Prior to the work of Machado et al. (Reference Machado, Muller, Hiltmann, Bhat, Půža, Malan, Castaneda-Alvarez, San-Blas, Duncan, Shapiro-Ilan, Karimi, Lalramliana, Lalramnghaki and Baimey2025b), genomic data for Heterorhabditis were also scarce, with nuclear genomes available only for H. bacteriophora (Bai et al., Reference Bai, Adams, Ciche, Clifton, Gaugler, Kim and Grewal2013; McLean et al., Reference McLean, Berger, Laetsch, Schwartz and Blaxter2018) and H. indica (Bhat et al., Reference Bhat, Somvanshi, Budhwar, Godwin and Rao2022), along with a mitochondrial sequence for H. bacteriophora (Regeai et al., Reference Regeai, Fitzpatrick, Burnell and Kakouli-Duarte2022). Based on the assembled nuclear and mitochondrial genomes of 47 Heterorhabditis isolates representing over 17 described species, Machado et al. (Reference Machado, Muller, Hiltmann, Bhat, Půža, Malan, Castaneda-Alvarez, San-Blas, Duncan, Shapiro-Ilan, Karimi, Lalramliana, Lalramnghaki and Baimey2025b) established a whole-genome-based phylogenetic framework that robustly delineated the taxonomy and systematics of the genus, elucidated its co-phylogenetic relationships with Photorhabdus symbionts, and reconstructed its biogeographic history – thereby strengthening the foundation for deploying these nematodes as biocontrol agents in sustainable agriculture.
In this study, we presented the draft genome assemblies and partial mitochondrial genome sequences of the newly described S. tarimense (Zhan et al., Reference Zhan, Tian, Li, Yang, Bao, Zhang, Zhang, Shi, Tomalak, Půža and Guo2025) and a putative novel Heterorhabditis species (strain XJ-55). These genomic resources provided valuable references for future research with significant evolutionary and agricultural implications. Through comprehensive annotation, we identified and characterised protein-coding genes in both genomes. These data will facilitate comparative genomic analyses to elucidate the evolutionary origins of insect parasitism in nematodes, genetic diversity among entomopathogenic species, and molecular adaptations underlying host-parasite interactions.
Materials and methods
Nematode collection and morphological observation
Specimens of S. tarimense isolate Z32 were collected from soil samples of Populus euphratica forest in Yuli County, Xinjiang, China (41°4′25.0824″N and 86°7′3.3708″E) and maintained in laboratory cultures using Galleria mellonella larvae as hosts. Similarly, specimens of Heterorhabditis sp. XJ-55 were collected from soil samples of a wild walnut forest in Gongliu County, Xinjiang, China (43°29′43.4″N, 82°04′01.2″E). Individual nematodes were heat-killed and temporarily mounted in distilled water on glass slides. Morphological identification and photography were performed using an Olympus BX51 differential interference contrast (DIC) microscope (Olympus Optical, Tokyo, Japan) equipped with an Olympus C5060Wz digital camera.
DNA extraction and sequencing
Genomic DNA was extracted from pools of ten individuals using the REPLI-g Single Cell Kit (Qiagen, Hilden, Germany). DNA concentration was quantified using a Qubit® fluorometer (Invitrogen, Carlsbad, CA, USA). Sequencing libraries were prepared using the TruSeq DNA Sample Prep Kit (Illumina Inc., CA, USA) following the manufacturer’s protocols. Sequencings were performed on the Illumina NovaSeq PE150 platform.
Genome assembly, annotation, and functional prediction
Raw Illumina reads were quality-filtered using FASTP (Chen et al., Reference Chen, Zhou, Chen and Gu2018) with the following parameters: removal of reads <50 bp, trimming of bases with Q-score <20, and exclusion of reads containing >3 ambiguous nucleotides (N). De novo assembly of genome was performed in SPAdes v3.15.5 (Bankevich et al., Reference Bankevich, Nurk, Antipov, Gurevich, Dvorkin, Kulikov and Pevzner2012). Possible non-nematode contaminations were examined using BlobToolKit v4.0.7 (Challis et al., Reference Challis, Richards, Rajan, Cochrane and Blaxter2020), followed by prokaryotic sequence removal with Tiara v1.0.3 (Karlicki et al., Reference Karlicki, Antonowicz and Karnkowska2022). Assembly metrics were evaluated using QUAST v5.3.0 (Gurevich et al., Reference Gurevich, Saveliev, Vyahhi and Tesler2013), with genome completeness assessed via BUSCO v5.8.3 (Seppey et al., Reference Seppey, Manni and Zdobnov2019) against the Nematoda odb12 lineage (accessed 2025-07-01). GenomeScope 2.0, Smudgeplot 0.2.1 (Ranallo-Benavidez et al., Reference Ranallo-Benavidez, Jaron and Schatz2020), and findGSE v1.94 (Sun et al., Reference Sun, Ding, Piednoël and Schneeberger2018) were employed for k-mer analysis (k = 21) to estimate genome size, heterozygosity, and repetitiveness. Redundant sequences were removed using RepeatMasker v4.2.2 (Tarailo-Graovac & Chen, Reference Tarailo-Graovac and Chen2009).
Gene prediction of assembled genome was performed in AUGUSTUS v3.5.0 (Stanke & Morgenstern, Reference Stanke and Morgenstern2005) and GeneMark-ES v4.73 (Brůna et al., Reference Brůna, Lomsadze and Borodovsky2020). Reference-based predictions were performed in GeMoMa v1.9 (Keilwagen et al., Reference Keilwagen, Hartung and Grau2019) using Caenorhabditis elegans, S. carpocapsa, and S. hermaphroditum as references. The predictions from different pipelines were summarised in EVidenceModeler v2.1.0 (Haas et al., Reference Haas, Salzberg, Zhu, Pertea, Allen, Orvis and Wortman2008). Functional annotation of predicted genes was performed using eggNOG-mapper v2.1.12 (Cantalapiedra et al., Reference Cantalapiedra, Hernández-Plaza, Letunic, Bork and Huerta-Cepas2021) based on eggNOG database v5.0.2, with functional assignments retrieved against KEGG, GO, and COG categories (Kanehisa et al., Reference Kanehisa, Sato, Kawashima, Furumichi and Tanabe2016). The assembled whole-genome sequences have been submitted to NCBI and are currently under processing (BioProject PRJNA1365490).
Mitochondrial genome assembly and phylogenetic analysis
Mitochondrial genomes were also used to infer the phylogeny placement of S. tarimense and Heterorhabditis sp. XJ-55. Briefly, initial genome assemblies were performed using MitoZ v3.3 (Meng et al., Reference Meng, Li, Yang and Liu2019), and assembly extension was carried out using NOVOPlasty v4.3.5 (Dierckxsens et al., Reference Dierckxsens, Mardulyn and Smits2017) with the cytochrome c oxidase subunit I (COI) gene as the seed sequence. The resulting contigs were manually examined and extended to maximise sequence continuity. Gene prediction and annotation were performed using the MITOS webserver (http://mitos.bioinf.uni-leipzig.de/index.py) (Bernt et al., Reference Bernt, Donath, Jühling, Externbrink, Florentz, Fritzsch, Pütz, Middendorf and Stadler2013) with the invertebrate genetic code. The annotations were further refined through alignment with reference mitogenomes of Steinernema (Kikuchi et al., Reference Kikuchi, Afrin and Yoshida2016) and Heterorhabditis (Regeai et al., Reference Regeai, Fitzpatrick, Burnell and Kakouli-Duarte2022) and were visualised using CGView (Stothard et al., Reference Stothard, Grant and Van Domselaar2019). The sequences of 12 protein-coding genes (PCGs) were aligned using MAFFT v7.526 (Katoh & Standley, Reference Katoh and Standley2013), and concatenated in the order cox1, cox2, nad3, nad5, nad6, nad4L, nad1, atp6, nad2, cob, cox3, and nad4 using Concatenator (Vences et al., Reference Vences, Patmanidis, Kharchev and Renner2022). Phylogeny based on concatenated amino acid sequences of PCGs was inferred by the Maximum Likelihood (ML) analysis as implemented in RAxML 8.2.12 (Stamatakis, Reference Stamatakis2014) under the MtZoa protein substitution model. An isopod-parasitic mermithid nematode, Thaumamermis cosgrovei, was selected as the outgroup. The final assembled mitogenomes have been deposited in GenBank, with the accession number PX591264 for S. tarimense, and PX591264 for Heterorhabditis sp. XJ-55.
Phylogenetic analysis of nematodes and symbiotic bacteria
To validate the morphological identifications, we extracted ribosomal RNA (rRNA) gene sequences from the whole-genome sequencing data of Steinernema tarimense and Heterorhabditis sp. XJ-55, respectively. The analytical pipeline included sequence alignment, read extraction and assembly, and rRNA sequence identification. Briefly, raw reads were aligned against reference rRNA sequences from Steinernema and Heterorhabditis species using NextGenMap v0.5.5 (Sedlazeck et al., Reference Sedlazeck, Rescheneder and Von Haeseler2013). Mapped reads were extracted with SAMtools v1.21 (Li et al., Reference Li, Handsaker, Wysoker, Fennell, Ruan and Homer2009), and rRNA sequences were subsequently assembled de novo using SPAdes v3.15.5 (Bankevich et al., Reference Bankevich, Nurk, Antipov, Gurevich, Dvorkin, Kulikov and Pevzner2012). The resulting contigs were aligned to reference rRNA databases to identify the 18S rRNA gene, 28S rRNA gene, and internal transcribed spacer (ITS) region. Meanwhile, the 16S rRNA gene sequences of the symbiotic bacteria Xenorhabdus and Photorhabdus associated with the two nematode species were screened and extracted from the genomic sequencing data using RiboTaxa v1.5 (Chakoory et al., Reference Chakoory, Comtet-Marre and Peyret2022). The similarity of the obtained rRNA gene sequences of nematodes was evaluated by BLAST analysis against the NCBI GenBank nr/nt and RefSeq databases.
To determine the phylogenetic positions of the nematode isolates and their associated symbiotic bacteria, the ML trees were reconstructed using IQ-TREE v2.4.0 (Minh et al., Reference Minh, Schmidt, Chernomor, Schrempf, Woodhams, von Haeseler and Lanfear2020). For the nematodes, phylogenetic analysis was performed based on the ITS sequences and the COI gene sequences extracted from mitogenomes. The obtained sequences were aligned with representative homologs from the genera Steinernema and Heterorhabditis retrieved from GenBank. For the bacterial symbionts, the 16S rRNA gene sequences were used to reconstruct the phylogeny within the genera Xenorhabdus and Photorhabdus. The extracted sequences have been deposited in GenBank under the following accession numbers: Steinernema tarimense Z32 ITS (PZ392409), Heterorhabditis sp. XJ-55 ITS (PZ392135), Xenorhabdus sp. Z32 16S rRNA (PZ392895), and Photorhabdus sp. XJ-55 16S rRNA (PZ392896).
Results
Morphological observation
Steinernema tarimense (Figure 1 A–G) adults have short stoma, pharynx robust with rounded basal bulb; males monorchid with ventrally curved spicules, gubernaculum fusiform in the first and second generations, tail conoid and slightly ventrally curved, with blunt terminus; females didelphic-amphidelphic with shorter conoid tail bearing a fine mucron in the first generation and longer conoid tail lacking mucron in the second generation; and infective juvenile (IJ) with short body, poorly developed pharynx, and tail elongated conoid.
Morphological characters of Steinernema tarimense (A–G) and Heterorhabditis sp. XJ-55 (H–M). S. tarimense, A: Stoma and pharynx region of first-generation female, excretory pore showing by arrow; B: Tail of first-generation female; C: Tail of second-generation female; D: Stoma and pharynx region of first-generation male; E: Tail region of first-generation male; F: IJ anterior body region; G: IJ tail region. Heterorhabditis sp. XJ-55; H: Stoma and pharynx region of hermaphroditic female; I: Stoma and pharynx region of amphimictic female; J: Amphimictic female tail; K: Male tail region; L: IJ anterior body region; M: IJ tail region (Scale bars, H = 100 μm; A, C, I–M = 50 μm; B, D–G = 10 μm).
Figure 1. Long description
Panel A shows the stoma and pharynx region of the first-generation female Steinernema tarimense, with an arrow marking the excretory pore. Panel B displays the tail of the first-generation female. Panel C presents the tail of the second-generation female. Panel D depicts the stoma and pharynx region of the first-generation male. Panel E shows the tail region of the first-generation male, highlighting spicules. Panel F illustrates the anterior body region of the infective juvenile (I J). Panel G shows the tail region of the I J. Panel H presents the stoma and pharynx region of the hermaphroditic female Heterorhabditis sp. XJ-55. Panel I shows the stoma and pharynx region of the amphimictic female. Panel J displays the tail of the amphimictic female. Panel K depicts the male tail region, with visible spicules. Panel L shows the anterior body region of the I J. Panel M presents the tail region of the I J. Scale bars are 100 micrometers for panel H, 50 micrometers for panels A, C, I to M, and 10 micrometers for panels B, D to G. Each panel focuses on specific morphological features for taxonomic comparison.
Heterorhabditis sp. XJ-55 (Figure 1 H–M) hermaphrodite female has six prominent labial papillae, funnel-shaped stoma with denticles, pharynx with cylindrical corpus, isthmus and pyriform basal bulb with reduced valve; amphimictic female similar to hermaphroditic female but smaller, vulval lips smooth, elliptical and slightly protruding, post anal swelling slightly developed, tail conoid with terminal tip rounded or pointed; male spicules ventrally curved and manubrium slightly rounded, tail with peloderan bursa supported by nine pairs (1–2–3–3) of papillae with enlarged tips; and IJ with slender and elongate body, longitudinal striae presented, labial region with prominent dorsal cuticular tooth, mouth and anus closed, stoma collapsed, pharynx similar to amphimictic female but longer and narrower, and tail long with pointed tip.
Whole genome assembly and annotation
The draft whole-genome assemblies exhibited considerable fragmentation, as expected for complex nematode genomes. A total of 5,099 contigs were assembled for S. tarimense with an N50 of 15,715 bp and a maximum contig length was 183,305 bp. Similarly, a total of 2,385 contigs were assembled for Heterorhabditis sp. XJ-55 with an N50 of 50,921 bp and a maximum contig length of 343,276 bp. Despite this fragmentation, BUSCO analysis against the Nematoda odb12 lineage revealed relatively high completeness for both assemblies: 84.06% for S. tarimense and 92.28% for Heterorhabditis sp. XJ-55 (Table 1).
Features and assemble statistics of the genomes generated in Steinernema tarimense and Heterorhabditis sp. XJ-55
Table 1. Long description
Beginning at the top row, the table lists features in the left column, with S. tarimense values in the middle and Heterorhabditis sp. X J dash 55 values on the right. Number of contigs is 5,099 for S. tarimense and 2,385 for Heterorhabditis sp. X J dash 55. Largest contigs are 183,305 and 343,276 respectively. G C content is 46.33 percent for S. tarimense and 32.74 percent for Heterorhabditis sp. X J dash 55. N 50 value is 15,715 base pairs for S. tarimense and 50,921 base pairs for Heterorhabditis sp. X J dash 55. For contigs greater than or equal to 0 base pairs and 1,000 base pairs, both species have the same counts as their total contigs. For contigs greater than or equal to 5,000 base pairs, S. tarimense has 3,292 and Heterorhabditis sp. X J dash 55 has 2,006. For contigs greater than or equal to 10,000 base pairs, counts are 1,576 and 1,527. For contigs greater than or equal to 25,000 base pairs, counts are 411 and 879. For contigs greater than or equal to 50,000 base pairs, counts are 88 and 402. Total length is 54,647,016 for S. tarimense and 66,911,292 for Heterorhabditis sp. X J dash 55. Complete B U S C O s are 501 (84.06 percent) and 550 (92.28 percent). Complete and single copy B U S C O is 483 (81.04 percent) and 546 (91.61 percent). Complete and duplicate B U S C O s are 18 (3.02 percent) and 4 (0.67 percent). Fragment B U S C O s are 28 (4.70 percent) and 26 (4.36 percent). Missing B U S C O s are 67 (11.24 percent) and 20 (3.36 percent).
Whole genome annotation
The sizes of assembled genomes estimated by k-mer frequencies are 84.27 Mbp for S. tarimense and 75.11 Mbp for Heterorhabditis sp. XJ-55 (Figure 2A, B). Ploidy analysis strongly supported diploid genomes (100% probability for S. tarimense; 93% for Heterorhabditis sp. XJ-55) (Figure 2C, D). Taxonomic annotation of assembled contigs revealed significant contamination from symbiotic and environmental sources. In S. tarimense assembly, 8.4 Mb were assigned to Proteobacteria (potential Xenorhabdus symbionts) and 8.8 Mb to Arthropoda (Figure 2E). Similarly, in Heterorhabditis sp. XJ-55 assembly, 64 Mb were assigned to Proteobacteria (likely symbiotic Photorhabdus), and 3.9 Mb to Arthropoda (probable host DNA) (Figure 2F). Minor contaminations from Bacteroidota, Chordata, and Streptophyta were also detected at low abundance (Figure 2E, F).
The evaluation of genome assembly, genome features, and contamination of S. tarimense (A, C, and E) and Heterorhabditis sp. XJ-55 (B, D, and F). A, B: Genome size estimated by k-mer counting; C, D: Estimation of genome ploidy; E, F: Contamination of assembly shown by phylum-level-annotated GC-coverage plots.
Figure 2. Long description
Panel A, top left, is a line graph titled Genome size estimation by find G S E for S. tarimense. X axis is frequency of k-mer, Y axis is number of k-mers. Two lines are shown: a solid purple line with a peak near frequency 40, and a dashed blue line. A legend at top right lists genome size estimates and parameters. Panel B, top center, is a similar line graph for Heterorhabditis sp. XJ-55, with a single peak near frequency 60 and a corresponding legend. Panel C, top right, is a heatmap titled estimated diploid, showing total coverage of the lower k-mer pair A plus B on the Y axis and normalized minor k-mer coverage B over A plus B on the X axis. A color scale at top right indicates k-mer pairs, with a dense cluster at low X and Y values labeled A B. Marginal histograms are shown above and to the right. Panel D, top far right, is a similar heatmap for Heterorhabditis sp. XJ-55, with a dense cluster at higher X and Y values labeled A A B. Panel E, bottom left, is a scatterplot with GC content on the X axis and mean per-base coverage on the Y axis, both with log scales. Colored points represent phylum-level annotation, with a legend at top right listing total, Nematoda, Proteobacteria, Arthropoda, Chordata, Ascomycota, Streptophyta, Mollusca, and other. Marginal histograms show sum length distributions by GC content and coverage. Panel F, bottom right, is a similar scatterplot for Heterorhabditis sp. XJ-55, with a different distribution of phyla and coverage. Legends and axes are consistent across panels.
Genome annotation predicted 14,658 protein-coding genes in S. tarimense and 12,274 in Heterorhabditis sp. XJ-55. Functional classification using the KEGG database demonstrated remarkably conserved pathway distributions between the two species (Figure 3A, B). In primary classification, the Metabolism pathway dominated both species (2,474 in S. tarimense; 2,640 genes in Heterorhabditis sp. XJ-55), and subsequent pathways exhibited consistent rankings: Organismal Systems (2,057 vs. 2,241); Cellular Processes (1,194 vs. 1,325); Environmental Information Processing (1,525 vs. 1,599). In secondary classification, identical pathway patterns were also observed between the two species. The Global and overview maps were the most abundant, followed by the Signal transduction, Endocrine system, Transport and catabolism, and Translation.
Gene function annotation of S. tarimense (A, C, and E) and Heterorhabditis sp. XJ-55 (B, D, and F). A, B: KEGG annotation; C, D: GO annotation; E, F: COG annotation.
Figure 3. Long description
Panel A and B: Horizontal stacked bar charts with x-axis labeled Counts and y-axis listing KEGG pathways. S. tarimense (A) and Heterorhabditis sp. XJ-55 (B) show highest counts in Metabolism and Cellular Processes, with color-coded categories: Organizational Systems (blue), Metabolism (yellow), Genetic Information Processing (green), Environmental Information Processing (light blue), Cellular Processes (orange). Panel C and D: Vertical bar charts with x-axis labeled Description and y-axis Number of Genes, showing GO annotation for S. tarimense (C) and Heterorhabditis sp. XJ-55 (D). Bars are grouped and colored by Ontology: Biological Process (green), Cellular Component (orange), Molecular Function (blue). Most genes are annotated under Biological Process. Panel E and F: Vertical bar charts with x-axis labeled by COG function letters and y-axis Number of genes, for S. tarimense (E) and Heterorhabditis sp. XJ-55 (F). Tallest bars are for category S (Function unknown), followed by J (Translation, ribosomal structure and biogenesis), and others. COG function legend is provided at the right of panel F.
The GO and COG functional annotations both demonstrated that S. tarimense and Heterorhabditis sp. XJ-55 exhibited highly similar gene distribution patterns. The GO analysis revealed three predominant functional categories: Regulation of cellular process (Biological Process), Intracellular organelle (Cellular Component), and Heterocyclic compound binding (Molecular Function) (Figure 3C, D). In the COG classification, the majority of genes were categorised as having unknown functions. Among the annotated genes, those involved in Signal transduction mechanisms were most prevalent, followed by genes associated with Posttranslational modification, protein turnover, and chaperones (Figure 3E, F).
Mitochondrial genome assembly and gene content
The complete mitochondrial genome of Steinernema tarimense (13,836 bp) (Figure 4A) and the partial mitochondrial genome of Heterorhabditis sp. XJ-55 (16,865 bp) (Figure 4B) were successfully assembled, both comprising a complete set of 12 PCGs (atp6, cob, cox1-3, nad1-6, and nad4L), 22 tRNA genes and two rRNA genes. In the S. tarimense mitogenome, the recovered 22 tRNAs ranged from 54 bp (trnS1 and trnS2) to 62 bp (trnK), and two rRNA genes were positioned as rrnL (949 bp) between trnH and nad3, and rrnS (694 bp) between trnE and trnS2. In the Heterorhabditis sp. XJ-55 mitogenome, the recovered 22 tRNA genes ranged from 54 bp (trnP, trnV, trnS1, and trnH) to 63 bp (trnK), and two rRNA genes were located as rrnL (964 bp) between trnH and nad3, and rrnS (703 bp) between trnE and trnS2.
Schematic representation of the mitochondrial genomes of S. tarimense (A) and Heterorhabditis sp. XJ-55 (B). Protein-coding genes, ribosomal RNA (rRNA) genes, and transfer RNA (tRNA) genes are denoted by orange, red, and blue blocks, respectively. All genes are encoded on the anticlockwise strand. Uncoloured segments represent non-coding regions.
Figure 4. Long description
Panel A on the left shows the complete mitochondrial genome of Steinernema tarimense as concentric circles. The innermost ring displays a scale from 1 kbp to 14 kbp, with the genome size labeled as 13836 bp. The next ring outward shows GC content in green and GC skew in black. The outermost ring contains color-coded gene blocks: orange for protein-coding genes (nad1, nad2, nad3, nad4, nad4L, nad5, nad6, atp6, cox1, cox2, cox3, cob), red for rRNA genes (rrnS, rrnL), and blue for tRNA genes (trnS2, trnN, trnY, trnE, trnW, trnV, trnP, trnT, trnL, trnK, trnC, trnM, trnG, trnF, trnQ, trnA, trnR, trnH, trnD, trnI, trnS1, trnL2, trnT, trnG, trnL1). Uncolored segments represent non-coding regions. All genes are annotated on the anticlockwise strand. Panel B on the right shows the partial mitochondrial genome of Heterorhabditis sp. XJ-55, with a similar layout. The innermost ring shows a scale from 1 kbp to 18 kbp, with the genome size labeled as 16865 bp. The next ring outward displays GC content in green and GC skew in black. The outermost ring contains orange protein-coding genes (nad1, nad2, nad3, nad4, nad4L, nad5, nad6, atp6, cox1, cox2, cox3, cob), red rRNA genes (rrnS, rrnL), and blue tRNA genes (trnS2, trnN, trnY, trnE, trnW, trnV, trnP, trnT, trnL, trnK, trnC, trnM, trnG, trnF, trnQ, trnA, trnR, trnH, trnD, trnI, trnS1, trnL2, trnT, trnG, trnL1). Uncolored segments indicate non-coding regions. All genes are encoded on the anticlockwise strand.
The nucleotide composition of both mitogenomes exhibited a strong AT bias, with 73.16% A+T content in S. tarimense and 74.69% in Heterorhabditis sp. XJ-55. Codon usage analysis revealed a striking T-richness. In S. tarimense, the most frequently used codons are TTT (11.8%), TTA (7.2%), ATT (5.1%), and TAA (4.8%); similarly, in Heterorhabditis sp. XJ-55, are TTT (12.2%), ATT (6.2%), TAT (6.1%), and GTT (4.3%).
The two species exhibited distinct preferences for mitochondrial start and stop codons. In S. tarimense, ATT served as the start codon for seven PCGs (apt6, cox1, cox2, nad1, nad2, nad4L, and nad5), TTG for four (cox3, cob, nad3, and nad4), and ATA for nad6. Regarding stop codons, TAA was used in seven genes (cob, cox1, cox2, nad1, nad3, nad4, and nad6), TAG in four (atp6, cox3, nad4L, and nad5), and TTA in nad2. In contrast, among the 12 PCGs of Heterorhabditis sp. XJ-55, the start codons were ATT (for cox1, cox3, nad3 and nad5), TTG (for cob, nad1, nad2, nad4L, and nad6), ATG (for cox2), and ATA(for nad4). The stop codons were exclusively TAA (in six genes: atp6, nad1, nad2, nad4, nad5, and nad6) and TAG (in the other six: cob, cox1-3, nad3, and nad4L).
Phylogenetic analysis based on PCGs of mitochondrial genome
Based on the amino acid sequences of 12 protein-coding genes (PCGs), a concatenated maximum likelihood phylogenetic tree was constructed, comprising 39 nematode taxa. These taxa represent multiple suborders – including Dorylaimina, Plectina, Rhabditina, Spirurina, and Tylenchina – and diverse families such as Aphelenchoididae, Neodiplogastridae, Rhabditidae, and Steinernematidae. In the resulting phylogenetic topology (Figure 5), both Steinernema and Heterorhabditis were recovered as highly supported monophyletic genera. Specifically, S. tarimense formed a strongly supported sister relationship with S. kushidai (bootstrap = 96), whereas Heterorhabditis sp. XJ-55 clustered with H. indica and a previously reported H. bacteriophora strain (bootstrap = 100). All taxa were consistently resolved within their respective evolutionary clades, as illustrated in the phylogenetic tree.
Maximum likelihood phylogenetic tree inferred from concatenated amino acid sequences of 12 protein-coding genes. Newly obtained sequence is indicated in bold with asterisk. Scale bars indicate the number of substitutions per site.
Figure 5. Long description
Starting at the bottom left, Thaumamermis cosgrovei forms the root. Moving upward, branches split into two main clades: Enoplea on the far right with Dorylaimina, Plectina, and Trichuridae, and Chromadorea on the left. Chromadorea further divides into Spirurina and Tylenchina, with colored vertical bars marking families such as Rhabditidae, Neodiplogasteridae, Ascarididae, Anisakidae, Rhigonematidae, Steinernematidae, Panagrolaimidae, Aphelenchoididae, Aphelenchidae, Oxyuridae, Heteroderidae, Pratylenchidae, Meloidogynidae, Onchocercidae, Cephalobidae, Longidoridae, and Mermithidae. Taxa names are listed at branch tips, with Heterorhabditis sp. XJ55 and Steinernema tarimense in bold, the latter marked with an asterisk to indicate the newly obtained sequence. Support values are shown at nodes. The scale bar at the bottom left represents 0.3 substitutions per site.
Phylogeny of nematodes and symbiont bacteria
To confirm the taxonomic identity of two species used in this study, the rRNA gene sequences extracted from the genome sequencing were analysed. BLAST results revealed high similarity to known sequences. For S. tarimense, the 18S rRNA gene exhibited 99.61% identity with S. kushidai (LC157426), the 28S rRNA gene showed 99.78% identity with S. tarimense strain Z32 (PQ590688), and the ITS region displayed 99.89% identity with the same strain (PQ590685). Regarding Heterorhabditis sp. XJ-55, sequence identities reached 99.65% for the 18S rRNA gene with H. bacteriophora (AF036593), and 100% for both the 28S rRNA gene and ITS region with isolates P5 (OR398615) and TP-LP12H (OQ211104) of H. bacteriophora, respectively. These molecular data, consistent with morphological observations, unequivocally confirm the genus designation.
Maximum Likelihood (ML) phylogenetic trees were constructed for species in the genera Steinernema and Heterorhabditis, focusing on S. tarimense and Heterorhabditis sp. XJ-55. The Steinernema dataset included 34 species representing 10 clades, with Caenorhabditis elegans and Oscheius myriophilus as outgroups. For Heterorhabditis, the dataset comprised 21 species across 3 clades (including seven strains of H. bacteriophora), with O. myriophilus as the outgroup. In the ITS and COI trees (Figure 6A, B), S. tarimense was resolved as sister to the clade containing S. akhursti and S. kushidai with strong support. In the ITS tree of Heterorhabditis (Figure 6C), Heterorhabditis sp. XJ-55 formed a sister relationship to H. bacteriophora isolates, which grouped as an independent clade. In contrast, in the COI tree (Figure 6D), Heterorhabditis sp. XJ-55 was placed close to H. casmirica (OQ517980), which was sister to the H. bacteriophora isolates. Both trees suggest that Heterorhabditis sp. XJ-55 may represent a novel species.
Maximum-likelihood phylogenies of nematodes based on ITS sequences (A, C) and COI gene sequences (B, D), and of their symbiotic bacteria based on 16S rRNA gene sequences (E, F). (A, B) Steinernema; (C, D) Heterorhabditis; (E) Xenorhabdus; (F) Photorhabdus. Newly obtained sequences are bold with an asterisk. Scale bars indicate the number of substitutions per site.
Figure 6. Long description
Panel A at top left shows a maximum-likelihood phylogenetic tree for Steinernema nematodes based on I T S sequences, with clades labeled S. feltiae, S. kushidai, S. monticolum, S. glaseri, S. karii, S. bicornutum, S. longicaudum, S. carpocapsae, and S. affine. Newly obtained sequences are bold with an asterisk. The scale bar at the bottom indicates substitutions per site. Panel B at top center presents a similar tree for Steinernema using C O I gene sequences, with the same clade structure and labeling. Panel C at top right displays a tree for Heterorhabditis nematodes based on I T S sequences, with clades H. bacteriophora, H. megidis, and H. indica. Panel D at bottom left shows Heterorhabditis based on C O I gene sequences, with the same clade structure. Panel E at bottom center is a phylogenetic tree for Xenorhabdus symbiotic bacteria based on 16 S r R N A gene sequences, with multiple species and strains labeled, and the outgroup Photobacterium leiognathi at the base. Panel F at bottom right shows a tree for Photorhabdus symbiotic bacteria, also based on 16 S r R N A sequences, with species and subspecies labeled, and Xenorhabdus bovienii as the outgroup. All scale bars indicate the number of substitutions per site. The trees show the evolutionary relationships among nematodes and their symbiotic bacteria, with newly obtained sequences highlighted.
The ML phylogenetic trees based on 16S rRNA gene sequences were reconstructed separately for the genera Xenorhabdus and Photorhabdus. The Xenorhabdus analysis included 40 species, with P. luminescens subsp. venezuelensis as the outgroup, whereas the Photorhabdus analysis comprised 33 species, with X. nematophila as the outgroup. In the Xenorhabdus tree (Figure 6E), the Xenorhabdus sp. Z32 was sister to the group of X. beddingii (NR 042822) and X. cabanillasi (NR 042945) (bootstrap = 74), and this branch together with the group of X. magdalenensis (HQ877464) and X. nematophila (NR 042821) formed an independent clade (bootstrap = 74). In the Photorhabdus tree (Figure 6F), the Photorhabdus sp. XJ-55 was sister to P. laumondii subsp. clarkei (MK039078) (bootstrap = 99), and both were grouped with P. africana (OR835571) and P. laumondii subsp. laumondii (NR 028870) as a highly supported independent branch (bootstrap = 96). The phylogenetic analyses confirmed that the Steinernema strain Z32 corresponded to S. tarimense and its symbiont to Xenorhabdus sp. Z32, while the Heterorhabditis strain XJ-55 represented a potential novel species associated with the symbiont Photorhabdus sp. XJ-55.
Discussion
Steinernema tarimense strain Z32 was collected from Populus euphratica riparian forests in the Tarim Basin, while Heterorhabditis sp. strain XJ-55 was isolated from a wild walnut forest in a valley of the Tianshan Mountains in Xinjiang, China. The natural insect hosts for both species remain unknown. As the first EPN species described from this vast region, S. tarimense exhibited several distinctive biological traits compared to its congeners, including an accelerated developmental rate and the production of light black to nearly colourless cadavers in infected Galleria mellonella larvae (Zhan et al., Reference Zhan, Tian, Li, Yang, Bao, Zhang, Zhang, Shi, Tomalak, Půža and Guo2025). These characteristics may represent evolutionary adaptations to the extreme environmental conditions of the Tarim region, which is characterised by aridity, low precipitation, nutrient-poor soils, and limited insect host availability (Zhang et al., Reference Zhang, Wang and Liu2018; Wang & Li, Reference Wang and Li2020; Chen et al., Reference Chen, Zhao and Ma2021). In contrast, the Heterorhabditis sp. XJ-55 examined here was collected from the Tianshan Mountains region, characterised by extensive forests and grasslands (Wang & Zhang, Reference Wang and Zhang2021). Morphologically, it is consistent with generic characters, whereas phylogenetic analyses indicated that it may represent a novel species. By conducting an integrated analysis of nuclear and mitochondrial genomes in these two nematode species, we can precisely resolve their phylogeny and taxonomy, elucidate their environmental adaptations and co-evolutionary history with symbiotic bacteria, and thereby provide theoretical support for their utilisation in sustainable agriculture.
Through comprehensive bioinformatic analysis, we successfully generated draft genome assemblies and partial mitochondrial genome sequences of S. tarimense and Heterorhabditis sp. XJ-55. The fragmented nature of our genome assemblies, characterised by numerous contigs and modest N50 values (15,715 bp for S. tarimense; 50,921 bp for Heterorhabditis sp. XJ-55), reflects well-documented challenges in de novo assembly of non-model eukaryotic genomes (Alkan et al., Reference Alkan, Sajjadian and Eichler2011; Bankevich et al., Reference Bankevich, Nurk, Antipov, Gurevich, Dvorkin, Kulikov and Pevzner2012). These results are consistent with recent Illumina-based assemblies of 60 nematode species, which ranged from 484 to 33,537 contigs (Qing et al., Reference Qing, Zhang, Sun, Ahmed, Lo, Bert, Holovachov and Li2025). Notably, despite the fragmentation, our assemblies achieved high completeness scores (BUSCO: 84.06% for S. tarimense; 92.28% for Heterorhabditis sp. XJ-55), suggesting they contain nearly complete gene complements suitable for comparative genomic analyses. Moreover, the assembled genome of Heterorhabditis sp. XJ-55 is 66.91 Mb with a GC content of 32.74%, both of which fall within the expected range for Heterorhabditis (genome size ~64–77 Mb, GC content ~33%; Machado et al., Reference Machado, Muller, Hiltmann, Bhat, Půža, Malan, Castaneda-Alvarez, San-Blas, Duncan, Shapiro-Ilan, Karimi, Lalramliana, Lalramnghaki and Baimey2025b), indicating that the contamination has been effectively resolved.
Taxonomic annotation of the assembled genomes revealed substantial contamination from symbiotic and environmental organisms. Specifically, we identified significant sequence contributions from Proteobacteria (8.4 Mb in S. tarimense; 64 Mb in Heterorhabditis sp. XJ-55) and Arthropoda (8.8 Mb and 3.9 Mb, respectively), a finding consistent with prior genomic studies on entomopathogenic nematodes (Dillman et al., Reference Dillman, Macchietto, Porter, Rogers, Williams, Antoshechkin and Mortazavi2015; McLean et al., Reference McLean, Berger, Laetsch, Schwartz and Blaxter2018). These contaminants likely arise from the nematodes’ obligate bacterial symbionts (Xenorhabdus/Photorhabdus) residing in the digestive tract, residual DNA from the Galleria mellonella host used in the laboratory culture, or environmental DNA co-extracted during sampling. While such biological contamination complicates genome assembly, it conversely validates key ecological and life-history traits of the nematodes. Moreover, when properly characterised, these sequences can provide valuable incidental data for investigating host-symbiont coevolution (Dillman et al., Reference Dillman, Macchietto, Porter, Rogers, Williams, Antoshechkin and Mortazavi2015).
We successfully extracted the 16S rRNA gene sequences of the symbiotic bacteria Xenorhabdus and Photorhabdus from the genomic sequencing data of associated nematode hosts. Phylogenetic analysis based on these sequences provided sufficient resolution to assign both symbionts to known genera and to identify their closest relatives. The distinct phylogenetic positions of Xenorhabdus sp. Z32 and Photorhabdus sp. XJ-55, along with the concurrent discovery of their novel nematode hosts, suggest the possibility of parallel host–symbiont diversification. Nevertheless, final taxonomic conclusions require more comprehensive genomic and phenotypic analyses.
Our study confirms that S. tarimense and Heterorhabditis sp. XJ-55 shared identical mitochondrial gene arrangements with other reported species within their respective genera, indicating a highly conserved gene order in these nematode lineages (Kumar et al., Reference Kumar, Singh and Perry2021; Zhang et al., Reference Zhang, Liu and Li2023). Both mitochondrial genomes exhibited a marked A+T bias, consistent with patterns observed in other nematodes (Li et al., Reference Li, Wang and Chen2022; Gendron et al., Reference Gendron, Qing, Sevigny, Li, Liu, Blaxter, Powers, Thomas and Porazinska2024), which may reflect evolutionary constraints on mitochondrial DNA structure. Furthermore, variations in the usage of start and stop codons in protein-coding genes suggest potential plasticity in the mitochondrial genome, possibly associated with adaptive evolution (Wang et al., Reference Wang, Zhao and Liu2024).
The genomic resources presented in this study constitute a valuable dataset for enhancing the application of EPNs in sustainable pest management (Bai et al., Reference Bai, Adams, Ciche, Clifton, Gaugler, Kim and Grewal2013; Dillman & Sternberg, Reference Dillman and Sternberg2012). By systematically characterising genes implicated in host infection, symbiosis, and environmental adaptation, this work will facilitate the development of targeted strategies to optimise EPN efficacy against key agricultural pests (Lu et al., Reference Lu, Baiocchi and Dillman2016). Additionally, the mitochondrial genomes can serve as valuable phylogenetic markers for reconstructing evolutionary history and clarifying the origins of parasitism within the Rhabditida order.
Despite expanding the genomic resources available for EPNs, this study has certain limitations. The current genome assemblies remain fragmented and contain detectable contamination, reflecting common challenges in non-model organism genomics (Bai et al., Reference Bai, Adams, Ciche, Clifton, Gaugler, Kim and Grewal2013; Dillman et al., Reference Dillman, Macchietto, Porter, Rogers, Williams, Antoshechkin and Mortazavi2015). Future work would benefit from employing long-read sequencing technologies (e.g., PacBio or Oxford Nanopore) to generate more contiguous and complete genomes while minimising the inclusion of exogenous DNA (Koren et al., Reference Koren, Walenz, Berlin, Miller, Bergman and Phillippy2017; Miller et al., Reference Miller, Staber, Zeitlinger and Hawley2023). Additionally, integrated functional analyses combining transcriptomic and proteomic approaches will be essential to validate the roles of candidate genes in nematode pathogenicity and symbiosis.
Conclusion
This study presented the first genomic resources for entomopathogenic nematodes from Xinjiang, China, featuring Steinernema tarimense and a novel Heterorhabditis candidate species (XJ-55). Despite fragmented nuclear assemblies, high BUSCO completeness (84.06% and 92.28%) confirmed their utility for comparative genomics. Mitochondrial genomes revealed conserved structure with distinct A+T bias. These resources support phylogenetic studies and genetic adaptation research, providing a foundation for developing enhanced biocontrol strategies in sustainable agriculture.
Author contribution
FZ, HL, and XQ conceived the concept of the study. FZ and WG collected the samples. FZ, IMM, and HL performed the morphological identification and species description. CS and YR did sequencing and bioinformatic analysis. FZ, CS, YR, and HL participated in the drafting of the manuscript. FZ, HL, WG, and XQ carried out the critical revision of the manuscript. All authors read and approved the final manuscript.
Acknowledgements
This study was supported by the National Natural Science Foundation of China (Grant No. 32160377), the National Science and Technology Fundamental Resources Investigation Program of China (2024FY100400), the Project of Fund for Stable Support to Agricultural Sci-Tech Renovation (xjnkywdzc-2025002-08), and the Public Welfare Project of Xinjiang Uygur Autonomous Region (No. KY2024024).
Competing interests
The authors declare no competing interests.
Ethics approval and consent to participate
Not applicable.
Consent for publication
Not applicable.
Open access
