Complete Mitochondrial Genomes and Phylogenetic Positions of Two Longicorn Beetles, Anoplophora glabripennis and Demonax pseudonotabilis (Coleoptera: Cerambycidae)

Anoplophora glabripennis (Motschulsky, 1854) and Demonax pseudonotabilis Gressitt & Rondon, 1970 are two commonly found longicorn beetles from China. However, the lack of sufficient molecular data hinders the understanding of their evolution and phylogenetic relationships with other species of Cerambycidae. This study sequenced and assembled the complete mitochondrial genomes of the two species using the next-generation sequencing method. The mitogenomes of A. glabripennis and D. pseudonotabilis are 15,622 bp and 15,527 bp in length, respectively. The mitochondrial gene content and gene order of A. glabripennis and D. pseudonotabilis are highly conserved with other sequenced longicorn beetles. The calculation of nonsynonymous (Ka) and synonymous (Ks) substitution rates in PCGs indicated the existence of purifying selection in the two longicorn beetles. The phylogenetic analysis was conducted using the protein-coding gene sequences from available mitogenomes of Cerambycidae. The two species sequenced in this study are, respectively, grouped with their relatives from the same subfamily. The monophyly of Cerambycinae, Dorcasominae, Lamiinae, and Necydalinae was well-supported, whereas Lepturinae, Prioninae, and Spondylidinae were recovered as paraphyletic.

The phylogeny and early evolution of Cerambycidae have been comprehensively reviewed by Haddad & Mckenna (2016) [13]. The phylogeny of longicorn beetles, especially the monophyly of Cerambycidae s.s. and s.l., as well as the subfamily and tribe-level relationship, remains debatable due to the high species richness and highly variable morphological characters [5,14,15]. Haddad et al. (2018) [5] reconstructed the higher-level phylogeny of Cerambycidae with anchored hybrid enrichment of nuclear genes. Their results recovered a monophyletic Cerambycidae s.s. in most analyses and a polyphyletic

Mitogenome Assembly, Annotation, and Analyses
Before the assembly, the high-quality reads were trimmed again using BBDuk with default settings implemented in Geneious Prime [26]. The high-quality reads of A. glabripennis and D. pseudonotabilis were, respectively, mapped to the reference mitogenome of the previously sequenced A. glabripennis (NC_008221) [23] and amplified bilaterally by Geneious Prime [26], with the parameters set as follows: 95% minimum overlap identity, 50 bp minimum overlap, and maximum ambiguity as 4. The completeness of each circular mitogenome was confirmed when both ends of the final assembled contigs overlapped (100% coverage). The assembled mitogenomes of A. glabripennis and D. pseudonotabilis were deposited in GenBank under the accession numbers OP096420 and OP096419, respectively.
The two mitogenomes were annotated in the MITOS web server [27]. The resultant gene boundaries of the protein-coding genes (PCGs) were checked manually by the NCBI's ORF Finder (https://www.ncbi.nlm.nih.gov/orffinder/, accessed on 9 August 2022). The location and secondary structures of the transfer RNA (tRNA) and ribosomal RNA (rRNA) genes were predicted and visualized by MITOS. The mitogenome structure and GC skews were visualized by the CGView Server [28]. The nucleotide composition, skews, codon usage, and relative synonymous codon usage (RSCU) were calculated by MEGA11 [29]. The synonymous substitution rate (Ks) and nonsynonymous substitution rate (Ka) were calculated using KaKs_Calculator v2.0, with the mitogenome of Aoria nigripes (Baly, 1860) (Chrysomelidae) as the outgroup [30,31]. The alignment file of each PCG was uploaded to the Datamonkey web server for a more thorough exploration of selective pressure in PCGs of the two newly sequenced mitogenomes. BUSTED (Branch-Site Unrestricted Statistical Test for Episodic Diversification) was used to test whether each PCG has experienced positive selection [32]. FEL (fixed-effects likelihood) was employed to infer site-specific Ks and Ka values and detect the following four types of sites in each PCG: diversifying sites, purifying sites, neutral sites, and invariable sites [33]. Tandem repeats in the control regions were identified using the Tandem Repeats Finder web server [34]. The stem-loop structures in the control region were predicted by the Mfold web server with default settings [35].

