Comparative Mitgenome Analysis of Anoplophora horsfieldi and Other Chrysomeloidea, Cucujiformia Insects Reveals Conserved Mitogenome Organization and Phylogeny

Anoplophora is an important genus in the longhorned beetle family Cerambycidae. They are widely distributed in the Oriental and eastern Palearctic regions, and contain 36 species (McDougall, 2001; Richard and Gregory, 2002). Several notorious pest insects belong to this genus (McDougall, 2001). Anoplophora horsfieldi has body length of about 35 mm. Camellia sinensis, Celtis sinensis, Quercus glauca and Ulmus pumila are main host plants. Adults appeared in summer, and lived in low and middle mountainous areas. This species can be found in India, Vietnam and China (Chen et al., 1959). Most of the Cerambycidae species, as a large family in the Coleoptera order, are forestry pests with above 25,000 species worldwide (Sama et al., 2010; Wang et al., 2012). Although Cerambycidae shows large taxonomic diversity, it is still limited knowledge about the Cerambycidae mitogenome.

Mitogenomes are important molecular markers, which can be used for exploring the phylogenetics, population genetics and phylogeographics of species due to low recombination rate, maternal inheritance, conserved gene organization and orientation and high mutation rate (Dowton et al., 2002; Jin et al., 2004; Wang and Tang, 2018). In general, the metazoan mitogenome forms a closed-circular structure, with a length of 15 to 18 kilobases (kb) (Dowton et al., 2002; Lu et al., 2013; Cameron, 2014), including 13 PCGs, 22 tRNAs, 2 rRNAs and a non-coding region (A+T-rich region) that in other Cucujiformia and in invertebrates is named as A+T-rich region because of its extremely high content in Adenines and Thymines (Wolstenholme, 1992; Cameron, 2014; Sun et al., 2017; Zhang et al., 2019). Recent high throughput sequencing development has enriched insect mitogenome datasets and accelerated the research on insect molecular evolution, too. At present, insect mitochondrial genome data has been used for species identification, comparative genomics, evolutionary genetics, reconsecration of phylogeny and phylogeography (Cameron and Whiting, 2008; Hao et al., 2012; Ma et al., 2012; Timmermans et al., 2014; Li et al., 2016).

Coleoptera is an ancient order and can adapt to a variety of habitats (Li et al., 2016). Coleoptera has approximately 400,000 described species and contains Adephaga, Archostemata, Myxophaga and Polyphaga (Crowson, 1960; Nie and Yang, 2013; Zhang, 2013). Based on morphological features, Crowson (1960) suggested Polyphaga can be divided into five infraorders and 18 superfamilies. However, Bouchard et al. (2011) think that Polyphaga can be divided into six infraorders and 16 superfamilies. Because of different phylogenetic positioning of some superfamilies, so these classification systems contain different numbers of infraorders (Nie and Yang, 2013). Cucujiformia, as the most highly diversified infraorders of Polyphaga (Crowson, 1960; ØDegaard, 2000; Bouchard et al., 2009; Wang and Tang, 2018) includes six superfamilies and most of two superfamilies species (Curculionoidea and Chrysomeloidea) are plant-feeding beetles (Grimaldi and Engel, 2005). Recently, Robertson et al. (2015) and Yuan et al. (2016) suggest that Cucujiformia can be divided into seven superfamilies (Cleroidea, Coccinelloidea, Lymexyloidea, Tenebrionoidea, Cucujoidea, Chrysomeloidea, Curculionoidea). Moreover, long-horned beetles and leaf beetles in Chrysomeloidea and weevils in Curculionoidea are considered as Coleoptera ‘Phytophaga’. Chrysomeloidea contains the families Chrysomelidae, Cerambycidae, Megalopodidae, Vesperidae, Oxypeltidae, Disteniidae and Orsodacnidae (Hunt et al., 2007).

