The First Eight Mitogenomes of Leaf-Mining Dactylispa Beetles (Coleoptera: Chrysomelidae: Cassidinae) Shed New Light on Subgenus Relationships

Simple Summary With 387 species, Dactylispa is a large genus of family Chrysomelidae, subfamily Cassidinae, and the most species-rich genus of leaf-mining hispine beetles. Among leaf-mining hispines, Dactylispa also feeds on the largest number of host plants, including 29 families and 80 genera, with the main hosts belonging to Poaceae, Rosaceae, and Fagaceae. Some species of Dactylispa are economic crop pests. However, the subgenus classification of Dactylispa species, which is currently based on morphology only, is problematic. Molecular phylogenetic studies may help solve this problem. In this study, we aimed to report the first eight mitochondrial genomes of Dactylispa and construct the first molecular phylogenetic trees of Dactylispa. The evolutionary relationships among three subgenera of Dactylispa were partially resolved by the molecular trees. However, more species should be included to solve the evolutionary issues for this genus of leaf-mining beetles. Abstract The taxonomic classification of Dactylispa, a large genus of leaf-mining beetles, is problematic because it is currently based on morphology alone. Here, the first eight mitochondrial genomes of Dactylispa species, which were used to construct the first molecular phylogenies of this genus, are reported. The lengths of the eight mitogenomes range from 17,189 bp to 20,363 bp. All of the mitochondrial genomes include 13 protein-coding genes (PCGs), 22 transfer RNA genes (tRNAs), 2 ribosomal RNA genes (rRNAs), and 1 A + T-rich region. According to the nonsynonymous/synonymous mutation ratio (Ka/Ks) of all PCGs, the highest and the lowest evolutionary rates were found for atp8 and cox1, respectively, which is a common phenomenon among animals. According to relative synonymous codon usage, UUA(L) has the highest frequency. With two Gonophorini species as the outgroup, mitogenome-based phylogenetic trees of the eight Dactylispa species were constructed using maximum likelihood (ML) and Bayesian inference (BI) methods based on the PCGs, tRNAs, and rRNAs. Two DNA-based phylogenomic inferences and one protein-based phylogenomic inference support the delimitation of the subgenera Dactylispa s. str. and Platypriella as proposed in the system of Chen et al. (1986). However, the subgenus Triplispa is not recovered as monophyletic. The placement of Triplispa species requires further verification and testing with more species. We also found that both adult body shape and host plant relationship might explain the subgeneric relationships among Dactylispa beetles to a certain degree.


