Insights into the Divergence of Chinese Ips Bark Beetles during Evolutionary Adaptation

Simple Summary Bark beetle species of the genus Ips are among the major pests of Chinese conifer forests. Based on mitochondrial genome and SNP, we investigated the phylogenetic relationships and evolutionary trends of 19 populations of six Ips species that had serious outbreaks in recent years. Our results demonstrated the relationships between Ips evolution and host plants, pheromones, and altitudinal differences, and provided new insights into the mechanism of adaptive evolution of Ips bark beetles. Abstract Many bark beetles of the genus Ips are economically important insect pests that cause severe damage to conifer forests worldwide. In this study, sequencing the mitochondrial genome and restriction site-associated DNA of Ips bark beetles helps us understand their phylogenetic relationships, biogeographic history, and evolution of ecological traits (e.g., pheromones and host plants). Our results show that the same topology in phylogenetic trees constructed in different ways (ML/MP/BI) and with different data (mtDNA/SNP) helps us to clarify the phylogenetic relationships between Chinese Ips bark beetle populations and Euramerican species and their higher order clades; Ips bark beetles are polyphyletic. The structure of the mitochondrial genome of Ips bark beetles is similar and conserved to some extent, especially in the sibling species Ips typographus and Ips nitidus. Genetic differences among Ips species are mainly related to their geographic distribution and different hosts. The evolutionary pattern of aggregation pheromones of Ips species reflects their adaptations to the environment and differences among hosts in their evolutionary process. The evolution of Ips species is closely related to the uplift of the Qinghai-Tibet Plateau and host switching. Our study addresses the evolutionary trend and phylogenetic relationships of Ips bark beetles in China, and also provides a new perspective on the evolution of bark beetles and their relationships with host plants and pheromones.


Introduction
Tree-destroying bark beetles (Coleoptera, Curculionidae, Scolytinae) are the most destructive pests of conifer forests worldwide and, therefore, are of considerable ecological importance [1]. The diversity of factors that determine the population dynamics of bark beetles in forest ecosystems is a challenge for many researchers. These include host phytochemistry [2,3], host range expansion [4], beetle-microbe interactions [5,6], influences of natural enemies, climatic changes, and microhabitat differences [7,8]. Bark beetles of the genus Ips are mainly distributed in conifer forests of the Northern Hemisphere. To plants [44]. Live adults were preserved in anhydrous ethanol and stored in a laboratory freezer at −20 • C until biological analysis. Detailed information on collection sites, hosts, and abbreviated population codes can be found in Table 1. Genomic DNA was extracted from thoracic and leg muscle tissues (five individuals per population) and purified using the Wizard Genomic DNA Purification Kit (Promega Corporation, Madison, WI, USA), according to the manufacturer's protocol. Then, DNA preparations were quantified using PicoGreen ® (Thermo Fisher Scientific, Wilmington, DE, USA) and a NanoDrop 2000 fluorospectrometer (Thermo Fisher Scientific, Wilmington, DE, USA), with DNA integrity detected by 1% agarose gel electrophoresis.

