The First Mitogenomes of the Subfamily Odontiinae (Lepidoptera, Crambidae) and Phylogenetic Analysis of Pyraloidea

Simple Summary The Odontiinae is a small group in the Pyraloidea comprised of 388 species in 88 genera, but externally, these moths are diverse, including heterogeneous maculation and a size range from 9 to 50 mm in total wingspan. The monophyly of Pyraloidea and the two families (Pyralidae and Crambidae) is well supported by phylogenetic analyses based on morphology and molecular data of multiple nuclear genes. However, only a few mito-phylogenetic analyses have been conducted and no mitogenome of Odontiinae species has been reported. Three complete mitogenomes of odontiine species were sequenced and analyzed for the first time herein. The results showed that Odontiinae mitogenomes shared similar genomic characters with other Pyraloidea. The phylogenetic analyses based on 13 PCGs of mitogenomes confirmed the monophyly of Odontiinae and its position within Crambidae. Abstract The complete mitochondrial genomes of three species of Odontiinae were newly sequenced: Dausara latiterminalis Yoshiyasu, Heortia vitessoides (Moore), and Pseudonoorda nigropunctalis (Hampson). These circular and double-stranded mitogenomes vary from 15,084 bp to 15,237 bp in size, including 13 protein-coding genes (PCGs), two ribosomal RNA genes (rRNAs), and 22 transfer RNA genes (tRNAs) and an A + T-rich region. The nucleotide composition indicated a strong A/T bias. Most PCGs are initiated with an ATN codon and terminated by a codon of TAR. All tRNAs could be folded into the clover-leaf structure with the exception of trnS1 (AGN), in which the dihydrouridine (DHU) arm formed a simple loop, and the motif ‘ATAG’ and ‘ATTTA’ in the A + T-rich region was also founded. The phylogenomic analyses covering Odontiinae + 11 subfamilies of Pyraloidea were conducted. Similar topologies were generated from both Bayesian inference (BI) and maximum likelihood (ML) analyses based on the nucleotide and amino acid sequence data. There was some discrepancy in the sister-group relationship of Odontiinae and Glaphyriinae, and the relationships among the subfamilies in the ‘CAMMSS clade’ of the Crambidae. The results of this study suggest that mitogenomic data are useful for resolving the deep-level relationships of Pyraloidea and the topologies generated from amino acid data might be more realistic and reliable. Moreover, more mitogenomic taxon sampling and larger scale analyses with more genes or a combination of mitogenomic and nuclear genes are needed to reconstruct a comprehensive framework of the pyraloid phylogeny.


Introduction
Within Lepidoptera, the Pyraloidea is one of the largest superfamilies that includes two families-Pyralidae and Crambidae. This superfamily, found on all continents except Antarctica, comprises about 16,000 named species worldwide [1,2]. As one of the small subfamilies of Pyraloidea, Odontiinae is comprised of 388 species in 88 genera, and it is numerous in eremic habitats of all major biogeographic regions except New Zealand. Many odontiine species are important pests, and their larvae are leaf miners and folders, flower and bud feeders, and fruit borers [3][4][5]. They usually cause serious damage to medicinal plants and fruit trees such as the incense trees (Aquilaria sinensis) and olive trees (Canarium album) in tropical and subtropical regions [6,7]. Historically, the relationships among Odontiinae, Glaphyriinae, Noordinae, Evergestinae and Cathariinae have been controversial for a long time until the landmark studies conducted by Regier et al. [8] and Léger et al. [4] which discussed the phylogeny of Pyraloidea including Odontiinae based on several nuclear genes combined with a single mitochondrial gene COI. However, a mitogenome-based investigation about the relationships among Odontiinae and the other subfamilies of Pyraloidea has never been discussed.
The typical arthropod mitochondrial genome (mitogenome) is a circular, doublestranded molecule which encodes 37 genes, including 13 protein-coding genes (PCGs), two ribosomal RNA genes (rRNAs), 22 transfer RNA genes (tRNAs), and an A + T-rich region [9,10]. Due to cellular abundance, an absence of introns, rapid evolutionary rate, and a lack of extensive recombination, mitogenome sequences can be easily amplified and have been proven to be a useful source that has been extensively employed in systematics, population genetics and evolutionary biology in the past decade [11][12][13][14][15]. In recent years, mitochondrial genomes from different subfamilies of Pyraloidea have been obtained, and several mitogenome-based investigations with respect to the phylogeny of this superfamily have been performed [16][17][18][19][20][21]. Unfortunately, a mitochondrial genome of the subfamily Odontiinae has remained unknown. Moreover, 50+ mitogenomes of Pyraloidea are a tiny part comparing with 16,000 named species of the superfamily. Therefore, more mitogenomes of pyraloids are needed to be sequenced and included in phylogenetic analyses.
To begin to rectify some of the above issues, three complete mitogenomes of odontiine species, including Dausara latiterminalis Yoshiyasu, Heortia vitessoides (Moore), and Pseudonoorda nigropunctalis (Hampson) were sequenced and annotated for the first time. In addition, a comparative analysis was conducted to reveal the genomic organization, nucleotide composition, codon usage, and tRNA secondary structure. Moreover, we used the newly sequenced mitochondrial genomes of this study and mitogenomes available online to reconstruct a phylogenetic tree of Pyraloidea to better understand their evolutionary history and relationships within Pyraloidea.

