Specimens of J. variolaris were collected on H. ammodendron in the vicinities of Turpan City, Xinjiang Uygur Autonomous Region, China, on 14 May 2022. Specimens of P. chinensis were collected from Dayaoshan Mountains in Guangxi Zhuang Autonomous Region, China, on 20 March 2021. Specimens of C. japonica were collected from Quanzhou City, Fujian Province, China, on 23 February 2021. The above specimens are preserved in 95% alcohol at -24 °C in specimen collection at China West Normal University, Nanchong, China. Next-generation sequencing and assembly were performed by Beijing Aoweisen Gene Technology Co. Ltd. (Beijing, China) to obtain the complete mitogenome sequences.

The raw data were processed using Trimmomatic v. 0.35 (Bolger et al. 2014) to remove low-quality reads and obtain a high-quality clean data. Finally, 4.8 Gb, 5.28 Gb, and 6.8 Gb clean data were obtained to assemble complete mitogenome of J. variolaris, P. chinensis, and C. japonica, respectively. Three mitogenome sequences were annotated using Geneious 11.0.2 (Kearse et al. 2012) based on the invertebrate mitochondrial genetic code. All tRNA genes were reconfirmed using the online tool MITOS Web Server (Bernt et al. 2013) and the second structures were further predicted using tRNAscan-SE server v. 1.21 (Lowe and Chan 2016). Two rRNA genes were identified by alignment with other buprestid rRNA sequences. Three mitogenome maps were drawn using Organellar Genome Draw v. 1.3.1 (Greiner et al. 2019). Strand asymmetry of mitogenome sequence was calculated using the formulae reported by Perna and Kocher (1995): AT-skew = (A – T)/(A + T), and GC-skew = (G – C)/(G + C). The base composition and relative synonymous codon usage (RSCU) values of three mitogenome sequences were determined using MEGA v. 12.0.0 (Kumar et al. 2016). The non-synonymous substitutions (Ka) and synonymous substitutions (Ks) of all PCG genes were calculated using DnaSP v. 5 (Librado and Rozas 2009). The tandem repeat elements of control region (CR, also known as A + T-rich region) were detected by the online tool Tandem Repeats Finder (Benson 1999). The heterogeneous analysis of nucleotide matrix (13 PCGs + 2 tRNAs) was performed using AliGROOVE v. 1.06 (Kück et al. 2014).

To investigate mitogenome arrangement patterns in Buprestidae, the gene orders of all known buprestid mitogenomes were compared with that of closely related taxa. A total of 22 buprestid mitogenomes (Table 1), including three newly generated sequences in this study, were subjected for phylogenetic analyses, using Heterocerus parallelus Gebler, 1830 (Heteroceridae) and Dryops ernesti Gozis, 1886 (Dryopidae) as outgroups (Xiao et al. 2019; Huang et al. 2022; Wei 2022). The test of substitution saturation for the dataset (13 PCGs + 2 rRNAs) was performed with DAMBE to test whether the sequence is suitable for constructing a phylogenetic tree (Xia 2017). Then, the phylogenetic trees were reconstructed using nucleotide matrix 13 PCGs + 2 rRNAs based on ML and BI methods. The nucleotide matrix was aligned using ClustalW (Thompson et al. 1994) and trimmed by trimAl v. 1.2 (Capella-Gutiérrez et al. 2009). In BI and ML analyses, the best-fit models were deduced by ModelFinder (Kalyaanamoorthy et al. 2017). The phylogenetic trees were reconstructed using IQ-tree v. 1.6.8 (Guindon et al. 2010) and MrBayes v. 3.2.6 (Ronquist et al. 2012) integrated into PhyloSuite v. 1.2.2 (Zhang et al. 2020). During this analyzing process, PhyloSuite was run with previous parameters (Wei 2022).

We sequenced and annotated the complete mitogenome of J. variolaris (GenBank No. OP390084), P. chinensis (No. OP388449), and C. japonica (No. OP388437). Overall, these mitogenome sequences were 15,759 to 16,227 bp in length, which are medium length in Buprestidae (Table 1). It is a circular, double-stranded ring that includes 37 insect mitochondrial genes (13 PCGs, 22 tRNAs, and 2 rRNAs) and an A + T-rich region (control region, CR).

