Comparative Analysis of Eight Mitogenomes of Bark Beetles and Their Phylogenetic Implications

Simple Summary Many bark beetles are destructive pests in coniferous forests and cause extensive ecological and economic losses worldwide. Comparative studies of the structural characteristics of mitogenomes and phylogenetic relationships of bark beetles can improve our understanding of mitogenome evolution. In this study, we sequenced eight mitogenomes of bark beetles. Our results show that the use of start and stop codons, the abundance of amino acids, and the relative frequency of codon use are conserved among the eight bark beetles. Different regions of tRNA exhibit different degrees of conservatism. Together with the analysis of evolutionary rates and genetic distance among bark beetle species, our results reveal phylogenetic relationships among bark beetles of the subfamily Scolytinae. Abstract Many bark beetles of the subfamily Scolytinae are the most economically important insect pests of coniferous forests worldwide. In this study, we sequenced the mitochondrial genomes of eight bark beetle species, including Dendroctonus micans, Orthotomicus erosus, Polygraphus poligraphus, Dryocoetes hectographus, Ips nitidus, Ips typographus, Ips subelongatus, and Ips hauseri, to examine their structural characteristics and determine their phylogenetic relationships. We also used previously published mitochondrial genome sequence data from other Scolytinae species to identify and localize the eight species studied within the bark beetle phylogeny. Their gene arrangement matched the presumed ancestral pattern of these bark beetles. Start and stop codon usage, amino acid abundance, and the relative codon usage frequencies were conserved among bark beetles. Genetic distances between species ranged from 0.037 to 0.418, and evolutionary rates of protein-coding genes ranged from 0.07 for COI to 0.69 for ND2. Our results shed light on the phylogenetic relationships and taxonomic status of several bark beetles in the subfamily Scolytinae and highlight the need for further sequencing analyses and taxonomic revisions in additional bark beetle species.


Introduction
Bark beetles (Coleoptera: Curculionidae: Scolytinae) feed mainly on the phloem and xylem of their hosts (e.g., pines, spruces, larches), severely affecting forest ecology and reducing forestry production [1,2]. Because they bore cryptically under the bark and into the wood, they are difficult to detect at the onset of the threat, resulting in many bark beetle species becoming established outside their original range [3]. In addition, international trade has led to the increasing establishment of exotics, including many invasive species from the subfamily Scolytinae [4]. A recent example is Gnathotrichus materiarius (Fitch, 1868) and Cyclorhipidion bodoanum (Reitter, 1913), which have been newly recorded Insects 2021, 12, 949 2 of 12 as alien species for the British Isles [5]. Because bark beetles are both highly damaging and difficult to control, many researchers have studied the developmental trends of bark beetle species/populations using molecular markers [6][7][8][9][10]. The mitochondrial genome of insects, which is characterized by low molecular weight, simple structure, maternal inheritance, and rapid evolutionary rate [11][12][13][14][15], is widely used to study insect phylogeny and population inheritance [16][17][18]. With the development of DNA sequencing technology, the mitochondrial genomes of more and more insect species are being sequenced [19][20][21]. Nevertheless, data on the mitochondrial genomes of bark beetles remain limited. Cognato et al. [22] speculated that bark beetles in China may have high genetic diversity. In addition, taxonomic identification of some bark beetles is very important due to their pest status. Although thousands of species in the subfamily Scolytinae are not known to have economic impacts, they are ecologically important because they are the first decomposers of woody materials [23].
In this study, we sequenced and annotated the mitochondrial genomes of eight bark beetles, including Dendroctonus micans (Kugelann, 1794), Orthotomicus erosus (Wollaston, 1857), Polygraphus poligraphus (Linnaeus, 1758), Dryocoetes hectographus (Reitter, 1913), Ips nitidus (Eggers, 1933), Ips typographus (Linnaeus, 1758), Ips subelongatus (Motschulsky, 1860), and Ips hauseri (Reitter, 1894). By comparing these mitochondrial genomes, we can not only understand their structural characteristics, but also explore events in the evolutionary process of bark beetles that involve rearrangements and variations in mitochondrial genes. In addition, we determined the phylogenetic relationships between the eight species analyzed here and 18 other bark beetle species and 22 undetermined species from the subfamily Scolytinae whose mitochondrial genomes were sequenced. Our analysis provides new data for studying the phylogenetic relationships of some important bark beetles.

