Mitochondrial genomes of four slug moths (Lepidoptera, Limacodidae): Genome description and phylogenetic implications

Abstract The family Limacodidae belongs to the superfamily Zygaenoidea, which includes 1672 species commonly referred to as slug moths. Limacodidae larvae are major pests for many economically important plant species and can cause human dermatitis. At present, the structure of the mitochondrial genome (mitogenome), phylogenetic position, and adaptive evolution of slug moths are poorly understood. Herein, the mitogenomes of Parasa lepida, Phlossa conjuncta, Thosea sinensis, and Setora sinensis were sequenced and compared with other available mitogenome sequences to better characterize the mitogenomic diversity and evolution of this moth family. The mitogenomes of P. lepida, P. conjuncta, T. sinensis, and S. sinensis were confirmed to be circular in structure with lengths of 15,575 bp, 15,553 bp, 15,535 bp, and 15,529 bp, respectively. The Limacodidae mitogenomes exhibited similar nucleotide composition, codon usage, RNA structure, and control region patterns, indicating the conservation of the mitogenome in the family Limacodidae. A sliding window, Ka/Ks, and genetic distance analyses revealed that the atp8 and nad6 genes exhibited the highest levels of variability and the most rapid evolutionary rates among the 13 protein‐coding genes (PCGs) encoded in these Limacodidae mitogenomes, suggesting that they may offer value as candidate DNA markers. The phylogenetic analysis recovered the overall relationship as Tortricoidea + (Sesiidae + (Zygaenoidea + (Cossoidea/+Choreutoidea + (others)))). Within Zygaenoidea, Limacodidae was recovered as monophyletic, and the phylogenetic relationships were recovered as (Phaudidae + Zyganidae) + Limacodidae in all six phylogenetic trees. The analysis indicated that P. lepida, P. conjuncta, T. sinensis, and S. sinensis are members of the Limacodidae.


| INTRODUC TI ON
Given their status as eukaryotic organelles that originally evolved from endogenous symbiotic bacteria, mitochondria have been major targets of research interest for over a century (Andersson et al., 1998;Jacobs, 2023;Sagan, 1967).The mitochondrial genome (mitogenome) in insects consists of a closed double-stranded circular DNA loop that offers substantial value to studies of genetic evolution, species origins, population biology, phylogenetic relationships, and taxonomic classification owing to the high mutation rates, rapid evolution, and low recombination rates to which mitogenomic DNA is subjected (Chen et al., 2022;Sankoff et al., 1992;Tao et al., 2023;Wolstenholme, 1992).
The slug moth family Limacodidae (1672 species) belongs to the superfamily Zygaenoidea (van Nieukerken et al., 2011).They are a class of plant pests in subtropical and tropical regions.The taxonomic status of Limacodidae within Lepidoptera has been controversial.Previously, Limacodidae was classified under Tineoidea or Psychoidea.However, Brock (1971) classified it as a member of the superfamily Cossoidea based on the morphological characteristics of the adult thorax and forewings.Other researchers further confirmed and supported this classification (Fletcher & Nye, 1982;Heppner, 1998;Holloway, 1986).
However, Common (1975) placed it in Zygaenoidea based on immaturestage characters.Epstein (1996) classified Limacodidae, Aididae, Dalceridae, Megalopygidae, and Somabrachyidae as monophyletic groups based on the characteristics of adults and larvae, and classified them into Zygaenoidea.Thus, the phylogenetic relationship of Zygaenoidea and its internal phylogeny remains unclear.
As of January 2023, more than 1000 Lepidoptera mitogenomes have been published.However, only 17 are available for species in the superfamily Zygaenoidea.Due to limitations of molecular data, phylogenetic studies of the superfamily Zygaenoidea have long remained at the morphological level.To increase the existing genomic information for these insects, the present work was conducted to sequence and annotate the mitogenomes of four slug moths belonging to the family Limacodidae (Figure 1).The mitogenomes of Phlossa conjuncta and Setora sinensis sequenced herein are the first time for the genera Phlossa and Setora.In contrast, Parasa lepida is the second species sequenced in the genus Parasa, and an additional resequencing of Thosea sinensis was performed for these analyses (Bian et al., 2020).The analysis and comparison of the mitogenomes of different slug moth species belonging to the family Limacodidae established a foundation for research focused on the genetic structure of this family.T. sinensis is one of the dominant species among the plant-damaging slug moth species (Xiao et al., 2020).Therefore, resequencing the T. sinensis mitogenome will provide a valuable reference for future research on the genetic diversity, population structure, and origin of this species.To gain a deeper understanding of the superfamily Zygaenoidea and its internal phylogenetic structure, 144 Apoditrysia mitogenomes from three datasets were utilized to reconstruct six phylogenetic trees.Two species from the superfamily Yponomeutoidea (Lyonetia clerkella and Plutella armoraciae) were used as outgroups, which served as a basis for future research on the phylogenetic relationships among members of the superfamily Zygaenoidea.