Phylogenetic Analyses Methods
The phylogenetic relationships were reconstructed based on the nucleotide sequences of 13 PCGs derived from 186 mitogenomes of Cerambycidae s.l. (Table 1). Overall, 30 of the 186 mitogenomes were originally unannotated in GenBank; they were re-annotated by MITOS and manual homology alignments in this study. Other mitogenomes from GenBank that had incomplete set of 13 PCGs or incorrect PCG sequences were omitted from the dataset. The mitogenome of A. nigripes (Chrysomelidae) was used as the outgroup [31]. The 13 PCGs were, respectively, aligned using MUSCLE with a codon mode [36], followed by the deletion of stop codons and the concatenation of sequences by SequenceMatrix v1.7.8 [37]. The best-fit partitioning schemes and substitution models for each PCG region were determined by PartitionFinder v2.1.1 using the Bayesian information criterion (BIC) and a greedy search algorithm of all available models [38]. Phylogenies were inferred using maximum-likelihood (ML) and Bayesian inference (BI) methods. The best-fit model was GTR+I+G for two partitioned subsets: one subset included ND1, ND4, ND4L, and ND5; the other subset included the remaining 9 PCGs. IQ-Tree was used to perform the ML analysis under the edge-unlinked partition model for 5000 ultrafast bootstraps as well as the Shimodaira-Hasegawa-like approximate likelihood-ratio test [39][40][41]. The BI analysis was conducted by MrBayes v3.2.7 [42] with four independent Markov chains for 30 million generations and sampled every 100 generations. The first 25% of the trees were discarded as burn-in. FigTree v1.4.4 was used to edit and visualize the phylogenetic trees [43].

Genome Structure and Composition
The assembled complete mitogenomes of A. glabripennis and D. pseudonotabilis are circular DNA molecules of 15,622 bp and 15,527 bp in length (Figure 1), respectively, which is within the range of the sequenced mitogenomes of Cerambycidae in GenBank (Table 1). Due to the presence of a shorter COX1 gene, the newly obtained A. glabripennis mitogenome is slightly shorter than the previously sequenced mitogenome (15,774 bp) based on samples from Hebei Province [31]. Both newly sequenced mitogenomes contain the standard set of 37 mitochondrial genes (13 PCGs, 22 tRNA genes, and 2 rRNA genes) as all other longicorn beetles. The gene order is identical to all other species of Cerambycidae as well as the ancestral mitogenome type of Drosophila yakuba Burla, 1954 [14,44,45]. Among the 37 genes, 23 (9 PCGs and 14 tRNAs) genes are on the majority strand (J-strand), while the remaining 4 PCGs, 8 tRNAs, and 2 rRNA genes are on the minority strand (N-strand).
A total of nine gene overlapping regions were found in the A. glabripennis mitogenome with a total of 29 bp in length, and the longest overlapping sequence (8 bp) was located between trnCys and trnTyr. In the D. pseudonotabilis mitogenome, there are 12 overlapping regions with a total of 21 bp in length, and the longest overlapping sequences were only 4 bp in length. The universally found 7 bp overlapping regions between ATP8 and ATP6, as well as NAD4 and NAD4L in Cerambycidae and many other insects [14,15], are restricted to the overlapping between NAD4 and NAD4L in the A. glabripennis mitogenome, which might be resulted from the different annotation methods. In addition to the overlapping regions, multiple intergenic spacers are scattered throughout both mitogenomes (Tables 2 and 3). The base composition is 38.8% A, 14.2% C, 9.2% G, and 37.8% T for the A. glabripennis mitogenome and 39.7% A, 14.5% C, 10.5% G, and 35.3% T for D. pseudonotabilis. The two mitogenomes are highly skewed towards A and T nucleotides, with an A + T content of 76.6% in A. glabripennis and 75.0% in D. pseudonotabilis (Table 1).
are restricted to the overlapping between NAD4 and NAD4L in the A. glabripennis mitogenome, which might be resulted from the different annotation methods. In addition to the overlapping regions, multiple intergenic spacers are scattered throughout both mitogenomes (Tables 2 and 3). The base composition is 38.8% A, 14.2% C, 9.2% G, and 37.8% T for the A. glabripennis mitogenome and 39.7% A, 14.5% C, 10.5% G, and 35.3% T for D. pseudonotabilis. The two mitogenomes are highly skewed towards A and T nucleotides, with an A + T content of 76.6% in A. glabripennis and 75.0% in D. pseudonotabilis (Table 1).