Sampling and DNA Extraction
Eight species of bark beetles, including Dendroctonus micans, O. erosus, P. poligraphus, Dryocoetes hectographus, I. nitidus, I. typographus, I. subelongatus, and I. hauseri, were collected in the field. Dr. Fang used a taxonomic search [24] to identify the species based on their body characteristics (e.g., teeth, body size, frontal tubercles, number of frontal hairs) along with the host tree species, and the detailed information is provided in Table 1. All collected specimens were stored in the Insect Museum of Chinese Academy of Forestry (Curator: Xiangbo Kong, Beijing, China). Live adult specimens were preserved in anhydrous ethanol and stored in a freezer at −20 • C. Total genomic DNA was extracted from pronotum and leg muscle tissue (five individuals per species) using a Wizard ® Genomic DNA Purification Kit (Promega Corporation, Madison, WI, USA) according to the manufacturer's protocol. In conjunction with the NanoDrop 2000 Fluorospectrometer (Thermo Fisher Scientific, Wilmington, DE, USA), the PicoGreen ® assay (Thermo Fisher Scientific, Wilmington, DE, USA) was used to quantify dsDNA, and dsDNA integrity was determined by 1% agarose gel electrophoresis.

Mitogenome De Novo Sequencing and Assembly
Genomic DNA was fragmented to 400-500 bp using a Covaris ® M220 focused ultrasonicator (Covaris, Woburn, MA, USA). A 460 bp paired-end library was generated from each sample separately and finally sequenced using an Illumina HiSeq X Ten platform (Majorbio Bio-pharm Technology, Shanghai, China) to obtain 4 Gb of data. Reads with low sequencing quality were filtered out using Trimmomatic software (v0.36, http: //www.usadellab.org/cms/?page=trimmomatic). Sequences were removed if they: contained splice fragments, had low mass values (Q < 25), or had lengths of low mass bases greater than half the total sequence length. The software Spades (Center for Algorithmic Biotechnology, St. Petersburg, Russia) [25] was used to assemble clean reads, and the assembly results were compared by mapping to obtain the longest segment. The MitoZ program (BGI-Shenzhen, Shenzhen, China) [26] was used to annotate the assembled genome, and the coding sequences (CDS), transfer RNA (tRNA), and ribosomal RNA (rRNA) were identified. The mitochondrial genome sequences of all eight species were deposited at GenBank (see Table 1 for deposit numbers).

Comparative Mitochontrial Genome Analysis
The nucleotide composition of the mitochondrial genome, including the content of A, T, C, and G bases, AT and GC skew, amino acid usage, and relative synonymous codon usage (RSCU) was calculated using MEGA 5.2.2 software (MEGA Software, Paris, France) [27]. Compositional skew was calculated using the formulas AT skew = (A − T)/(A + T) and GC skew = (G − C)/(G + C) [28]. tRNA genes were determined using the tRNAscan-SE Search Server (http://lowelab.ucsc.edu/tRNAscan-SE) (v2.0) and the MITOS [29] web server with invertebrate mitochondrial genetic codes [30].

Genetic Distance and Selection Pressure Analysis
To measure genetic distances among bark beetles of the subfamily Scolytinae, we used 26 mitogenomes, 8 of which were sequenced in our study and 18 of which were from the NCBI database to perform the corresponding analyses. Sequences were aligned using the Clustal X program and pairwise genetic distances were calculated using the MEGA v5.2.2 (MEGA Software, Paris, France) program based on Kimura 2-parameter model [31,32]. The software ClustalX was used to align the nucleic acid sequences of 13 protein-coding genes (PCGs). The nonsynonymous substitution rate (Ka), synonymous substitution rate (Ks), and Ka/Ks of the PCGs were determined using DnaSP v6.12.03 software [33].

Phylogenetic Analysis
To illustrate the phylogenetic relationships of the eight species of the subfamily Scolytinae in a broader evolutionary context, we constructed phylogenetic trees by analyzing the 50 mitochondrial genomes of the species of the subfamily Scolytinae (Supplementary  Table S1), using Sitophilus zeamais (Motschulsky, 1855) and Sitophilus oryzae (Linnaeus, 1763) as outgroups [34]. PCGs, rRNA, and tRNA were aligned and concatenated in each mitochondrial genome using standard ClustalW parameters. Phylogenetic trees were constructed using the maximum likelihood method (ML) and Bayesian inference method (BI). The software MEGA v5.2 (MEGA software, Paris, France) was used to analyze the alternative models (GTR+G) and construct the tree ML. The confidence values of each branch node of the phylogenetic tree were bootstrapped with 1000 replicates [35]. The software MrBayes 3.2.2 [36] was used to construct the tree BI and MrMTgui software [37] to compare and analyze the alternative models of the sequences, and the optimal alternative model GTR+G was used to construct the evolutionary tree. The MCMC method was used to simulate 10,000,000 generations, sampling once every 1000 generations to ensure sampling independence. The first 25% of the simulations were discarded as burn-in. Stationarity was considered achieved when the average standard deviation of the split frequencies was less than 0.01 [38].

Nucleotide Composition and Condon Use
The AT content of the whole mitochondrial genome was above 70% in all eight species studied here, with Dendroctonus micans having the highest value (75.72%) and O. erosus the lowest (70.69%). In the sequences of PCGs, tRNA, and rRNA, the AT content was also significantly higher than that of CG, with the AT content being highest in rRNA and lowest in PCGs (Supplementary Table S10). All eight Scolytinae mitogenomes showed a positive AT skew (0.335 to 0.414) and a negative GC skew (−0.2722 to −0.2421).
Except for the ND1 gene in Ips bark beetle which used TTG, PCGs of all species used ATN as the start codon (Supplementary Tables S2-S9). The use of stop codons also varied widely within and between species. In this study, we found two types of stop codons: the complete stop codons TAA and TAG and the incomplete stop codon T-. The ND1 and ND5 genes of I. hauseri used the incomplete stop codon T- (Supplementary Table S6), while the complete stop codon TAG was mainly found in the ATP8, ND4L, ND3, ND5, ND1, and CYTB genes of all eight species, while other genes used TAA as the stop codon.

Comparison of Codon Usage in Protein-Coding Genes and tRNA Secondary Structure
Analysis of the amino acid content and relative codon usage frequency in the PCGs of the eight species revealed the presence of all possible synonymous codons of the 22 amino acids ( Figure 1). Among these, Phe, Leu1, and Ile were the most abundant with frequencies greater than 300, and Cys was the least abundant with frequencies less than 50. Analysis of the relative use of synonymous codons showed that RSCU was highest for NNA and NNU codons, generally greater than 1, and lowest for NNG and NNC codons ( Figure 2). Codons with relatively high G and C content were rare, confirming a finding in other insect species. Overall, UUA, UCU, and GUU were the three most common relatively synonymous codons.
The tRNA length in the eight species analyzed here ranged from 62 bp to 72 bp. Most base pairs in the stem region conformed to the Watson-Crick pairing principle, and a few followed the G-U vibrational pairing principle. There were 166 G-U pairs, with I. hauseri having the most mismatched pairs and Dendroctonus micans having the fewest. In our data, mismatched bases in the secondary structure of tRNA occurred mainly in the amino acid acceptor arm and least frequently in the TΨC arm (Table 2). Thus, the different regions showed different degrees of conservatism, with the TΨC arm being the most conservative. Differences in tRNA secondary structure were small among species belonging to the same genus or closely related species. Insects 2021, 12, x 5 of 12  The tRNA length in the eight species analyzed here ranged from 62 bp to 72 bp. Most base pairs in the stem region conformed to the Watson-Crick pairing principle, and a few followed the G-U vibrational pairing principle. There were 166 G-U pairs, with I. hauseri having the most mismatched pairs and Dendroctonus micans having the fewest. In our data, mismatched bases in the secondary structure of tRNA occurred mainly in the amino acid acceptor arm and least frequently in the TΨC arm (Table 2). Thus, the different regions showed different degrees of conservatism, with the TΨC arm being the most conservative. Differences in tRNA secondary structure were small among species belonging to the same genus or closely related species.   The tRNA length in the eight species analyzed here ranged from 62 bp to 72 bp. Most base pairs in the stem region conformed to the Watson-Crick pairing principle, and a few followed the G-U vibrational pairing principle. There were 166 G-U pairs, with I. hauseri having the most mismatched pairs and Dendroctonus micans having the fewest. In our data, mismatched bases in the secondary structure of tRNA occurred mainly in the amino acid acceptor arm and least frequently in the TΨC arm (Table 2). Thus, the different regions showed different degrees of conservatism, with the TΨC arm being the most conservative. Differences in tRNA secondary structure were small among species belonging to the same genus or closely related species.

Phylogenetic Analysis
Bayesian phylogenetic trees constructed from the complete mitochondrial genome had high confidence values, and the monophyly of the tribe and genera of the species studied Insects 2021, 12, 949 6 of 12 was well supported ( Figure 3A). The phylogenetic relationships of the subfamily Scolytinae by tribe were determined as follows: (((Ipini + (Xyleborini + Dryocoetini)) + (Trypophloeini + Corthylini)) + (Polygraphini + Xyloterini)) + (Hylurgini + Hylastini)). The topologies of the phylogenetic trees constructed based on the methods of ML are similar to those of the BI tree ( Figure 3B), and the monophyly of the tribe and most genera was also well supported in the ML tree. However, the phylogenetic relationships between some species are still confusing. For example, both the ML tree and the BI tree show that the genus Dryocoetes is poorly defined due to the position of the species Dryocoetes villosus (Fabricius, 1792) in the phylogenetic tree. In addition, Xylosandrus crassiusculus (Motschulsky, 1866) and Anisandrus dispar (Fabricius, 1792) have different positions in the ML tree and BI tree. The genus Xylosandrus is monophyletic in the BI tree but not in the ML tree. We also constructed the phylogenetic tree for the undetermined species of the subfamily Scolytinae using the data from the NCBI database (Supplementary Figure S9)

Phylogenetic Analysis
Bayesian phylogenetic trees constructed from the complete mitochondrial genome had high confidence values, and the monophyly of the tribe and genera of the species studied was well supported ( Figure 3A). The phylogenetic relationships of the subfamily Scolytinae by tribe were determined as follows: (((Ipini + (Xyleborini + Dryocoetini)) + (Trypophloeini + Corthylini)) + (Polygraphini + Xyloterini)) + (Hylurgini + Hylastini)). The topologies of the phylogenetic trees constructed based on the methods of ML are similar to those of the BI tree ( Figure 3B), and the monophyly of the tribe and most genera was also well supported in the ML tree. However, the phylogenetic relationships between some species are still confusing. For example, both the ML tree and the BI tree show that the genus Dryocoetes is poorly defined due to the position of the species Dryocoetes villosus (Fabricius, 1792) in the phylogenetic tree. In addition, Xylosandrus crassiusculus (Motschulsky, 1866) and Anisandrus dispar (Fabricius, 1792) have different positions in the ML tree and BI tree. The genus Xylosandrus is monophyletic in the BI tree but not in the ML tree. We also constructed the phylogenetic tree for the undetermined species of the subfamily Scolytinae using the data from the NCBI database (Supplementary Figure S9)

