The voucher specimens of 5 adult leaf beetles were preserved in anhydrous ethanol and stored in a refrigerator at -20°C before DNA extraction. All specimens were collected from tobacco fields and detailed collection information is provided in Supplementary Table S1 . After being identified based on morphological characters, the genomic DNA was extracted using DNeasy Blood and Tissue kit (Qiagen, Germany) on the basis of the manufacturer’s protocol, and the voucher specimen was stored at Entomological Museum of China Agricultural University.

The Illumina TruSeq library was prepared with an average insert size of 350 bp and sequenced with the paired-end reads length of 150 bp on Illumina NovaSeq 6000 platform (Berry Genomic, Beijing, China). A total of 6 Gb raw data was obtained. We used Prinseq version 0.20.4 (Schmieder and Edwards, 2011) to remove short and low-quality reads with the poly-Ns > 15 bp, or > 75 bp bases with quality score < 3. The remaining reads were de novo assembled using IDBA-UD (Peng et al., 2012), with minimum and maximum k values of 41 and 141 bp, respectively. The DNA barcode regions of each species (~658 bp) were amplified by the following primer pair LCO1490/HCO2198. PCR cycling condition was: 95°C for 1 min, 40 cycles of 95°C for 20 s, 50°C for 50 s, and 68°C for 1.5 min, followed by 72°C for 5 min. The PCR products were sequenced by Sanger sequencing at Tsingke Biotechnology (Beijing, China). To identify the corresponding mitogenome assemblies, the assembled contigs were searched with COX1 sequence using BLAST with at least 98% similarity. To investigate the assembled accuracy and sequencing depth, clean reads were mapped using Geneious version 10.1.3 (Kearse et al., 2012). Finally, we obtained five complete circular mitogenomes.