| Sample collection and DNA extraction
The larvae samples of P. lepida, P. conjuncta, T. sinensis, and S. sinensis were collected from different provinces in China (Zhejiang, Henan, and Anhui) in 2021 (Table 1).A thorax muscle tissue sample from each fresh specimen was stored in 100% ethanol at −20°C for future use.Morphological characteristics were used to identify these specimens (Epstein, 1996;Wang & Wang, 2009), and voucher specimens were kept in the Department of Medical Parasitology, Wannan Medical College, Anhui Province, China.Based on the provided protocols, the DNeasy Blood & Tissue Kit (Qiagen, Hilden, Germany) was used to extract total DNA from these samples.The extracted DNA was stored at −20°C.

| Mitogenomic annotation and analysis
The annotation of assembled mitogenomic sequences was performed with the MITOS tool (Bernt et al., 2013)  mpimp -golm.mpg.de/ OGDraw.html).In addition, a reannotation of nine mitogenomes was performed for phylogenetic and comparative genomic analyses.Their annotated information is detailed in Table S1.

| Primer design, PCR amplification, and sequencing
To validate the start codon of the cox1 gene of P. lepida, we designed two pairs of primers using Primer Premier v5.0 (Premier Biosoft, Palo Alto, CA, USA): 1L_F: 5′-TCGCTTAATAACTCAGCCATTT-3′, 1L_R: A total of 25 μL of the PCR reaction system was used (12.5 μL of Prime STAR Max Premix (2×) (TaKaRa, Beijing, China), 2 μL of template DNA, 1 μL each of 10 μM forward and reverse primers, and the rest was made up with double distilled water).PCR amplification was performed under the following conditions: 95°C for 5 min, 35 cycles of 95°C for 30 s, 50-52°C for 60 s, and 72°C for 60 s.PCR products were detected by 1% agarose gel electrophoresis and sent to General Biosystems (Anhui) Co., Ltd for Sanger sequencing.Then, the sequencing results were compared to P. lepida mitogenome using DNAMAN version 6.0 (Lynnon Biosoft, USA).
BI trees were calculated inferred with MrBayes 3.2.6 (Ronquist et al., 2012), with 10,000,000 generations and a partition model using two parallel runs and four independent Markov chains, performing sampling every 1000 generations.After reaching an average standard deviation for split frequencies <0.01, the first 25% of the sampled data were discarded, with the remaining trees then being used to represent posterior probability (PP) values.The best-fit partition model (Edge-linked) was selected using ModelFinder v2.2.0 (Kalyaanamoorthy et al., 2017) by a distribution-free rate model using the Bayesian information criterion for ML analyses.ML trees were calculated with IQ-TREE v1.6.8 (Nguyen et al., 2015).Using an edge-linked partition model, 5000 ultrafast bootstrap replicates were used to assess bootstrap support (BS) (Minh et al., 2013).The best-fit partition model details for each tree are provided in Table S3.

| Mitogenomic organization and base composition
After sequencing was completed, the mitogenomes of P. lepida, P. conjuncta, T. sinensis, and S. sinensis were confirmed to be circular in structure with respective lengths of 15,575 bp, 15,553 bp, 15,535 bp, and 15,529 bp (Figure 2).These results were largely consistent with those for other published Limacodidae mitogenomes, which are 15.2-15.7 kb (Table S4).These four mitogenomes encode 37 total genes, including 22 tRNAs, 13 PCGs, and two rRNAs (rrnL and rrnS).
Of these, 14 tRNAs and 9 PCGs were found to be encoded on the majority strand (J-strand), while the remaining genes were encoded on the minority strand (N-strand) (Table 2).Each of these four mitog-  Negative AT skew and GC skew were observed for the full mitogenomes of Limacodidae species other than Monema flavescens (AT skew = 0.013; Figure S1).

