﻿First mitochondrial genome of subfamily Julodinae (Coleoptera, Buprestidae) with its phylogenetic implications

﻿Abstract Complete mitochondrial genomes of three species of the family Buprestidae were sequenced, annotated, and analyzed in this study. To explore the mitogenome features of the subfamily Julodinae and verify its phylogenetic position, the complete mitogenome of Julodisvariolaris was sequenced and annotated. The complete mitogenomes of Ptosimachinensis and Chalcophorajaponica were also provided for the phylogenetic analyses within Buprestidae. Compared to the known mitogenomes of Buprestidae species varied from 15,499 bp to 16,771 bp in length, three newly sequenced mitogenomes were medium length (15,759–16,227 bp). These mitogenomes were encoded 37 typical mitochondrial genes. Among the three studied mitogenomes, Leu2 (L2), Ser2 (S2), and Pro (P) were the three most frequently encoded amino acids. Within the Buprestidae, the heterogeneity in sequence divergences of Agrilinae was highest, whereas the sequence homogeneity of Chrysochroinae was highest. Moreover, phylogenetic analyses were performed based on nucleotide matrix (13 PCGs + 2 rRNAs) among the available sequenced species of Buprestidae using Bayesian Inference and Maximum Likelihood methods. The results showed that the Julodinae was closely related to the subfamily Polycestinae. Meanwhile, the genera Melanophila, Dicerca, and Coomaniella were included in Buprestinae, which was inconsistent with the current classification system of Buprestidae. These results could contribute to further studies on genetic diversity and phylogeny of Buprestidae.


Introduction
The family Buprestidae is one of the largest families in Coleoptera, including six subfamilies, 521 genera, and more than 15,000 species distributed worldwide (Bellamy 2008;Kubáň et al. 2016). In this family, all species are phytophagous. The adults are feeders on flowers, leaves and stems, whereas the larvae are internal feeders in roots and stems, or feed on the foliage of woody and herbaceous plants, the larvae of Julodinae are soil habitants feeding externally by the roots (Bellamy and Volkovitsh 2016). Different groups have different functions covered ecological, social and economic functions, such as: most larvae of Buprestinae and Chrysochroinae are important decomposers of woody plants; with most species being ornamental beetles with attractive metallic luster; many species of Agrilinae are forest and agricultural pests; and some species of the tribes Stigmoderini, Acmaeoderini, and Anthaxiini are pollinator taxa. Although some buprestid taxonomists have made important contributions to the classification based on morphological analyses (Cobos 1980(Cobos , 1986Tôyama 1987;Hołyński 1988Hołyński , 1993Hołyński , 2009Kolibáè 2000;Bellamy 2003), the problems of the overall classification of Buprestidae remain.
In the past two decades, the mitochondrial genome emerged as important molecular data for higher-level phylogenetic analyses (Saccone et al. 1999;Timmermans et al. 2010Timmermans et al. , 2016Cameron 2014;Li et al. 2015;Qin et al. 2015;Nie et al. 2020Nie et al. , 2021Motyka et al. 2022;Zheng et al. 2022), evolutionary strategies (Krzywinski et al. 2011;Nie et al. 2019;Motyka et al. 2022;Zhang et al. 2022), and genetic diversity analyses (Lim et al. 2021). The buprestid mitogenome also caught the attention of taxonomists. In Buprestidae, the first complete mitogenome of Chrysochroa fulgidissima (Schönherr, 1817) was reported by Hong et al. (2009). In the same year, the mitogenome of Acmaeodera sp. was used to analyze the nonstationary evolution and compositional heterogeneity of Coleoptera. To date, only 22 buprestid mitogenomes (Table 1) have been reported worldwide, including three newly generated in this study.
To date, the mitogenome of the subfamily Julodinae has not been reported. The lack of the data on complete mitogenome of Julodinae species has limited our understanding of the real phylogenetic relationships within jewel beetles. The single molecular phylogenetic analysis, including Julodinae, showed that Julodinae is monophyletic group and close to Polycestinae (Evans et al. 2015). The subfamily Julodinae includes one tribe and six genera (Hołyński 2014). The described Julodinae species are mainly distributed in the arid and semiarid zones of the Ethiopian and Palaearctic regions, except for the species of the genus Sternocera Eschscholtz, 1829 distributed in humid tropical zones of Asia and Africa (Bellamy 2008;Hołyński 2014).
In the present study, three complete mitogenomes are sequenced and annotated, of which that of Julodis variolaris (Pallas, 1771) is the first complete mitogenome sequence to be reported in the subfamily Julodinae. In China, this species is widely distributed in Xinjiang Uygur Autonomous Region. The adults, appearing in May and June, feeder on the leaves of Haloxylon ammodendron (Meyer, 1829) and the larvae feeder on the roots of this plant. Additionally, the complete mitogenomes of Chalcophora japonica (Gory, 1840) (Chrysochroinae: Chalcophorini) and Ptosima The total length of the mitogenome in C. japonica was consistent with the results of Weng et al. (2022). In order to explore the phylogenetic position of the subfamily Julodinae, phylogenetic analyses of the family Buprestidae were performed based on a nucleotide matrix (13 PCGs + 2 rRNAs) among buprestid species using Bayesian Inference (BI) and Maximum Likelihood (ML) methods.

Sequence assembly, annotation, and analysis
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).