In these three mitogenome, the N-strand encoded the sense-strand of 14 genes (nad1, nad4L, nad4, nad5, trnQ, trnV, trnL1, trnP, trnH, trnF, trnY, trnC, rrnL, and rrnS), while the J-strand encoded the sense-strand of the remaining 23 genes (Table 2), which was consistent with the known buprestid species (Cao and Wang 2019a, b; Xiao et al. 2019; Chen et al. 2021; Peng et al. 2021; Huang et al. 2022; Wei 2022; Weng et al. 2022).

These three mitogenome sequences had a high A + T content, with an average of 68.47%, showing a strong AT bias (Suppl. material 1: table S1). Among them, the A + T content of J. variolaris (70.43%) was higher than of both C. japonica (67.97%) and P. chinensis (67.00%). These three mitogenome sequences showed a positive AT skew (0.12–0.13) and negative GC skew (-0.22), which is consistent with the known buprestid species. In this study, there were 21 gaps in three mitogenome sequences, which varied from 1 bp to 57 bp. The longest intergenic spacer (bp) was located between trnD and atp8 genes in C. japonica. There were 41 overlapping gene regions in total, ranging from 1 bp to 27 bp in length.

In Julodinae, the concatenated length of 13 PCGs of J. variolaris (Julodinae) was 11,170 bp, which encoded 3715 amino acid residues. In P. chinensis (Polycestinae), the total length of 13 PCGs was 11,162 bp, which encoded 3710 amino acid residues. In C. japonica (Chrysochroinae), the total length of 13 PCGs was 11,161 bp, which encoded 3710 amino acid residues. Compared with the other known buprestid species (Chen et al. 2021; Peng et al. 2021; Huang et al. 2022; Wei 2022; Weng et al. 2022), the concatenated length of 13 PCGs and the number of amino acid-coding codons of Julodinae is slightly higher than in other subfamilies.

The majority of PCGs directly used ATN as the start codon, but the exceptions were nad1 (J. variolaris, P. chinensis, and C. japonica), nad4L (C. japonica), and nad5 (C. japonica) genes which started with TTG, GTG, and GTG, respectively. The unusual start codon TTG was also reported in Agrilinae (Wei 2022) and Buprestinae (Huang et al. 2022). The start codon of the cox1 gene in these three mitogenomes was not determined, which may use non-canonical start codons (Friedrich and Muquim 2003; Fenn et al. 2007; Yang et al. 2013; Wang et al. 2021; Wu et al. 2022). There were three types of stop codons, TAA, TAG, and an incomplete stop codon T, which was completed by the addition of 3’ A residues to the mRNA.

To investigate further, the frequency of synonymous codon usage and relative synonymous codon usage (RSCU) values were calculated and presented. Taken together, the three most frequently used amino acids were L2, S2, and P (Fig. 1A, B), and the most frequently used codons were TTA (L2), TCT (S2), and CCT (P) (Fig. 2).

Figure 1. 

Numbers of different amino acids in the three new mitogenome sequences A and the percentages of the top ten amino acids B the stop codon is not included in these graphs.

Figure 2. 

Relative synonymous codon usage (RSCU) of the three newly sequenced mitogenomes.

The Ka/Ks ratio can be used to estimate whether a sequence is undergoing negative, neutral, or positive selection (Hurst 2002; Mori and Matsunami 2018). The ratio of Ka/Ks for each mitogenome sequence was calculated using Anthaxia chinensis Kerremans, 1898 as the reference sequence (Fig. 3A). In three mitogenome sequences, values of Ka, Ks, and Ka/Ks were all less than 1, suggesting the presence of purifying selection in these three species.

Figure 3. 

Evolutionary rates of mitochondrial genomes in three new mitogenome sequence (A) and the heterogeneity of two dataset in Buprestidae (B).

The rRNA genes were located between the A + T-rich region and trnL1, and separated by trnV, which is consistent with previous studies (Duan et al. 2017; Cao and Wang 2019a, b; Xiao et al. 2019; Sun et al. 2020; Chen et al. 2021; Peng et al. 2021; Huang et al. 2022; Wei 2022; Weng et al. 2022). The total length of rRNA genes ranged from 2035 bp (C. japonica) to 2093 bp (J. variolaris), of which the length of 16S gene ranged from 1299 bp (C. japonica and P. chinensis) to 1301 bp (J. variolaris). The A + T content of rRNA genes ranged from 71.50% (C. japonica) to 74.30% (J. variolaris).

