Adult Monochamus sutor (Linnaeus), Monochamus guerryi Pic, Monochamus galloprovincialis (Olivier), and Monochamus latefasciatus (Breuning) species collected information was shown in Table 1, and preserved in anhydrous ethanol (purity ≥ 99.7%) at −20°C. Specimens were first identified according to external morphological characters before the DNA extraction (Sun 2019; Silva-Pinheiro et al. 2022). A small amount of muscle tissue was removed from the chest of each four longhorn beetles for total DNA extraction. Total genomic DNA was extracted using TIANamp Genomic DNA Kit (TIANGEN, Beijing, China) following the manufacturer’s instructions. The cox1 sequences were amplified with the primers LCO149/HCO2198 and then retrieved from the National Center for Biotechnology Information database to identify from the molecular level (Folmer 1994).

DNA samples of four longhorn beetles were sequenced by next-generation sequencing at Biomarker Technologies Co., Ltd. (Beijing, China) to obtain complete mitochondrial genes. An Illumina TruSeq library was generated, with average lengths of 350 bp and 250 bp paired-end reads sequenced on Illumina Hiseq 2500 platform (San Diego, CA). After filtering using the NGS QC Toolkit v2.3.3 (Patel et al. 2012), clean reads with high quality were used for subsequent genome assembly and identified through MitoZ 3.6 (Meng et al. 2019), with the mitochondrial genomes of M. saltuarius and Monochamus urussovii (Fischer-Waldheim) (OP169419 and OP169420) being used as references.