Mitogenome De Novo Sequencing and Assembly
Each genomic sample was sheared using a Covaris M220 ultrasonicator (Covaris, Woburn, MA, USA) to obtain 400-500 bp DNA fragmentation. A 460 bp paired-end library was generated from each sample and 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 (http://www. usadellab.org/cms/?page=trimmomatic (accessed on 1 July 2020)). Raw data were sheared as follows: to delete sequences containing splice fragments and to remove sequences with low mass values (Q < 25). The software Spades [45] was used to assemble the cleaned data and the results were compared by mapping to obtain the longest segments. The mitoZ program [46] was used to annotate the genome of the assembly results. The coding sequence (CDS), transfer RNA (tRNA), and ribosomal RNA (rRNA) regions were filtered out. The mitochondrial genome sequences for each species were submitted to GenBank and the accession numbers are listed in Table 1.

Phylogenetic Analysis
Based on the complete mitochondrial genome sequence, we constructed the maximum likelihood tree (ML) and the maximum parsimony tree (MP) using MEGA v5.2.2 software for the 19 Ips populations with Dendroctonus micans (Curculionidae, Scolytinae) as an outgroup. The confidence values of each branch node of the phylogenetic tree were tested by bootstrapping with 1000 replicates [47]. We also constructed the Bayesian inference tree (BI) using MrBayes v3.2.2 software (Uppsala, Sweden) [48] and used the optimal alternative model (GTR+G) to build the evolutionary tree. The MCMC method was used to calculate 2,000,000 generations, sampling every 100 generations to ensure sampling independence, and the original 25% of the data were discarded as burn-in. Stationarity was considered achieved when the average standard deviation of split frequencies was less than 0.01 [49]. To illustrate the phylogenetic relationships of Ips species within the Scolytinae, we also constructed the ML and BI trees by combining the mitochondrial genomes of 26 published species of Scolytinae (Table S1) and our studied six species with Curculio elephas (Curculionidae, Curculioninae) as an outgroup. In addition, a total of 43 Ips species with 105 COI sequences (19 from our study and 86 from the published literature [9,12,13]) were selected to construct a phylogenetic tree using the BI method with D. micans, Pseudips concinnus, and Pseudips mexicanus as outgroups (Table S2).

Estimation of Divergence Time
We used the sequences of 13 protein-coding genes to estimate the divergence time for the Ips species using the software BEAST v2.3.0 [50]. MrMTgui software was used to select the most appropriate sequence replacement model (GTR) for comparative analyses, and the site heterogeneity model was set to gamma; clock was set to strict clock. The Yule procedure calculated 20,000,000 generations using the MCMC method, sampling every 1000 generations and discarding the first 20% of the data as burn-in. In the absence of relevant fossil records as time markers, the application of the molecular clock was chosen based on the insect mitochondrial genome, which contained the protein-coding genes of COXI and had an average substitution rate of 0.00175 substitutions/site/Myr (s/s/Myr) [51]. The substitution rates of other mitochondrial genes were scaled to the mean COXI rate, and a weak lognormal distribution with a large standard deviation was used to allow COXI calibration to control the analysis. We imported the results generated by the BEAUti 2 software in .xml format into BEAST v2.3.0 to perform the analysis. Then, the Tracer v1.5 software (http://tree.bio.ed.ac.uk/software/tracer (accessed on 1 November 2020)) was used to check the results and obtain a standard distribution based on the effective sample size (ESS) of each estimated parameter, with ESS requiring at least 200 for the posterior and prior. The Tree Annotator V2.3.0 software was used to calculate the consensus tree and annotate the divergence time, and the FigTree V1.3.1 software (http://tree.bio.ed.ac.uk/software/figtree (accessed on 1 December 2020)) was used to open it and check the analysis results.

Comparative Mitogenome Analyses of Bark Beetles
Circle maps of the mitogenome of six Ips species (e.g., ItyJG, IsuJM, IshQL, IhaTS, In-iXZ, and IduBY) and one Dendroctonus species (DmiMX) were generated and displayed using the CGView software (http://stothard.afns.ualberta.ca/downloads/CCT/tutorials.html (accessed on 1 July 2020)). The nucleotide composition, codon usage (except for stop codons), and relative synonymous codon usage (RSCU) of the mitochondrial genomes of the analyzed Ips species were calculated using the MEGA 5.2.2 software (Paris, France). The compositional skew was calculated using the following formulas: AT skew = (A − T)/(A + T) and GC skew = (G − C)/(G + C) [52]. To test the interspecific distinctions of the Ips species, the numbers of base differences per site between sequences were presented. The analyses included seven nucleotide sequences and the codon positions were 1 + 2 + 3. All positions with less than 95% coverage were eliminated, i.e., less than 5% alignment gaps, missing data, and ambiguous

Genetic Distance, Selection Pressure and Isolation by Distance (IBD) Analysis
Sequences were preliminarily aligned using the Clustal X program (lllkirch, France) [54]. Pairwise genetic distances were calculated using the MEGA v5.2.2 software based on the Kimura-2 parameter model [55]. The CodeML software implemented in the PAML package was used with the free ratio model to calculate the non-synonymous (Ka) and synonymous (Ks) substitution rates along with each branch of the tree [56]. The selection pressure of the mtDNA genome was assessed by calculating the Ka/Ks ratios of 13 mitochondrially encoded protein genes in different species. The haplotype network of Ips species was analyzed using a median-joining algorithm in the program Network 4.6 [57]. We calculated the geographic distance of the sampling sites using the longitude and latitude data and evaluated the correlations between the genetic distance and the geographic distance of the species using the Mantel test of IBD 1.53 [58].

SNP Sequencing and Analysis
In this study, we also sequenced the SNP of the Ips populations. The DNA of each sample was extracted using a TIANamp Genomic DNA Kit (Osce Biological Technology, Shanghai, China), and its quality and concentration were determined using the agar-agar assay and a Nanodrop2000 fluorospectrometer. Paired-end sequencing was performed after qualified DNA extraction on the Illumina Hiseq Xten platform (Osce Biological Technology, Shanghai, China). Paired-end sequencing was also performed using 2b-RAD technology, with all samples attached to the enzyme digestion label using the 5 -NNN-3 linker. Upon completion of sequencing, reads containing restriction sites were extracted from the sequencing data. The Stacks software v 1.34 (Eugene, OR, USA) [59] was used for clustering and reference sequences were generated. The sequencing data were aligned to the reference sequence using the software SOAP v 2.21 (Shenzhen, China) [60]. The methods of haplotype analysis and RAD sequence-based ML tree construction are the same as for mitochondria. The software STRUCTURE v 2.3.4 (Oxford, UK) [61] was used to analyze the population structure from K = 2 to K = 10, and 10 different seeds were selected for 10 repeated analyses. According to the optimal K value determined by the cross-validation error, a box shape plot was drawn for each K value, and the CV value was repeated N times.

Phylogenetic Relationships, Divergence Times, and Pheromone Differences in Ips Bark Beetles
The phylogenetic trees ML, MP, and BI, which were constructed based on the PCGs of the Ips species/populations, showed similar topology with high nodal support in most of the central nodes ( Figure 1A). The phylogenetic relationships of the Ips species are as follows: (((I. typographus + I. nitidus) + I. subelongatus) + I. shangrila) + (I. duplicatus + I. hauseri). Each species belongs to a separate branch, although I. typographus and I. nitidus are phylogenetically closely related, and I. duplicatus and I. hauseri are also closely related. In addition, using Curculio elephas as an outgroup, we created ML and BI trees by combining 26 species of Scolytinae (whose mitochondrial genomes were sequenced) with mtDNA sequences of the Ips species analyzed in this study ( Figure 1B). Initially, the phylogenetic trees based on the mtDNA and the 13 PCGs of the Ips species showed a similar topology as in Figure 1A. The phylogenetic relationship within the trees ML and BI was constructed as follows: ((((Ipini + Polygraphini) + (Drycoetini + Xyleborini)) + (Trypophloeini + Corthylini)) + Xyloterini) + (Hylurgini + Hylastini), and the phylogenetic tree exhibited high bootstrap values. The phylogenetic analysis using SNP also showed that the phylogenetic relationships of Ips species were similar to those established using the mitochondrial genome ( Figure S1). relationships of Ips species were similar to those established using the mitochondrial genome ( Figure S1). The results of BEAST dating analysis showed that the age of the youngest common ancestor of Ips species was 12.02 Mya, indicating the earliest divergence time of Ips species in China. Among them, the divergence time of I. subelongatus was estimated to be 6.49 Mya; I. duplicatus and I. hauseri began to diverge in 4.89 Mya, and I. typographus and I. nitidus in 1.54 Mya. The Ips species have gradually evolved to their present status over the past 400,000 years ( Figure S2).
The major and minor pheromone components of Ips bark beetles were superimposed on the phylogenetic tree, and a significant difference was found between pheromone components in closely related Ips species ( Figure 1A). The pheromone components 2-methyl-3-butten-2-ol, (S)-(+)-ipsdienol, (R)-(-)-ipsdienol, and (S)-(-)-ipsenol, which are biosynthesized de novo in the anterior midgut of bark beetles, did not differ gradually among Ips species. A similar pattern was observed for the pheromone component (S)-cis-verbenol, which is derived from a host precursor. It should be noted that the major and minor components of aggregation pheromones of a given species are not tied to a particular synthetic pathway and that interspecific changes in pheromone components do not parallel their evolutionary relationships.

Phylogenetic Relationships among Ips Bark Beetles Worldwide
We used the Bayesian method to reconstruct phylogenetic trees based on 105 COI sequences for 43 Ips species worldwide ( Figure 2). The results showed that Ips bark beetles were polyphyletic. When we combine the phylogenetic relationships of Ips species with The results of BEAST dating analysis showed that the age of the youngest common ancestor of Ips species was 12.02 Mya, indicating the earliest divergence time of Ips species in China. Among them, the divergence time of I. subelongatus was estimated to be 6.49 Mya; I. duplicatus and I. hauseri began to diverge in 4.89 Mya, and I. typographus and I. nitidus in 1.54 Mya. The Ips species have gradually evolved to their present status over the past 400,000 years ( Figure S2).
The major and minor pheromone components of Ips bark beetles were superimposed on the phylogenetic tree, and a significant difference was found between pheromone components in closely related Ips species ( Figure 1A). The pheromone components 2methyl-3-butten-2-ol, (S)-(+)-ipsdienol, (R)-(-)-ipsdienol, and (S)-(-)-ipsenol, which are biosynthesized de novo in the anterior midgut of bark beetles, did not differ gradually among Ips species. A similar pattern was observed for the pheromone component (S)-cisverbenol, which is derived from a host precursor. It should be noted that the major and minor components of aggregation pheromones of a given species are not tied to a particular synthetic pathway and that interspecific changes in pheromone components do not parallel their evolutionary relationships.

Phylogenetic Relationships among Ips Bark Beetles Worldwide
We used the Bayesian method to reconstruct phylogenetic trees based on 105 COI sequences for 43 Ips species worldwide ( Figure 2). The results showed that Ips bark beetles were polyphyletic. When we combine the phylogenetic relationships of Ips species with the correlations of their hosts, we find that the species are more tightly linked to a specific range and the ranges of the hosts are also more tightly linked. The hosts of spruce, larch, and pine do not overlap, and 23 Ips species colonize pine, 18 species feed mainly on spruce, and I. subelongatus and Ips cembrae Heer 1836 mainly damage larch. We also performed a haplotype analysis based on the 105 COI sequences of Ips species (Figure 3). A total of 81 haplotypes were identified, including 13 common haplotypes and 48 exclusive haplotypes. The sequences within each common haplotype belong to the same species. The haplotype network results also showed that different haplotypes of the same species were relatively close to each other in the diagram, which intuitively illustrated the relationship between Ips species. the correlations of their hosts, we find that the species are more tightly linked to a specific range and the ranges of the hosts are also more tightly linked. The hosts of spruce, larch, and pine do not overlap, and 23 Ips species colonize pine, 18 species feed mainly on spruce, and I. subelongatus and Ips cembrae Heer 1836 mainly damage larch. We also performed a haplotype analysis based on the 105 COI sequences of Ips species (Figure 3). A total of 81 haplotypes were identified, including 13 common haplotypes and 48 exclusive haplotypes. The sequences within each common haplotype belong to the same species. The haplotype network results also showed that different haplotypes of the same species were relatively close to each other in the diagram, which intuitively illustrated the relationship between Ips species.

Analysis of Selection Pressure, Genetic Distance, and IBD Test
In this study, the six Ips species analyzed are mainly distributed in five mountains in China, which have a wide distribution range and high altitudinal variation (Figure 4). Analysis of selection pressure with 13 PCGs in six Ips species revealed a separate Ka/Ks ratio for each terminal branch in the tree ML (Table 2). Ips nitidus with the highest altitude distribution had the highest mean (0.0547), but I. subelongatus with the lowest altitude distribution retained the lowest mean (0.0186). The same change trends were also observed for ATP8, COII, COIII, ND1, ND3, ND4, ND4L, and ND5 genes. The genetic distance analyses were consistent with the results of the phylogenetic analyses ( Figure 5). The genetic distance value between I. nitidus and I. typographus was the lowest among Ips species, indicating their closest genetic relationship; however, the value of I. shangrila and I. duplicatus was the highest, indicating the most distantly related species. The genetic distances between species were greater than 0.1, except for I. nitidus and I. typographus, which were greater Biology 2022, 11, 384 8 of 18 between species than between populations. The results of the Mantel test showed that there were no significant correlation relationships between genetic distance and geographic distance among species of the genus Ips. However, only significant correlations were found between geographic distance and genetic distance between populations of I. typographus as well as populations of I. nitidus based on different algorithms (Table 3).

Analysis of Selection Pressure, Genetic Distance, and IBD Test
In this study, the six Ips species analyzed are mainly distributed in five mountains in China, which have a wide distribution range and high altitudinal variation ( Figure 4). Analysis of selection pressure with 13 PCGs in six Ips species revealed a separate Ka/Ks ratio for each terminal branch in the tree ML (Table 2). Ips nitidus with the highest altitude distribution had the highest mean (0.0547), but I. subelongatus with the lowest altitude distribution retained the lowest mean (0.0186). The same change trends were also observed for ATP8, COII, COIII, ND1, ND3, ND4, ND4L, and ND5 genes. The genetic distance analyses were consistent with the results of the phylogenetic analyses ( Figure 5). The genetic distance value between I. nitidus and I. typographus was the lowest among Ips species, indicating their closest genetic relationship; however, the value of I. shangrila and I. duplicatus was the highest, indicating the most distantly related species. The genetic distances between species were greater than 0.1, except for I. nitidus and I. typographus, which were greater between species than between populations. The results of the Mantel test showed that there were no significant correlation relationships between genetic distance and geographic distance among species of the genus Ips. However, only significant correlations were found between geographic distance and genetic distance between populations of I. typographus as well as populations of I. nitidus based on different algorithms (Table 3).

Comparative Mitochondrial Genome Analyses of Ips Bark Beetles
We obtained mitogenomic sequences from 19 populations of six Ips species with all 36 functional mitochondrial genes except the control region. Sequence lengths vary from 15,259 bp to 15,641 bp and all of these sequences have the same gene arrangement as Drosophila yakuba Burla 1954 [62] (Figure S3). We compared the similarity of mitochondrial genome sequences of Ips bark beetles, including PCGs, tRNA, and rRNA regions, with those of I. typographus ( Figure 6). Among them, the sequence similarity of I. nitidus was the highest at over 94%; the sequence similarity of I. hauseri, I. duplicatus, I. shangrila, and I. subelongatus was mostly between 82% and 88% as compared with that of I. typographus; the sequence similarity between D. micans and I. typographus was practically less than 82%.
Among genes in the entire coding region of the mitochondrial genome of six Ips species, the number of spacer regions between adjacent genes ranged from 23 to 27, with the highest number in I. typographus and I. subelongatus. The number of gene overlap regions ranged from 3 to 5. In the sequence of PCGs, tRNA and rRNA regions of the mitochondrial genome of six Ips species, the AT-content was highest in I. subelongatus. The values of AT-skew and GC-skew among Ips species differed only slightly, and the average AT content in each region of mtDNA was significantly lower than that of D. micans (Table S3).            The ND1 protein-coding genes of I. subelongatus, I. typographus, I. hauseri, and I. nitidus all used TTG as the start codon, whereas the others used ATN as the start codon. The ATP8 gene of six Ips species used TAG as the termination codon, whereas most of the other TAA and a few genes used T as the termination codon (Tables S4-S10). We also analyzed the protein-coding genes in the mitochondrial genome, with the amino acids Leu used most frequently and Cys used least frequently. The most frequently used condons were TTT, ATT, and TTA, and AGC with the lowest frequency, which also reflected the AT preference of their nucleotide composition ( Figure S4). Among genes in the entire coding region of the mitochondrial genome of six Ips species, the number of spacer regions between adjacent genes ranged from 23 to 27, with the highest number in I. typographus and I. subelongatus. The number of gene overlap regions ranged from 3 to 5. In the sequence of PCGs, tRNA and rRNA regions of the mitochondrial genome of six Ips species, the AT-content was highest in I. subelongatus. The values of ATskew and GC-skew among Ips species differed only slightly, and the average AT content in each region of mtDNA was significantly lower than that of D. micans (Table S3).
The ND1 protein-coding genes of I. subelongatus, I. typographus, I. hauseri, and I. nitidus all used TTG as the start codon, whereas the others used ATN as the start codon. The ATP8 gene of six Ips species used TAG as the termination codon, whereas most of the other TAA and a few genes used T as the termination codon (Tables S4-S10). We also analyzed the protein-coding genes in the mitochondrial genome, with the amino acids Leu used most frequently and Cys used least frequently. The most frequently used condons were TTT, ATT, and TTA, and AGC with the lowest frequency, which also reflected the AT preference of their nucleotide composition ( Figure S4).

Phylogenetic Relationships of Ips Species in Relation to Their Biological Characteristics
The phylogenetic relationships and evolutionary trend of Ips species in China are well established in our study, further confirming the taxonomic status of I. shangrila in the genus Ips identified by Cognato et al. [63]. We also analyzed the phylogenetic relationship between the six species of Ips in China and the other bark beetles in Scolytinae, which was essentially consistent with the reports of Du et al. [17]. The phylogenetic results of Ips species are positively correlated with their host ranges. Feng et al. [64] constructed a deep lineage of the conifer genus Picea, the host of most Ips bark beetles. Here, we compared the phylogenetic relationship between Ips bark beetles and their hosts. The main hosts of I. typographus and I. nitudus are Picea koraiensis Nakai 1919 and Picea crassifolia Kom. 1923, respectively, which are closely related in the phylogenetic analysis of the hosts. Similarly, the main host of I. hauseri is Picea schrenkiana Fisch. et Mey., which is genetically distant from Korean spruce and thick-leaved spruce. The host of I. duplicatus is Picea meyeri Rehd. et Wils. 1914 [32], a species endemic to China, distributed only on the eastern edge of the Hunshandake Sandy Land in Inner Mongolia, and whose systematic status has not yet been determined. The host of I. subelongatus is Larix sp. (e.g., Larix gmelinii Rupr., Larix kaempferi Carr., Larix principis-rupprechtii Mayr, and Larix sibirica Ledeb.), indicating its special taxonomic status within the genus Ips. Moreover, the phylogenetic relationship is related to the biological characteristics of Ips bark beetles. For example, the body length of I. hauseri, I. duplitatus, and I. shangrila is about 3 mm; in contrast, the body length of I. typographus and I. nitidus is about 7 mm, and I. subelongatus is slightly larger than 7 mm. In addition, the corresponding characteristics in the elevations of the frontal tubercles and the number of frontal hairs (I. typographus and I. nitidus), and the shape and spacing of the teeth (I. subelongatus, I. hauseri, I. duplitatus, and I. shangrila) also help to distinguish the sexes [65,66].

Evolutionary Lineages of the Ips Species in Relation to the Geographic Environment
After analyzing the divergence times, we thought that the evolutionary trend of Ips species might be closely related to the geological events in the mountain ranges where they are distributed, as well as to their hosts. The bark beetles of the genus Ips evolved in the Miocene when the Qinghai-Tibet Plateau and the Tienshan Mountains were uplifted. Geological events caused their ancestors to seek refuge and gradually split into new species as they adapted to the new habitat. The evolution and dispersal of their hosts was also influenced by geological events at an appropriate time [77]. The uplift of the Qinghai-Tibet Plateau and the Tienshan Mountains from the Miocene to the Pliocene accelerated the climatic drought and Asian monsoon within the Asian continent [78], which played an important role in the spread of spruce/larch [79]. Similar results were found for other highland species in China, such as the Asian honeybee [80] and locust [81].
Since most bark beetle species in China are distributed in a wide range of elevations, some of them at high or even very high elevations, the genetic adaptation mechanisms of Ips bark beetles to high elevation environments were also analyzed from the perspective of mitogenomics. We found that selection pressure is significantly higher in I. nitidus and I. shangrila, which are common at higher elevations, than in I. typographus and I. subelongatus, which are common at lower elevations. This suggests that the low oxygen and cold environment at higher elevations causes selection on the mitochondrial gene of Ips species. More than 95% of the total cellular energy in eukaryotic cells is generated by mitochondrial OXPHOS [82]. Proteins encoded by 13 mitochondrial genes play a key role in adaptation to extreme, high-altitude environments [83,84]. The results of selection pressure reflect the adaptability of mitochondrial genes of Ips species to different microhabitats. Similar results have been found in high-altitude species such as Tibetan wild yak [85], Tibetan human [86], and schizothoracine fish [87].
The results of genetic distance of Ips species were consistent with those of phylogenetic trees, with significantly less genetic distance within species than between species. There were no significant correlations between geographic distance and genetic distance among Ips species, but they are closely related to their ranges, hosts, and biological characteristics. Ips typographus and I. subelongatus have stronger flight ability and a wider host range than other Ips species, resulting in a relatively small genetic distance between them and other species; in contrast, Ips duplicatus, I. hauseri, and I. shangrila have poorer flight ability and breed only in species of the genus Picea. The barriers between the Tarim Basin, the Qilian Mountains, and the Plateau increase the genetic distance between them and other species. In contrast, the genetic distance between I. nitidus and I. typographus was small, suggesting that this is the result of the evolution of I. nitidus to adapt to the extreme environment. The fact that the divergence time of I. nitidus and I. typographus not sympatric species, is the youngest in our analyses confirms their relatedness.

Adaptive Evolution Inferred from Comparative Mitochondrial Genome ANALYSIS
The base content of the entire mitochondrial genome of six Ips species follows the same pattern: A > T > C > G, with the AT content being significantly higher than the GC content, similar to the species I. grandicllis and I. sexdentatus [88]. The sorting order of gene sequences and the use of start and end codons are also consistent with the reported species (Ips grandicollis, Ips typographus, and Ips sexdentatus) [88], suggesting evolutionary conservatism of the mitochondrial genome in Ips species [82]. The results of sequence consistency comparisons among Ips species are identical to those of phylogenetic trees. For example, the sequence consistency of I. nitidus and I. typographus is over 94% in each region of the mitochondrial genome, the similarity between I. duplitatus and I. hauseri is relatively high, and I. shangrila and I. subelongatus also regularly agree with phylogenetic relationships. In addition, the frequency of amino acid codon use in I. subelongatus appears to differ from that of other species, suggesting that it may be related to its larch host.

Phylogenetic Status of Ips Species from China Based on COI Gene Analysis
In studying the phylogenetic and taxonomic relationships between Chinese and European and American Ips species based on COI sequences, we found that the species groups in the current classification for Ips are not monophyletic. The Chinese Ips species do not form a monophyletic group, which are distributed among the European and North American species in the phylogenetic tree. The same conclusions were drawn for the European Ips bark beetles [9]. Therefore, the phylogenetic relationships of Ips species on the same continent are not more distant from each other than those on other continents, i.e., geographic distance is not a determining factor for phylogenetic relationships. From the early Miocene to the Pliocene, East Asia and western North America were connected by the Bering land bridge, which was covered with temperate conifer forests [89]. Conifer-related insect groups such as the bark beetles Ips may have established widespread populations on these continents at this time. Multiple invasions between American and Eurasian Ips species, as well as the fact that most Ips species also interbred, could account for polyphyletic groups [90]. The phylogenetic relationships of Ips species are related to their hosts, with sibling host species being phylogenetically closely related. Of course, different Ips species may share the same host, which may be related to adaptive evolution between Ips beetles and hosts. Ips borealis Swaine 1911, I. hunteri Swaine 1917, I. pilifrons Swaine 1912, and I. tridens Mannerheim 1852 are morphologically very similar and difficult to distinguish, and our results suggest that these beetles are closely related and can be resolved by phylogenetic analyses. We hypothesize that I. knausi Swaine 1915, I. emarginatus LeConte 1876, and I. sexdentatus Boerner 1767 are basal species of the genus Ips. To date, mitochondrial genome sequences of Ips species are relatively sparse and incomplete, limiting phylogenetic analysis to distinguish monophyletic relationships and determine phylogenetic relationships. Stauffer et al. [13] hypothesized that speciation of I. acuminatus was accompanied by host switching from Picea abies to Pinus silvestris L.; however, we hypothesize that Pinus trees are the original hosts of bark beetles and that they were also the hosts of the basal species Ips latidens, Ips spinifer, Ips mannsfeldi, and Ips nobilis, which is consistent with the conclusion of Cognato et al. [9]. The genus Ips split from the genus Dendroctonus 60 Mya ago [91], and the genus Picea gradually separated more recently, resulting in different Ips species adapting to host switching and gradually forming the present status quo in the process of host dispersal and evolution. The host switch of Dendroctonus spp. from Pinus-associated ancestors to Larix and Picea is thought to have occurred simultaneously in the western Nearctic regions after the late Oligocene/early Miocene [92], which is also consistent with the timing of the Ips divergence. The ancestors of some Ips species did not breed in species of the genus Picea [9], as Ips latidens, Ips spinifer, Ips mannsfeldi, and Ips nobilis are basal to Ips and all breed in the Pinus species. Further studies are needed to refine estimates of the age of these events. Independent studies of different insect species associated with conifers tend to have relatively similar biogeographic and host evolutionary scenarios, which adds to their credibility [93].

Conclusions
In our study, we sequenced the mitogenomes of Ips bark beetles from 19 populations, enriched the baseline data, and examined their phylogenetic and evolutionary relationships. Ips species show adaptive evolution under the influence of the uplift of the Qinghai-Tibet Plateau and the Tienshan Mountains, consistent with the evolution of their hosts, as indicated by interspecific divergence time, genetic distance, selection pressure, and comparison of mitochondrial genomes. The phylogenetic relationship between species and the evolutionary relationship of pheromones are also discussed from a new perspective, namely the source of pheromone synthesis, which provides a new idea to study their relationships. Our results provide a better understanding of the mechanisms of adaptive evolution in relation to the environment, hosts, and phenotypes of bark beetles.

Supplementary Materials:
The following supplementary information can be downloaded from https://www.mdpi.com/article/10.3390/biology11030384/s1. Figure S1: Phylogenetic tree (ML) of Ips geographic populations based on SNP. Numbers above or below branches indicate bootstrap values, Figure S2: Tree of divergence times of Ips species/populations. The numbers above the branches indicate the divergence time (in millions of years) of the branches. The middle solid line indicates the mean population size, and the gray area below and above the solid line indicates the 95% confidence interval, Figure S3: Circular maps of six newly sequenced bark beetles. Genes that were transcribed clockwise are shown outside the circle, while genes that were transcribed counterclockwise are shown inside the circle, Figure S4: Amino acid usage of protein-coding genes of the mitochondrial genome in the genus Ips, Table S1: The list of sample information of mitochondrial genomes of 26 published species of Scolytidae, Table S2: Ips species used in this study, Table S3: Base composition in mitochondrial genomes of bark beetles, Table S4: Organization of the mitochondrial genome of I. typographus, Table S5: Organization of the mitochondrial genome of I. subelongatus, Table S6: Organization of the mitochondrial genome of I. shangrila, Table S7: Organization of the mitochondrial genome of I. duplicatus, Table S8: Organization of the mitochondrial genome of