Complete mitochondrial genomes of four entomopathogenic nematode species of the genus Steinernema

Nematodes belonging to the genus Steinernema are insect parasites and are used as effective biological agents against soil-dwelling insect pests. Although the full nuclear genomes of multiple Steinernema species have become available recently, mitochondrial genome information for the genus is limited. In this study, we sequenced the complete mitochondrial genomes of four species of Steinernema and analysed their structure, codon usage and phylogenetic relationships. Mitochondrial genomes of Steinernema carpocapsae, S. glaseri, S. kushidai and S. litorale comprised 13,924, 13,851, 15,182 and 21,403 bp, respectively, with highly AT-rich nucleotide contents (AT ratio of 71.05–76.76 %). All the expected genes, including 12 protein-coding genes (encoding ATP6, CYTB, COX1-3, ND1-6 and ND4L), two rRNA genes and 22 tRNA genes were identified in the four genomes. Phylogenetic analyses based on the amino acid sequences of the 12 protein-coding genes identified the Steinernema species as monophyletic, representing a sister clade of Rhabditina and Ascaridida. In addition, they were more closely positioned to other Clade 10 nematodes, including Bursaphelenchus xylophilus, Aphelenchoides besseyi and Panagrellus redivivus, than to Strongyloides species. Gene arrangements and codon usage analyses supported this relationship. Mitochondrial genome comparison of two distinct strains of S. carpocapsae detected high intra-specific diversity. The mitochondrial genomes of four species of Steinernema determined in this study revealed inter- and intra-species divergences/diversities of mitochondrial genomes in this genus. This information provides useful insights into the phylogenetic position of the genus Steinernema within the Nematoda and represents a useful resource for selecting molecular markers for diagnosis and population studies. These data will increase our understanding of the interesting biology of insect parasites.