Protein-Coding Genes
The PCGs have identical arrangement and similar size between the two mitogenomes and also other cerambycids. Most PCGs of the two species start with the standard ATN start codons (ATA, ATC, ATG, and ATT), whereas ND1 of both mitogenomes begins with the special codon TTG (Tables 2 and 3), which was similar to all other published Cerambycidae mitogenomes [14,15]. Most PCGs of each mitogenome have the complete termination codon TAN (TAA, TAT, or TAG), whereas four PCGs (COX1, COX2, ND4, and ND5) of A. glabripennis and four PCGs (COX1, COX3, ND3, and ND5) of D. pseudonotabilis end with an incomplete stop codon T. These incomplete stop codons are considered to be caused by the post-transcriptional polyadenylation [46] and can be completed by the addition of 3 nucleotide residues to the neighboring mitochondrial genes.
The relative synonymous codon usage (RSCU) values indicate the most frequently used codon is TTA (Leu) for both mitogenomes (Figure 2), which appears to be a common feature of other sequenced longicorn beetles [14]. ATP8 of both mitogenomes has the highest A + T content among the 13 PCGs (Tables 2 and 3). The Ka/Ks ratios for each PCG of each mitogenome are calculated to assess the selective pressure of the two cerambycid species ( Figure 3A). The evolutionary rate of ND6 was the highest among the 13 PCGs. The Ka/Ks ratios of all the 13 PCGs calculated by KaKs_Calculator v2.0 were below 1, which suggests the existence of purifying selection in the two species ( Figure 3A). The results of Ka/Ks calculation were similar to a recent mitogenomic work [47], which used DnaSP for the calculation. The gene-wide BUSTED analysis based on the likelihood-ratio test found no evidence of episodic diversifying selection in the PCGs. The site-specific FEL analysis detected ND4 and ND4L each had one codon site under diversifying positive selection at p ≤ 0.1 ( Figure 3B). Nearly one-third of each PCG's codon sites were under purifying selection at p ≤ 0.1. The calculation of KaKs_Calculator v2.0 was consistent with the results of FEL analysis that the PCGs with lower Ka/Ks ratios tended to have more purifying codon sites (Figure 3).