| Analyses of protein-coding genes and codon usage
The total respective sizes of the 13 PCGs identified in the P.  S4).While the majority of these PCGs utilized ATN (ATT, ATG, and ATA) start codons, the cox1 gene in the P. conjuncta, T. sinensis, and S. sinensis mitogenomes was found to utilize a CGA start codon.In contrast, ATT was used in P. lepida (Figure 3).The results of the two primer pairs sequenced by Sanger were in complete agreement with the result of high-throughput sequencing (Figure S2).

| Analyses of tRNAs and rRNAs
These four mitogenome sequences were all found to contain 22 tRNA genes, of which 14 and 8 were encoded on the J-strand and N-strand, respectively.The concatenated tRNA genes for these 11 Limacodidae genomes exhibited a positive AT skew and GC skew (Table S4).Except for trnS1 (AGN), most tRNAs exhibited typical cloverleaf secondary structures (Figure S5).In addition to appropriate base pairing, some mismatched base pairs (U-U, A-C, G-U) were also observed in the secondary structure of these tRNAs.In these four mitogenomes, the two identified rRNA genes were separated by trnV, and were located between trnL1 and the control region (CR).The length of rrnL ranged from 1330 bp (P.lepida) to 1434 bp (T.sinensis), while the length of rrnS ranged from 778 bp (P.conjuncta) to 785 bp (P.lepida) (Table 2).

| Putative control region analyses
The predicted length of the control regions (CRs) in these four mitogenomes, which were positioned between rrnS and trnM, ranged from 393 bp (S. sinensis) to 425 bp (P.conjuncta) with the A + T content ranging from 91.3% (S. sinensis) to 94.8% (P.conjuncta), which is in line with similar findings for other Limacodidae species (Table S4).
In general, four conserved CR structures were present in the mitogenomes, including (Figure S6): (1) an ATAGA motif adjacent to rrnS followed by a poly-T sequence that serves as the origin of minority or light-strand replication; (2) microsatellite-like repeat regions (AT)n or (TA)n beginning with an ATTTA motif that is widely used for analyses of genetic diversity, kinship, and species origins; (3) a poly-A region upstream of trnM; and (4) tandem repeats throughout the CR (Table S6).

| Phylogenetic relationships
We discussed the position of the Zygaenoidea in the clade Apoditrysia and the phylogenetic relationships within it based on the 20 superfamilies, 144 mitogenomes from 40 families, and the two Yponomeutoidea outgroup species available on NCBI.Both analysis methods (BI and ML) based on three datasets (PCGs_rRNA, PCGs_rRNA_codPosOneTwo, and PCGs_AA) produced six phylograms that showed similar topological features with high statistical support (Figure 5; Figure S7).In five trees, Sesiidae within Cossoidea formed the base of Zygaenoidea (Figure 5; Figure S7a,c-e), whereas in the BI analysis using the PCGs_rRNA dataset, Sesiidae was positioned as the base of Tortricoidea (Figure S7b).Urodoidea and Copromorphoidea were sister branches in four trees (Figure 5; Figure S7a-c), while in the other two trees, Urodoidea was placed at the base of Tortricoidea (Figure S7d,e).Here, we had chosen to show a tree with more of the same topology and high support for the Zygaenoidea (Figure 5).The overall relationship is roughly as follows: Tortricoidea + (Sesiidae + (Zygaenoidea + (Coss oidea /+Choreutidae + (others)))).The monophyly of Zygaenoidea (PP = 1, BS = 100) was strongly supported in all six phylogenetic trees (Figure 5; Figure S7).In Zygaenoidea, the family-level relationships were determined as (Phaudidae + Zyganidae) + Limacodi dae.The family Limacodidae (PP = 1, BS = 100) was recovered as monophyletic in all six phylogenetic trees.These results provide strong support for the attribution of P. lepida, P. conjuncta, T. sinensis, and S. sinensis to the Limacodidae, as well as the phylogenetic relationships within the Zygaenoidea.

