A Mitochondrial Genome Phylogeny of Cleridae (Coleoptera, Cleroidea)

Simple Summary The family Cleridae is a cosmopolitan group with approximately 4000 species and 320 genera. Within the family, the phylogenetic relationships among the subfamilies, and the timing of divergence, remain not yet fully resolved. Mitochondrial genomes have been widely used to reconstruct phylogenies of various insect groups, but never introduced to Cleridae until now. In this study, we newly generated 18 complete or nearly complete mitochondrial genomes, which are conserved in the organization and structure, as well as exhibit typical high A+T-bias and a preference of nucleotides A and G over T and C, as other insects. Further based on these sequences, a phylogeny of this family is reconstructed of different datasets by both maximum likelihood (ML) and Bayesian inference (BI) methods. The results are congruent and support the monophylies of the family and each subfamily, and the subfamilial relationships are recovered as Korynetinae + (Tillinae + (Clerinae + Hydnocerinae)). Moreover, a molecular clock analysis estimated the divergence time of Korynetinae from others no later than 160.18Mya (95% HPD: 158.18–162.06Mya). The current study presents the first mitochondrial genome-based phylogeny of Cleridae, which provides new evidence in reconstructing the phylogenetic relationships among the subfamilies and understanding the mitochondrial features of this family. Abstract The predaceous beetle family Cleridae includes a large and widely distributed rapid radiation, which is vital for the ecosystem. Despite its important role, a number of problems remain to be solved regarding the phylogenetic inter-relationships, the timing of divergence, and the mitochondrial biology. Mitochondrial genomes have been widely used to reconstruct phylogenies of various insect groups, but never introduced to Cleridae until now. Here, we generated 18 mitochondrial genomes to address these issues, which are all novel to the family. In addition to phylogenomic analysis, we have leveraged our new sources to study the mitochondrial biology in terms of nucleotide composition, codon usage and substitutional rate, to understand how these vital cellular components may have contributed to the divergence of the Cleridae. Our results recovered Korynetinae sister to the remaining clerids, and the calde of Clerinae+Hydnocerinae is indicated more related to Tillinae. A time-calibrated phylogeny estimated the earliest divergence time of Cleridae was soon after the origin of the family, not later than 160.18 Mya (95% HPD: 158.18–162.07 Mya) during the mid-Jurassic. This is the first mitochondrial genome-based phylogenetic study of the Cleridae that covers nearly all subfamily members, which provides an alternative evidence for reconstructing the phylogenetic relationships.