Phylogenetic analysis
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.  (Zhang et al. 2020). During this analyzing process, PhyloSuite was run with previous parameters (Wei 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.

Protein-coding genes, codon usage, and nucleotide diversity
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 . 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).   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.

Ribosomal and transfer RNA genes, and heterogeneity
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).
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 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.

Phylogenetic analysis
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.
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.

Discussion
The gene composition and arrangement of these three mitogenomes are the same as other known buprestid mitogenomes (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 had a positive AT skew, which was similar to most known buprestid mitogenomes (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 genes nad1 (J. variolaris, P. chinensis, and C. japonica), nad4L, and nad5 (C. japonica) which started with TTG, GTG, and GTG, respectively, was also reported by previous studies in Buprestidae Wei 2022). The Julodinae are closest to Polycestinae with high support values, which is consistent with the results of a previous study (Evans et al. 2015). The monophyly of Buprestidae has been corroborated once more, as all the buprestid species converge together as an independent clade (Evans et al. 2015;Huang et al. 2022;Wei 2022). In this study, the Coraebini was also found to be polyphyletic with the genera Meliboeus Deyrolle, 1864 and Coraebus Gory & Laporte, 1839 in different clades, also consistent with the previous studies (Evans et al. 2015;Huang et al. 2022;Wei 2022). Compared to Melanophilini, Coomaniellini is more closely related to Dicercini, which is in line with previous studies (Volkovitsh 2001;Evans et al. 2015;Huang et al. 2022).
In the present study, the sampling might be too limited to address the comprehensive phylogeny of Buprestidae. In the future, classification problems could be solved when enough mitogenomes are accumulated for more buprestid species, which requires the cooperation of taxonomists around the world.

Conclusions
In this study, the complete mitogenomes of Julodis variolaris, Chalcophora japonica, and Ptosima chinensis were annotated and analyzed, of which the mitogenome of J. variolaris was the first complete mitogenome representative of the subfamily Julodinae. The three mitogenome sequences were of medium length (15,759-16,227 bp) in Buprestidae. These three mitogenomes shared the same gene order, which was consistent with those of known buprestid species. These three mitogenome sequences all had a high A + T content, and strong AT bias. All PCGs of the three species began with the typical ATN codon except nad1 (J. variolaris, P. chinensis, and C. japonica), nad4L (C. japonica), and nad5 (C. japonica) which were initiated with TTG, GTG, and GTG, respectively. In the present study, the BI and ML trees had exact same topologies with high-value support. The results of phylogenetic analyses also show that Julodinae is close to Polycestinae, the clade composed of Julodinae and Polycestinae is close to that of Buprestinae, and the Agrilinae clade is sister to that of (Chrysochroniae + ((Julodinae + Polycestinae) + Buprestinae)), and all the subfamilies are grouped in a monophyletic group with high support.