Specimen Collection and Genomic DNA Extraction
Three odontiine species were collected by light traps in China, Malaysia and Laos (Supplementary Table S1), respectively, and were preserved in 99.5% ethanol during fieldwork. They were then stored at −50 centigrade degree environment in the Insect Collection of Nankai University (NKU), Tianjin, China.
All the above specimens were identified based on the morphological characters. Thorax and legs were used to extract genomic DNA with the DNeasy tissue kit (Qiagen, Hilden, Germany) according to the manufacturer's protocol. DNA concentrations were measured with a DeNovix DS-11 Spectrophotometer, DNA integrity was examined with agarose gel electrophoresis by 0.5 × TBE (Tris base, Boric acid and EDTA) buffer with 3 Volt/centromere for 45 min.

High throughout Sequencing
All three odontiine genomic DNA were qualified for high throughput sequencing (HTS). The genomic DNA was fragmented to 350-500 bp by Covaris S220 Focused Ultrasonicator (Covaris, MA, USA). The sequence libraries were constructed using TruSeq DNA LT Sample Preparation Kit (Illumina, Inc., San Diego, CA, USA). After repairing the blunt ends, adenylating 3 ends and ligating adapters, the fragmented DNA were enriched.
Those libraries were pooled and sequenced by an Illumina Hiseq X10 platform. The raw data were filtered as follows: (1) the adapters were removed; (2) the reads that contained more than five Ns were removed; (3) 4 base slide windows were performed, and the reads or the averaged Phred value of less than Q20 were removed; (4) after the above steps, the reads shorter than 75 bp or the Phred values less than Q15 were removed. Finally, approximately 40 Gb clean data of paired-end reads of 150 bp length were generated.