Protein-coding genes (PCGs) and ribosomal RNA (rRNA) genes were annotated by aligning with known homologous Cerambycidae species mitogenome sequences from GenBank. Additionally, the secondary structures of transfer RNA (tRNA) genes were identified using tRNAscan-SE Search Server v2.0 (http://lowelab.ucsc.edu/tRNAscan-SE/) (Lowe et al. 2016). The limits of PCGs and rRNA genes were determined by comparing them with analogous genes found in other Cerambycidae species and by referencing the positions of transfer RNA (tRNA) genes. To affirm the precision of the deduced locations, the nucleotide sequences of 13 PCGs were translated into amino acids through the utilization of MEGA X (Kumar et al. 2018). The mitogenomes were visualized using the proksee website (https://proksee.ca/).

The mitogenomes of 89 species obtained from the NCBI (https://www.ncbi.nlm.nih.gov) and 4 newly sequenced species were used for analysis (Table S1). Subsequently, an online software tool called Interactive Tree Of Life (iTOL) v 6.8 (Letunic et al. 2019) was employed to visualize the gene order of all species with different gene represented in distinct colors for enhanced clarity and comparison. Mitogenome rearrangement events were completed through the manual comparison, with Drosophila yakuba Burla as a reference (Clary et al. 1985). Pairwise comparisons of DNAs between M. sutor and other longhorn beetles mitochondrial genomes were performed using BLASTN searches in the CGView Comparison Tool (CCT) (https://github.com/paulstothard/cgview_comparison_tool). The base composition, amino acid usage, codon usage, and relative synonymous codon usage (RSCU) were analyzed through MEGA11/ PhyloSuite v1.2.3 (Tamura et al. 2021; Xiang et al 2023). The formula: AT-skew = (A−T) / (A+T) and GC-skew = (G−C) / (G+C) was calculated manually to evaluate nucleotide composition bias (Perna et al. 1995). The MEGA 11 (Tamura et al. 2021) software was utilized to compute the average genetic distances employing the Kimura-2-parameter model across 13 PCGs. Synonymous (Ks) and non-synonymous (Ka) substitution rates among PCGs were calculated by DNASP v5.10 (Librado et al. 2009). The effective number of codons (ENC) and codon bias index (CBI) for PCGs were also calculated by the software DNASP 5.10 (Librado et al. 2009).

The mitochondrial genome data were obtained from the NCBI nucleotide database. A total of 86 species in 10 groups were selected from the longhorn beetles. Additionally, the dataset incorporated the mitochondrial genomes of the four newly sequenced species.

After filtering the redundant sequence, 13 PCGs and 2 rRNAs were extracted by PhyloSuite v1.2.3 (Xiang et al. 2023) for downstream analysis. Multiple sequence alignments were conducted for amino acids, PCGs, and rRNA gene using MAFFT v7.313 (Katoh et al. 2013). For the alignment of PCGs, the L-INS-i algorithm was employed, while the G-INS-i algorithm was utilized for the alignment of rRNA genes (Katoh et al. 2013). The alignments of PCGs were imported into the MACSE v2.03 (Ranwez et al. 2018) for optimization. Sequences were concatenated using PhyloSuite version 1.2.3 to create four distinct datasets: (i) PCG123: PCGs with all three codon positions; (ii) PCG123R: PCG123 dataset plus two rRNA genes; (iii) PCG12R: PCG12 dataset plus two rRNA genes; (iv) PCGAA: amino acid sequences translated from PCGs. ModelFinder (Kalyaanamoorthy et al. 2017) and PartitionFinder2 (Lanfear et al. 2017) were used to select optimal partitioning strategies and evolutionary models for PCGsRNA and AA datasets. The heterogeneity and substitution saturation of each dataset were assessed using the AliGROOVE v1.08 (https://github.com/PatrickKueck/AliGROOVE) (Kück et al. 2014) and DAMBE v5.5.29 (Xia 2013), respectively.

Maximum likelihood (ML) phylogenies and Phylobayes MPI 1.5a were used to construct phylogenetic trees mainly about Cerambycidae on PhyloSuite v1.2.3 and CIPRES Science Gateway (https://www.phylo.org). Seven species with complete sequences from other families that similar to Cerambycidae were selected as outgroups (Eucryptorrhynchus scrobiculatus Motschulsky, Eucryptorrhynchus brandti (Harold), Curculio davidi Fairmaire, Sympiezomias velatus (Chevrolat), Naupactus xanthographus Germaer, Aegorhinus superciliosus (Guérin), Sitophilus zeamais Motschulsky). For ML analysis, the IQ-TREE v. 1.6.8 (Minh et al. 2020) was used to reconstruct the ML phylogenetic tree under the Partition Mode (Table S2). A standard bootstrap approximation approach was selected, and 1000 replicates were performed to assess the tree’s robustness. Phylobayes MPI on XSEDE (version 1.8c) was utilized through the CIPRES Science Gateway (https://www.phylo.org) to construct Bayesian inference (BI) phylogenies. The site-heterogeneous mixture model CAT+GTR was employed for both nucleotide and amino acid datasets (Lartillot et al. 2013). Two independent chains, a Markov chain, and a Monte Carlo chain, were run after removing a fixed point from a sequence. The chains were stopped after two runs of 30000 generations converged satisfactorily, indicated by a maxdiff value of less than 0.1. The initial 25% of trees generated were discarded as burn-in.

Complete mitogenomes of four Lamiinae species were newly sequenced: M. sutor, M. guerryi, M. galloprovincialis and M. latefasciatus, which ranged in size from 15759 (M. guerryi) to 15929 (M. latefasciatus) base pairs and each contained 37 genes (13PCGs, 22 tRNA and 2 rRNA) (Fig. 1). These mitogenomes exhibit a circular and double-stranded structure. The majority strand (J-strand) encodes 23 genes, including 9 PCGs (atp6, atp8, cytb, cox1, cox2, cox3, nad2, nad3 and nad6), and 14 tRNA genes. The minority strand (N-strand) encodes 14 genes, comprising 4 PCGs (nad1, nad4, nad4L and nad5), 8 tRNA genes and 2 rRNA genes (rrnL and rrnS). The gene arrangement in these four new species is consistent with that of other published Cerambycidae mitogenomes and follows the ancestral insect mitogenome arrangement. However, Vesperus conicicollis Fairmaire & Coquerel exhibited gene rearrangement, with nad6 positioned ahead of trnP. And the translocation of trnL1 and trnL2 and a local inversion of rrnS were observed in Apomecyna saltato (Fabricius) and Arhopalus unicolor (Gahan), respectively (Fig. 10). These gene rearrangements mainly occur in the replication origin, contributing to species evolution acceleration. Overall, Cerambycidae mitochondrial genomes display a high degree of conservation, indicating relative stability in mitochondrial function (Liu et al. 2023).

Figure 1. 

Complete mitochondrial genomes of M. sutor (A), M. guerryi (B), M. latefasciatus (C), and M. galloprovincialis (D). The transcriptional direction is indicated by arrows.

The mitochondrial genomes of the longhorn beetles display a noticeable AT bias, with an A+T content ranging from 67.9% (Prioninae, Agrianome spinicollis (MacLeay)) to 81.2% (Lepturinae, Sachalinobia koltzei (Heyden)). Interestingly, only four species exhibit an A+T content exceeding 80%. Lepturinae stands out with the highest A+T content, averaging 79.5%, while Prioninae has the lowest, averaging 70.8%. The AT-skew values are predominantly positive, ranging from –0.016 (Lamiinae, Pterolophia sp. ZJY-2019) to 0.125 (Prioninae, Priotyrannus closteroides Thomson), while GC-skew values are mainly negative, spanning from –0.286 (Prioninae, P. closteroides) to –0.136 (Cerambycinae, Xystrocera globose (Oliver)). Notably, the newly sequenced Lamiinae species all have A+T contents higher than the average, with positive AT skew. Although the nucleotide compositions across the 10 groups are generally similar, there are differences in A+T content, AT-skew, and GC-skew (Table S3; Fig. 2).

Figure 2. 

Three-dimensional scatter plots of the AT-skew, GC-skew, and A + T content of 86 mitochondrial genomes. Balls of different colors correspond to different groups.

To explore the base content characteristic in the coding region, we specifically analyzed the PCGs of mitogenomes. The results revealed a significant AT bias in the entire PCGs, with A+T content ranging from 65.42% (Prioninae, A. spinicollis) to 79.73% (Lepturinae, S. koltzei). The Lepturinae again exhibits the highest A+T content, averaging 77.9%, while the Prioninae has the lowest, averaging 67.9%, the same trends observed in the entire mitochondrial genome. The AT-skew values are negative across all PCGs, while most GC-skew values are also negative. This AT bias is particularly prominent at the 3^rd codon position, consistent with a high mutational pressure toward AT nucleotides at this position (Meng et al. 2019). Notably, the 1^st condon position shows a typical T, G bias, while the 2^nd and 3^rd codon positions show a typical T, C bias. Among the four newly sequenced species, both the entire PCGs and the second positions of these genes exhibit higher A+T content compared to the corresponding average values. Lepturinae maintains the highest A+T content, while Prioninae retains the lowest (Table S3). Sequences with high A+T content are generally less stable and more mutation-prone, potentially indicating adaptation to specific environments, such as extreme conditions, high altitudes, or high latitude areas (Song et al. 2016). This suggests that species within the Lepturinae may have a competitive advantage in adapting to particular habitats compared to those within the Prioninae. The AT content in the genus Monochamus is generally similar, and similar AT bias may indicate shared evolutionary trajectories, especially in the case of the M. alternatus, M. galloprovincialis, and M. urussovii, where A+T content has reached a consistent level. This information enhances our understanding of mitochondrial genome functionality, evolution, and adaptive relationships within different environments.

Furthermore, there is a clear AT bias in RNA genes, with AT-skew and GC-skew values mainly positive in tRNA but predominantly negative in rRNA. Among the four newly sequenced beetles, the AT content in RNA genes is higher than the average, while the GC content is lower than the average (Table S3).

As depicted in Figure 3, the species obtained from the BLASTN analysis using the CCT and the mitochondrial genome of M. sutor were used as the reference species. The sequence similarity between M. sutor and other mitochondrial genomes ranged from 75.80% to 97.6%. Notably, M. sutor exhibited its closest relationship with M. galloprovincialis, followed by M. urussovii, M. saltuarius, and M. alternatus. Furthermore, the control region displayed the highest degree of variation, and there were varying levels of variation observed in genes such as ND2, ND4, ND4L, ND5 and ND6. These differences pointed to the diversity among species and variations in evolutionary rates. In contrast, the cox1 gene exhibited a lower degree of variation and was relatively more conservative.

Figure 3. 

Graphical map of the BLASTN results showing the nucleotide identity between the M. sutor mitochondrial genome and that of 85 other species listed in Table S1. CCT arranges the BLASTN results such that the sequence most similar to the reference (M. sutor) is placed closer to the outer edge of the map.

After analyzing the PCGs of 86 mitogenomes, we observed that the starting codon usage bias was quite similar at the species level when ATN codons were used as the start codons. Most species predominantly employed ATT or ATG as the start codon, while some species used ATA, ATC, and a few even utilized TTG and GGA. As for stop codons, TAA, TAG, and TTA were used in these mitogenomes, with TAA being the most frequently employed. At the gene level, most PCGs preferred ATT and ATG as start codons and still commonly used TAA as a stop codon. However, there was diversity in the types of start codons for NADH dehydrogenase genes, with some less common ones like TTG. Additionally, 8 genes exhibited incomplete stop codons in PCGs, with only atp6, atp8, nad6, and nad1 being exceptions. In most cases, the incomplete stop codon was “T”, and very rarely “TA” (Table S5). These incomplete stop codons are relatively common in insect mitochondrial genomes (Wu et al. 2023). Overall, the start codon patterns align with the frequent use of NNU and NNA in each codon.

Relative synonymous codon usage (RSCU) values of all 10 groups and the four new sequenced species were plotted in Fig. 4. The codon usage of PCGs exhibited a strong bias towards NNU and NNA codons. UUA (Leu2) had the highest RSCU in all samples, indicating that it had undergone strong selection during evolution. Conversely, the RSCU values for CGC (Arg), ACG (Thr), and CUG (Leu1), which end with G and C, were relatively low. Notably, Leu2 and Ser2 consistently possessed the highest RSCU values in all groups. These results demonstrated a connection between RSCU values and codon base selection bias (Table S6). Commonly used codons consisted of A and T, contributing to a higher A+T content in the mitochondrial genome. Furthermore, almost all frequently used codons ended with A or T, reinforcing the bias toward AT content in the third codon position (Song et al. 2016). In Dorcasominae, Prioninae, Cerambycinae and Disteniinae, the RSCU value of Leu2, particularly UUA in Leu2, distinctly decreased, while that of Leu1 increased. This could be a contributing factor to the lower AT content in these four groups (Zheng 2022). Overall, the differences in codon usage among different groups were minimal, suggesting that they might share a similar genetic evolutionary background.

Figure 4. 

Relative synonymous codon usage (RSCU) in protein-coding genes of 10 groups and four new sequenced species. Codon families are indicated on the x-axis, along with different combinations of synonymous codons that code for specific amino acids. RSCU is shown on the y-axis.

In the 10 groups, the frequency of amino acid usage was relatively consistent. Overall, Leu was the most frequently used amino acid, followed by Ile, Phe, Met, and Ser in that order. However, some variations were observed in Spondylidinae, Prioninae, Cerambycinae and Disteniinae. In these groups, there was a general trend of Ser being more frequently used than Met (Table S6; Fig. 5). To sum it up, the frequency of amino acid usage in PCGs of Cerambycidae exhibited stability, with only minor differences in specific subfamilies.

Figure 5. 

Amino acid usage in the mitochondrial genome. The different colored dashed lines represent different groups, the x-axis represents different amino acid types, and the y-axis represents the usage of amino acids.

The effective number of codons (ENC) and the codon bias index (CBI) in these species ranged from 37.03 (Lepturinae, S. koltzei) – 51.2 (Disteniinae, Typodryas sp. N143) for ENC, and from 0.275 (Prioninae, P. closteroides) – 0.643 (Lepturinae, S. koltzei) for CBI. Notably, Spondylidinae and Orsodacinade had the lower average ENC value of 43.79 and 43.61, while Prioninae had the highest average ENC value of 50.32. As for CBI, Prioninae had the lowest average value of 0.33, while Spondylidinae, Lepturinae and Orsodacinade had the higher average CBI value of 0.50 (Table S4). The ENC and CBI results showed that Spondylidinae and Orsodacinade species exhibited a higher codon bias, suggesting that synonymous codons were more frequently used in their genomes with a lower codon selection, implying relatively lower genomic complexity. These Spondylidinae and Orsodacinade species appear to possess strong adaptability to specific ecological environments (Li et al. 2023; Wu et al. 2023). Among the four newly sequenced species, M. sutor showed a relatively low ENC value of 38.95 and a relatively high CBI value of 0.59. These values indicate that M. sutor has strong environmental adaptability and more uniform codon usage (Li et al. 2023; Wu et al. 2023). Furthermore, the correlation analysis depicted in Fig. 6 revealed that ENC and CBI were positively and negatively correlated with GC3s, respectively. This suggests that the genomic G+C content plays a crucial role in determining codon preference across different species (Xu et al. 2020; Yan et al. 2021). The correlations between GC3 and both ENC and CBI are high across different groups except Prioninae. Higher GC3 levels in the genome may exert stronger selective pressure, leading to more uniform codon usage (Wei et al. 2014). Moreover, elevated GC3 levels could result in specific functional genes displaying pronounced codon usage preferences with the principles of neutral mutational theories (Andersson et al. 1996).

Figure 6. 

Nucleotide compositions of the mitogenomes. A Scatter plot of the GC content of 3^rd codon sites versus ENC; B Scatter plot of the GC content of 3^rd codon sites versus CBI.

The number of intergenic regions in the mitochondrial genomes of the entire 10 groups ranged from 2 (Cerambycinae, Ceresium sinicum ornaticolle Pic) to 14 (Lamiinae, Lamiinae sp. 1 ACP-2013 and Prioninae, P. closteroides). On the other hand, gene overlap regions varied from 5 (Disteniinae, Distenia punctulatoides Hubweber) to 25 (Lamiinae, Oberea diversipes Pic). These gene overlaps were generally small, often spanning from 1 to 7 bp. The largest overlap regions were frequently observed between atp6 and atp8 (Wei 2009; Zheng et al. 2022), as well as between nad4 and nad4L. The lengths of the intergenic regions displayed significant variation, with the longest intergenic region reaching up to 1813 bp (Lamiinae, Jamesia sp. KM-2017). Among the mitochondrial genomes analyzed in the beetle family, a total of 69 instances of a 7 bp intergenic region between atp6 and atp8 were observed, and 62 instances of a 7 bp gene intergenic region between nad4 and nad4L were identified (Table 2; Table S5).

The genes atp6, atp8, nad4, and nad4L are likely to play crucial roles in essential cellular processes like energy production and oxidative phosphorylation (Pogoda et al. 2018; Yang et al. 2019). The presence of special intergenic regions between these genes may reflect their functional relationships and coordination (Zhang et al. 2013). Additionally, these distinctive intergenic regions could have implications for the adaptation of insect species to specific environments or ecological conditions (Yévenes et al. 2022).

As depicted in Fig. 7, within these groups , certain mitochondrial genes displayed distinct rates of evolution. Specifically, nad6 (0.354), atp8 (0.335), nad3 (0.312), and nad2 (0.300) were identified as fast-evolving genes, whereas cox1 (0.191), cox2 (0.202), nad1 (0.207) and atp6 (0.207) exhibited a slower rate of evolution. Genes characterized by smaller genetic distances are indicative of relatively stable functions and fewer accumulated mutations within their sequences (Nei 1978). In groups, it was noted that nad6 consistently exhibited the higher genetic distances, while cox1 consistently displayed the lower genetic distances, except for Spondylidinae. In Cerambycidae, all calculated values of ka/ks were less than 1, ranging from 0.067 (cox1) to 0.441 (atp8). This indicates a higher likelihood of retaining synonymous mutations within these genes. Notably, genes such as atp8 (0.441), nad4L (0.318), and nad6 (0.300) exhibited a high rate of evolution, whereas cox1 (0.067), cox2 (0.094) and cytb (0.110) showed a low rate of evolution. Further analysis of ka/ks revealed that atp8 and nad6 had higher ka/ks ratios, suggesting that these genes have experienced relaxed selection pressures. In contrast, cox1 and cox2 displayed lower ka/ks ratios, indicating a history of strong purifying selection acting upon them (Sun et al. 2020). Combining the findings from both Ka/Ks and genetic distance analyses, it became evident that cox1 and cox2 are exceptionally conserved genes, relatively unaffected by natural selection (Lynch et al. 2000; Wertheim et al. 2015). The cox1 gene, in particular, is employed as a DNA barcode for species identification, leveraging the high level of conservatism in its sequence (Yévenes et al. 2022; Fonseca et al. 2022). This utility is particularly valuable in biodiversity surveys and species identification, aiding in distinguishing species that are otherwise challenging to discern. Conversely, genes such as atp8, nad4L, and nad6 play crucial roles in cellular respiration and energy metabolism pathways (Hui et al. 2018; Yang et al. 2021). They necessitate ongoing variation and evolution to adapt to diverse ecological environments, resulting in greater sequence variability.

Figure 7. 

Mean Ka/Ks (A, B) and Genetic distances (C, D) (Kimura-2-parameter) of mitochondrial gene sequences in groupes. ±SD, hollow dots; Median, middle line inside each box; IQR (interquartile range), the box containing 50% of the data;

All four newly sequenced species exhibited a complete set of 22 tRNA genes, with entire tRNA lengths 1441 (M. latefasciatus), 1444 (M. guerryi), 1452 (M. sutor), 1452 (M. galloprovincialis). When individually examining these tRNAs, their lengths ranged from 61 to 70 in M. latefasciatus, from 63 to 70 in M. guerryi, from 63 to 69 in M. sutor, and from 63 to 69 in M. galloprovincialis, respectively. Except for trnS1, all tRNAs exhibited the characteristic cloverleaf structure, where the dihydrouridine arm of trnS1 formed as a single loop. This peculiar structure of trnS1 may arise from abnormalities during transcription or post-transcriptional processes, which could subsequently affect the development of its typical cloverleaf structure (Feng et al. 2020) (Table S5; Fig. 8). Such an occurrence has been frequently observed within the longhorn beetle genus Monochamus. The absence of the standard cloverleaf configuration implies functional versatility for mitochondrial trnS1, potentially enabling diverse roles in various biological processes (Garey et al. 1989; Wolstenholme et al. 1992). Nevertheless, this also indicates a deviation from the traditional functionality associated with a typical tRNA.

The tRNA structures in all four newly sequenced beetles exhibited a cloverleaf-like secondary structure overall; the phenomenon of tRNA-Ser1 lacking the DHU arm was present in all four newly sequenced beetles. Additionally, some degree of base mismatch occurred, with the most common type being the G-U mismatch. Specifically, there were 25 (21 G-U, 3 U-U, 1 A-G) mismatched base pairs in M. latefasciatus, 24 (20 G-U, 3 U-U, 1 A-G) mismatched base pairs in M. guerryi, 20 (15 G-U, 4 U-U, 1 A-G) mismatched base pairs in M. sutor, 21 (16 G-U, 4 U-U, 1 A-G) mismatched base pairs in M. galloprovincialis (Fig. 8). The absence of the DHU arm and the occurrence of base mismatches in mitochondrial tRNA may result from various factors, including adaptive changes during evolution, genetic mutations, and selective pressures (Ji et al. 2023). The DHU arm played a role in the stability and function of the tRNA molecule, and its absence could potentially impact the stability and proper folding of tRNA-Ser1, thereby affecting its role in protein synthesis. Base mismatches may lead to mutations in genes encoding proteins, accumulating across generations and influencing the genetic traits of the population (Acuna-Hidalgo et al. 2016). Overall, the structural variations observed in these four newly sequenced beetles were consistent and largely consistent with the tRNA structures found in most species.

Figure 8. 

Putative secondary structures of tRNAs from the M. sutor, M. guerryi, M. galloprovincialis, and M. latefasciatus mitochondrial genome. Dashes indicate the Watson–Crick base pairs; dots indicate wobble GU pairs; and the other non-canonical pairs are not marked.

The rRNA components in the mitochondrial genomes of the four newly sequenced species encompass both lrRNA and srRNA. Specifically, lrRNA is situated between tRNA-Leu and tRNA-Val, while srRNA is positioned between tRNA-Val and the control region. The length of lrRNA genes ranged from 1329 bp in M. guerryi to 1334 bp in M. latefasciatus, while the sizes of srRNA genes varied from 775 bp in M. latefasciatus to 781 bp in M. guerryi, M. galloprovincialis, and M. sutor (Table S5).

The substitution saturation index (Iss) for the six nucleotide sequence datasets was notably lower than the critical Iss values (Iss.cSym and Iss.cAym) for both symmetrical and asymmetrical trees. This observation strongly indicated the absence of substitution saturation, providing clear evidence of a robust phylogenetic signal (all Iss < Iss.cSym or Iss.cAsym, P < 0.05) (Table 3, Fig. S1) (Yan et al. 2021).

To assess nucleotide divergence heterogeneity, pairwise comparisons were conducted within a multiple sequence alignment. The results revealed that datasets involving the third codon position (PCG123, PCG123R) exhibited lower heterogeneity in sequence composition. In comparison, datasets with the third codon position removed (PCG12, PCG12R) displayed intermediate heterogeneity levels. The AA dataset exhibited the highest heterogeneity. It’s worth noting that some species, such as C. claviger, Dynamostes audax Pascoe, H. nigromaculata, Oberea formosana Pic, and V. conicicollis, exhibited greater heterogeneity due to incomplete sequences caused by missing data (Fig. 9) (Li et al. 2020).

Figure 9. 

Heterogeneity of the sequence composition of the mitochondrial genomes in different datasets. The pairwise Aliscore values are indicated by colored squares. Darker colors indicate full random similarity, and lighter colors indicate the opposite. All taxon names are listed on the left side of the heat map.

We selected a total of 86 species of longhorn beetles and analyzed four datasets using IQ-TREE and Phylobayes two methods. The results indicated that most nodes within the phylogenetic trees displayed strong support values, indicating robust branching patterns (Fig. 10, Fig. S2). The overall topology of the phylogenetic trees remained largely consistent, with only a few branches showing instability.

Figure 10. 

Phylogenetic tree inferred from the PCG123 and AA dataset using maximum likelihood (ML) and arrangements of genes in the mitochondrial genomes. Numbers beside each node are support values. The values indicate ML bootstrap values (BPs).

In these phylogenetic trees, results obtained from the AA dataset and the PCG123 dataset both constructed using IQ-TREE are notably similar with relatively high node support (Fig. 10, Fig. S2). Through analyzing all the phylogenetic trees, we observed that Cerambycinae and Prioninae are sister groups in all of the tree. And the position of V. conicicollis within Vesperidae exhibited notable instability across all phylogenetic trees. Specifically, in the results generated by IQ-TREE, V. conicicollis appeared to be distantly related to the other sequenced longhorn beetles. In AA dataset, the relationships among Disteniinae, Cerambycinae and Prioninae remained consistent across both methods: (Disteniinae+(Cerambycinae+Prioninae)) (BPs = 99, PPs = 0.99). The Orsodacnidae, Oxypeltinae and Vesperidae relationship is unstable. In PCG12RNA, Dorcasominae were highly supported, but the phylogenetic positions of them were unstable. This may be attributed to the extremely limited amount of data available. In PCG123 dataset and PCG123RNA dataset, Orsodacnidae and Oxypeltinae both form sister group but with different node support. In a comparative analysis of the results from eight phylogenetic trees, we observed that the tree constructed by the IQ-TREE method using the AA dataset is better, demonstrating relatively high node support and strong consistency of phylogenetic relationships. The results in this tree revealed four major clades, which are consistent with Nie et al., except for the ‘lamiinae’ clade in our study is ((Lepturinae + Spondylidinae) + Lamiinae), while it is ((Lamiinae + Spondylidinae) + Lepturinae) in Nie et al. (Nie et al. 2020) (Fig. 10, Fig. S2).

Through constructing phylogenetic tree and the analysis of mitochondrial genome data presented in this study, we have presented following perspectives on some contentious issues. The consistent clustering of Spondylidinae and Aseminae in the phylogenetic trees supports the idea that these two groups may be treated as a single entity, Spondylidinae, with some cases of synonymous nomenclature (Linsley et al. 1961; Gressitt et al. 1970; Nakamura et al. 1981; Napp et al. 1994). In addition, some other scholars believe that Oxypeltidae is a family at the same taxonomic level as Cerambycidae (Monne 2012). The morphological identification of Orsodacnidae suggests that it is an independent family within the order Coleoptera (Li 2018). In the phylogenetic results of current study, Oxypeltidae and Orsodacnidae consistently form a single branch. This indicates that they are relatively closely related and their taxonomic status may be the same. However, in our phylogenetic tree analysis, there are still some unresolved issues. For instance, the phylogenetic relationships of certain groups remain unresolved (Vesperidae, Oxypeltinae, Orsodacnidae, and Dorcasominae), and the support values at key nodes are not notably high. Despite employing various methods and constructing different datasets, the results are consistently similar. These challenges are also observed in the study by Nie et al. (Nie et al., 2020). Therefore, we propose that future research should focus on expanding mitochondrial genome sequencing to include a more extensive range of species, thereby enriching our sample size. Additionally, exploring alternative tree-building methods may be beneficial to identify an approach better suited for the phylogenetic analysis of the longhorn beetles. It is essential to note that, although we, like Nie et al. have related to research on phylogeny of the longhorn beetles, our methods for constructing the phylogenetic tree and the focal points of our research are different. Nie et al. chose a more extensive set of species, offering new insights into the phylogeny of longhorn beetles. In contrast, our objective is to employ diverse methods to analyze mitochondrial genome data from multiple perspectives. Our aim is to conduct comparative analyses of these datasets, providing additional support for the phylogeny of the longhorn beetles. And the results we obtained in the phylogenetic section are generally consistent.

In two methods employed for the analysis, the four newly sequenced beetles consistently grouped within the Lamiinae. This clustering demonstrates their close genetic affinity with other species within the Lamiinae. The proximity of closely related species within the same subfamily is evident in the phylogenetic trees, indicating a shared evolutionary history and genetic similarity among these species (Fig. 10, Fig. S2). In 2022, Shi et al. conducted a comparative analysis of mitochondrial genomes within Lamiinae. Although both studies focus on the mitochondrial genome analysis of longhorn beetles, it’s crucial to note that our research targets differently. While the study by Shi et al. concentrates on Lamiinae, the scope of this article encompasses the entire longhorn beetles families or subfamilies. Additionally, new content has been incorporated into the analysis methods, including analyses such as CBI and ENC. In our study, the research results regarding Lamiinae align with those of Shi et al. Apart from a few individual species that are categorized under Lamiinae in the Shi et al., but our research indicates that they belong to Spondylidinae, such as A. unicolor. Furthermore, compared to other groups, we observed that the phylogenetic relationships within Lamiinae are relatively stable, and they exhibit a close affinity with Spondylidinae and Lepturinae.