Background
Nematodes belonging to the genus Steinernema are insect parasites that are used as biological agents to control soil-dwelling insect pests [1,2]. Steinernema nematodes form symbiotic associations with enterobacteria of the genus Xenorhabdus and the ability to kill insects is formed by a complex association of the nematodes and the bacteria, which makes the nematodes and the bacteria an attractive model for studying animalmicrobe symbiosis [3,4].
Steinernema nematodes belong to infraorder Panagrolaimomorpha, which includes, among others, the vertebrateparasitic genus Strongyloides and free-living bacteriovores of the genera Panagrolaimus and Panagrellus [5]. Phylogenetic reconstruction based on the small subunit ribosomal RNA gene placed the Steinernema spp. in Clade IV [6] and, in a later study, in Clade 10 [7] with the Panagrolaimomorpha nematodes and the plant-parasitic/fungivorous order Aphelenchida, which includes Bursaphelenchus xylophilus and Aphelenchoides besseyi.
The number of species in the genus has increased rapidly in recent years because of extensive surveys around the world that aimed to identify more efficient biological control agents (for examples, see [8][9][10]). Currently, about 80 described species can be found in the genus in the NCBI taxonomy database (http://www.ncbi.nlm.nih.gov/Taxonomy/). Considering the great diversity of hosts (insects), further specific diversity in this nematode group is expected.
Mitochondrial genome sequences have been widely used in population or evolutionary studies of nematodes. Mitochondrial genomes provide a rich source of molecular variation and widespread utility in population genetics and evolutionary biology [11]. Furthermore, mitochondria are deeply involved in many biological processes, from energy production to programmed cell death and ageing [12]. Studying mitochondrial genomes is, therefore, essential for our understanding of fundamental biology of eukaryotes.
Recently, the full nuclear genomes of five Steinernema species were sequenced [13], which paved the way to investigate the population structures and genetic factors involved in their key biological processes and parasitism using genome-wide analyses. By contrast, mitochondrial genome information was limited to only one species, S. carpocapsae [14]. Complete mitochondrial genome sequences from other Steinernema species will certainly be an important resource to study their ecology and parasitic biology. In this study, we sequenced the complete mitochondrial genomes of four species of Steinernema and analysed their structure, codon usage and phylogenetic relationships.

Biological materials
Steinernema carpocapsae (strain All), Steinernema glaseri (strain Sds102), Steinernema kushidai (strain Hamakita) and Steinernema litorale (strain IbKt142) were maintained using insect larvae (Galleria mellonella or Anomala cuprea) at the laboratory of plant nematology of Meiji University or Kyushu Okinawa Agricultural Research Center, National Agriculture and Food Research Organization. Infective third-stage juveniles of the nematodes were isolated from the infected insect cadavers using a White trap and further purified by the Berman funnel method, as described previously [15,16].

Library preparation and sequencing
Genomic DNA was extracted from the nematodes using a QIAamp DNA Mini Kit (Qiagen, Tokyo, Japan) according to the manufacturer's instructions. Illumina libraries were constructed using a Nextera DNA Sample Prep Kit (Illumina, San Diego, USA) with 100 ng of DNA, according to the manufacturer's instructions. The libraries were sequenced on an Illumina MiSeq using a v3 Reagent kit (600 cycles), according to the manufacturer's recommended protocol (https://icom.illumina.com/) to produce 300-bp paired-end reads.

Mitochondrial genome assembly
Mitochondrial genomes were reconstructed from the Illumina reads using MITObim ver. 1.6 [17]. Initial assemblies were generated using SGA assembler [18] and mitochondrial fragments in the assembly were identified by BlastX using Caenorhabditis elegans mitochondrial genes as queries. Those fragments were extended by iterative mapping of Illumina short reads using MITObim [14]. Gap regions in the assemblies were PCR amplified using Tks Gflex DNA Polymerase (Takara, Shiga, Japan) and sequenced on an ABI 3130 sequencer with a BigDye Terminator v3.1 Cycle Sequencing Kit (Applied Biosystems, Foster City, USA) or Illumina MiSeq, as mentioned above, to obtain a complete mitochondrial genome sequence. The assembled mitochondrial genomes were annotated for protein-coding, tRNA and rRNA genes using the MITOS web server [19] and by manual curation, with support from sequence similarity to other published mitochondrial genomes using Artemis genome annotation tool [20]. Repeats in the assemblies were detected using Tandem Repeat Finder ver.4.09 [21]. General statistical values were calculated using R (ver3.1.1) and in-house Python scripts.

Phylogenetic analysis
Thirty-nine nematode mitochondrial genomes were selected and included in the analysis, with an emphasis on species from Holterman's Clade 10, to which Steinernema nematodes belong. Amino acid alignments for 12 protein-coding genes were generated separately using MAFFT [22], with options (−L-INS-i), and trimmed by Gblocks [23] with less stringent options to remove nonwell-aligned sites. The best substitution model for each alignment was estimated using ProtTest [24]. All the protein alignments were concatenated and maximum likelihood trees were constructed using RaxML ver. 7.2.8 [25], under the MtArt model, for each gene partition. The tree topologies obtained from ML analyses were evaluated with 1,000 bootstrap pseudoreplications. Bayesian inference was performed using MrBayes 3.2.2 [26] from 0.5 million Markov Chain Monte Carlo generations, under a strict clock model using the Mtrev model for each gene partition. Markov chains were sampled at intervals of 100 generations. The first 0.1 million generations were discarded as 'burn-in' and three independent Markov Chain Monte Carlo runs converged to the same posterior probability.
Sequences of the small subunit ribosomal RNA gene were aligned using MAFFT version 7 [22], with option E-INS-I, and a maximum likelihood tree was obtained using RaxML ver. 7.2.8 [25] under the GTR-gamma model with 500 bootstrap pseudoreplications.

Mitochondrial genome assembly
We sequenced mitochondrial genomes of four Steinernema species: S. carpocapsae, S. glaseri, S. kushidai and S. litorale. The assembly resulted in single circular molecules for all four species with genome sizes ranging from 14 to 21 kb (13,924, 13,851, 15,182 and 21,403 bp for S. carpocapsae, S. kushidai, S. glaseri and S. litorale, respectively) ( Table 1). This size difference mainly reflected the length of the non-coding regions in the genomes. Steinernema litorale and S. glaseri had large non-coding regions (8,137 and 1,962 bp, respectively), while the non-coding regions of S. carpocapsae and S. kushidai were shorter than 1 kb (Fig. 1, Table 1). Similar to mitochondrial genomes of other species, the nucleotide compositions were highly AT rich (AT ratio of 71.05-76.76 %). Steinernema kushidai and S. litorale had a slightly lower AT % than the other two species. Non-coding regions showed higher AT %s than coding regions (Table 1). All four mitochondrial genomes contain 12 protein-coding genes (encoding ATP6, CYTB, COX1-3, ND1-6 and ND4L), two rRNA genes and 22 tRNA genes (Fig. 1). Like most other nematodes, an atp8 gene was not found in the genomes. All 36 genes were encoded on the same strand.

Codon usage
Codon usages of the four mitochondrial genomes are shown in Additional file 1: Table S1. In S. glaseri, S. litorale and S. kushidai, ATT was the most frequently used start codon (seven, five and eight cases, respectively) followed by TTG, ATG and ATA. Two GTTs in S. kushidai were the only start codons that were exceptions to the translation table (transl_table 5; invertebrate mitochondrial code). By contrast, in S. carpocapsae, exceptional start codons were used frequently (TTT for six genes and TTA for two genes). ATT and ATA start codons were observed in only two genes each in S. carpocapsae. TAA and TAG were used as termination codons in many of the genes. A truncated stop codon T was also observed in a small number of genes in all four species. The most frequently used codon in all four species was TTT (phe) followed by ATT (Ile) and TTA (leu). Principal component analysis (PCA) suggested that the overall codon usage patterns were similar in the four species compared with other Clade10 nematodes (Additional file 2: Figure S1). Within the Steinernema, the pattern of S. carpocapsae was similar to that of S. glaseri, whereas S. litorale showed a similar pattern to S. kushidai (Additional file 2: Figure S1).

Non-coding regions
The S. glaseri mitochondrial genome had a long noncoding region (1,525 bp). The non-coding region was highly AT rich (AT ratio of 85.71 %) and contained repeat sequences. Tandem repeat finders identified a short tandem repeat (112 copies of "AT") and a long repeat (two copies of the 458-bp unit) in the sequence (Fig. 1, Additional file 1: Table S2). Two long non-coding sequences were found in the S. litorale mitochondrial genome: one between tRNA-Ser2 and tRNA-Asn (1,249 bp), and the other between tRNA-Asn and tRNA-Tyr (6,260 bp) (Fig. 1). The latter region contained many repeat sequences with variable unit sizes, with the longest unit being 383 bp (Additional file 1: Table S2).

Phylogenetic analysis
A maximum likelihood tree based on 12 protein-coding genes from 39 nematode mitochondrial genomes, with enoplean nematode sequences as outgroups, is shown in Fig. 2  Within the Steinernema clade, S. carpocapsae occupied the basal position and S. glaseri plus S. kushidai was at the most derived position. The relationships of the four Steinernema species inferred from the nuclear 18S rRNA gene were slightly different (Additional file 2: Figure S2). In the 18S rRNA gene tree, S. litorale and S. kushidai were positioned most internally, although S. carpocapsae was placed at the base of the four species.

Gene arrangements
All the genes in the Steinernema mitochondrial genomes were transcribed in one direction, as in other Chromadorea nematodes (Fig. 3). The gene arrangements were also well conserved within the Steinernema species and Rhabditina species, including C. elegans and Pristionchus pacificus, and Ascaridida species, including A. suum and Toxocara canis. Steinernema glaseri and S. litorale showed perfectly identical gene arrangements to C. elegans, whereas S. kushidai and S. carpocapsae were different only in the position of tRNA-His (Fig. 3). Other Clade 10 nematodes, including B. xylophilus, A. besseyi, P. redivivus and H. gingivalis, showed similar gene arrangements to the Steinernema spp. with only few rearrangements of tRNA genes. Strongyloides clade species (including Strongyloides spp. and Parastrongyloides trichosuri) were the exception, showing highly rearranged mitochondrial genomes [27].

Intraspecific variation
A mitochondrial genome of a different strain of S. carpocapsae was previously reported (strain Breton, originating from France) [14]. Comparison of the French sequence with our S. carpocapsae (strain All, originating from Japan) revealed intraspecific diversity. A nucleotide alignment identified 201 single nucleotide variants (SNV) and nine indel positions, which accounted to 1.5 % of the total genome. Among the 201 SNVs, 177 were located in protein-coding regions, 17 in the tRNA region, four in the rRNA region and three in non-coding regions (Table 3). Within the protein-coding regions, 146 variants were synonymous and 31 were nonsynonymous. Indels were identified only in non-coding regions (six indels) and in the rRNA region (three indels) ( Table 3).

Discussion
Although the full nuclear genomes of five Steinernema species were published recently [13], information on the mitochondrial genomes of this genus was limited. In this study, we sequenced four complete mitochondrial genomes using next-generation sequencing techniques. To the best of our knowledge, this is the first study employing multi-species comparison of mitochondrial genomes in this genus.
Inconsistencies between mitochondrial and nuclear rRNA gene phylogenies of Nematoda have been reported previously [27,28]. The mitochondrial gene tree obtained in this study showed a similar topology to those in previous studies. For example, the non-monophyly of Clade III (sensu Blaxter et al. 1998 [6]) and the nonmonophyly of Tylenchomorpha (sensu De Ley & Blaxter, 2002 [5]) were also observed in our tree. In our mitochondrial gene phylogeny, the four Steinernema species were clustered as a sister clade of Rhabditina plus Ascaridida, with strong support (Fig. 2). Although Steinernema and Strongyloides have been placed in the same superfamily (Strongyloidoidea), the phylogeny did not support closer relationships between the two genera. This finding agrees with recent studies [14,[28][29][30]. In this study, we included more Steinernema species and Clade 10 species in the phylogenetic analysis, which confirmed the distant relationship of Steinernema spp. to Strongyloides spp. Comparisons of codon usage or gene arrangements are useful to obtain insights into phylogenetic relationships. The codon usage in Steinernema species was distinguishable from that of other related species in the PCA plot (Additional file 2: Figure S1). The gene arrangements suggested that Steinernema is more closely related to Rhabditina and Ascaridida species than to Strongyloides species (Fig. 3). These results are consistent with the phylogeny based on amino acid data. However, it is not straightforward to interpret the relationships of the four species in the genus. The codon usage analysis placed S. glaseri and S. carpocapsae in a sub-group (Additional file 2: Figure   S1); however, the gene arrangements of S. glaseri were more similar to S. litorale than to S. carpocapsae (Fig. 3). Either of these results are not completely consistent with mitochondrial amino acid phylogeny and 18S phylogeny ( Fig. 2; Additional file 2: Figure S2).
Comparison of two strains of S. carpocapsae of different origins revealed 1.5 % variations in the mitochondrial genomes at the nucleotide level, indicating that the mitochondrial genome can provide good markers for population studies. We observed a location-biased distribution of SNVs and indels, suggesting that marker selection will be important for population studies. Interestingly, the rRNA gene regions have fewer SNVs than the protein-coding region, while indels were more frequently represented in the rRNA gene and were not Fig. 2 Maximum-likelihood tree inferred from 12 mitochondrial proteins. Amino acid sequences were aligned before concatenation, and the phylogenetic analysis was performed using RAxML v7.2.8 with 500 bootstrap resampling replicates. Almost the same topology was obtained from the same amino acid alignment using Bayesian analysis with MrBayes v3.2.2. Numbers above the branches indicate bootstrap values for maximumlikelihood and Bayesian posterior probabilities, respectively. The scale-bar shows the number of amino acid substitutions per site found in protein-coding regions. This is probably because of functional differences between protein-coding and rRNAs. Indels in protein-coding region can be more problematic because they could cause large protein changes by insertions of stop or frame-shift mutations, whereas mismatches by SNVs, rather than indels, could be critical in rRNA genes because they could produce more erratic RNA secondary structures. These are only the results from a comparison of a small number of strains. For a deep understanding of their evolution patterns, it will be interesting to re-sequence and compare more mitochondrial genomes from a larger numbers of strains using high-throughput sequencers. The genome information provided in this study will be useful as references in those kinds of studies.

Conclusions
In this study, the complete mitochondrial genome sequences of four Steinernema species were determined. The complete mitochondrial genomes of S. carpocapsae, S. glaseri, S. kushidai and S. litorale were 13,924, 13,851, 15,182 and 21,403 bp in length, respectively, and all contained the 36 expected genes (12 protein-coding genes, two ribosomal RNA genes and 22 transfer RNA genes), all transcribed in the same direction. Phylogenetic analyses using the amino acid sequences of the 12 protein-coding genes showed that the Steinernema species formed a monophylic clade as a sister clade of Rhabditina and Asacridida. They were also more closely positioned to other Clade 10 nematodes, including B. xylophilus, A. besseyi and P. redivivus, than to Strongyloides species, which belong to the same infraorder as the Steinernema species. The data on codon usage and gene arrangement also supported this result. We also compared the mitochondrial genomes of two strains of S. carpocapsae and detected high intraspecific diversity. This mitochondrial genome information will form a useful resource for evolutionary and population studies of Steinernema nematodes.