The concatenated lengths of all tRNA genes ranged from 1437 bp (C. japonica) to 1456 bp (J. variolaris), whereas individual tRNA genes ranged from 61 bp (trnR) to 71 bp (trnK), of which eight tRNA genes were encoded on the N-strand and the remaining 14 genes encoded on the J-strand. The predicted secondary structure of tRNAs showed a standard clover-leaf structure (Suppl. material 1: figs S2–S4), except for trnS1 (Fig. 4A), which lacked the dihydrouridine arm, and formed a loop commonly found in other insects (Xiao et al. 2011; Park et al. 2012; Yu et al. 2016; Yan et al. 2017; Yu and Liang 2018; Li et al. 2019). The UG mismatches were detected in some tRNAs (Suppl. material 1: figs S2–S4), which also appeared in other buprestid species (Sun et al. 2020; Chen et al. 2021; Huang et al. 2022; Wei 2022; Weng et al. 2022).

Figure 4. 

The predicted secondary cloverleaf structure for the trnS1 of three new mitogenomes (A) and the gene order of known buprestid mitogenomes (B).

The degree of heterogeneity of the PCGs + RNAs dataset was higher than that of the PCGs dataset (Fig. 3B). Additionally, the heterogeneity in sequence divergences was slightly stronger for Agrilinae than for other families (Fig. 3B). The heterogeneity in sequence homogeneity was higher for Chrysochroinae than other families.

The A + T-rich region was the largest non-coding region in mitogenome, located between trnI and rrnS. This region, containing regulatory elements correlated with the regulation of replication and transcription (Zhang et al. 1995), plays a very important role in molecular evolution (Zhang and Hewitt 1997). The length of A + T-rich region ranged from 1080 bp (C. japonica) to 1510 bp (J. variolaris), which are of medium length in the Buprestidae (Sun et al. 2020; Huang et al. 2022; Wei 2022). The A + T content of the A + T-rich region of C. japonica (75.93%) and P. chinensis (78.38%) was found to be higher than that of the whole genome (67.97%, 67.00%), PCGs (66.46%, 64.55%), rRNAs (71.50%, 72.51%), and tRNAs (68.82%, 71.46%), whereas the A + T content of J. variolaris (72.85%) was lower than that of whole genome (70.43%), PCGs (68.82%), rRNAs (74.30%), and tRNAs (74.79%).

The tandem repeat regions of three species were detected in this study. The repeat regions in each of the three new mitogenomes differ from each other in length and copy number of tandem repeat units. The repeat region of J. variolaris was 43 bp in length, comprising a 17 bp and a 26 bp tandem repeat element. In contrast, in P. chinensis, the total length of the repeat sequence was 111 bp, consisting of three incomplete repeat units. These tandem repeat elements are slightly shorter than those of Agrilinae (Wei 2022).

The gene rearrangements were regarded as important molecular markers for exploring the evolution and phylogeny of insects (Dowton et al. 2002; Cameron 2014). All the buprestid mitogenomes released in GenBank were compared and analyzed, with one mitogenome arrangement pattern exhibited in Buprestidae (Fig. 4B). The mitochondrial gene order of these three species was consistent with other known buprestid mitogenomes.

For the concatenated sequences, the test of substitution saturation showed that the value of I_ss = 0.3910 was significantly smaller than I_ss.c = 0.8537 and p (0.0000) < 0.01, suggesting the sequences suitable for phylogenetic analysis. In the present study, both ML and BI trees using a nucleotide matrix (13 PCGs + 2 rRNAs) produced identical topologies (Fig. 5, Suppl. material 1: fig. S5), (Chrysochroniae + ((Julodinae + Polycestinae) + Buprestinae) + Agrilinae), in terms of subfamily-level relationship.

Figure 5. 

Phylogenetic relationships of studied species of Buprestidae using BI analyses based on 13 PCGs + 2 rRNAs of mitogenomes. The numbers on the branches show posterior probabilities.

The target species J. variolaris, representing Julodinae, formed an independent clade close to Polycestinae with high support values (BI: 1; ML: 94), which supported the results of a previous study (Evans et al. 2015). The target species P. chinensis and Acmaeodera sp. are grouped together as an independent clade with high support values (BI: 1; ML: 100), representing Polycestinae. The Julodinae and Polycestinae formed a clade which was sister to Buprestinae with high support values (BI: 1; ML: 84). The target species C. japonica was clustered with other chrysochroine species as a clade, representing Chrysochroinae, with high support values (BI: 1; ML: 100). All the species of Agrilinae were clustered on one branch with high support values (BI: 1; ML: 100) and close to other buprestid clades, while the Coraebini was polyphyletic.