The phylogeny and evolution of Chrysomeloidea have received considerable attention of systematists. Nonetheless, the interrelationships (families and subfamilies) are unclear and many chrysomeloid evolutionary questions still persist. This is especially true for Cerambycidae, for which there is a few molecular phylogenies, however, the numbers and species used in the construction of phylogenetic trees are limited. Several studies have explored the phylogenetic relationships of Chrysomeloidea by using the combined data including morphological data (Farrell and Sequeira, 2004; Gómez-Zurita et al., 2007) and molecular data (18S rRNA Hunt et al., 2007; 28S rDNA Marvaldi et al., 2009, and partial mitochondrial genes Bocak et al., 2014; Li et al., 2016; Wang and Tang, 2018). In addition, there are very few molecular studies about Chrysomeloidea taxa. Currently, some researchers have found some new evidence to understand the Chrysomeloidea phylogenetic relationships, but some conflicting results are often achieved due to the differences of datasets and analytical methods. Therefore, the phylogenetic relationships of Cucujiformia and Chrysomeloidea need to be further evaluated.

In this study, we achieved the complete sequence of the mitogenome of A. horsfieldi using next-generation sequencing and compared it with other cucujiform insects. We aimed to analyze the structural characteristics of A. horsfieldi mitogenome and compared it with previously sequenced Anoplophora mitogenomes. Moreover, based on all 117 mitogenome sequences also used BI and ML methods to reconstruct a tentative phylogeny to evaluate relationships among Cucujiformia and Chrysomeloidea. Our review is structured taxonomies, commencing with the phylogenetic neighborhood of Chrysomeloidea and their relatives, followed by sections on Chrysomelidae, Cerambycidae and then Cucujiformia and its relatives. The sequence characteristics and annotation of A. horsfieldi mitogenome will be a significant increase in further research of Cucujiformia and Chrysomeloidea mitogenome structures and phylogenetics.

The specimen of the A. horsfieldi, was sampled from an adult vulture captured in the Yongxing Town, Mianyang City, Sichuan Province, China in June 2020 (104°39′5.43″E, 31°27′39.06″N, 475 m.a.s.l). The specimen was preserved in 95% ethanol and stored -76°C until genomic DNA extraction. Total genomic DNA was extracted from the thorax muscle of a single specimen using an E.Z.N.A.® Tissue DNA Kit (Omega, Norcross, GA) abiding by the manufacturer’s instructions. The extracted DNA quality was examined by 0.9% agarose gel electrophoresis (w/v) and used to sequence the whole mitogenome of A. horsfieldi with next-generation sequencing.

Genome sequencing libraries with about 400 bp of insertion fragment were constructed with a NEXTflex™ Rapid DNA-Seq Kit (Illumina, San Diego, CA) according to the manufacture’s protocols. The library was sequenced on Illumina Hiseq X Ten platform to produce 150 bp paired end reads (300 cycles).