| DISCUSS ION
In this study, four mitogenomes of P. lepida, P. conjuncta, T. sinensis, and S. sinensis were sequenced, annotated, and compared with other Limacodidae mitogenomes.By comparison, the organization and base composition of the Limacodidae mitogenomes were highly consistent, indicating the conservation of the mitogenome within the family Limacodidae.For PCG analysis, the cox1 CGA start codon is a common finding in many of the published lepidopteran mitogenomes (Bian et al., 2020;Cameron, 2014;Djoumad et al., 2017;Jiang et al., 2022;Liu et al., 2016;Ma et al., 2016;Zhu et al., 2018).
However, the ATT start codon was predicted to be utilized by the cox1 gene of P. lepida in this study.To validate this, we designed two pairs of primers for Sanger sequencing and found that the results matched exactly with the high-throughput sequencing results (Figure S7).In addition, no possible CGA codon or other ATN start codons were found near the putative start codon ATT.Thus, it is reasonable to assume that the start codon of the cox1 gene of P. lepida is ATT.This is the first time this codon has been identified in the Limacodidae family (Figure 3).
Nucleotide diversity is a fundamental measure employed to evaluate genetic diversity, and it is crucial in studies centered around genetic diversity (Clark et al., 2007;Olson et al., 2010).
A sliding window was used to analyze the nucleotide diversity (π) of 13 PCGs in the Limacodidae species (Figure 4a).Of these 13 PCGs, atp8 (π = 0.197) and nad6 (π = 0.19) exhibited the greatest degree of variability.The ω value is used to evaluate whether a specific gene is influenced by purifying selection (0 < ω < 1), neutral evolution (ω = 1), or positive selection (ω > 1) (Meiklejohn et al., 2007).The Ka/Ks ratio (ω) for each PCG identified in these mitogenomes was found to be less than 1, indicating that all of these genes are under purifying selection.Of these 13 PCGs, cox1 (ω = 0.08) exhibited the strongest purifying selection and the lowest evolutionary rate, whereas atp8 (ω = 0.522) and nad6 (ω = 0.419) were subjected to relatively weak purifying selection and thus had relatively rapid evolutionary rates.Similar results were also obtained through pairwise calculations of genetic distances, with higher genetic distance values for atp8 (0.251) and nad6 (0.224) genes, indicating the comparatively rapid evolution of these genes.This is in stark contrast to the lower values observed for cox1 (0.124) and cox2 (0.136), which are consistent with slower relative evolution (Figure 4b).Thus, the atp8 and nad6 genes were also identified as promising candidate DNA markers that may enable species delimitation by clarifying phylogenetic relationships among slug moth species.According to the latest phylogenetic hypothesis of Lepidoptera, Yponomeutoidea is a sister to Apoditrysia (Kawahara et al., 2019), so two Yponomeutoidea species (L.clerkella and P. armoraciae) were chosen as the outgroups in this study.Long-branched attraction (LBA) is a major obstacle to phylogenetic reconstruction, and excluding the third codon reduces the effect of LBA (Qu et al., 2017;Sanderson et al., 2000).Therefore, in this study, two different nucleotide datasets (PCGs_rRNA and PCGs_rRNA_codPosOneTwo) and PCGs_AA datasets for phylogenetic analysis were used.The topological features of the six phylogenetic trees obtained exhibited minor variations.The majority of topological discrepancies were observed on branches with low support, possibly due to variations in the datasets used and differences between the BI and ML phylogenetic methods (Liu et al., 2021).In the future, expanding genetic data may assist in resolving this problem.
Cossoidea and Zygaenoidea were sister groups supported by some previous morphological and molecular evidence (Bao et al., 2019;Bian et al., 2020;Minet, 1991;Scott, 1986).Phylogenetic analysis of 13 PCGs based on mitogenomes showed a close phylogenetic relationship among Zygaenoidea, Papilionoidea, and Tortricidea by Liu et al. (2016).Some studies consider Phaudidae as a subfamily of Zygaenidae (Fänger, 1999;Minet, 1991).Niehuis et al. (2006)  ).In addition, Zygaenidae was shown to be the sister group of this big branch.Zhang, Gao, et al. (2020) reported that Phaudidae is a monophyletic family sister to Zygaenidae.However, Phauda flamman of the Phaudidae belongs to the Zygaenidae, as determined in our study of six trees (Figure 5; Figure S7).Phaudidae may belong to a subfamily of Zygaenidae.Increasing the dataset may help to solve such problems.Within the family, Limacodidae, Monema, Parasa, and Latoia formed a branch, as did Setora, Phlossa, Iragoides, Thosea, and Quasithosea.These two sister branches and the sister clade Narosa and Apoda exhibited the strongestsupport (PP = 1, BS = 100) among the six trees constructed in this study (Figure 5; Figure S7).However, Parasa consocia was found to be more closely related to Latoia hilarata than to P. lepida, while T. sinensis Jiangsu (NC_061059.1) was more closely related to Quasithosea sythoffi than to T. sinensis Henan (OP132388) (Figure 5; Figure S7b-e).However, given the sample limitations for these analyses, the phylogenetic connections established in these studies cannot be considered definitive.
While mitogenomic sequences are commonly used when inferring phylogenetic relationships among insect species, the sole reliance on these data is nonetheless subject to certain limitations.Host mitochondrial DNA variations at the intra-specific or inter-specific levels can be influenced by Wolbachia bacteria, leading to divergence during mitogenome-based phylogenetic analysis observed in butterflies (Jiang et al., 2018;Kodandaramaiah et al., 2013;Sucháčková Bartoňová et al., 2021;Whitworth et al., 2007).At the intra-specific level, Wolbachia generates selective sweep through selective pressure, which can cause the rapid expansion of host mitochondrial haplotypes associated with Wolbachia through the hitch-hike effect, thereby increasing the proportion of these mitochondrial haplotypes in the host population and consequently reducing the mitochondrial DNA diversity in the host population (Deng et al., 2021;Jiang et al., 2018;Morrow & Riegler, 2021).
Wolbachia may lead to mitochondrial introgression between closely related species at inter-specific levels (Gompert et al., 2008;Jackel et al., 2013).Large integrated datasets are needed in the future to improve phylogenetic resolution among lepidopteran taxa, such as greater integrated datasets of mitogenomes, nuclear genes, and morphological characters.