Data Assemble and Annotation
The raw data were assembled by MitoZ v2.4 [22] with the 'all' option, which finalized data assembly, de novo assembly, genome annotation, and visualization only within one step. The assembled circular mitogenomes were reordered trnM as a start gene with the script 'Mitogenome_reorder.py' in MitoZ. All calculation was performed in the highperformance computing platform in Capital Normal University.
As for the sequence length of the control region, less than 600 bp is not annotated by mitoZ, and the known length of the A + T-rich region in most lepidopteran species ranges from 280 to 500 bp [20,21,[23][24][25]. The annotation of the three mitochondrial genomes was also performed by the MITOS2 online server with default parameters (http://mitos2.bioinf. uni-leipzig.de/index.py, accessed on 30 March 2021). The A + T-rich region was ensured according to the results of MITOS2 and the other available mitochondrial genomes of Crambidae [20].

Statistics of the Odontiine Mitochondrial Genomes
Nucleotide composition and relative synonymous codon usage (RSCU) of the odontiine mitogenomes were calculated in MEGA 5 [26]. Nucleotide compositional skew was calculated according to the formulas: [27].

Phylogenetic Analyses
To investigate the phylogenetic implications of the mitogenomes of the three species in Pyraloidea, we reconstructed the subfamily-level relationships within Pyraloidea using three different datasets of the 13 protein coding genes (PCG) and two inference methods.
The mitogenomic phylogeny of Pyraloidea was reconstructed with 40 ingroups (37 online data and 3 newly produced data in this study) and 5 outgroups ( Table 1). The three datasets are PCG123 (13 PCGs including all codon positions), PCG12 (13 PCGs without third codon positions) and AA (amino acid of 13 PCGs). Bayesian inference (BI) and maximum likelihood (ML) methods were used to reconstruct phylogenetic trees. For PCG123 and PCG12 datasets, the best DNA model based on Akaike information criterion (AIC) was performed using jModeltest 2.1.7 [52] (Table S2), and those selected models were used by BI with software MrBayes 3.2.6 [53]. To ensure that the average standard deviation of split frequencies was less than 0.01, eight million generations were run with sampling every 1000 generations. Node support was assessed by posterior probabilities (PPs). The ML analyses were performed using IQ-TREE v2.1.2 [54], selecting the best model and constructing phylogenetic trees automatically using 1000 non-parametric bootstrap replicates (BS) and SH-aLRT test with an unpartitioned strategy ('-m MFP -b 1000 alrt 1000 ), whilst other settings were default.
To mitigate the possible effects of base-composition bias and among-site rate heterogeneity, the AA dataset was analyzed using Phylobayes v4.1c [55], and the MtZoa model, which fit the mitogenomic phylogenetic analyses in metazoan groups, was chosen to perform the BI analyses [56]. Two independent MCMC chains in Phylobayes were run until convergence (maxdiff < 0.1 and minimum effective size > 50). As for the ML analysis, the dataset AA was performed by IQ-TREE v2.1.2 with the parameters '-m MtZoa + F + I + G4 -b 1000 -alrt 1000 . Tracer v1.6 [57] was used to check the likelihoods of all parameters of BI analyses of the three datasets to ensure the effective sample size (ESS) values greater than 200. The consensus tree was calculated by discarding the first 25% trees. To verify the consistencies of the topologies, both BI and ML analyses were repeated three times, and the phylogenetic trees were visualized by Figtree v.1.4.3 [58].

Mitogenome Organization and Nucleotide Composition of Odontiinae
The annotations for the mitogenomes of the three odontiine species are shown in Table 2, and the circular maps of the mitogenomes of the three species are shown in Figure  1. The complete mitogenomes of the three species were investigated here, and all were found to be composed of circular double stranded molecules with mildly varying sizes. Each mitogenome contains the typical set of 37 genes, including 13 typical protein-coding genes (PCGs), 22 transfer RNA genes (tRNAs), two ribosomal RNA genes (rRNAs) and an A + T rich area. The majority strand (J-strand) encoded 23 genes (9 PCGs, 14 tRNAs), while the remaining genes were located on the minority strand (N-strand) (four PCGs, eight tRNAs and two rRNAs) ( Figure 1, Table 2). The sizes were as follows: Dausara latiterminalis 15,147 bp, Heortia vitessoides 15,237 bp, and Pseudonoorda nigropunctalis 15,084 bp. Gene orders of these three Odontiinae species were identical to those of other species reported in the Crambidae [35,37], specifically in reference to the tRNA gene cluster trnI-trnQ-trnM that was rearranged to trnM-trnI-trnQ.
In the whole mitogenomes of the three Odontiinae species, the nucleotide composition is as shown in Table S3. It indicated a strong A and T bias, and the A + T% content ranged from 80.6% (in Dausara latiterminalis and Heortia vitessoides) to 81.0% (in Pseudonoorda nigropunctalis). Comparing the AT content of the whole mitogenome, control region, PCGs, tRNAs, and rRNAs, the A + T-rich region was the highest while the PCGs was the lowest for all the three species of Odontiinae (Table S3). The AT skew of the whole mitogenome within the three odontiine species ranged from −0.012 to −0.003, and GC skew ranged from −0.201 to −0.172, which was consistent with that of several previously reported pyraloid species [35,49,59,60].

Protein-Coding Genes of Odontiinae
The PCGs within species of Odontiinae ranged from 162 bp (ATP8) to 1752 bp (ND5) in size, and the total PCGs lengths ranged from 11,268 bp to 11,313 bp. The three mitogenomes of Odontiinae exhibited similar start and stop codons as follows ( Table 2): all the initiation codons of PCGs were ATN, and ATT was the most frequently used start codon for ND2, ATP8, ND3, and ND5; ATG was the most frequently used for COII, ATP6, COIII, ND4, ND4L, Cob and ND1, except for the COI which started with CGA. Except for COII which terminated with a single T residue, twelve PCGs were terminated by the standard stop codon TAR, and TAA was the most frequently used codon. Truncated termination codons are commonly used in metazoan mitogenomes and are modified by the post-transcriptional poly-adenylation to a complete TAA stop codon [61]. The RSCU values of the three odontiine species are shown in Figure 2. The codons UUA-Leu2, UCU-Ser2, GGA-Gly and UCA-Ser2 were the most commonly used in the Odontiinae mitogenome, while AGC-Ser1 was not present in these PCGs.

RNA Genes of Odontiinae
The 22 transport RNA (tRNA) genes of the three odontiine species were discovered (Figure 3), and the entire lengths of the three mitogenomes was 1452 bp in Pseudonoorda nigropunctalis, 1463 bp in Dausara latiterminalis, and 1471 bp in Heortia vitessoides. The size of the 22 tRNA genes ranged from 63 to 71 bp ( Table 2). Among them, 21 tRNAs could be folded into the typical clover-leaf structure, except for trnS (AGN), which lost a dihydrouridine (DHU) arm (Figure 3). This phenomenon is common in Pyraloidea and in other insect mitogenomes [50,[62][63][64][65]. The secondary structures comprised of the anticodon loop (7 nt), anticodon stem (5 bp), and the acceptor stem (7 bp) were conserved in length, while the length of DHU (3-4 bp) and TψC (4-5 bp) stems was variable, except for trnS1. Additionally, the identified unmatched base pairs in different stems of tRNAs are shown in Figure 3, and these mismatched nucleotides might be restored during the post-transcriptional editing processes [66]. As for the ribosomal RNA (rRNA) of the three species of Odontiinae, both of l-rRNA (rrnL) and s-rRNA (rrnS) genes were encoded on the N-strand, and the rrnL ranged from 1364 bp (in Heortia vitessoides) to 1369 bp (in Dausara latiterminalis and Pseudonoorda nigropunctalis) in length, while the rrnS ranged from 780 bp (in Heortia vitessoides) to 784 bp (in Pseudonoorda nigropunctalis). Both of the ribosomal RNAs showed a positive AT skew and negative GC skew, and the A + T% was about 84% in rrnL while about 86% in rrnS among the three odontiine species.

A + T-Rich Region of Odontiinae
In mitogenome, the largest non-coding region is normally the A + T-rich region (also called the control region). The A + T-rich region of odontiine mitogenomes are located between the rrnS and trnM genes, and the length was 280 bp in Pseudonoorda nigropunctalis, 335 bp in Dausara latiterminalis, and 347 bp in Heortia vitessoides. These lengths were similar to those in other mitogenomes of Pyraloidea (about 340 bp on average) [18]. The A + T% content was 95.2% in Dausara latiterminalis, 96.8% in Heortia vitessoides, and 96.4% in Pseudonoorda nigropunctalis. The comparison of the control regions of the representatives of 10 subfamilies of Pyraloidea is shown in Figure 4. Although the length of the control regions varied in different subfamilies, several conserved elements were discovered: the motif ATAG and the following large poly-T stretch, the motif ATTTA and the following microsatellite structures (AT)n, and the poly-A stretch. These conserved blocks were considered to play a key role in controlling the replication and transcription of the mitogenome [67].

Phylogenetic Analysis
The phylogenetic trees of 40 pyraloid mitogenomes were fully resolved with identical topology, except for the 'non-PS clade' (Regier et al. [8]) that includes Odontiinae based on PCG123, PCG12 and AA for both BI and ML analyses (Figures 5-7; Tables S4 and S5). In all of the phylogenetic trees obtained, we think that the topology of the AA tree is optimal, which is congruent with the previous results of nuclear genes. In all trees based on three datasets, the monophyly of the two families, Pyralidae and Crambidae, was well supported, as has been indicated with the morphology-based and the molecular-based results [1,4,8,68].   Within the Pyralidae, all analyses consistently supported its monophyly of three of the subfamilies, i.e., Pyralinae, Galleriinae and Phycitinae, and the inclusion of the subfamilies, i.e., Pyralinae, Galleriinae, Epipaschiinae and Phycitinae, with the exception of the Chrysauginae, that was not included for a lack of online mitogenomic data. The relationship of the four subfamilies was recovered as Galleriinae + (Phycitinae + (Pyralinae + Epipaschiinae)), which is compatible with previous studies based on mitogenomic data or multiple gene markers [4,8,20]. The subfamily Galleriinae was placed as the sister group to the other three subfamilies with very strong support (PP = 1, BS = 100) (Figures 5-7).
As for the Crambidae, our results show a basal split into two sister lineages, one consisting of Pyraustinae and Spilomelinae, and another comprising the remaining six subfamilies of Crambidae included in this study (Acentropinae, Crambinae, Glaphyriinae, Odontiinae, Schoenobiinae and Scopariinae). The two lineages generally correspond to the 'PS clade' and 'non-PS clade' as defined by Regier et al. [8] based on nuclear genes. Within the 'PS clade', our analyses confirm the sister relationship of Pyraustinae and Spilomelinae with very strong support (PP = 1, BS = 100), which is compatible with previous studies [4,20,59,69].
In this study, the 'non-PS clade' included six subfamilies, and a subset of the currently recognized subfamilies in this taxon: Acentropinae, Crambinae, Glaphyriinae, Odontiinae, Schoenobiinae, Scopariinae [4,8]. In our study, the 'non-PS clade' was clustered into two branches, which corresponds to the 'OG clade' and 'CAMMSS clade' according to Regier et al. [8], respectively. In the 'OG clade', all the results based on PCG123, PCG12 and AA for BI and ML analyses supported the monophyly of the Odontiinae, however, the Glaphyriinae was recovered as paraphyletic based on PCG12 and PCG123, and as monophyly with moderate (PP = 0.5) or high support values (BS = 87) based on AA, which was consistent with the previous studies [4,8].
Within the 'non-PS clade', our results based on three datasets exhibited different topologies, and the inconsistence mainly focused on the 'CAMMSS clade'. The four subfamilies in this study were recovered as Acentropinae + (Scopariinae + (Crambinae + Schoenobiinae)) based on the PCG123 in both BI and ML analyses ( Figure 5). This confirmed the mitogenome-based results of Zhu et al. [18], which were recovered as Acentropinae + (Crambinae + Schoenobiinae) based on the PCG123 and PCG12 in BI and ML analyses, despite the fact that Scopariinae was not involved. Moreover, in our result, the relationships Scopariinae + (Schoenobiinae + (Acentropinae + Crambinae)) and Schoenobiinae + (Scopariinae + (Acentropinae + Crambinae)) were recovered in BI and ML analyses based on PCG12 with low-to-moderate support, respectively ( Figure 6).
As Cameron [10] mentioned, the inclusion of third codon positions may result in serious artifacts due to the faster rate of evolution, and it seems that the results of the inconsistent relationships within the 'CAMMSS clade' based on PCG123 and PCG12 are an example of this phenomenon. Moreover, the long branch of Schoenobiinae was very distinctive in all topologies, and it might be another reason for the discrepant topologies in different analyses. As only one species of this subfamily was sequenced, its position within the 'CAMMSS clade' needs extra samplings to confirm. In addition, in the ML analysis of AA (Figure 7), the relationship (Scopariinae + Crambinae) + (Acentropinae + Schoenobiinae) was recovered. It was identical to the findings of Regier et al. [8] and Léger et al. [4], which was based on several nuclear genes combined with a single mitochondrial gene COI. The BI analysis of AA (Figure 7) showed an identical topology for the 'CAMMSS clade' with the ML analysis, but the 'OG clade' was formed as polyphyletic.
On the basis of the above analyses, a close relationship between Scopariinae and Crambinae and between Schoenobiinae and Acentropinae for the 'CAMMSS clade' and the monophyly of Odontiinae and Glaphyriinae might be more realistic. Perhaps, more mitogenomic data of pyraloid taxa, a combination of mitogenomic and nuclear genes or even the genome-scaled analysis would help to confirm or understand the phylogenetic relationships within Pyraloidea.

Conclusions
In this study, we reported the mitogenome sequences of Dausara latiterminalis, Heortia vitessoides, and Pseudonoorda nigropunctalis, which are the first complete mitogenomes in the subfamily Odontiinae. Compared to other previously reported mitogenomes of Pyraloidea, the newly sequenced odontiine complete mitogenomes are conserved in gene organization, base composition, codon usage of PCGs and secondary structures of tRNAs. The phylogenetic analyses inferred from mitogenomes (PCG123, PCG12 and AA) produced a well-resolved framework for the relationships of Pyraloidea, and the monophyly and position of Odontiinae were also confirmed. Our results were largely congruent with the previous results based on nuclear and mitogenomic genes, except for the relationship within the 'CAMMSS clade'. As for the 'non-PS clade', the ML analysis of the AA dataset was inconsistent with the previous studies based on mitogenomic genes but was identical with the former results based on multiple nuclear genes. The results of this study provided a new or alternative insight to improve the current phylogeny of the Pyraloidea. Moreover, it suggested that mitogenome data were useful for resolving the phylogenetic issues of Pyraloidea and the topologies generated from amino acids might be more realistic and reliable. In addition, extensive samples of multiple taxa and larger scale analyses with more genes, or even a whole genome may be helpful to reconstruct a comprehensive framework of the pyraloid phylogeny.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/insects12060486/s1. Table S1. Collecting information of specimens in present study. Table S2. The best fit DNA model of each locus in datasets PCG123 and PCG12 selected by jModeltest 2.1.7. Table S3. Nucleotide composition of mitochondrial genomes of three odontiine species. Table S4. Phylogenetic relationships among pyraloid species by BI analysis based on the PCG12 dataset. Table S5. Phylogenetic relationships among pyraloid species by BI analysis based on the AA dataset.