SOAPdenovo v2.0 (Luo et al., 2012), MITObim v1.8 (Christoph et al., 2013) and NOVOPlasty v2.7.1 (Dierckxsens et al., 2017) were alternate application to assembly the mitochondrial genome. The assembled mitochondrial fragments were identified by BlastX using A. chinensis (Li et al., 2015) (NC_029230) and A. glabripennis (Fang et al., 2016) (NC_008221) mitochondrial genes as queries. Prediction and annotation of 13 PCGs, 22 tRNA and 2 rRNA genes were performed with DOGMA (http://dogma.ccbb.utexas.edu/) or MITOS (http://mitos.bioinf.uni-leipzig.de/index.py) using the support of annotation from reference mitogenome.

The mitogenome sequence is further verified and verified below. 12S and 16S rRNA genes were determined by comparison using homologous sequences of mitochondrial DNA from other Anoplophora species with ClustalX version 2.0 (Larkin et al., 2007) 13 PCG sequences were translated into putative proteins according to the Invertebrate Mitochondrial Genetic Code. The base composition of nucleotide sequences was calculated by skewness on the basis of the following formulas: AT skew = [A−T]/[A+T] and GC skew = [G−C]/[G+C] (Perna and Kocher, 1995). The content of A+T and relative synonymous codon usage (RSCU) values were described with MEGA 7.0 (Kumar et al., 2016). The tRNA genes were resolved with the tRNAscan-SE 1.21 (http://lowelab.ucsc.edu/tRNAscan-SE/) (Lowe and Chan, 2016) and the MITOS Web Server (Perna and Kocher, 1995; Bernt et al., 2013). The tRNAs secondary structure were achieved with RNAviz v2.0 (De Rijk et al., 2003). The tandem repeats of the A+T-rich region were explored by the tandem repeat finder (http://tandem.bu.edu/trf/trf.html) (Benson, 1999; Timmermans et al., 2014).

The newly sequenced mitogenomes of A. horsfieldi was aligned with 116 mitogenomes of Cucujiformia and Chrysomeloidea available in GenBank, with Prosopocoilus gracilis (Coleoptera: Scarabaeiformia: Lucanidae: KP735805), Cheirotonus jansoni (Coleoptera: Scarabaeiformia: Scarabaeidae: NC_023246) and Necrophila americana (Coleoptera: Staphyliniformia: Silphidae: NC_018352) as outgroups (Supplementary Table S1). The multiple alignments (13 PCGs concatenated nucleotide sequence datasets) were performed with ClustalX soft (v2.0 Larkin et al., 2007) with the default settings. Accession numbers of all mitogenomes are listed in the phylogenetic trees (Supplementary Table S1). We employed the nucleotide sequences of the 13 PCGs as the dataset to construct the maximum likelihood (ML) and Bayesian inference (BI) phylogenetic trees. The nucleotide sequences of 13 PCGs from all 117 Cucujiformia mitogenome were aligned separately in MEGA 7.0 (Kumar, et al., 2016). The final alignments were refined and the conserved sequences were identified with Gblocks v0.91b (Talavera and Castresana, 2007). Geneious 8.1.6 was used to concatenate the resulting alignments (Kearse et al., 2012). The optimal nucleotide substitution model was selected on the basis of the Akaike Information Criterion (Posada and Buckley, 2004) with J Model test v.0.1.1 (Posada, 2008). ML and BI phylogenetic trees show different algorithms. ML phylogenetic tree was constructed by PhyML v3.0 (Guindon et al., 2010) based on the best-fit model (as above) with 1,000 replicates. BI tree was conducted using MrBayes 3.2 with a GTR + I + G model, each of four chains (three hot and one cold), with run length of 20 million generations and sampling every 1,000 generations (Ronquist et al., 2012). Convergence was assessed with Tracer 1.5 (Rambaut and Drummond, 2007) and trees from the first 25% of the samples were removed as burn-in. Based the value of Bayesian posterior probabilities (BPP), node support was assessed. Figtree v1.4.3 software was used to view and edit the consensus trees (Rambaut, 2009).

We deposited the complete mitogenome of A. horsfieldi in the GenBank database (MN248534) and the mitogenome sequence of A. horsfieldi is 15,796 bp in size (Table I and Fig. 1), which is located in the range of whole sequenced Cucujiformia species with the length ranging from 15,064 bp in Naupactus xanthographus (Curculionidae) to 20,124 bp in Ceutorhynchus obstrictus (Curculionidae) (Supplementary Table S1). Alignment results revealed 38 mitogenome regions, including 13 PCGs, two rRNA, 22 tRNA and a non-coding region with high A+T-rich composition, which is common in most animal mtDNAs (Table I). The majority strand (J-strand) of the mitogenome of A. horsfieldi encodes nine PCGs and fourteen tRNAs, while the remaining genes are encoded on the minority strand (N-strand) (Table I).

The gene arrangement and orientation of A. horsfieldi mitogenome are similar to the mitogenome of typical Cerambycidae species. The differences in size of mitogenomes among Cucujiformia insects can be explained by the numbers and types of repetitive sequences in the A+T-rich regions.

Three species of the genus Anoplophora possess the same characteristics in the mitochondrial genome. The total base composition of A. horsfieldi was A (39.74%), T (39.57%), C (12.47%) and G (8.22%); A (39.48%), T (38.17%), C (13.57%) and G (8.78%) in A. chinensis and A (39.62%), T (38.71%), C (13.07%) and G (8.59%) in A. glabripennis, respectively. A heavy AT bias was found in the three species of mitochondrial genomes (79.31% in A. horsfieldi, 77.65% in A. chinensis and 78.34% in A. glabripennis, Supplementary Table S2-S4), which is similar to other sequenced Cerambycidae species (Wang et al., 2013; Fang et al., 2016; Li et al., 2016). All three Anoplophora species present a positive AT skew and n et al egative GC skew in the whole mitogenome and had a higher A + T content in rRNAs than in tRNAs (Table II). In fact, the region rich in A + T has been viewed as the main source of the variations in the length of the whole mitogenome, and in the process of polynucleotide operation, insertion/deletion, the number of nucleotides is variable, and the copy numbers and types of tandem repeat elements are also very different between different species (Zhang and Hewitt, 1997).

The length of the total PCGs is 11,133 bp for A. horsfieldi, 11,130 bp for A. chinensis and 11,092 bp for A. glabripennis, accounting for 70.48%, 70.42% and 70.25% of their total mitogenomes, respectively (Table II). Four PCGs (ND1, ND4, ND4L and ND5) are from N strand, and the remaining genes are found on J strand (Fig. 1 and Table I). Thirteen PCGs showed variable range from 156 bp (ATP8) to 1711/1714 bp (ND5) in A. horsfieldi, A. chinensis and A. glabripennis (Table I). All three mitogenomes showed similar characteristics with the smallest size of ATP8 and the largest that of ND5 among PCGs (Fang et al., 2016; Li et al., 2016). In addition, their AT skewness and GC skewness are negative and positive, respectively which are not consistent with those of complete mitogenome (Table II).

Almost all PCGs of three species start with the ATN codon, only ND5 in A. chinensis starts with the GTG codon (Table II). In A. horsfieldi, an incomplete stop codon “T–” is used for COXI, COXII, COXIII, ND3, ND4 and ND5, while ND1 used TAG as stop codon, and the remaining five PCGs used TAA as stop codon (Table I). In A. chinensis, an incomplete stop codon “T” is used for COXI, COXII, ND3, ND4 and ND5, while ND1 and Cytb used TAG as a complete stop codon, and the remaining six PCGs used TAA as complete stop codon (Table I). Moreover, in A. glabripennis an incomplete stop codon “T–” is used for COXI, COXII, COXII, ND1, ND3, ND4 and ND5, while TAG is used as a complete stop codon for Cytb, and TAA is used as a complete stop codon for the remaining five PCGs (Table I). In all three mitochondrial genomes, the termination codon TAA appeared higher than TAG and at least five incomplete stop codons “T--” were present (Table I) (Fang et al., 2016; Li et al., 2016). The incomplete termination codons could be added as TAA by post-transcriptional polyadenylation during the mRNA process maturation (Ojala et al., 1981; Schuster and Stern, 2009; Guindon et al., 2010; Ronquist et al., 2012). And this phenomenon has also been discovered in other insects as well (James and Andrew, 2006; Larkin et al., 2007; Wang et al., 2013; Wu et al., 2014; Huang et al., 2015; Liu et al., 2018).

The similar amino acid base composition and the relative synonymous codon usage (RSCU) of the three Anoplophora species are found (Fig. 2 and Table III).The genome-wide bias towards AT was viewed to reflect in the codon usage by the PCGs. The total number of codons of the PCGs ranges from 3695 to 3709 (Table II). UUA (Leu2), AUU (Ile), UUU (Phe) and AUA (Met) are the most frequently utilized amino acids (Fig. 3), all high frequency codons are composed of A or U. The third codon frequency of A/T is significantly higher than that of G/C, reflecting nucleotide A + T bias in the mitochondrial PCGs among Cucujiformia. The composition of most commonly used amino acids, like Leu, Ile, and Phe, varied from 42.47% (A. glabripennis) to 43.12% (A. horsfieldi) (Fig. 3). This pattern with rich with A or T nucleotides in all PCGs was also similar to other Cucujiformia insects (Kim et al., 2009; Du et al., 2016; Wang and Tang, 2018).

Table II. Nucleotide composition and skew values in Anoplophora horsfieldi, Anoplophora chinensis and Anoplophora glabripennis.

Twenty-two tRNAs of A. horsfieldi, A. chinensis and A. glabripennis mitogenomes scattered discontinuously across the complete mitogenome (Table I). The tRNAs length of these three species were 1,437, 1,446 and 1,431 bp, respectively, and accounted for 9.10%, 9.15% and 9.08% of the total mitogenomes, respectively (Table I). The three mitogenomes possess 22 tRNA genes, eight transcribed from N strand and 14 from J strand (Fig. 1, Table I). The length of these 22 tRNAs range from 60 (trnR) to 76 bp (trnD) in A. horsfieldi, from 62 (trnC) to 84 bp (trnD) in A. chinensis and from 62 (trnH) to 74 bp (trnP) in A. glabripennis (Table I). 21 of 22 tRNAs can form typical cloverleaf secondary structure, while trnS1 (AGN) formed a simple loop due to lacking the DHU arm (Fig. 4, Supplementary Figs. S1-S2), as discovered in other coleopteran species (Stewart and Beckenbach, 2003; Song et al., 2010; Cabrera-Brandt and Gaitan-Espitia, 2015). The lack of DHU stem in trnS1 is a common phenomenon in insect mitogenomes (Beckenbach, 2011; Cameron, 2014; Li et al., 2015), and has been proven as a typical characteristic of metazoan mitogenomes (Wang et al., 2019). Many nucleotide substitutions are spotted in four different stems and the anticodon loop is highly conserved (Fig. 4). Besides the orthodox AU and CG pairs, some mismatched base pairs are also discovered in different stems. A total of 15 GU mismatches, 3 UU mismatches, and 1 AG mismatch are found in A. horsfieldi (Fig. 4), 17 GU mismatches, 2 UU mismatches, 2 AG mismatches, and 1 AA mismatch are found in A. glabripennis while 17 GU mismatches, 2 UU mismatches, 1 AG mismatch, and 1 AC mismatches are found in A. chinensis (Fig. 4, Supplementary Figs. S1-S2, Table I).

Two rRNA genes (rrnL and rrnS) are encoded by the N-strand in A. horsfieldi, A. chinensis and A. glabripennis and rrnL is distributed between trnL1 and trnV, while rrnS lies between trnV and CR. The genes rrnL and rrnS are 1273 and 804 bp long in A. horsfieldi and A. chinensis, and 1236 and 806 bp long in A. glabripennis, respectively (Table I). The lengths range from 1236 to 1273 bp in rrnL, and from 804 to 806 bp in rrnS in these three mitogenomes of genus Anoplophora (Fang et al., 2016; Li et al., 2016). The A+T content of the two rRNAs (80.88%) in A. horsfieldi was higher than 80.01% in A. chinensis and 79.79% in A. glabripennis, with positive AT skew (0.010) and negative GC skew (-0.385), however, in the other two species, the AT and GC skews are negative and positive, respectively (Table III).

A total of 15 gene overlap in the A. horsfieldi mitogenome from 1 to 13 bp, adding up to 59 bp, 14 gene overlaps exist in the A. chinensis mitogenome with the length from 1 to 8 bp, amounting to 45 bp, while 18 gene overlaps occur in the A. glabripennis mitogenome with the length from 1 to 9 bp, adding up to 64 bp (Table I). The longest overlap region (13 bp) of the three mitogenomes is found between trnK and trnD. All three Anoplophora species have eight identical overlap regions, including trnQ-trnM (1 bp), ATP6-COXIII (1 bp), trnE-trnF (1 bp), ND6-Cytb (1 bp), Cytb-trnS (2 bp), ND4-ND4L (4 bp), ATP8-ATP6 (7 bp), trnY-COXI (8 bp) and trnW-trnC (8 bp) (Coates, 2014; Fang et al., 2016; Li et al., 2016; Yao et al., 2017). Eight intergenic spacers appear in A. horsfieldi, ranging in length from 2 to 17 bp and amounting to 42 bp. Nine intergenic spacers occur in A. chinensis, ranging in size from 1 to 17 bp and adding up to 41 bp. A total of 115 bp in A. glabripennis distributed in 12 intergenic spacers, ranging in size from 1 to 35 bp. All three mitochondrial genomes share four identical intergenic spacers (trnR-trnN (2 bp), ND4L-trnT (2 bp), ND5-trnS (3 bp), trnS-ND1 (17 bp)) (Table I). The size of the intergenic spacers is more variable than overlaps.

Two 7-bp long overlaps (ATGATAA) were discovered in the Anoplophora species, which were also observed in other Polyphaga insects (Fig. 5), however, 7-bp long overlaps (ATGATTA) was discovered in Epicauta chinensis (KP692789) and 7-bp long overlaps (ATGATAG) was found in Tribolium castaneum (AJ312413) (Friedrich and Muqim, 2003) and Cryptolestes pusillus (NC_028204) (Li et al., 2016) between ATP8 and ATP6. The overlaps lie between ATP8 and ATP6 on the J-strand and between ND4L and ND4 on the N-strand, respectively. The overlapped sequences were translated as a bicstron (Stewart and Beckenbach, 2005; Wang and Tang, 2018). The other two 5 bp (TTAAT) and 7 bp (TTTAGT) long motif were detected between trnSer (UCN) and ND1 in mitogenomes of the three Anoplophora species, which was also present in other cucujiform beetle (Fig. 6). However, another 5 bp long motif (TAGTA) was found at this location (Wang and Tang, 2018). This has been explained as the possible

binding site of mtTERM because it is located at the end of the H-strand region in the circular mitochondrial genome (Taanman, 1999). Based on mitogenomic comparisons, the motifs between ATP8 and ATP6, and between trnSer and ND1 were relatively conserved in four suborders (Polyphaga, Adephaga, Archostemata and Myxophaga). However, the motif “ATGATAA” was only found between ND4 and ND4L of the mitogenomes of Polyphaga insects, while “ATGTTAA” was found in the other three suborders’ mitogenomes (Adephaga, Archostemata and Myxophaga) (Wang and Tang, 2018).

The A+T-rich region in insect mitochondrial genome, similar to the control region of vertebrate mitochondrial genome, possess the transcription and replication origin sites (Andrews et al., 1999; Yukuhiro et al., 2002; Zhu et al., 2013). The largest A+T-rich region lied between the genes rrnS and trnI. This region plays important roles in regulating transcription and replication (Zhang and Hewitt, 1997). The lengths of A+T-rich region are 1,143 bp in A. horsfieldi. In all three Anoplophora species, A. glabripennis has the smallest control region of 1,115 bp, while A. horsfieldi has the largest that of 1,143 bp (Fig. 7, Table I). The A + T contents are 86.44% in A. horsfieldi, 85.59% in A. chinensis and 87.98% in A. glabripennis. The motif of tandem repeats in A+T-rich region has been recorded in many insect orders (Zhang and Hewitt, 1997; Du et al., 2017; Wang et al., 2018; Huang et al., 2019; Li et al., 2019), but is found in A. horsfieldi for the first time. A. horsfieldi and A. glabripennis have three types of long repeat tandem units and one types of microsatellite-like repeat, while there are only two different long repeats in A. chinensis (Fig. 7B). The repeat units presented obvious differences in size and copy number among Coleoptera species, thus resulting in various lengths of the A+T-rich region.

Tandem repeats (due to difference in type and number) can be used as microsatellites for phylogeography and population genetics researches and also as a molecular marker in evolutionary biology (Wang et al., 2018; Zhang et al., 2019). In addition, in some insect species (Kim et al., 2006; Yin et al., 2012), there was no long tandem repetitive sequence in this region, but it did have some microsatellite-like repeats (such as, (A)9, (T)15, (C)8, (TA)6, (TA)8). The A+T-rich region in A. horsfieldi possess a 16 bp long poly-T and 14 bp long poly-A stretches which has been considered as a possible recognition site for the initiation of replication of the mtDNA minor strand (Fig. 7A, Andrews et al., 1999; Kim et al., 2009; Du et al., 2016). When compared with the A+T-rich region of three Anoplophora specie, three conserved sequence blocks (CSBs) were recognized in those species. These CBSs ranged in length from 23 to 111 bp, and their sequence similarity among species was generally more than 50%. These CSBs were also found in other stoneflies mitogenomes (Qian et al., 2014; Chen and Du, 2015, 2017; Cao et al., 2019), but the functions and characteristics of these conserved blocks need to be further studied. In A. horsfieldi and A. glabripennis, the AT and GC skewness of the A+T-rich region are negative and positive, respectively, while they are opposite in A. chinensis (Table II).

The phylogenetic analyses were carried out using the concatenated nucleotide sequences of 13 PCGs from 117 Cucujiformia mitogenomes (Supplementary Table S1). Similar topologies were achieved by using Bayesian (BI) and maximum likelihood (ML) analyses (Fig. 8). Phylogenetic analysis showed that A. horsfieldi grouped with A. glabripennis and A. chinensis with high nodal supports. In Anoplophora, A. horsfieldi is at the base of the genus. In our phylogenetic analyses, Anoplophora was clustered with Monochamus, rather than genus Psacothea and it was consistent with the results of previous phylogenetic analyses based on mitochondrial data (Coates, 2014; Yuan et al., 2016; Liu et al., 2018; Wang et al., 2019). In addition to Orsodacne lineola (Vesperidae), Prioninae, Cerambycinae, Disteniinae, Cassidinae and

Philinae species, Cerambycidae is closely related with family Chrysomelidae, and Anastrangalia sequensi was positioned as the basal lineage of Cerambycidae (Fig. 8). At the same time, the phylogenetic topological structure supports the two families as monophyletic. However, comparing our results to previous morphological and molecular analyses of Cerambycidae and Chrysomelidae relationships was harder to determine (Marvaldi et al., 2009; Yuan et al., 2016; Liu et al., 2018). Previous phylogenetic analyses of molecular data resulted in a monophyletic Cerambycidae and Chrysomelidae (Zhou et al., 2016; Liu et al., 2017; Yao et al., 2017). This may be because of the limited mitogenomic data, especially for Chrysomelidae species. Based on phylogenetic tree, the monophyly of Chrysomelidae and Cerambycidae cannot be supported in Chrysomeloidea, and their evolutionary relationships are relatively complex (Fig. 8) and need to be confirmed by additional samples and species and phylogenetic information, including nuclear genes.

Previously, the series Cucujiformia insects contain six superfamilies (Bouchard et al., 2011; Lawrence et al., 2011). However, our results support that Cucujiformia is divided into seven superfamilies, which is in agreement with the views of Robertson et al. (2015) and Yuan et al. (2016). Although the topological structures of phylogenetic trees are different among different data strings and analysis methods, these studies do not support the previous classification of Cucujoidea as monophyletic (Hunt et al., 2007; Marvaldi et al., 2009; Lawrence et al., 2011; Bocak et al., 2014; Robertson et al., 2015; Timmermans et al., 2016), which is consistent with our phylogenetic results (Fig. 8), and found that Elateroidea is embedded between Coccinelloidea and Cleroidea classification. Recently, Robertson et al. (2015) constructed a high-level phylogenetic analysis for Cucujoidea according to four mitochondrial and four nuclear genes, where two superfamilies (Cleroidea and Cucujoidea) were redefined and a new superfamily Coccinelloidea (cerylonid series) was proposed, and our reconstructed phylogenetic results also support the Coccinelloidea as a new superfamily (Fig. 8). In these seven superfamilies, except for Dastarcus helophoroides and Cyllorhynchites ursulus, five superfamilies (Lymexyloidea, Cucujoidea, Coccinelloidea, Cleroidea and Curculionoidea) are consistently supported as monophyletic with ML and BI methods, which is in accord with the results of Timmermans, et al. (2016) and Yuan et al. (2016) (Fig. 8). Although the monophyly of other two superfamilies (Chrysomeloidea and Coccinelloidea) are recovered in PhyloBayes analyses with the P123R dataset, our results support that Chrysomeloidea and Tenebrionoidea belong to paraphyletic relationships based on 13 PCGs (Fig. 8). Yuan et al. (2016) carried out the PhyloBayes analyses with the P123R dataset, and formed such a cucujiform relationship of (Cleroidea + (Coccinelloidea + ((Lymexyloidea + Tenebrionoidea) + (Cucujoidea + (Chrysomeloidea + Curculionoidea))))) (Fig. 8), however, our topology tree supports such phylogenetic relationships of ((Lymexyloidea + Tenebrionoidea) + (Curculionoidea + (Chrysomeloidea + (Cucujoidea + (Cleroidea + Coccinelloidea))))) (Fig. 8). This is different from the studies of (((Chrysomeloidea + Cucujoidea) + Tenebrionoidea) + Cleroidea) obtained by Kim et al. (2009), ((((Curculionoidea + Chrysomeloidea) + Cucujoidea) + (Cucujoidea + Cleroidea)) Tenebrionoidea) in the PCG-BI tree obtained by Li et al. (2016), (((Curculionoidea + Chrysomeloidea) + Cucujoidea) + (Cleroidea + (cerylonid lineages + (Lymexyloidea + Tenebrionoidea))) constructed by Timmermans et al. (2016) and (Cleroidea + (Erotylid series + (Tenebrionoidea + (Cucujoidea + (Chrysomeloidea + Curculionoidea))))) of Crampton-Platt et al. (2015) where these superfamily relationships within Cucujiformia cannot be well resolved. Lymexyloidea + Tenebrionoidea were positioned as the basal lineage of the series Cucujiformia insects in the PCG-BI/ML tree (Fig. 8), which was inconsistent with previous reports (Coates, 2014; Li et al., 2016), whereas superfamily Tenebrionoidea is located at the basal position in the ML and BI trees, and is also different from the basal position of Cleroidea in Cucujiformia (Kim et al., 2009; Marvaldi et al., 2009). These differences may be due to the use of different numbers of samples, molecular markers and analytical methods.

From the view of topology tree, the layout of Chrysomeloidea species is complicated and paraphyletic, which is consistent with Bocak et al. (2014). Relationships of some families in Chrysomeloidea remain to be poorly supported, especially for the small families Orsodacnidae and Vesperidae. Besides Prioninae, Cerambycinae, Disteniinae, Cassidinae and Philinae, the phylogenetic relationship of the family Cerambycidae is monophyletic. Unsatisfactory results were also found for internal relationships of Cerambycidae for which the sampling number and density of species were low and especially the representativeness of conservative loci was poor (Bocak et al., 2014). The Prioninae, Cerambycinae, Disteniinae, Cassidinae and Philinae species are not clustered in the family Chrysomelidae, which leads to the formation of paraphyly of the family. Besides the five subfamilies, Chrysomelidae were well supported and this is in accord with previous study results from Gomez-Zurita et al. (2007) and Bocak et al. (2014). The major clades were resolved (Fig. 8), containing the basal splits of a group of Donaciinae, Criocerinae and Bruchinae, the monophyly of Cassidinae. However, our results support the paraphyly of Galerucinae, Chrysomelinae and Cerambycinae (Fig. 8). The phylogenetic relationships at subfamily Prioninae, Cerambycinae, Disteniinae, Cassidinae and Philinae, especially in family Vesperidae and Orsodacnidae, still have not been resolved due to the lack of sufficient phylogenetic information and many more new species.

In word, as far as Cucujiformia, Chrysomeloidea and Chrysomelidae is concerned, the conflict results about phylogeny might be due to the different types of molecular markers and analytical methods (Kim et al., 2009; Marvaldi et al., 2009; Coates, 2014; Crampton-Platt et al., 2015; Meiklejohn et al., 2014). Therefore, the phylogenetic relationship and taxonomic status within the Cucujiformia insects require a large number of additional samples and more molecular markers to elucidate their evolutionary relationship. Undoubtedly, this study will contribute to explore the phylogeny, systematics and taxonomy of the Cucujiformia insects.

In this study, the entire mitogenome of A. horsfieldi (16,759 bp) was sequenced and annotated. The mitochondrial gene order and orientation of A. horsfieldi were similar to Chrysomeloidea species. The motifs, ‘ATGATAA’ between ND4L and ND4, was more conserved than that between trnS (UCN) and ND1 and between ATP8 and ATP6 in Cucujiformia. Phylogenetic tree showed that A. horsfieldi grouped with A. glabripennis and A. chinensis with high nodal supports. Within Cucujiformia, such phylogenetic relationships of (Lymexyloidea + (Tenebrionoidea + (Curculionoidea + (Chrysomeloidea + (Cucujoidea + (Cleroidea + Coccinelloidea)))))) were rebuilt, and the paraphyletic relationships that Chrysomeloidea, Coccinelloidea and Curculionoidea were supported. We believe that the mitogenomes of A. horsfieldi will helps people to further the studies of molecular evolution, and Cucujiformia mitogenome architecture and phylogenetics.

We would like to express our gratitude to all those who helped us during the writing of this manuscript.

This research was supported by the Natural Science Foundation Project of Science and Technology, Department of Sichuan Probince (No. 2023NSFSC0206), the National Natural Science Foundation of China (No. 32270457), the Research Project of Ecological Security and Protection Key Laboratory of Sichuan Province (No. ESP2003).

The animal handling procedures conformed with the China Animal Welfare guidelines and have been approved by the Ethic and Animal Welfare Committee and the Animal Protection and Use Commission of Mianyang Normal University (MYSY2017JC02).

This study was carried out in accordance with the Animal Care and Use Committee at the Mianyang Normal University. Efforts were taken to minimize suffering and included administering anesthesia. The study did not involve endangered or protected species.

There is supplementary material associated with this article. Access the material online at: https://dx.doi.org/10.17582/journal.pjz/20220824100800

The authors have declared no conflict of interest.