All five mitogenome sequences were preliminarily annotated by MitoZ (Meng et al., 2019) and accurately corrected in Geneious. The locations and secondary structures of tRNA genes were determined by tRNAscan-SE search server (Lowe and Eddy, 1997) and ARWEN version 1.2 (Laslett and Canbäck, 2008). PCGs and rRNA genes were identified by alignment with homologous genes of other Galerucinae species, and the secondary structure of rRNAs was predicted by RNAfold WebServer Version 2.4.18 (Lorenz et al., 2011) online platform. The control region of each species was confirmed by mapping the read sequences to the mitogenome with a depth > 1,000×. The circular mitogenome maps with A + T content and G + C content were generated using Geneious (Kearse et al., 2012). The nucleotide composition of mitogenome and relative synonymous codon usage (RSCU) of PCGs were analyzed by MEGA version 7.0 (Kumar et al., 2016). Nucleotide compositional differences were calculated by the formulae: AT skew = (A – T)/(A + T) and GC skew = (G - C)/(G + C) (Perna and Kocher, 1995). The tandem repeat units of the control region were identified using the Tandem Repeats Finder server (https://tandem.bu.edu/trf/home, accessed on 5 July 2024) (Benson, 1999).

The phylogenetic analysis of Galerucinae was performed by merging 92 mitogenomes published in GenBank and our newly determined 5 mitogenomes. Chrysomela vigintipunctata and Gonioctena intermedia from Chrysomelinae, as well as Physosmaragdina nigrifrons and Cryptocephalus dimidiatipennis from Cryptocephalinae were selected as outgroups, according to a prior study (Douglas et al., 2023) ( Supplementary Table S2 ).

All 13 PCGs of each mitogenome were aligned separately using the L-INS-I strategy of the MAFFT algorithm (Katoh and Standley, 2013) implemented in TranslatorX (Abascal et al., 2010). Two rRNA genes were aligned individually using the G-INS-I strategy of MAFFT version 7.0 online server (Katoh and Standley, 2013). All alignments were checked manually in MEGA (Kumar et al., 2016). Gene fragments were imported into Geneious and concatenated into two datasets: (1) the PCGRNA matrix with 13,499 nucleotides, corresponding to 13 PCGs and two rRNA genes; and (2) the PCG12RNA matrix with 9,881 nucleotides, corresponding to the first and second codon positions of the 13 PCGs and two rRNA genes.

Both datasets were analyzed under the strategy of Bayesian inference (BI) and maximum-likelihood (ML). The BI phylogenetic trees were recovered using PhyloBayes MPI version 1.5a (Lartillot et al., 2013) under the site-heterogeneous mixture CAT + GTR model (Tavare, 1986; Lartillot and Philippe, 2004). The Markov Chain Monte Carlo chains were run independently after removing constant sites from the alignment and were stopped after the two runs had satisfactorily converged (maxdiff < 0.1). A consensus tree was yielded from the remaining trees after discarding initial 25% trees of each run as “burn-in”. The ML phylogenetic trees were recovered using IQ-TREE web server (Trifinopoulos et al., 2016), with automatic model prediction and 1000 ultrafast bootstrap replicates.

A Bayesian algorithm of MCMCTree in PAML 4.9j (Yang, 2007) package was implemented to estimate divergence times within Galerucinae with the molecular clock under the auto-correlated rates model. We used the topology generated by PCRRNA and ML method as the fixed tree for subsequent analysis. The Dirichlet-gamma prior to the overall substitution rate (rgene gamma) was set to G (1, 1.72), calculated by baseml.

Divergence time estimation was performed based on three vetted internal calibrations nodes. The root age for the divergence between Alticini and Galerucini was set to 102.21 million years ago (Ma) according to the results of Nie et al. (2021). The fossil Theonenosa nigripennis served as the calibration node for the tribe Galerucini with the upper- and lower-bound was assigned to 99.4 and 98.2 Ma following the recommendation of Nie et al. (2021). The maximum- and minimum-bound values of fossil record (15.97-20.44 Ma), Altica tholimorpha provided by the Paleobiology Database (https://paleobiodb.org/) was assigned to the crown node of Altica. MCMCTree analyses used two independent MCMC chains with a 1,000,000 generation burnin, subsequently sampling every 100 generations until 10,000 samples were collected. To further ensure convergence, effective sample size (ESS) values of parameters were calculated for the combined results in Tracer v1.7.1 (Rambaut et al., 2018).

We obtained complete mitogenomes of five leaf beetles (GenBank accession number PQ038342-PQ038346) ( Figure 1 ; Supplementary Tables S3 – S7 ). The mitogenomes of these leaf beetles ranged in size from 15,737 (O. maculatus) to 16,484 bp (Euphitrea sp.) and included conventional 37 coding genes (13 PCGs, two rRNAs, and 22 tRNAs) and a control region (CR). Four PCGs (ND5, ND4, ND4L and ND1), eight tRNA genes (trnQ, trnC, trnY, trnF, trnH, trnP, trnL^(CUN) and trnV) and two rRNA genes (srRNA and lrRNA) were encoded on the minority strand (N-strand) while the other 23 genes were encoded on the majority strand (J-strand). The gene orders of mitogenomes of five leaf beetles were consistent with the ancestral gene order of Drosophila yakuba, which is thought to be the putative ancestral pattern of insect mitogenomes (Clary and Wolstenholme, 1985; Boore, 1999). The AT content of the five complete mitogenomes of Galerucinae varied from 77.7% in Altica sp. to 80.9% in O. maculatus ( Table 1 ), suggesting a strong AT bias, similar to that observed in other mitogenomes of leaf beetles (Wang et al., 2019; Li et al., 2023). The five mitogenomes all exhibited a positive AT skew, ranging from 0.063 (O. bowringii) to 0.033 (Euphitrea sp.), and a negative GC skew, varying from -0.234 (Euphitrea sp.) to -0.107 (O. maculatus) ( Table 1 ), indicating that the nucleotide composition was significantly biased toward A and C.

All newly sequenced Galerucinae mitogenomes comprised 13 PCGs (ND4, ND4L, ND5, and ND1 coded on the N-strand and COX1, COX2, COX3, CYTB, ATP6, ATP8, ND2, ND3, and ND6 coded on the J-strand). The total length of these 13 PCGs ranged from 11,076 bp (O. maculatus) to 11,189 bp (Euphitrea sp.), and the A + T content of PCGs ranged from 75.9% (Altica sp.) to 79.5% (O. maculatus) ( Table 1 ). The mean A + T content of PCGs (77.1%) were slightly lower than that of the whole mitogenome (78.9%). The AT skew (-0.147 to -0.060) indicated a bias towards T in the nucleotide composition of PCGs, which was different from that of the whole genome ( Table 1 ). Additionally, the GC skew of PCGs in the five Galerucinae species ranged from -0.156 to 0.054 ( Table 1 ). The negative AT skew was consistent in all three codons, while the positive GC skew was consistent in the first codon and the negative GC skew in the third codon of PCGs across five leaf beetles ( Supplementary Table S8 ).

We calculated the RSCU values, and found the six most prevalent codons UUA, GUU, GGA, UCU, CGA and GCU were mainly composed by T and/or A ( Figure 2 ). At the third codon position, A or T were overrepresented in PCGs. The biased usage of A + T nucleotides was reflected in codon frequencies. The codon usage pattern of five newly sequenced mitogenomes was similar to the pattern found in previously reported Galerucinae species (Li et al., 2023). In five leaf beetles, all 13 PCGs were initiated with ATN as the start codon (eight with ATA, six with ATC, 27 with ATG and 21 with ATT), except for ND1 in G. nigropicta, O. bowringii and O. maculatus, which used TTG as the start codon ( Supplementary Tables S3 – S7 ). The complete stop codon TAA and TAG was most widely assigned to ten PCGs, while single T and TA residue were found in COX2 and ND4 respectively as incomplete stop codons ( Supplementary Tables S3 – S7 ), and the stop codons are presumed to be completed through post-transcriptional polyadenylation (Stewart and Beckenbach, 2009; Zhao et al., 2023).

To analyze the evolutionary rate of PCGs, we calculated the nonsynonymous substitution rate (Ka), synonymous substitution rate (Ks), and the ratio of Ka/Ks for each PCG in five Galerucinae species ( Figure 3 ). The highest average Ka value was observed in ATP8, whereas ND2 exhibited the highest Ks value. The average Ka/Ks value ranged from 0.092 (COX1) to 1.065 (ATP8). All Ka/Ks ratios were lower than 1 except for ATP8, indicating that mutations in these PCGs were primarily exchanged through synonymous substitutions and evolved under purifying selection. On average, the Ka/Ks value of ATP8 was greater than 1, indicating that positive selection played a predominant role in the evolution of ATP8 gene. In contrast, COX1 exhibited the lowest evolutionary rate, and underwent the strongest purifying selection. The analysis revealed that almost all the genes mitochondrial functional genes in Galerucinae have been influenced by purifying selection and only one by positive selection, with none evolving neutrally. The phenomenon was also reported in other Galerucinae species or insect groups (Huang et al., 2025; Li et al., 2023; Liu et al., 2025).

Twenty-two typical tRNA genes were identified in the five mitogenomes of Galerucinae. The mean length of tRNAs varied from 64.4 bp (O. maculatus) to 65.7 bp (Euphitrea sp.) ( Supplementary Tables S3 – S7 ). The longest tRNA is trnK of Altica sp. and Euphitrea sp., whereas the shortest tRNA is trnR of O. maculatus. All tRNAs could be folded into a typical clover-leaf structure, except for trnS^(AGN) , where the dihydrouridine (DHU) arm forms a simple loop ( Supplementary Figure S1 ). The structure is commonly found in other beetles’ mitogenomes (Li et al., 2023; Huang et al., 2023b; Zhao et al., 2025). The lrRNAs of Galerucinae mitogenomes were located between trnV and trnL^(CUN) , and srRNAs were located between trnV and the CR. The lengths of lrRNA genes ranged from 1,265 (O. bowringii) to 1,284 bp (Euphitrea sp.), and those of srRNAs ranged from 793 (O. bowringii) to 826 bp (Euphitrea sp.). The mean AT content of lrRNAs is 82.9%, and the mean AT content of srRNAs is slightly lower at 81.9%.

As an indispensable noncoding fragment, the CR of mitogenome initiates the transcription and replication (Zhang et al., 1995; Ji et al., 2019). The CR of five mitogenomes of Galerucinae was located between srRNA and trnI, ranging in size from 1,220 bp (O. maculatus) to 1,836 bp (Euphitrea sp.) ( Supplementary Tables S3 – S7 ). The A + T content was ranged from 84.3% (Altica sp.) to 88.7% (O. maculatus), with mean AT skew of 0.025 and mean CG skew of -0.309. The CRs of insect mitogenomes comprise several structural elements (tandem repeat sequences, poly-T regions, GA-rich regions, and replication origin sequences etc.), that lack evident patterns of conservation. We analyzed the CRs of five Galerucinae species and identified their structural elements. The tandem repeats and poly-T regions were all found in five CRs, while the poly-A regions were only found in O. maculatus ( Figure 4 ). Two types of tandem repeat regions were present in Altica sp., with lengths of 314 bp and 77 bp, and three poly-T regions were found between these tandem repeat regions. In Euphitrea sp., three tandem repeat regions were identified, with repeat units of 72 bp, 35 bp, and 22 bp, each containing three units. Three poly-T regions were located between the first and second tandem repeat regions, and another poly-T region was situated between the third repeat region and trnI. The CR of G. nigropicta contained only one 40 bp repeat region and one 16 bp poly-T region. In O. bowringii, three tandem repeat regions with lengths of 185 bp, 95 bp and 77 bp were identified, and there were four poly-T regions interspersed between them. In O. maculatus, four tandem repeats were found, along with three poly-T regions. Additionally, there were two poly-A regions each 10 bp in size, located between the third and fourth tandem repeats, and between the fourth tandem repeat and trnI, respectively. The tandem repeat regions, poly-T regions, and poly-A regions are commonly found in the CRs of insect mitogenomes (Huang et al., 2023a; Xu et al., 2021; Li et al., 2023). The variability in the number and length of structural elements, particularly tandem repeats, results in significant differences in the CR among species. These characteristics were then proposed as novel and reliable markers for studies in phylogenetics and evolution (Zhang et al., 1995; Li and Liang, 2018).

The phylogenetic studies of Galerucinae were conducted based on the two datasets using BI and ML approaches yield broadly consistent topology with high node support values in Galerucinae ( Figure 5 ). The subfamily Galerucinae has been split into two tribes, the Alticini and Galerucini (Bouchard et al., 2011). Both Alticini and Galerucini were inferred as monophyletic groups, exhibiting strong sister-group relationships supported by high Bayesian posterior probabilities (BPP = 1/1) and ML bootstrap values (BSV = 83) inferred from PCG12RNA, but moderately supported by ML analysis based on PCGRNA dataset (BSV = 65). The phylogenetic relationships were congruent with those of previous studies (Nie et al., 2018; Douglas et al., 2023; Li et al., 2023). Within the tribe Alticini, the high consistency of topological structures obtained from various methods and data sets indicates the robustness and reliability of the results. Specifically, the cluster (Altica group + (Blepharida group + Nisotra group)) was observed to form sister groups with other groups of Alticini. This finding is consistent with previous studies based on 13 PCGs of mitogenome and nuclear genes, including 18S, 28S-D2, and 28S-D3 regions (Nie et al., 2018). Furthermore, the sister relationship between the Sphaeroderma group and the Chabria group recovered in our results was also supported by research conducted by Nie et al. (2018). However, in an earlier study, the Sphaeroderma group formed a sister group with the Nisotra group by analyzing the lrRNA and nuclear genes 18S and 28S (Ge et al., 2012).

Within the tribe Galerucini, all analyses recovered the subtribes Oidina, Diabroticina, Galerucina and Monoleptina as monophyletic, whereas the subtribes Metacyclina, Aulacophorina and Luperina were paraphyletic. Both of genera Nonarthra and Acrocrypta were nested in the tribe Hylaspina, which is consisted with previous study (Nie et al., 2018). Metacyclina consistently grouped with Galerucina despite the limited taxon sampling (Nie et al., 2018; Ge et al., 2011; Gillespie et al., 2004). The subtribe Oidina formed a monophyletic clade, branching at the base of Galerucini and constituting sister groups with the rest of the Galerucini, which was similar to previous findings (Li et al., 2023). In addition, within the six-subtribe classification system, the subtribe Luperina was found to be paraphyletic (Gillespie et al., 2004; Ge et al., 2011), and it was subsequently split into three subtribes: Aulacophorina, Diabroticina and Monoleptina based on mitochondrial protein-coding and nuclear genes (Ge et al., 2011). In our analysis, only the monophyly of Diabroticina and Monoleptina was supported, the subtribe Luperina and Aulacophorina were still recovered as paraphyletic. Therefore, we considered that the phylogenetic relationships within the subtribe Luperina still need to be discussed prudently with the addition of more samples and molecular markers. In addition, regarding the classification of the Alticini and Galerucini, some genera are considered “problematic” (see green clades in Figure 5 ) as their placement is not easily determined based on the presence or absence of the metafemoral extensor tendon (MET) (Seeno and Wilcox, 1982; Furth and Suzuki, 1994; Ge et al., 2012). All analyses in our study supported that “problematic genera” Sangariola and Phygasia were included in Alticini, and Luperomorpha, Hespera, Mandarella, Laotzeus, Nonarthra and Acrocrypta were grouped in Galerucini with high support. The result was in line with the previous research based on different molecular data (Ge et al., 2011; Nie et al., 2018).

In order to further understand the evolutionary modes of the CR in Galerucinae, we calculated the AT content and the CR length across 82 species of Galerucinae with complete CR ( Supplementary Table S9 ). Additionally, we investigated structural features in CRs such as tandem repeats, poly-A, and poly-T regions. These findings were then mapped onto phylogenetic trees to provide insights into their evolutionary patterns ( Figure 5 ). The CR in Galerucinae exhibited significantly high AT contents, varying from 77.6% in Laotzeus sp. to 90.3% in Paleosepharia posticata ( Supplementary Table S9 ). The average AT content of CR in Alticini (88.4%) was higher than Galerucini (87.8%), while the variation in AT content was lower in Alticini than in Galerucini. The length of CR exhibited obvious variation among species within subfamily Galerucinae, ranging from 512 bp in Arthrotus chinensis to 2,909 bp in Sphenoraia micans ( Supplementary Table S9 ). The average length of the CR in Alticini (1,348.6 bp) was shorter than that in Galerucini (1,409.4 bp), and the variation in CR length was also less in Alticini.

Many CRs contain tandem repeats, which are associated with the initiation and termination of transcription in mitogenomes (Ji et al., 2019). Previous studies have found that tandem repeats are present in the mitogenomes of various insect groups (Huang et al., 2023a; Li et al., 2023), and that variations in the size and copy number of these tandem repeat units lead to differences in CR sequence length, resulting in heteroplasmy of mitochondrial molecules (Zhang et al., 1995; Oliveira et al., 2007; Huang et al., 2023a). Tandem repeats were observed in 79 out of 82 complete CR sequences, excluding Chaetocnema pelagica from Alticini and two species from Galerucini: Acrocrypta assamensis and Monolepta pallidula ( Figure 5 ). The length of tandem repeats ranged from 40 bp in G. nigrolineata to 1,897 bp in Isotes sp., and the average length of tandem repeats in Alticini (336.3 bp) was shorter than that in Galerucini (429.8 bp), which aligned with the general trend observed in CR lengths ( Figure 5 ). The structural features, poly-A and poly-T regions, were found in 69 and 35 complete CRs, respectively. Both poly-A and poly-T elements simultaneously detected in the CRs of 31 Galerucinae species, including 21 species from the tribe Galerucini and 10 from the tribe Alticini ( Figure 5 ). The poly-A and poly-T are considered to be the signal motifs for initiating mitogenome transcription and are commonly observed in insects (Liang et al., 2022).

In this study, we analyzed the evolutionary patterns of CR of the subfamily Galerucinae for the first time. As a novel molecular marker, the CR exhibits characteristics that offer significant potential for resolving the phylogenetic problems within the subfamily Galerucinae. The “problematic” genera Sangariola and Phygasia exhibited higher AT content compared to Luperomorpha, Hespera, Mandarella, Laotzeus, Nonarthra and Acrocrypta. This finding aligns with the observed trend of AT content between the Alticini and Galerucini tribes, thereby supporting the classification of these “problematic” genera. Furthermore, our results showed that the CRs of most Galerucinae species contained tandem repeats and the poly-A motifs, whereas fewer species contained poly-T motifs. These structural features were unique to each Galerucinae species, and varied significantly between tribes and even within genera. This result was also found in the CRs of other insects, with differences in the number and size of tandem repeats, even for different taxa that were closely related (Li and Liang, 2018; Huang et al., 2023a). By investigating the structural features of the CRs within the phylogenetic framework of Galerucinae, we found that species with TRs > 500 bp were more numerous in Galerucini than in Alticini. No tandem repeats greater than 100 bp in length were detected in the Blepharida group, whereas these tandem repeats were observed in the Nisotra group, which is a closely related sister clade to the Blepharida group. Similarly, poly-A motifs were found in CRs of Chabria and Sphaeroderma groups, but only poly-T motifs were found in the Lanka group, which forms a sister clade to the Chabria + Sphaeroderma groups. The diversity in length and structure organization of the CR was attributed to variations in specific structural elements, which enriched our understanding of the differences in CR structural organization within the economically significant insect groups. Despite our research currently being limited to the genetic level in examining the evolutionary features of CR, our discoveries still underscore the rich molecular diversity present within the subfamily Galerucinae. To gain a more comprehensive understanding, additional experiments are essential to elucidate the specific biological roles of the different elements within the CR.

Divergence time estimates from MCMCTree are summarized in Figure 6 . Based on three calibration nodes, the initial divergences within the Galerucinae were predicted to have occurred during the Middle Cretaceous period (102.69 Ma, 95%HPD: 100.34-105.23 Ma), which is very close to a previous analysis conducted by Nie et al. (2021). Our fossil-calibrated Chronogram also suggests a similar Middle Cretaceous origin of Alticini (96.56 Ma, 95%HPD: 92.65-99.88 Ma) and Galerucini (98.72 Ma, 95%HPD: 98.18-99.37 Ma), which are in line with a prior study (Nie et al., 2021). We found the main groups of Galerucinae diverged and occurred during the period from Middle Cretaceous to Middle Paleogene, indicating this period may serve as an important temporal interval for the diversification of Galerucinae. Flowing plants (angiosperms) diversified quickly during the Cretaceous period and the coordination of diversification between angiosperms and phytophagous beetles have been extensive proved (Nie et al., 2021; Zhang et al., 2018). Therefore, we speculate the extraordinary diversity of Galerucinae may be attributed to co-evolution with angiosperms. However, the studies regarding divergence time estimation within the Galerucinae are still limited hinders further comparison with other results. Based on mitogenomic data, we have provided an initial framework for exploring the internal evolutionary history of Galerucinae.