Most checkered beetles occur on plants and tree trunks, both larvae and adults are evidently predators on other insects, especially wood-infesting beetles and their larvae [4].
Although the taxonomic work can be traced to Linnaeus [5], the systematics of Cleridae has been in a hot debate in the history. Many scholars have made their contributions to this issue, including Spinola [6,7], Lacordaire [8], Schenkling [9,10], Gahan [11], Chapin [12] and Böving and Craighead [13]. The eight subfamilies system proposed by Crowson [14] has been widely accepted until recently, which is largely based on the classification of the latter two works. Later, Winkler [15][16][17] established another two subfamilies and proposed a system of neutral terms as an interim measure to show relatedness among higher taxa, but none of his concepts were followed. The morphology-based classification of Kolibáč [18,19] reduced the number of the subfamilies into four by employing Transformation Series Analysis method, including Tillinae, Hydnocerinae, Clerinae and Korynetinae, which was followed by Leschen [20]. On the contrary, Opitz [1,21] erected another three subfamilies and recognized a total of 12 subfamilies, but no one agreed with his concept. More recently, a molecular phylogeny based on four gene markers was constructed by Gunter et al. [3], who proposed another subfamily, Epiclininae (separated from Clerinae), and meanwhile suggested Thaneroclerinae be reassigned to Cleridae. However, the latter was retained as a separate family by Gimmel et al. [22] as conducted by Kolibáč [19], on basis of a four-gene based phylogeny. In a most recent work [23], Hydnocerinae was downgraded to be a tribe of Clerinae.
On a more macro-scale, some groups within Cleridae are also interesting for their contribution to the fossil records, and debate continues over when this group first appeared. Even the timing of the emergence of clerid beetles is still debate, with the oldest unambiguous fossil appeared in the Middle Jurassic. A robust timing of the emergence of clerids via phylogeny would be an independent way to test any of these hypotheses, which has profound consequences for our interpretation of the origin and evolution of Cleridae. A number of previous attempts have been made to combine fossil data with molecular evidence in the Coleoptera and gain reliable estimates for the emergence of major clades in this order [28][29][30][31]. Particularly, a detailed time-scaled phylogeny of the superfamily Cleroidea was inferred by Kolibáč [25] most recently, and the origin and divergence time of Cleridae was dated. However, the divergence-time of the clerid subfamilies needs to be reassessed, due to inconsistent results in different studies [25,[28][29][30][31].
The paucity of phylogenomic-scale datasets available to address phylogeny prevents large-scale analysis of the Cleridae. Analysis using complete mitochondrial data provides an alternative means of addressing this issue. Mitochondrial genomes contain a range of useful information for phylogenetic investigations, as well as for understanding the basic biology of the organisms that contain them. Their mix of well-conserved and more variable regions render them perfect for understanding the inter-relationships of organism at a range of scales. Furthermore, the crucial roles of mitochondria in providing energy to the cell often means that changes at the genetic level are influenced by the environments inhabited by their species [32]. This direct reflection of evolutionary pressures mean that mitochondrial genome sequences provide a unique window into the biology of the species.
The mitochondrial genome (mitogenome) has been widely used to construct the phylogeny of a number of insect groups [33][34][35][36][37][38][39][40]; however, no studies were carried out on the Cleridae due to a lack of mitochondrial resources. To date, only one complete mitochondrial genome is available in the GenBank [41]. Alongside raw sequence data for phylogenetic comparison, complete mitochondrial genome also present a range of data for other investigations. The mitochondrial biological features, such as compositional bias, codon usage, and substitutional rate variation in mitochondrial genomes provide critical information for molecular evolution [42][43][44][45]. Moreover, phylogenomic analysis with all the 37 mitochondrial genes included has result in improved nodal confidence as compared with single-or multi-locus phylogenetics [33,43,46] but much remains left to be to be tested using such a dataset.
In the present study, we sequenced mitochondrial genomes from 18 species of Cleridae. We used these novel genomes to investigate both the phylogenetic relationships within Cleridae and general biological features of this group. The results will improve our understanding of the phylogeny, divergence time and mitochondrial biology within Cleridae.

Sampling and DNA Extraction
Adult specimens of 18 clerid species (Table 1) were collected and preserved in 100% ethanol at −20 • C before molecular experiments. The specimens were identified following Yang [47]. Total genomic DNAs were extracted using a DNeasy Blood& Tissue kit (QIAGEN, Beijing, China), according to the manufacturer's instructions. DNAs were stored at −20 • C for long-term storage and further molecular analyses, each species of which was attached with a voucher number and deposited in the Museum of Hebei University, Baoding, China (MHBU).

DNA Sequencing and Assembly
Whole mitochondrial genome sequencing was performed using an Illumina Novaseq 6000 platform (Illumina, Alameda, CA, USA) with 150 bp paired end reads at BerryGenomics, Beijing, China. The sequence reads were first filtered following Zhou et al. [48] and then the remaining high-quality reads were assembled using IDBA-UD [46] under a 98% similarity threshold and k values of a minimum of 40 and a maximum of 160 bp. The gene cox1 was amplified by polymerase chain reaction (PCR) using universal primers as "reference sequences" to target mitochondrial scaffolds by IDBA-UD [49] to acquire the best-fit, which is under at least 98% similarity. The PCR cycling conditions comprised a predenaturation at 94 • C for 5 min and 35 cycles of denaturation at 94 • C for 50 s, annealing at 48 • C for 45 s and elongation at 72 • C for 8 min at the end of all cycles. Geneious 2019.2 [50] software was used to manually map the clean readings on the obtained mitochondrial scaffolds to check the accuracy of the assembly.

Sequence Annotation and Analyses
Gene annotation was carried out by Geneious 2019.2 [50] software and the MITOS2 webserver (Available at http://mitos2.bioinf.uni-leipzig.de/index.py (accessed on 1 Feb. 2021) [51]. The circular map of the mitochondrial genome was produced using a visualization tool OrganellarGenomeDRAW (http://ogdraw.mpimp-golm.mpg.de/index. shtml (accessed on 3 April 2021) [52]. Base composition, component skew, codon usage, and relative synonymous codon usage (RSCU) were analyzed by PhyloSuitev1.2.2 [53]. DnaSPv5.10.01 [54] was used to estimate the nucleotide diversity (Pi) in a sliding window analysis (a sliding window of 200 bp and a step size of 20 bp) and non-synonymous (Ka)/synonymous (Ks) substitution rates among the 13 protein-coding genes (PCGs). The genetic distances were computed using MEGA 7.0 [55] with the Kimura-2-parameter model. Substitution saturation of each codon position of PCGs was measured based on Xia's test [56] implemented in DAMBE program v6.4.81 [57]. SymTest v2.0.47 with Bowker's matching pair symmetry test [58] was used to analyze the differences of heterogeneous sequences in the datasets, and the heat maps were generated according to the inferred p-values. The tandem repeats finder program (http://tandem.bu.edu/trf/trf.html (accessed on 1 April 2021) was used to predict the tandem repeat elements in the A+T-rich region [59].

Phylogenetic Analyses
Mitochondrial genomes of 18 species representing all four subfamilies of Cleridae were selected, which are all novel in this family ( Table 1). The one sequence available in the GenBank was not included in the analyses due to possible misidentification [41]. Two species of Dasytinae in Melyridae sensu lato were chosen as the outgroups [60,61]. To test the impact of the third codon position of the PCGs and gene combination types on the phylogenetic analysis, four datasets were concatenated: (i) the first and second codon positions of 13 PCGs (PCG12); (ii) all three codon positions of PCGs (PCG); (iii) all PCGs and rRNAs (PCGrRNA); (iv) all PCGs, rRNAs and tRNAs (PCGRNA).
Alignment of PCGs, tRNAs and rRNAs was performed by using Mafftv7.313 [62] in PhyloSuite v1.2.2 (alignment strategy: auto) [53]. Intergeneric gaps and ambiguous sites were removed using Gblocks v 0.91b [63], and individual alignments were concatenated using PhyloSuite. All matrices were analyzed using maximum likelihood (ML) with IQ-TREE v1.6.12 [64] [65] in PhyloSuite, respectively. A 1000 replicate bootstrapping was performed for ML analyses using the "ultrafast" option [66] implemented in IQ-TREE, with the SH-alerts test used to assess branch supporting values. The best model (Tables S1 and S2) was inferred by Partition-Finder (v2.1.1) [67]. Four simultaneous Markov chain Monte Carlo (MCMC) runs of 1 million generations twice, with trees sampled every one thousand generations, and the first 25% of 1 million generations twice, with trees sampled every one thousand generations, and the first 25% of steps were discarded as burn-in. Stationarity was considered to be reached when the average standard deviation of split frequencies was below 0.01. Trees produced from all analyses were visualized and edited using iTOL (https://itol.embl.de (accessed on 1 July 2021) [68].

Divergence Time Estimate
Divergence times among subfamilies were estimated using the nucleotide sequences of 13 PCGs with a relaxed clock log normal model in BEAST1.10.4 [69,70]. We adopted the Calibrated Yule model for the prior tree, and used the GTR+I+G for concatenation by Phylosuite v1.2.2. For estimating divergence time calibration, Protoclerus korynetoides, the oldest reported fossil of Cleridae from the Middle Jurassic in NE China (mean value of normal prior distribution c.160.2 Mya, SD = 1.0) [71] was used to assign age calibration. The final Markov chain was run twice for each 1 × 10 8 generations, sampling every 10,000 generations with the first 25% of generations discarded as burn-in, after confirming the convergence of chains with Tracer v.1.7.2 [72]. The effective sample size of the majority of parameters was >200. We summarized the subsamples trees in a maximum clade credibility tree with mean heights using Tree Annotator v1.10.4, and then the mean heights and 95% highest probability density (95%HPD) were displayed in Figtree v1.4.3 [73].

Phylogenetic Analyses
Heterogeneity of nucleotide divergence was examined under pairwise comparisons in a multiple sequence alignment. The heterogeneous sequence divergence of PCG12 dataset (Figure 1a) is much lower than that of the other three datasets (Figure 1b-d), indicating that the third codon positions are more rate-heterogeneous than the first and second ones. In all phylogenetic analyses, the monophyly of the family and each subfamily is well supported (PPs = 1, BSs = 100) (Figures 2 and S3). Korynetinae is recovered next to the remaining clerids with high support value (PPs = 1, BSs = 100). Hydnocerinae and Clerinae are recovered as sister groups, which is greatly supported (PPs = 1, BSs = 100). The clade of Hydnocerinae + Clerinae is sister to Tillinae, which is highly supported (PPs = 1, BSs = 100).  In all phylogenetic analyses, the monophyly of the family and each subfamily is well supported (PPs = 1, BSs = 100) (Figures 2 and S3). Korynetinae is recovered next to the remaining clerids with high support value (PPs = 1, BSs = 100). Hydnocerinae and Clerinae are recovered as sister groups, which is greatly supported (PPs = 1, BSs = 100). The clade of Hydnocerinae + Clerinae is sister to Tillinae, which is highly supported (PPs = 1, BSs = 100).

Divergence-Time Estimation
The age estimates (average and 95% PHD) of each subfamily based on the topology recovered from BEAST analysis were summarized in Figure 3. The BEAST analysis indicated that the divergence events of Cleridae occurred during the mid-Jurassic and mid-Cretaceous period. The

Divergence-Time Estimation
The age estimates (average and 95% PHD) of each subfamily based on the topology recovered from BEAST analysis were summarized in Figure 3. The BEAST analysis indicated that the divergence events of Cleridae occurred during the mid-Jurassic and mid-Cretaceous period.
The subfamily Korynetinae was divergent from all others soon after the origin of the family during the mid-Jurassic, approximately at 160.18 Mya (95% HPD: 158.18-162.07 Mya). After a period of evolution, the Tillinae split from the rest during the early Cretaceous, at 138.58 Mya (95% HPD: 123.00-147.69 Mya), and finally Clerinae and Hydnocerinae became divided at 109.03 Mya (95% HPD: 92.80-123.02 Mya).

General Features of Mitochondrial Genome
The complete or nearly complete mitochondrial genomes of 18 clerid species were successfully sequenced. It is a double-strand circular molecule, which is made up of 37 genes, including 13 PCGs, 22 tRNA genes, 2 rRNA genes, and an A+T-rich region (or control region) ( Figures S1 and S2), of which 14 genes (8 tRNAs, 4 PCGs and 2 rRNAs) were transcribed on the minority strand (N-strand), whereas the rest (14 tRNAs and 9PCGs) on the majority strand (J-strand).

General Features of Mitochondrial Genome
The complete or nearly complete mitochondrial genomes of 18 clerid species were successfully sequenced. It is a double-strand circular molecule, which is made up of 37 genes, including 13 PCGs, 22 tRNA genes, 2 rRNA genes, and an A+T-rich region (or control region) ( Figures S1 and S2), of which 14 genes (8 tRNAs, 4 PCGs and 2 rRNAs) were transcribed on the minority strand (N-strand), whereas the rest (14 tRNAs and 9PCGs) on the majority strand (J-strand).
The complete mitogenomes of clerids range from 15,638 bp to 17,127 bp in size (Table  S3). The sizes of the control region vary greatly among different species (ranging from 987bp to 2,401 bp), whereas the PCGs, tRNAs, and rRNAs show little variation in length (Table S3).

Nucleotide Composition
In the full genomes, the A+T contents range from 77.1% to 80%, those of the rRNAs range from 81.2% to 83.4%, and those of the control region from 84.5% to 90.3% in Cleridae (Table S3).
Comparison among the PCGs (Figure 4b, Table S4) shows that the average A+T content of cox1 is the lowest (68.5%) in all clerids, followed by cox3 (72.5%) and cytb (73.5%), whereas those of apt8 and nad6 are the highest (85.0% and 84.7%, respectively). Additionally, the A+T contents of the third codon position of PCGs, ranging from 86.04% to 93.78%, are much higher than those of the other two codon positions, which range from 68.19% to 72.42% (Figure 4a, Table S3). Moreover, the nucleotide skew analysis (Figure 4c, d, Table  S5) shows that the AT skews are positive for most species of Cleridae, whereas the GC skews are all negative. The correlations of Cleridae's mitogenomes were calculated between A+T content versus AT skew (y = −0.0076x + 0.5923, R 2 =0.2771), and G+C content  (Table S3). The sizes of the control region vary greatly among different species (ranging from 987bp to 2,401 bp), whereas the PCGs, tRNAs, and rRNAs show little variation in length (Table S3).

Nucleotide Composition
In the full genomes, the A+T contents range from 77.1% to 80%, those of the rRNAs range from 81.2% to 83.4%, and those of the control region from 84.5% to 90.3% in Cleridae (Table S3).
Comparison among the PCGs (Figure 4b, Table S4) shows that the average A+T content of cox1 is the lowest (68.5%) in all clerids, followed by cox3 (72.5%) and cytb (73.5%), whereas those of apt8 and nad6 are the highest (85.0% and 84.7%, respectively). Additionally, the A+T contents of the third codon position of PCGs, ranging from 86.04% to 93.78%, are much higher than those of the other two codon positions, which range from 68.19% to 72.42% (Figure 4a, Table S3). Moreover, the nucleotide skew analysis (Figure 4c,d, Table S5) shows that the AT skews are positive for most species of Cleridae, whereas the GC skews are all negative. The correlations of Cleridae's mitogenomes were calculated between A+T content versus AT skew (y = −0.0076x + 0.5923, R 2 = 0.2771), and G+C content versus GC skew (y = −0.0184x + 0.2618, R 2 = 0.6213), respectively. Both of them showed negative linear correlations, implying that the quantity of A+T becomes more equivalent with the increase in A+T content, whereas G+C show a larger quantity gap with the increase in G+C content.
versus GC skew (y=−0.0184x+0.2618, R 2 =0.6213), respectively. Both of them showed negative linear correlations, implying that the quantity of A+T becomes more equivalent with the increase in A+T content, whereas G+C show a larger quantity gap with the increase in G+C content.

Codon Usage and Evolutionary Rates
In total, six initiation codons (ATA, ATT, ATG, ATC, TTG, GTG) were used in encoding the PCGs of Cleridae, which are terminated with TAG or TAA codons or truncated T codons (Tables S6). Comparison of these codons among PCGs shows that atp6 always

Codon Usage and Evolutionary Rates
In total, six initiation codons (ATA, ATT, ATG, ATC, TTG, GTG) were used in encoding the PCGs of Cleridae, which are terminated with TAG or TAA codons or truncated T codons (Table S6). Comparison of these codons among PCGs shows that atp6 always uses ATA, and cox3 and cytb use ATG as start codons, respectively. For the stop codons, atp6, atp8, and nad6 all use TAA, and nad4 uses T-, respectively. The most frequently used start codons are ATT and ATA, and the stop codon is TAA.
Relative synonymous codon usage (RSCU) shows that all synonymous codons of 22 amino acids are present in Cleridae. Among these codons, UUA-Leu2 and UCU-Ser2 are the first two frequently used codons, followed by CCU-Pro, GCU-Ala, and CGA-Arg (Table  S7). The RSCU values of the PCGs reveal that there is a higher frequency in the usage of AT than that of GC in the third codon positions. In Cleridae (Table S8), Leu2, Ile, and Phe are the most frequently encoded amino acids (over 10%), followed by Met, Ser2, and Asn.
Sliding window analysis was implemented to study the nucleotide diversity of 13 PCGs exhibited in Figure 5b. Among the genes, nad6 (Pi = 0.26) has the highest variability, followed by atp8 and nad2 (both Pi = 0.24); cox2, cytb, and nad4l (all Pi = 0.15) have the lowest variability. This result is roughly congruent with that of the above pairwise genetic distance calculation. atp6, atp8, and nad6 all use TAA, and nad4 uses T-, respectively. The most frequently used start codons are ATT and ATA, and the stop codon is TAA. Relative synonymous codon usage (RSCU) shows that all synonymous codons of 22 amino acids are present in Cleridae. Among these codons, UUA-Leu2 and UCU-Ser2 are the first two frequently used codons, followed by CCU-Pro, GCU-Ala, and CGA-Arg (Table S7).
The RSCU values of the PCGs reveal that there is a higher frequency in the usage of AT than that of GC in the third codon positions. In Cleridae (Table S8), Leu2, Ile, and Phe are the most frequently encoded amino acids (over 10%), followed by Met, Ser2, and Asn.
Sliding window analysis was implemented to study the nucleotide diversity of 13 PCGs exhibited in Figure 5b. Among the genes, nad6 (Pi = 0.26) has the highest variability, followed by atp8 and nad2 (both Pi = 0.24); cox2, cytb, and nad4l (all Pi = 0.15) have the lowest variability. This result is roughly congruent with that of the above pairwise genetic distance calculation.

Phylogeny and Divergence-Time Estimation
Comparing all the topologies produced by different datasets by both ML and BI analyses, the nodal support values were improved based on the PCG12 dataset when the third codon positions were excluded. This may be a result of the high heterogeneity ( Figure 1) and saturated substitution (Figure 2c) of the third codon positions, which is relatively free in evolution [74], and the change of nucleotide in this position rarely change the amino acid product, consistent with their codon usage (Table S7). When the third codon positions are highly heterogeneous or saturated, they should be excluded from the phylogeny reconstruction, generally because they will affect the reliability of the phylogenetic analysis results or be less informative, as suggested by other studies [75][76][77][78].
The monophyly of the family and each subfamily is well supported, which is consistent with some phylogenetic studies [2,19,28]. Korynetinae is recovered next to the remaining clerids with high support value, which is congruent with the morphology-based phylogeny carried out by Kolibáč [19]. Korynetinae was defined by a synapomophy (reduction in size of the fourth tarsomeres) by Kolibáč [19]. However, when the fossil evidence was implemented, Thaneroclerinae was recovered as the basal clade of Cleridae [27]. Otherwise, the molecular phylogenies suggested Tillinae were at the base of the clerid clade [3,23,25,28], regardless of whether Thaneroclerinae was sampled or not. Tillinae is the only clerid subfamily in which the procoxal cavities are closed internally, which is a synapomorphy of this subfamily [1] but may be apomorphic within Cleridae, since the Cleroid family Rentoniidae is also equipped with this feature [3].
Hydnocerinae and Clerinae are recovered as sister groups, which is consistent with the results of some studies [2,28], but they were suggested to be paraphyletic by others [3,23,25,78]. Based on the latter results, Hydnocerinae was synonymized with Clerinae by Bartlett [23]. Although our taxa sampling is too limited to test the monophylies of these two subfamilies, Opitz's [1] findings indicated that Hydnocerinae and Clerinae both possess two secondary stomodaeal valve lobes (four in Tillinae and Thaneroclerinae or completely reduced in Korynetinae), which is synapomorphic for supporting their monophyly or sister relationship.
The clade of Hydnocerinae + Clerinae is sister to Tillinae, which is highly supported and shows some similarity to the phylogeny when the fossil data were implemented [78]. However, Epiclininae, which was separated from Clerinae, together with Korynetinae recovered a sister relationship with this clade [3,25]. In terms of our results, Korynetinae is relatively distant from this clade in the affinity, although no Epiclininae was available for this study.
The estimated divergence time is earlier than those predicted by previous studies [28,29,31], while later than others [25,30]. Protoclerus korynetoides, one of the two earliest clerid fossil was not attributed to any subfamily [24] but was used to assign age calibration prior for the family in the present study, so the estimated time is later than in theory at least. Therefore, the origin of Cleridae is probably from the early Jurassic period, as predicted by Toussaint et al. [30] and Kolibáč et al. [25], and accordingly, the divergence time of the subfamilies was also earlier than that estimated in this study.

The Characteristics of Mitochondrial Genomes
All typical animal mitochondrial genes and control regions were identified in the 18 mitogenomes, and they are arranged in the same order as the hypothetical ancestral insect [79], indicating that the mitogenomes are highly conserved in Cleridae.
The sizes of the whole mitogenome among the species are comparable, in which the length of the control region varies greatly, whereas other components show little variation. This suggests that the mitogenome size of different clerids is largely determined by the size of control regions, such as other insects [80].
As with all other insects [36], the mitochondrial genomes of clerids exhibit the typical high A+T-bias either in the full genome or different components or different positions of PCGs (Figure 2a), with the A+T contents all higher than 68.19% (Table S3). Additionally, Cleridae usually have a preference of A and G over T and C in the mitogenomes. The causes for such base composition bias are multifactorial, but most of the hypotheses suggest that the asymmetric nucleotide composition is the result of mutations and selection pressures [81], and the value of the GC-skew of the insect mitochondrial genomes seems to be associated with the replication orientation [82].
The ratio (ω) of Ka/Ks can be used to estimate whether a sequence is undergoing purifying, neutral, or positive selection [83]. We found that all PCGs are evolving under a purifying selection (ω < 1), and cox2 exhibited the strongest purifying selection, whereas the nad family genes and atp8 relaxed. Furthermore, the results of pairwise genetic distance calculation and the nucleotide diversity analysis are consistent, indicating that cox2, cytb, and nad4l evolve comparatively slowly or have lower variability, whereas atp8, nad2, and nad6 are evolving faster or have higher variability. Nucleotide diversity analyses are critical for designing species-specific markers useful in taxa where morphological identification is difficult and ambiguous [84][85][86]. Usually, cox1 is the last variable and can be a potential marker for species identification and has been widely used in the taxonomic work of insects [87,88]. However, our analysis reveals that cox2 may be more suitable for Cleridae. Given this result, we suggest that the gene markers should be designed for different families or even for different subfamilies in the taxonomy if necessary.

Conclusions
The previous work attempted to resolve the phylogenetic relationship of Cleridae exclusively based on morphological characters or short nucleotide fragments. With this study, we documented the phylogenetic relationships of Cleridae based on the complete mitochondrial genomes. We chose the Chinese species as a start point, because they were the most accessible for us, with 18 species representing four subfamilies providing a diverse, but not unmanageable, number of taxa for analysis. Nevertheless, the molecular phylogenies, including this study, were analyzed on the basis of a minority part of species or limited molecular data in comparison with an estimated 4000 species of checkered beetles worldwide. Therefore, many more species need to be included in future analysis to establish a solid and dependable classification of Cleridae, especially for some taxa (i.e., Thanerocleri-dae/-nae, Epiclininae) whose status is controversial. In particular, the complete mitochondrial genomes should be encouraged to accumulate more for Cleridae, in view of their high value in investigating phylogenetic relationships of the insects.

Supplementary Materials:
The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/insects13020118/s1, Figure S1: Phylogenetic tree of Cleridae inferred from the BI and ML analyses of the PCGRNA (left), PCGrRNA (middle) and PCG (right) datasets, Figure S2: Maps of the nearly complete mitogenomes of four subfamilies of Cleridae. Different color represents different types of genes and regions, Figure S3: Phylogenetic tree of Cleridae inferred from the BI and ML analyses of the PCGRNA (left), PCGrRNA (middle) and PCG (right) datasets. Numbers at the branches are bootstrap values (upper) or posterior probabilities (lower), Table S1: Best models calculated by BI of four mitochondrial datasets, Table S2: Best models calculated by IQ of four mitochondrial datasets, Table S3: The length and A+T% of different components of mitogenomes in Cleridae, Table S4: A+T% of each PCG of mitogenomes in Cleridae, Table S5: A+T% and A+T skew, C+G% and G+C skew in the 13 PCGs of mitogenomens in Cleridae, Table S6: Codon usage of PCGs of mitogenomes in Cleridae; Table S7: RSCU of the mitochondrial PCGs of Cleridae, Table S8: Amino acid constituents among four subfamilies of Cleridae.