Transfer RNAs, Ribosomal RNAs, and Control Region
The two mitogenomes both contain the complete set of 22 tRNA genes typical of metazoan mitogenomes. These tRNAs range in size from 62 to 69 bp, which was consistent with previously sequenced mitogenomes of Cerambycidae [15]. The highest A + T content is found in trnGlu of both mitogenomes (Tables 2 and 3). Most of the tRNAs have typical cloverleaf secondary structures, whereas the dihydrouridine (DHU) arm of trnSer1 is shortened in both mitogenomes (Figure 4), which is a common phenomenon in hexapods and metazoan mitogenomes [48]. Numerous mismatched base pairs are found in the secondary structures of tRNA genes, and all of them are G-U pairs.
The large ribosomal RNA (rrnL) gene and small ribosomal RNA (rrnS) gene are found in the conserved location between trnLeu1 and the control region (Tables 2 and 3). The rrnL gene is 1272 bp long in A. glabripennis and 1266 bp long in D. pseudonotabilis, with an A + T content of 80.1% and 78.8%, respectively. The rrnS gene is 779 bp long in A. glabripennis and 774 bp long in D. pseudonotabilis, with an A + T content of 78.6% and 76.5%, respectively.
The control region (CR) is the longest non-coding area in the two mitogenomes ( Figure 1) and is functional in the regulation, transcription, and replication processes of the mitogenomes [49]. The CR of A. glabripennis is 1104 bp long and has an A + T content of 79.3%; the CR of D. pseudonotabilis is 907 bp long and has an A + T content of 82.1% (Tables 2 and 3). In the CR of A. glabripennis, 5.2 copies of 57 bp long tandem repeat "AAAATTTCATCAGCTAGCTCCGCTATATAAAATCGCCTACCTTTCAAATTTCC-CCTA" are detected near the 5 end of this region. A total of 22 standard (single stem with single loop) and another 7 more complicated stem-loop structures are predicted in the CR of A. glabripennis ( Figure S1). There are 17 standard and 4 complicated stem-loop structures in the CR of D. pseudonotabilis ( Figure S2). However, no tandem repeats are found in the CR of D. pseudonotabilis. Functions of these secondary structures are unclear.

Transfer RNAs, Ribosomal RNAs, and Control Region
The two mitogenomes both contain the complete set of 22 tRNA genes typical of metazoan mitogenomes. These tRNAs range in size from 62 to 69 bp, which was consistent with previously sequenced mitogenomes of Cerambycidae [15]. The highest A + T content is found in trnGlu of both mitogenomes (Tables 2 and 3). Most of the tRNAs have typical cloverleaf secondary structures, whereas the dihydrouridine (DHU) arm of trnSer1 is shortened in both mitogenomes (Figure 4), which is a common phenomenon in hexapods and metazoan mitogenomes [48]. Numerous mismatched base pairs are found in the secondary structures of tRNA genes, and all of them are G-U pairs.