Introduction
With 387 species, Dactylispa Weise is a large genus of leaf-mining hispine beetles in Chrysomelidae [1][2][3][4][5]. Dactylispa also feeds on more host plant species than any other Generally, there are only a few publications on Dactylispa, especially in recent years. For example, in a search for "allintitle: Dactylispa" on Google Scholar, there are only 33 hits, with only one article since 2015. Most recent publications focus on morphological descriptions, biological notes, and economic importance [6,7,[20][21][22][23][24], without enough attention devoted to phylogenetic relationships. There are no studies on the mitogenomes and molecular phylogeny of Dactylispa. In this study, we aim to report the first eight mitochondrial genomes of Dactylispa and construct the first molecular phylogeny of Dactylispa to explore the evolutionary relationships among the three subgenera (Dactylispa s. str., Triplispa, and Platypriella).

Sampling and DNA Extraction
The specimens were collected from the following places: Xiangshan, Xunwu County, Jiangxi Province (24 •

Genomic Sequencing, Assembly, and Annotation
Total mitogenomes were obtained by next-generation sequencing (NGS) with the whole-genome shotgun (WGS) strategy on the Illumina MiSeq platform [19]. The reads of sequences were cleaned to remove the adaptor sequences and low-quality reads with ambiguous sequences [19]. The high-quality second-generation sequences were de novo assembled to obtain contigs and scaffolds using A5-miseq v20150522 [25] and SPAdes v3.9.0 [26]. Sequences were extracted based on the sequencing depth, and the sequences with high sequencing depths were then identified by the NCBI NT library using BLASTN (BLAST v2.2.31+) [27]. MUMmer v3.1 [28] was used to perform collinearity analysis, confirm the contig positions, and fill the gaps between contigs. Pilon v1.18 [29] was applied to correct the results and obtain the final mitochondrial sequences. Each mitogenome was annotated using Geneious R11 software, and coding regions were found through ORFfinder (https://www.ncbi.nlm.nih.gov/orffinder/, accessed on 12 November 2020) and manually verified by the BLAST tool on the NCBI website. tRNA structures were predicted and determined by MITOS2 (http://mitos.bioinf.uni-leipzig.de/index.py, accessed on 12 November 2020). The rRNA gene boundaries were interpreted as the end of a bounding tRNA gene and the alignment of sequences with homologous regions of known coleopteran mitogenomes. Circular mitogenome maps were generated using the CG View server V 1.0 [30,31].

Bioinformatic Analysis
Geneious R11 was employed to calculate the A + T and G + C contents, as well as the AT skew and GC skew. Here, AT skew = (A − T)/(A + T) and GC skew = (G − C)/(G + C). DnaSP v5 [32] was used to calculate the nonsynonymous substitution rate to synonymous substitution rate (Ka/Ks) ratio, and PhyloSuite [33] was used to calculate the relative synonymous codon usage (RSCU) for PCG analysis. Annotated sequences of the eight mitogenomes of Dactylispa beetles mentioned above were deposited in GenBank with the accession numbers specified in Table 1. We used DnaSP [32] to calculate the nucleotide divergence of 13 PCGs, 22 tRNA genes and two tRNA genes among eight Dactylispa species and compared the nucleotide diversity among the three subgenera of Dactylispa.
Each mitogenomic gene was aligned in MAFFT [34] with the normal alignment mode for the 24 RNAs and with the codon alignment mode for the 13 PCGs. The code table was '5 Invertebrate mitochondrial', and the strategy was 'auto'. The multiple alignments were then concatenated into one FASTA file. The best partitioning scheme and evolutionary models for 37 predefined partitions were then selected using PartitionFinder2 [35] with the greedy algorithm and AICc criterion (File S1). Maximum likelihood (ML) inference was performed in IQ-TREE [36] with the optimal model of GTR + I + G for 1000 bootstrap replications. Bayesian inference (BI) was performed using MrBayes 3.2.6 [37] under the GTR + I + G + F model with 2 million MCMC generations and 2 parallel runs in which the initial 25% of sampled data were discarded as burn-in. All of the above phylogenetic analyses were performed step by step in PhyloSuite [33].
The amino acid sequences of eight Dactylispa species and two outgroup species were translated from the nucleotide sequences of the 13 PCGs. The sequences were extracted, aligned, and concatenated using the above procedures in PhyloSuite [33]. Using the model of CAT + GTR, phylo-Bayesian (PB) inference was performed in Phylobayes MPI [38,39] on XSEDE (1.8c) with default parameters in the CIPRES Science Gateway [40].

Genome Organization and Nucleotide Composition
The eight sequenced mitogenomes of Dactylispa species have lengths ranging from 17,189 bp to 20,363 bp (Table 1; Figure 1 and Figure S1). Each mitogenome contains the typical set of 37 genes (13 PCGs, 22 tRNAs, and 2 rRNAs) and an A + T-rich region (Table S1). In total, 23 genes (9 PCGs and 14 tRNAs) are located on the positive strand (N-strand); the other 14 genes (4 PCGs, 8 tRNAs and 2 rRNAs) are located on the reverse strand (J-strand) ( Figure S1). The nucleotide compositions of Dactylispa mitogenomes are biased toward A/T A + T contents ranging from 71.4% to 76.6%-values well within the range found in ously sequenced beetles (65.66%-80.72%) [43]. Among the mitogenomes, those of D. discalis and D. albopilosa exhibit the highest and lowest A + T contents, respectively. T skew and GC skew present similar patterns in all Dactylispa mitogenomes, with po AT skews (from 0.093 to 0.154) and negative GC skews (from −0.262 to −0.110) (Tab

Protein-Coding Genes
The length of a single PCG varies from 153 bp (atp8) in all species to 1709 bp ( in D. latispina. The total length of all 13 PCGs ranges from 10,959 bp in D. approxim 10,998 bp in D. paucispina without stop codons (Table S1).
The A + T content of the 13 PCGs ranges from 70.6% to 75.5%. Among the spec approximata shows the highest A + T content, while D. albopilosa presents the lowest content. The PCGs exhibit positive AT skews (from 0.107 to 0.144) and negative GC The nucleotide compositions of Dactylispa mitogenomes are biased toward A/T, with A + T contents ranging from 71.4% to 76.6%-values well within the range found in previously sequenced beetles (65.66-80.72%) [43]. Among the mitogenomes, those of D. nigrodiscalis and D. albopilosa exhibit the highest and lowest A + T contents, respectively. The AT skew and GC skew present similar patterns in all Dactylispa mitogenomes, with positive AT skews (from 0.093 to 0.154) and negative GC skews (from −0.262 to −0.110) (Table S2).

Protein-Coding Genes
The length of a single PCG varies from 153 bp (atp8) in all species to 1709 bp (nad5) in D. latispina. The total length of all 13 PCGs ranges from 10,959 bp in D. approximata to 10,998 bp in D. paucispina without stop codons (Table S1).
The A + T content of the 13 PCGs ranges from 70.6% to 75.5%. Among the species, D. approximata shows the highest A + T content, while D. albopilosa presents the lowest A + T content. The PCGs exhibit positive AT skews (from 0.107 to 0.144) and negative GC skews (from −0.254 to −0.212) (Table S2). Among the 13 PCGs of the eight Dactylispa species, cox1 has the lowest Ka/Ks ratio, while atp8 has the highest ( Figure 2).  (Table S2). Among the 13 PCGs of the eight Dactylispa species, cox1 has the lowest Ka/Ks ratio, while atp8 has the highest ( Figure 2).

Figure 2.
Ka/Ks ratios of 13 protein-coding genes. Ka is the nonsynonymous substitution rate, and Ks is the synonymous substitution rate.

Codon Usage
The most frequently used codons in the eight Dactylispa species are UUA, AUU, UUU, AUA, and AAU; in addition, GCG, CGC, AGC, CCG, and UGC are rarely used. However, RSCU indicated strong codon usage and amino acid composition biases. The codon usage comparison for the amino acids indicates that the codon UCN (Ser) is more frequently used than the codon AGN. For Leu, the UUN codon is more frequently used than CUN.

tRNAs and rRNAs
The 22 typical tRNAs were detected in all eight species (Figure 3), as in other beetles with published mitogenomes. All anticodons are highly conserved, as in other beetle species [48]. The total length of all tRNAs in the eight Dactylispa species mitogenomes range from 1331 bp to 1336 bp. Twenty-one tRNAs display the classic clover-leaf secondary structure, whereas trnS1 lacks the dihydrouridine (DHU) arm and forms a simple loop, which is common in other metazoan mitogenomes. Although somewhat less effective than conventional tRNAs, this abnormal tRNA is functional.
The A + T content of tRNAs ranges from 74.7% to 77.3%. Among the species, D. approximata and D. latispina exhibit the highest A + T content, while D. albopilosa presents the lowest A + T content. The tRNAs exhibit positive AT skews (from 0.048 to 0.100) and negative GC skews (from −0.154 to −0.080) ( Table S2).
As in other insect mitogenomes, the 16S and 12S rRNA genes in Dactylispa species are encoded on the J-strand and located at conserved positions between trnL1 and trnV and between trnV and the control region, respectively [48]. The length of 16S ranges from 1035 bp in D. longispina to 1290 bp in D. albopilosa, and that of 12S varies from 719 bp in D. planispina to 793 bp in D. nigrodiscalis. The A + T content of 12S ranges from 75.50% in D.

Codon Usage
The most frequently used codons in the eight Dactylispa species are UUA, AUU, UUU, AUA, and AAU; in addition, GCG, CGC, AGC, CCG, and UGC are rarely used. However, RSCU indicated strong codon usage and amino acid composition biases. The codon usage comparison for the amino acids indicates that the codon UCN (Ser) is more frequently used than the codon AGN. For Leu, the UUN codon is more frequently used than CUN.

tRNAs and rRNAs
The 22 typical tRNAs were detected in all eight species (Figure 3), as in other beetles with published mitogenomes. All anticodons are highly conserved, as in other beetle species [48]. The total length of all tRNAs in the eight Dactylispa species mitogenomes range from 1331 bp to 1336 bp. Twenty-one tRNAs display the classic clover-leaf secondary structure, whereas trnS1 lacks the dihydrouridine (DHU) arm and forms a simple loop, which is common in other metazoan mitogenomes. Although somewhat less effective than conventional tRNAs, this abnormal tRNA is functional.
The A + T content of tRNAs ranges from 74.7% to 77.3%. Among the species, D. approximata and D. latispina exhibit the highest A + T content, while D. albopilosa presents the lowest A + T content. The tRNAs exhibit positive AT skews (from 0.048 to 0.100) and negative GC skews (from −0.154 to −0.080) (Table S2). As in other insect mitogenomes, the 16S and 12S rRNA genes in Dactylispa species are encoded on the J-strand and located at conserved positions between trnL1 and trnV and between trnV and the control region, respectively [48]. The length of 16S ranges from 1035 bp in D. longispina to 1290 bp in D. albopilosa, and that of 12S varies from 719 bp in D. planispina to 793 bp in D. nigrodiscalis. The A + T content of 12S ranges from 75.50% in D. longispina to 79.20% in D. approximata. Hence, there is no substantial size variation in 12S among the mitogenomes of the eight Dactylispa species (Table S2).

Gene Overlap, Intergenic Spacers and Noncoding A + T-Rich Regions
The size of overlapping regions was examined, varying from 1 bp to 24

Genetic Diversity among Mitogenomes
Mitogenome diversity shows similar trends among the three subgenera (Dactylispa s. str., Triplispa, and Platypriella). In general, Dactylispa s. str. has the lowest nucleotide diversity, and Triplispa and Platypriella have similar patterns of nucleotide diversity ( Figure 5).

Phylogenetic Analysis
Both ML and BI trees based on all mitochondrial genes ( Figure 6) showed similar topologies. The genus Dactylispa was strongly supported as a monophyletic group (BP = 100, PP = 1), which was further divided into two major clades: the first clade includes the monophyletic subgenus  The largest noncoding region in the mitogenomes of the eight Dactylispa beetles as well as most insects is usually located between rrnS and trnI. It might be identified as the putative control region according to the conserved region aligned with other Coleoptera mitogenomes [49]. The control regions of the eight Dactylispa mitogenomes vary widely, ranging in length from 1819 bp in D. approximata to 3069 bp in D. albopilosa, with an average A + T content of 77.61% (Table S2).

Genetic Diversity among Mitogenomes
Mitogenome diversity shows similar trends among the three subgenera (Dactylispa s. str., Triplispa, and Platypriella). In general, Dactylispa s. str. has the lowest nucleotide diversity, and Triplispa and Platypriella have similar patterns of nucleotide diversity ( Figure 5).

Genetic Diversity among Mitogenomes
Mitogenome diversity shows similar trends among the three subgenera (Dactylispa s. str., Triplispa, and Platypriella). In general, Dactylispa s. str. has the lowest nucleotide diversity, and Triplispa and Platypriella have similar patterns of nucleotide diversity ( Figure 5).

Phylogenetic Analysis
Both ML and BI trees based on all mitochondrial genes ( Figure 6) showed similar topologies. The genus Dactylispa was strongly supported as a monophyletic group (BP = 100, PP = 1), which was further divided into two major clades: the first clade includes the monophyletic subgenus Dactylispa s. str. (D. approximata, D. albopilosa and D. longispina) and the species D. paucispina belonging to Triplispa with high support in the BI tree (PP =

Phylogenetic Analysis
Both ML and BI trees based on all mitochondrial genes ( Figure 6) showed similar topologies. The genus Dactylispa was strongly supported as a monophyletic group (BP = 100, PP = 1), which was further divided into two major clades: the first clade includes the monophyletic subgenus Dactylispa s. str. (D. approximata, D. albopilosa and D. longispina) and the species D. paucispina belonging to Triplispa with high support in the BI tree (PP = 0.99) but moderate support in the ML tree (BP = 76). The second clade presented a nested form and consisted of the two Triplispa species (D. xanthopuse and D. nigrodiscalis) and the monophyletic subgenus Platypriella (D. latispina and D. planispina) with similar support as the former clade (BP = 88, PP = 1). In addition, the PB tree based on the amino acids of the 13 PCGs (Figure 7) was similar to the two nucleotide-based trees (Figure 6), except for the fact that D. paucispina was not included on the Dactylispa s. str. branch. The monophyly of both the subgenera Dactylispa s. str. and Platypriella was strongly supported by the phylogenetic trees, while the subgenus Triplispa appeared to be a paraphyletic group.  (Figure 7) was similar to the two nucleotide-based trees (Figure 6), except for the fact that D. paucispina was not included on the Dactylispa s. str. branch. The monophyly of both the subgenera Dactylispa s. str. and Platypriella was strongly supported by the phylogenetic trees, while the subgenus Triplispa appeared to be a paraphyletic group. The Dactylispa species were separated according to their body shapes: slim species were on one branch, and broad species were on the other branch ( Figure 6). Insect-host plant relationships might also help explain part of the Dactylispa phylogeny: all three species of Dactylispa s. str. are Poaceae-feeding, and both Platypriella species are Fagaceaefeeding, but the position of Triplispa could not be explained by host relationships alone [1].

Discussion
In this study, we sequenced and analyzed the first eight mitochondrial genomes of Dactylispa species. All of the sequenced Dactylispa mitogenomes exhibited the same gene arrangement, which was also consistent with the putative ancestral pattern in insects [50].
All Dactylispa sequences have the same 7 bp overlap "ATGATAA" between atp8 and atp6 and 8 bp overlap "AAGCCTTG" between trnW and trnC ( Figure 4). Additionally, the atp6 and atp8 of the eight Dactylispa species have the same respective lengths of 663 bp and 153 bp. Similar intergenic spacers were also found in other Coleoptera, as well as in other insect orders [51]. Although these intergenic spacers have been hypothesized to be binding sites for translation termination [52], more research is still required to ascertain the functions of conserved noncoding regions in the mitogenomes of Coleoptera.
The Ka/Ks ratio is a powerful method of diagnosing the strength and mode of evolutionary selection on genes [53][54][55]. The Ka/Ks ratios for all PCGs of the eight Dactylispa species showed that cox1 had the lowest evolutionary rate, with atp8 having the highest (Figure 2). The cox1 gene was highly conserved in nearly all animals [19,56] and is thus the best DNA barcode in animal taxonomy and species phylogenetics [57]. The atp8 gene evolves rapidly in many animals [58,59], indicating that it is usually under weak selection pressure [19,60].
Our phylogenetic analysis indicates that the Dactylispa genus is separated into two clades, although these clades were strongly supported in the BI tree but only moderately supported in the ML tree ( Figure 6). All three phylogenomic trees (Figures 6 and 7) support the delimitation of the subgenera Dactylispa s. str. and Platypriella as proposed in Chen's system [1]. However, Triplispa seems to be an invalid paraphyletic group. One reason is that Chen's concept of Triplispa is broad, and he did not try to define any other The Dactylispa species were separated according to their body shapes: slim species were on one branch, and broad species were on the other branch ( Figure 6). Insect-host plant relationships might also help explain part of the Dactylispa phylogeny: all three species of Dactylispa s. str. are Poaceae-feeding, and both Platypriella species are Fagaceaefeeding, but the position of Triplispa could not be explained by host relationships alone [1].

Discussion
In this study, we sequenced and analyzed the first eight mitochondrial genomes of Dactylispa species. All of the sequenced Dactylispa mitogenomes exhibited the same gene arrangement, which was also consistent with the putative ancestral pattern in insects [50].
All Dactylispa sequences have the same 7 bp overlap "ATGATAA" between atp8 and atp6 and 8 bp overlap "AAGCCTTG" between trnW and trnC ( Figure 4). Additionally, the atp6 and atp8 of the eight Dactylispa species have the same respective lengths of 663 bp and 153 bp. Similar intergenic spacers were also found in other Coleoptera, as well as in other insect orders [51]. Although these intergenic spacers have been hypothesized to be binding sites for translation termination [52], more research is still required to ascertain the functions of conserved noncoding regions in the mitogenomes of Coleoptera.
The Ka/Ks ratio is a powerful method of diagnosing the strength and mode of evolutionary selection on genes [53][54][55]. The Ka/Ks ratios for all PCGs of the eight Dactylispa species showed that cox1 had the lowest evolutionary rate, with atp8 having the highest (Figure 2). The cox1 gene was highly conserved in nearly all animals [19,56] and is thus the best DNA barcode in animal taxonomy and species phylogenetics [57]. The atp8 gene evolves rapidly in many animals [58,59], indicating that it is usually under weak selection pressure [19,60].
Our phylogenetic analysis indicates that the Dactylispa genus is separated into two clades, although these clades were strongly supported in the BI tree but only moderately supported in the ML tree ( Figure 6). All three phylogenomic trees (Figures 6 and 7) support the delimitation of the subgenera Dactylispa s. str. and Platypriella as proposed in Chen's system [1]. However, Triplispa seems to be an invalid paraphyletic group. One reason is that Chen's concept of Triplispa is broad, and he did not try to define any other subgenera because all "fit" within his concept [1]. In general, it is difficult to find true synapomorphies of Dactylispa as they have so many homoplasies that each having various states is not conservative. In theory, one can propose as many subgenera as they wish using peculiar characters. However, this is not a good approach in practice. In our study, not all species classified within the subgenus Triplispa belonged and thus rendered the subgenus polyand paraphyletic. However, this does not contradict its validity; we simply do not have enough data to resolve this with certainty.
The most interesting result is the separation of D. paucispina from the rest of the Triplispa-Platypriella clade (Figures 6 and 7), rendering Triplispa polyphyletic. The entire D. angulosa group is odd looking, and our results indicate that it might (1) warrant its own subgenus or (2) belong to Dactylispa s. str. However, this branching has the lowest bootstrap value among the three that we are presenting (Figures 6 and 7). Setting D. paucispina aside does not invalidate Triplispa, instead only showing that the phylogenetic history of the group is especially complicated. For the rest of Triplispa-Platypriella, there are at least two possible hypotheses: (1) the separation of Platypriella renders Triplispa paraphyletic; or (2) both are monophyletic with some other internal groups that would have to be separated. However, we lack enough data to make a clear decision, and our study is more of a first insight into the phylogeny of this genus.
According to the morphological characters of Chen's system [1], the subgenus Triplispa has many intermediate or transitional characters between those of the other two subgenera. For example, the moderate body shape of Triplispa species is in transition from the slim shape of Dactylispa s. str. and the broad shape of Platypriella ( Figure 6). Mitogenome nucleotide diversity indicated that the differences between Triplispa and Platypriella are small but the differences between the two subgenera and Dactylispa s. str. are large ( Figure 5). Insect-host plant associations might somewhat explain the relationships of such subgenus ( Figure 6). Most Dactylispa s. str. species are specialized on monocotyledonous plants such as Poaceae [1][2][3][4][5], the same as many primitive groups of hispine beetles [61,62]. In contrast, most Triplispa species and most Platypriella species feed on dicotyledonous plants. However, Platypriella seems to prefer Fagaceae plants, while Triplispa seems to prefer other dicotyledonous families [1][2][3][4][5].
In the inference of high-level phylogeny in Coleoptera, the PB tree with the heterogeneoussite model performs better than both BI and ML trees with homogeneous-site models [43], and the amino acid-based trees have a more reliable topology than the nucleotide-based trees [43]. In our case, all the three inferences (PB, BI, and ML) are very similar in terms of tree topology except that D. paucispina is clustered with the Dactylispa s. str. branch in both BI and ML trees ( Figure 6) but not so in the PB tree (Figure 7). Both body shape and host plant relation support that D. paucispina should be isolated from the Dactylispa s. str. species. Compared to both BI and ML, PB might also infer better low-level phylogeny in Coleoptera such as the subgenus relationships among the Dactylispa genus.
The mitogenome is a good resource for separating species that are extremely similar in morphology, such as Dactylispa spp. However, our sample is still rather small. In the future, more species should be included to solve the evolutionary issues in this genus of leaf-mining beetles.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/insects12111005/s1. Figure S1: Circular diagrams of eight Dactylispa mitogenomes. Table S1: Summary of eight Dactylispa mitochondrial genomes. Table S2: Nucleotide composition of mitochondrial genomes of eight Dactylispa species. File S1: PartitionFinder output files. Data Availability Statement: All the mitochondrial genome sequences were submitted to Gen-Bank (accessions MN016958-MN016966, as in Table 1), and they will be accessible when the article is published.