Analysis of the Genetic Distance and the Evolutionary Rate
Genetic distances between sequenced bark beetles based on the whole mitogenome supported their phylogenetic relationships. The distances between species within a genus were shorter than the distances between genera. The shortest distance, and thus the closest genetic relationship, was between I. typographus and I. nitidus (0.037), whereas the greatest distance was between P. poligraphus and G. materiarius (0.418) (Figure 4). Among the bark beetles of the Scolytinae, most of the genetic distances between species were greater than 0.1, except for the distances between Dryocoetes hectographus and Dryocoetes autographus, Trypodendron signatum (Fabricius, 1792) and Trypodendron domesticum (Linnaeus, 1758), and I. typographus and I. nitidus, suggesting that they are sibling species.
The Ka/Ks ratio, a diagnostic statistical method for detecting molecular adaptations, is used to estimate the evolutionary rate among insects. Evolutionary rates of PCGs among Scolytinae species, as measured by rates of nonsynonymous (Ka) and synonymous (Ks) substitutions and Ka/Ks ratio, are shown in Figure 5. Average Ka/Ks ratios were consistently less than 1 and ranged from 0.07 for COI to 0.69 for ND2, indicating that these mitogenomes evolve under purifying selection. All mitochondrial PCGs could be used to examine phylogenetic relationships within the Scolytinae.  The Ka/Ks ratio, a diagnostic statistical method for detecting molecular adaptations, is used to estimate the evolutionary rate among insects. Evolutionary rates of PCGs among Scolytinae species, as measured by rates of nonsynonymous (Ka) and synonymous (Ks) substitutions and Ka/Ks ratio, are shown in Figure 5. Average Ka/Ks ratios were consistently less than 1 and ranged from 0.07 for COI to 0.69 for ND2, indicating that these mitogenomes evolve under purifying selection. All mitochondrial PCGs could be used to examine phylogenetic relationships within the Scolytinae.