Phylogenetic Analyses
The phylogenetic positions of A. glabripennis and D. pseudonotabilis are reconstructed based on the combined mitochondrial gene set of 13 PCGs. The ML and BI analyses generated similar tree topology ( Figures 5 and S3). The phylogenetic results are largely congruent with the recent comprehensive mitogenomic phylogenetic study of Nie et al. (2020) [15]. The monophyly of Cerambycidae s.s. is not well-supported in both ML and BI trees due to the inclusion of other families of Chrysomeloidea ( Figure 5), which is similar to the results of Haddad et al. (2018) [5] and Nie et al. (2020) [15]. The positions of Disteniidae and Oxypeltidae are variable, and Oxypeltidae is recovered as the sister group to all other taxa in the BI tree ( Figure S3). The phylogenetic position of Disteniidae remains uncertain, and this family has been recovered as the sister group to various other members of Cerambycidae s.l. based on either molecular or morphological datasets [5,9,15,[50][51][52][53][54][55][56][57]. The monophyly of Cerambycidae s.s. is still one of the most debatable subjects in the phylogeny and evolution of Chrysomeloidea [5].
The monophyly of Cerambycinae, Dorcasominae, Lamiinae, and Necydalinae is wellsupported in both ML and BI analyses (Figures 5 and S3). The subfamily Parandrinae is placed within Prioninae and should be treated as a tribe of Prioninae, as suggested in previous studies [15,58]. Similarly, Necydalinae is nested in Lepturinae and should be regarded as a tribe of Lepturinae [15,50]. Spondylidinae is rendered paraphyletic by the species of Vesperidae, which differs from the monophyletic condition in Haddad et al. (2018) [5] and Nie et al. (2020) [15]. The non-monophyletic condition of Vesperidae has also been recovered based on morphological and molecular characters [54,55,[59][60][61]. The two species sequenced in this study are, respectively, grouped with their relatives from the same subfamily.
Although numerous contributions have been made to explore the higher-level phylogeny of longhorn beetles, there are still some debatable points to be solved: the monophyly of Cerambycidae s.l. and s.s.; the relative relationship between Cerambycidae s.s., Disteniidae, Oxypeltidae, and Vesperidae; and the monophyly and relationship of subfamilies in Cerambycidae s.l., especially within Cerambycidae s.s. The incongruence between different molecular phylogenetic studies could be attributed to the usage of different molecular types, sample sizes, and analytical methods. The taxonomic misidentification of sequenced samples in online databases such as GenBank could also lead to bizarre tree topology, especially for those clades with few taxa. Most main clades of Cerambycidae s.l. still lack sufficient molecular data to clarify their phylogenetic positions. The sequencing of more mitogenomes, optimization of datasets and substitution models, and the supplement of nuclear genes are expected to improve the resolution of mitochondrial phylogenetic reconstruction of Cerambycidae s.l. in future works.  The large ribosomal RNA (rrnL) gene and small ribosomal RNA (rrnS) gene are found in the conserved location between trnLeu1 and the control region (Tables 2 and 3). The rrnL gene is 1272 bp long in A. glabripennis and 1266 bp long in D. pseudonotabilis, with an A + T content of 80.1% and 78.8%, respectively. The rrnS gene is 779 bp long in A. glabripennis and 774 bp long in D. pseudonotabilis, with an A + T content of 78.6% and 76.5%, respectively.
The control region (CR) is the longest non-coding area in the two mitogenomes (Figure 1) and is functional in the regulation, transcription, and replication processes of the mitogenomes [49]. The CR of A. glabripennis is 1104 bp long and has an A + T content of 79.3%; the CR of D. pseudonotabilis is 907 bp long and has an A + T content of 82.1% (Tables  2 and 3). In the CR of A. glabripennis, 5.2 copies of 57 bp long tandem repeat "AAAATTTCATCAGCTAGCTCCGCTATATAAAATCGCCTAC-CTTTCAAATTTCCCCTA" are detected near the 5′ end of this region. A total of 22 standard (single stem with single loop) and another 7 more complicated stem-loop structures are predicted in the CR of A. glabripennis ( Figure S1). There are 17 standard and 4 complicated stem-loop structures in the CR of D. pseudonotabilis ( Figure S2). However, no tandem repeats are found in the CR of D. pseudonotabilis. Functions of these secondary structures are unclear.

Conclusions
In this study, we sequenced and analyzed the mitogenomes of two longicorn beetles, which are important pests of cultivated ecosystems in China. The structure and content of the two mitogenomes are conserved in comparison to other sequenced mitogenomes of

Conclusions
In this study, we sequenced and analyzed the mitogenomes of two longicorn beetles, which are important pests of cultivated ecosystems in China. The structure and content of the two mitogenomes are conserved in comparison to other sequenced mitogenomes of Cerambycidae, but the intraspecific mitogenomic variation is also detected. The monophyly of four subfamilies was supported by the phylogenetic analysis based on the nucleotide sequence of PCGs. The results provided basic genetic information for understanding the phylogeny and evolution of longicorn beetles.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/genes13101881/s1, Figure S1: Predicted stem-loop structures in the control region of A. glabripennis; Figure S2: Predicted stem-loop structures in the control region of D. pseudonotabilis; Figure S3: Bayesian inference phylogeny of Cerambycidae s.l. inferred from mitogenomic data. Numbers at the nodes are posterior probabilities.  Institutional Review Board Statement: No special permits were required to retrieve and process the samples because the study did not involve any live vertebrates or regulated invertebrates.

Informed Consent Statement: Not applicable.
Data Availability Statement: The data presented in this study are available in NCBI GenBank (Accession numbers: OP096420 and OP096419).