F
The cox1 start codons from Limacodidae species.Genes located on the N-strand are shown in light gray, and the cox1 start codons are circled in red.F I G U R E 4 Nucleotide diversity (π) and selection pressures for 13 protein-coding genes (PCGs) in 11 Limacodidae species.(a) Thirteen PCGs were aligned and subjected to a sliding window analysis.Average π values are indicated above the corresponding genes.(b) Nonsynonymous (Ka) to synonymous (Ks) substitution ratios and genetic distance analyses using a Kimura 2-parameter model for these 13 PCGs from 11 Limacodidae species.
conducted a study that found Phaudidae to be closely related to Lacturidae and the sister group of Limacodidae + (Somabrachyidae + Himantopteridae F I G U R E 5 Phylogenetic relationships of 20 superfamilies of 144 Apoditrysia mitogenomes inferred from Bayesian inference (BI) using the PCGs_AA dataset.Lyonetia clerkella (NC_037944.1) and Plutella armoraciae (NC_064059.1) of the superfamily Yponomeutoidea were used as outgroups.*The sequences for species with names in bold were generated in this study.

TA B L E 1 Sampling details for four specimens used in this study.
(Katoh & Standley, 2013)no acid sequences for the 13 PCGs and the nucleotide sequences for the RNAs were aligned with MAFFT v7.313(Katoh & Standley, 2013)(Code (Capella- Gutierrez et al., 2009)al).Nucleotide sequence alignment results for these 13 PCGs were optimized with MACSE v2.03(Ranwez et al., 2018), and ambiguously aligned fragments from these 13 alignments were removed in batches with Gblocks 0.91b(Talavera & Castresana, 2007)in Phylosuite v1.2.3.Any gaps in RNA and AA sequences were removed with TrimAl v1.2rev57(Capella- Gutierrez et al., 2009), and the best-fit partition model was selected using a greedy algorithm and the Akaike information criterion (AIC)