Discussion
In this study, we performed a comparative mitogenome analysis and revealed the conservatism of the eight bark beetle mitogenomes. The sequences of the control region

Discussion
In this study, we performed a comparative mitogenome analysis and revealed the conservatism of the eight bark beetle mitogenomes. The sequences of the control region of mitogenome were not assembled because their distribution and length are uncertain. Most of the sequences in the control region are located between the rrnS and trnI genes, and some of them are divided into two segments by the trnI gene [40,41]. The structure of PCGs in the bark beetle mitogenome followed the same pattern as in Iberobaenia beetles; gene rearrangement events, especially those related to PCGs, are rare in mitogenomes of Coleoptera [42]. Bark beetles exhibit the typical AT-biased composition of insect mitogenomes [43,44]. In them, AT content was significantly lower in PCGs and significantly higher in rRNA and tRNA, which might be due to PCGs being more functionally restricted than other regions. The biological reasons for such AT-biased composition have been extensively studied [45,46], and it has been hypothesized that AT-biased composition is more energy efficient [47,48]. We have also found that base shifts in the same region of the mitogenome are similar between species, with PCGs showing a strong preference for T and C, rRNA for T and G, and tRNA for A and G. Although the mechanisms for such skewed composition are complex, most hypotheses suggest that this phenomenon is due to mutation and selection pressure [49]. In insect mitogenomes, the extent of GC-biased composition appears to be related to gene replication and selection pressure during transcription [45], which is not consistent with the direction of mutations [49,50]. Regarding the use of start and stop codons in PCGs, PCGs of most species used ATN as the start codon. The start codons for ND1 are often non-standard in holometabolous insects, with the non-standard TTG being used in almost half of the known mitogenomes of species in the Curculionidae [51]. Stop codons include TAA, TAG, and T-, and the incomplete stop codon (T-) can become the complete stop codon by polyadenylation during posttranscriptional mRNA processing TAA [52], which is also common in the insect mitogenome [53]. Moreover, the start and stop codons also had the property of being AT-biased.
According to the secondary structure model of mitochondrial tRNA genes proposed by Kumazawa and Nishida [54], tRNA sequences can be arranged in stem-loop structures. With the exception of the tRNASer (AGN), which lacks the DHU (dihydrouridine) arm, all others are classical cloverleaf structures. In addition, we found that base mismatches are abundant in the secondary structure of bark beetles, which is common in the secondary structure of insect tRNA. For example, the ghost moth Thitarodes yunnanensis Yang has 18 mismatches [55], while the nymphalid butterfly Parathyma sulpitia Cramer has 22 mismatches [56]. Yokobori and Paabo [57] speculated that these mismatches can be explained by the lack of recombination in the insect mitogenome, making such mutations difficult to remove. These mismatches occurred mainly on the amino acid acceptor arm of the secondary structure of eight bark beetles, and least on the TΨC arm. The closer the relationship, the smaller the structural difference and the more similar the number of mismatches. This result, combined with the uneven distribution of mismatched bases in different tRNA regions, highlights the need for a large-scale comparison of tRNA secondary structure between species. Such an analysis could help reveal the evolutionary pathways of insect mitochondrial tRNA and the functions of individual regions and base mutations, providing a basis for further research into the evolution of mitogenome structure.
The constructed phylogenetic trees strongly supported the monophyly of the tribe and most of the genera in our analysis. However, the position of Dryocoetes villosus in the trees of ML and BI is different, which needs to be redefined and revised. Jordal et al. [58] also found that the genus Dryocoetes is not monophyletic by examining the 28S gene, which cannot be separated from the genus Taphrorychus in the phylogenetic tree. The genus Xylosandrus currently contains species with highly variable morphology, some of which resemble those of other genera, which has led to confusion regarding the generic boundaries of Xylosandrus and its relationship to and separation from the genus Anisandrus [59]. Therefore, many scholars have made great efforts to classify this genus [60][61][62][63]. For example, Dole and Cognato have recently revised Xylosandrus [64] and preliminary phylogenies suggest that Xylosandrus is monophyletic [61]. In addition, we have also investigated the phylogenetic relationships of some unidentified species in the subfamily Scolytinae using data from the NCBI database. However, these species are not placed in a well supported lineage, and much information and mitochondrial data are needed to define their phylogenetic relationships. We also analyzed the genetic distances among bark beetles, which supports the phylogenetic relationships of the studied species. An analysis of vertebrate and invertebrate COI gene sequences by Hebert et al. [65] revealed an average interspecific difference in genetic distance of approximately 0.113 in 98% of species. Within Ips species, Chang et al. [66] found an interspecific genetic distance between 0.056 and 0.431, similar to our results.
In general, we have described the characteristics of the mitogenomes of eight species from the subfamily Scolytinae and performed a comparative analysis of the sequence of gene arrangement, nucleotide composition, use of PCGs codons, tRNA secondary structure, genetic distance, and evolutionary rate among bark beetle species. The typical structural characteristics and conservatism of the bark beetle mitogenome suggest that it can be used to clarify phylogenetic relationships at higher taxonomic levels. Based on the mitogenome data, the constructed phylogenetic relationship among species of Scolytinae has produced an in-depth study of the genetic evolution of bark beetles. However, the sequenced species are still quite few, which directly limits a definitive species delimitation to an undisputed lineage in bark beetles.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/insects12100949/s1, Table S1: The list of sample information of bark beetles. Table S2: Annotation of the mitogenome of Orthotomicus erosus. Table S3: Annotation of the mitogenome of Dryocoetes hectographus. Table S4: Annotation of the mitogenome of Polygraphus poligraphus. Table S5: Annotation of the mitogenome of Dendroctonus micans. Table S6: Annotation of the mitogenome of Ips hauseri. Table S7: Annotation of the complete mitogenome of Ips subelongatus. Table S8: Annotation of the mitogenome of Ips typographus. Table S9: Annotation of the mitogenome of Ips nitidus. Table S10. Investigation of nucleotide composition in the mitochondrial genomes of eight bark beetles. Figure S1: Secondary structures for tRNA genes from the mtDNA of Dryocoetes hectographus. Figure S2: Secondary structures for tRNA genes from the mtDNA of Orthotomicus erosus. Figure S3: Secondary structures for tRNA genes from the mtDNA of Ips typographus. Figure S4: Secondary structures for tRNA genes from the mtDNA of Ips nitidus. Figure S5: Secondary structures for tRNA genes from the mtDNA of Ips hauseri. Figure S6: Secondary structures for tRNA genes from the mtDNA of Ips subelongatus. Figure S7: Secondary structures for tRNA genes from the mtDNA of Dendroctonus micans. Figure S8: Secondary structures for tRNA genes from the mtDNA of Polygraphus poligraphus. Figure Table 1.

Conflicts of Interest:
The authors declare no conflict of interest.