Mitochondrial genomes of blister beetles (Coleoptera, Meloidae) and two large intergenic spacers in Hycleus genera

Abstract

Background

Insect mitochondrial genomes (mitogenomes) exhibit high diversity in some lineages. The gene rearrangement and large intergenic spacer (IGS) have been reported in several Coleopteran species, although very little is known about mitogenomes of Meloidae.

Results

We determined complete or nearly complete mitogenomes of seven meloid species. The circular genomes encode 13 protein-coding genes (PCGs), 22 transfer RNAs (tRNAs) and two ribosomal RNAs (rRNAs), and contain a control region, with gene arrangement identical to the ancestral type for insects. The evolutionary rates of all PCGs indicate that their evolution is based on purifying selection. The comparison of tRNA secondary structures indicates diverse substitution patterns in Meloidae. Remarkably, all mitogenomes of the three studied Hycleus species contain two large intergenic spacers (IGSs). IGS1 is located between trnW and trnC, including a 9 bp consensus motif. IGS2 is located between trnS2 (UCN) and nad1, containing discontinuous repeats of a pentanucleotide motif and two 18-bp repeat units in both ends. To date, IGS2 is found only in genera Hycleus across all published Coleopteran mitogenomes. The duplication/random loss model and slipped-strand mispairing are proposed as evolutionary mechanisms for the two IGSs (IGS1, IGS2). The phylogenetic analyses using MrBayes, RAxML, and PhyloBayes methods based on nucleotide and amino acid datasets of 13 PCGs from all published mitogenomes of Tenebrionoids, consistently recover the monophylies of Meloidae and Tenebrionidae. Within Meloidae, the genus Lytta clusters with Epicauta rather than with Mylabris. Although data collected thus far could not resolve the phylogenetic relationships within Meloidae, this study will assist in future mapping of the Meloidae phylogeny.

Conclusions

This study presents mitogenomes of seven meloid beetles. New mitogenomes retain the genomic architecture of the Coleopteran ancestor, but contain two IGSs in the three studied Hycleus species. Comparative analyses of two IGSs suggest that their evolutionary mechanisms are duplication/random loss model and slipped-strand mispairing.

Background

The animal mitochondrial genome (mitogenome) is an informative model for phylogenetics, molecular evolution, and comparative genomic research, due to its simple genetic structure, maternal inheritance, and high evolutionary rate properties [1,2,3]. The mitogenome has been increasingly used to analyze the phylogeny and evolution of the highly diverse and rapidly radiating insects [4,5,6,7,8]. The typical mitogenome of metazoans is a circular molecule encoding a conserved set of 37 genes for 13 protein-coding genes (PCGs), two ribosomal RNA (rRNA) genes, and 22 transfer RNA (tRNA) genes, and comprising a non-coding control region [3]. The mitogenome commonly displays exceptional economy of organization evidenced by lacking introns, few intergenic spacers, incomplete stop codons, and even overlapping genes [9]. However, the extremely diverse Insecta also exhibits high diversity in their mitochondrial genomes, such as gene rearrangements and/or long non-coding regions (except the control region) in some lineages within Hymenoptera [10,11,12,13,14], Hemiptera [15,16,17], Dictyoptera [18, 19], Diptera [20], Orthoptera [21], Thysanoptera [22], Psocoptera [23], and Phthiraptera [24]. Within the order Coleoptera [25], most beetles retain the same gene content and gene organization as the hypothesized ancestral mitogenome for Insecta [26], while a few possess gene rearrangements of tRNAs [7, 27, 28] and PCGs [29] in their mitogenomes. For example, the mitogenome of the firefly, Pyrocoelia rufa contains a 1724-bp large intergenic spacer (IGS) composed of tandem repeat units [5]. Non-coding intergenic spacers were also found in the mitogenomes of other beetles (Trachypachus holmbergi, Priasilpha obscura, and Cyphon sp.) [30]. Although a recent study reported that Hycleus chodschenticus has a large intergenic spacer [31], very little is known about the rearrangement and large non-coding regions in Meloidae mitogenomes.

Meloidae is a medium sized family within Coleoptera Tenebrionoidea, containing more than 3000 species within approximately 125 genera [32]. Meloids are commonly referred to as blister beetles due to a defensive secretion, cantharidin. Cantharidin is an intoxicant that can be used for the removal of warts and may be effective in the treatment of primary liver cancer, leucocytopenia, chronic liver disease, neurodermatitis and other major illnesses [33, 34]. The medicinal properties of cantharidin, hypermetamorphosis development, and parasitoid habit of Meloidae species have been extensively researched [35, 36]. However, the mitogenomes of Meloidae species is less well researched. Two of the few studies of Meloidae mitogenomes found that Epicauta chinensis and Lytta caraganae retained the ancestral model of the insect mitogenome, without any gene rearrangement or long non-coding regions [31, 37]. Of the approximately 3000 species within Meloidae, only three complete mitogenomes (of Epicauta chinensis, Lytta caraganae, and Hycleus chodschenticus) have been described (excluding another without any description; Mylabris sp., JX412732.1). The lack of research considerably limits the genomic comparisons and molecular phylogenetic studies of Meloidae. Thus we believed there was an urgent need to explore the mitogenome evolution in the diverse family of Meloidae.

Consequently, we determined the mitogenomes of seven blister beetle species, representing four meloid genera. These species were E. gorhami, E. tibialis, L. caraganae, Mylabris aulica, H. phaleratus, H.marcipoli, and H. cichorii. We described the general features of the newly sequenced mitogenomes from the seven species and analyzed two IGSs in all Hycleus species to explore their evolutionary mechanisms. In addition, we attempted to assess the possibility of the IGS2 to be a molecular marker and indicators of phylogenetic relationships within Meloidea, based on mitogenomic datasets. The mitogenomes of the seven meloids will significantly add to the knowledge of Meloidae taxonomy, phylogeny, and evolution.

Results

Genome content and gene organization

This study presents six complete mitogenomes and one nearly complete mitogenome (H. cichorii) with the absence of the control region and three tRNAs (trnI, trnQ and trnM). The total lengths of complete mitogenomes range from 15,633 to 16,003 bp. All seven new sequences were submitted to GenBank under the accession numbers listed in Table 1. Every mitogenome of the seven meloid species is a circular DNA molecule, encoding the typical 37 genes including 13 PCGs, 22 tRNAs, two rRNAs, and a putative control region. The major strand (J strand) carries most of the genes (9 PCGs and 14 tRNAs), while the remaining genes are encoded on the minor strand (N strand) (Fig. 1).

Table 1 Locality information and accession numbers of meloid species employed in this study

Fig. 1

Circular map of the meloid mitochondrial genomes. tRNA genes are abbreviated to the capital letter of their coding amino acid, L, L2, S and S2 indicate trnL1 (CUN), trnL2 (UUR), trnS1 (AGN), and trnS2 (UCN). Arrows indicate the direction of transcription

All PCGs of the seven mitogenomes use typical ATN start codons. Conventional stop codons TAA/TAG are assigned to most of the PCGs, but the cox1, cox2, nad5 and nad4 genes of all meloids terminate with the incomplete stop codon T. This terminator is adopted by cox3 genes of Epicauta and Lytta species (Additional files 1, 2, 3, 4, 5, 6 and 7: Tables S1–S7). The A + T - rich regions of meloid mitogenomes range from 1015 to 1201, with the location between rrnS and trnI (Table 2). The poly-T stretch (15 bp) was detected in control regions of all meloids, but without tandem repeats.

Table 2 Annotations of the seven meloid mitochondrial genomes

A + T content and codon usage

The overall A + T contents of the seven meloid mitogenomes range from 67.53% to 71.94% (Table 1), and such an A + T bias is reflected in the codon frequencies of these mitogenomes (Fig. 2). Relatively synonymous codon usages (RSCU) were calculated over all seven meloid species, excluding stop codons (Additional file 8: Table S8). The RSCU demonstrate that codons with A or T in the third position are always overused as compared to other synonymous codons. Additionally, codons TTT (Phe), TTA (Leu), ATT (Ile), and ATA (Met) are the four most frequently used codons in these meloid mitogenomes. These codons are all comprised of A or T nucleotides, which is indicative of the biased usage of A and T nucleotides in the meloid mitochondrial PCGs.

Fig. 2

Relative synonymous codon usage (RSCU) in mitochondrial genomes of seven meloids. Average 3701 codons for meloids were analyzed, excluding stop codons. The RSCUs are averages over all seven meloid taxa. Codon families are provided on the x-axis. Leu, Leu*, Ser, and Ser* indicate trnL1 (CUN), trnL2 (UUR), trnS1 (AGN), and trnS2 (UCN), respectively

Evolutionary rates of PCGs

All available mitogenomes were used to assess the evolutionary rate of PCGs for Meloidae. The variable sites, nucleotide diversity (π), and the ratio of non-synonymous substitution (Ka) to synonymous substitution (Ks) were calculated for each PCG (Table 3). Nad6 and atp8 genes exhibit the highest level of nucleotide diversities, whereas cox1 gene is the most conserved. The Ka/Ks value is correlated with the nucleotide diversity, with the highest level in nad6 and atp8 and the lowest in cox1. Notably, the Ka/Ks ratio for every PCG is lower than 1, indicating that all PCGs are evolving under the purifying selection.

Table 3 Evolutionary rates of mitochondrial PCGs among meloid species

Comparison of tRNA secondary structures

All seven meloid mitogenomes encode 14 tRNAs on the J strand and 8 on the N strand (Fig. 1, Table 2). The comparative results of tRNA secondary structures are provided in Figs. 3 and 4. All tRNAs could be folded into the typical clover-leaf structure, except trnS1 (AGN) lacks a dihydrouridine (DHU) arm, which is replaced by a simple loop (Fig. 4).

Fig. 3

Secondary structure of tRNAs (trnA-trnL1) in meloid mitogenomes. The nucleotide substitution pattern for each tRNA was modeled using as reference the structure determined for Hycleus phaleratus

Fig. 4

Secondary structure of tRNAs (trnL2-trnV) in meloid mitogenomes. The nucleotide substitution pattern for each tRNA was modeled using as reference the structure determined for Hycleus phaleratus

The nucleotide conservation of tRNAs is markedly J strand-biased, trnA, trnE, trnG and trnK, which have the highest percentage of identical nucleotides, are all located on the J strand. The tRNAs on the N strand (i.e., trnQ, trnH, trnF) and close to the control regions (i.e., trnI and trnQ) have high level of nucleotide variation (Additional file 9: Table S9).

Intergenic spacers

The mitogenomes of E. gorhami, E. tibialis, and L. caraganae contain four intergenic spacers with a total length of 23, 25, 26 bp, respectively, while the M. aulica has six intergenic spacers with a total length of 28 bp (Table 2). Unexpectedly, H. phaleratus, H. marcipoli, and H. cichorii contain 7, 8, 8 intergenic spacers with a total length of 266, 203, 175 bp, respectively (Table 2), including two large intergenic spacers (>30 bp). IGS1 is located between trnW and trnC with a length of 62, 56, 39 bp in H. phaleratus, H. marcipoli, and H. cichorii, respectively (Table 2). The IGS1 sequences exhibit relatively high similarity among the three Hycleus species, and a 9 bp long congruent motif AAATTATGG was detected in the three Hycleus species (Fig. 5a). IGS2 is located between trnS2 and nad1 with a length of 181, 123, and 110 bp in H. phaleratus, H. marcipoli, and H. cichorii, respectively (Table 2). The alignment of IGS2s among all sequenced Hycleus mitogenomes, including the recently published sequence of H. chodschenticus [31], shows that a pentanucleotide motif (TACTA) exists in these Hycleus species. Furthermore, H. chodschenticus, H. phaleratus, H. marcipoli, and H. cichorii include five, four, three, and two repeats (respectively) of this motif (Fig. 5b). The organization of IGS2 indicates that these four Hycleus species contain two copies of an 18 bp conserved sequence (ATACTAAAYTTTRTTAAC) in both ends of IGS2 (Fig. 5b), but other meloid beetles have only one (Fig. 6a).

Fig. 5

Alignments of the intergenic spacers in Hycleus mitochondrial genomes. a IGS1 located between trnW and trnC. The nucleotides in red refer to consistent ones in three Hycleus species. The shaded nucleotides indicate the congruent motif (AAATTATGG) found in Hycleus insects. - refers to gap (b) IGS2 located between trnS2 and nad1. The shaded nucleotides indicate the conserved pentanucleotide motif (TACTA). The nucleotides in red indicate the 18 bp conservative sequence ATACTAAAYTTTRTTAAC. - refer to gaps, dots refer to omitted nucleotides

Fig. 6

a Organizations of the intergenic spacers between trnS2 (UCN) and nad1 in mitogenomes of all meloids. b Putative mispairing involves replicating DNA strands of Hycleus mitogenomes. c The proposed mechanism of IGS2 in Hycleus mitogenomes. The CS indicates the 18 bp conservative sequence ATACTAAAYTTTRTTAAC

Phylogenetic relationships

We carried out MrBayes, RAxML, and PhyloBayes analyses based on nucleotide and amino acid datasets to determine the influence of different datasets and analytical methods on tree topology and node reliability. Bayesian Inference (BI) and Maximum Likelihood (ML) analyses used the same datasets to generate congruent tree topologies. BI trees had higher node support values than ML trees (Fig. 7). PhyloBayes analyses generated different tree topologies with polytomies (Additional file 10: Figure S1). Four tree topologies derived from our six phylogenetic trees are consistent with the monophylies of Meloidae and Tenebrionidae, and the basal position of Mordellidae. Differences are present in the inter-family relationships of Aderidae, Ciidae, Oedemeridae, and Prostomidae.

Fig. 7

The phylogenetic tree of 16 species from superfamily Tenebrionoidea based on the nucleotide dataset (a) and the amino acids dataset (b) of 13 mitochondrial protein-coding genes, inferred from Bayesian inference and maximum likelihood. The numbers abutting branches refer to Bayesian posterior probabilities (left) and ML bootstraps (right), − not recovered. Branch lengths and topology are from BI analyses. The Diabrotica barberi (Coleoptera: Chrysomelidae) was employed to root the trees as outgroup

Within Meloidae, L. caraganae is sister taxon to the species belonging to the genus Epicauta. The Meloidae family results monophyletic and receives maximum statistical support. In Tenebrionidae, the Tenebrioninae and Diaperinae are never recovered as monophyletic groups (Fig. 7, Additional file 10: Figure S1).

Discussion

General features

All seven mitochondrial genomic arrangements share the ancestral type for insects [26, 38], as is reported in the published meloid mitogenomes [31, 37]. The nucleotide composition of mitochondrial genomes for meloids also corresponds well to the A + T bias generally observed in insect mitogenomes.

All mitogenomes of the seven meloids have incomplete stop codons, which have been described in many other insect species [39, 40]. It has been demonstrated that incomplete stop codons can produce functional stop codons in polycistronic transcription cleavage and polyadenylation processes [9]. Remarkably, the cox3 genes of Lytta and Epicauta species possess the same incomplete stop codon, while Hycleus and Mylabris beetles utilize complete terminators (Additional files 1, 2, 3, 4, 5, 6 and 7: Tables S1–S7). The similar preference for the adoption of stop codons seems to suggest that the genus Lytta is more closely related to Epicauta than the two other genera, and this relationship was confirmed by phylogenetic results.

The evolutionary rates of all mitochondrial PCGs indicate that their evolution is based on purifying selection (Table 3), as is reported in other insects [41, 42]. However, the cytochrome oxidase subunits (cox1, cox2, and cox3) and cytochrome b (cob) have lower Ka/Ks ratios than ATPase subunits (atp8 and atp6) and NADH dehydrogenase subunits (nad1–6 and 4 L). The nucleotide diversity also shows cox and cob genes are obviously more conserved than atp and nad genes. This phenomenon indicates that the various functional genes in the mitochondria of meloids underwent different selection pressures during evolution. Furthermore, the cox1 has the slowest evolutionary rate, demonstrating that functional constraints are more powerful for this gene than positive selection.

The absence of a DHU arm of trnS1 commonly exists in many metazoan mitogenomes, including insects [30, 43, 44]. However, this tRNA (missing DHU arm) was often suggested by authors to be functional [45, 46]. Another unusual feature is the use of TCT as the trnS1 anticodon in meloids, although most arthropods adopt a GCT anticodon in trnS1. This exceptional trnS1 anticodon was also found in many other beetles [5, 30, 47].

Some mismatched pairs in stems of tRNAs (e.g., T-T in the DHU stem of trnQ and in anticodon stem of trnK; C-T in the TΨC stem of trnI; A-C in anticodon stem of trnL2 (UUR); A-G in acceptor stem of trnW), are common in insect mitogenomes and can be corrected through editing processes or may represent unusual pairings [44]. It was not possible to model the substitution pattern of the TΨC loop in trnH (Fig. 3) because of the high variation among orthologous sequences. The increasing variation usually accompanies more compensatory base changes in stems, resulting in the tRNA more or less not conserved (Additional file 9: Table S9).

The ends of rRNA genes of meloid mitogenomes were assumed to extend to the boundaries of flanking genes because it is impossible to accurately determine by DNA sequence alone [48]. Consequently, rrnL was assumed to fill up the blank between trnV and trnL1 (CUN) (Fig. 1), but the boundary between the rrnS and the putative control region was defined by the alignment with homologous sequences [49].

The control region in the insect mitogenome is equivalent to the control region of vertebrate mitogenomes, which contains the origin sites for transcription and replication [50, 51]. The six complete mitogenomes include a poly-T stretch (15 bp) that was suggested to function as a possible recognition site for the initiation of replication of the mitochondrial DNA N strand [50]. Like other Coleopteran mitogenomes, the control regions of meloids also exhibit the highest A + T content in the whole mitogenome. This region is unlikely to be more variable than protein-coding genes due to such high A + T content and consequently limits its usefulness as a molecular marker [52].

Intergenic spacers

All newly sequenced mitochondrial genomic arrangements share the ancestral type for insects without rearrangement, but possess large non-coding regions (except the control region) in some lineages. The intergenic spacers in the mitogenomes of E. tibialis, E. gorhami and L. caraganae are similar to those in E. chinensis mitogenome [37]. The total length of M. aulica’s intergenic spacers is not significantly different from the former four meloids, but its genome does contain more intergenic spacers (Table 2). Unexpectedly, the whole lengths of intergenic spacers in the three Hycleus mitogenomes are much longer than those of other meloids. The most remarkable feature is the presence of two IGSs in the mitogenomes of three Hycleus species. A 494-bp long intergenic spacer was also reported in the recently published mitogenome of H. chodschenticus [31]. Consequently, the total lengths of known Hycleus mitogenomes are longer than those of other meloid mitogenomes, but predominantly due to the presence of IGSs rather than the lengths of genes or control regions.

The mitochondrial genome typically displays exceptional economy of organization, evidenced by lack of introns, few intergenic spacers, incomplete stop codons and even overlapping genes [9]. However, the large IGSs in mitochondrial genomes were observed in some Hymenopteran [11, 13], Hemipteran [15, 17] Dictyopteran [18] and Coleopteran [5, 30] insects. Previously reported IGSs contain tandem repeat units (in Pyrocoelia rufa and Evania appendigaster) [5, 11], or additional origin of replication [10] and similar to it [15]. The IGS in Reduviidae bugs have unidentified open reading frames encoding 103 or 104 amino acids but without blast similarity [15, 17]. In contrast, the two IGSs in Hycleus species have no significant similarity with other genes of their mitogenomes and lack open reading frames or tandem repeats. Nevertheless, the IGS2 of all studied Hycleus species contain discontinuous repeats of a 5 bp consensus motif (TACTA) (Fig. 5b). This motif was also found in many other Coleopteran insects [30], similar to the 7 bp motif ‘ATACTAA’ conserved in Lepidoptera [6] and the hexanucleotide motif ‘THACWW’ in Hymenoptera [11]. The pentanucleotide motif was suggested as the possible binding site of a transcription termination peptide (mtTERM), as its position signifies the end of the J strand coding region in the circular mitochondrial DNA [38]. However, we do not know the function of this discontinuous repeat.

The discontinuous repeats and the 18 bp long consensus sequence in both ends of IGS2 within all studied Hycleus species (Fig. 6a) suggest that the slipped-strand mispairing [53] may be the evolutionary mechanism of this IGS. According to this theory, mispairing involves dissociation of replicating DNA strands and then misaligned reassociation (Fig. 6b), following replication or repair lead to insertions of several repeat units. Formed tandem repeat experiences random loss and/or point mutation, only the repeat units in both ends are completely retained and the residues form the IGS2 (Fig. 6c). Although we are not absolutely certain about our assumption due to the highly divergent region, the fragmented repeat units in highly divergent region (Fig. 5b) and complete repeat units at both ends (Fig. 6a) suggest that the slipped-strand mispairing is the most convincing mechanism for IGS2.

The IGS between trnW and trnC was only previously reported in Trachypachus holmbergi [30], while no IGS (>30 bp) at this position has yet been found in the mitogenomes for Meloidae. Although H. chodschenticus has 3-bp intergenic spacer at the same position [31], it is too usual as many short intergenic spacers to be considered as IGS1. Therefore, the evolutionary mechanism of IGS1 in the three Hycleus species in present study is different from H. chodschenticus. The 9-bp consistent motif and the relatively high similarity among the three Hycleus species in the present study (Fig. 5) suggest that they have the mutual mechanism of IGS1. The duplication/random loss model [1, 54] may account for the IGS1. We speculate that tandem duplication of trnW-trnC-trnY is caused by some error in DNA replication, followed by random loss of partial duplicated genes, and then the residues form the IGS1 (Fig. 8). This model was also proposed as the mechanism for rearrangements in some Hemipteran [17] and Mantodea mitogenomes [19], and the IGSs in Blaptica dubia [18]. It is possible that the duplication/random loss event of Hycleus beetles occurred relatively early and many nucleotides were deleted during the random loss period. Consequently, the residual IGS1 has low similarity with the original tRNAs, and the IGS1 is somewhat not conserved in Hycleus.

Fig. 8

Putative mechanism of IGS1 in mitogenomes of Hycleus species under the duplication/random loss model. The random losses of partial genes are marked with *

Similar to the three Hycleus species in the present study, H. chodschenticus also possesses the IGS between trnS2 and nad1 [31]. Although intergenic spacers with approximately 20 bp at this position are common in Coleopteran mitogenomes, no large intergenic spacers have been found within other Coleopteran lineages to date. Based on the current knowledge of IGS2 across Coleoptera, it is reasonable to assume that the emergence of IGS2 in Hycleus mitogenomes might be a special evolutionary mutation. It is difficult to identify the closely related genera Hycleus and Mylabris since they share very similar morphological characteristics. However, IGS2 might be a potential marker to distinguish Hycleus from its closely related and indistinguishable genera with respect to the sizeable intergenic spacer that exists in all studied Hycleus species but is absent in other genera,. Although we are unable to adequately confirm that IGS2 exists in all Hycleus species due to limited samples, it provides a new candidate for molecular identification of this genus. Variations in the quantity, location and sequence of intergenic spacers might also be a valuable marker for phylogenetic and evolutionary studies at lower taxonomic levels, if these details of more taxa were obtained in future.

Phylogeny

Phylogenetic studies indicated that different datasets and inference methods influence the tree topology of Coleoptera [31, 55]. Our phylogenetic results showed that tree topologies are sensitive to datasets rather than inference methods, since the different inference methods with the same datasets generated consistent tree topologies. The heterogeneous-site model in PhyloBayes was suggested as being more reliable for phylogenetic inferences within Coleoptera [31, 55]. Our Bayesian analyses under the heterogeneous-site model are unable to resolve phylogenetic relationships within Tenebrionoidea, but perhaps this is due to insufficient taxa. Although the nucleotide dataset of PCGs was proposed as better for phylogenetic analyses at superfamily level [31], we could not assess the quality of different datasets because some key lineages share the same tree topologies. However, nodal support values are sensitive to inference methods. For the same datasets, ML trees had significantly lower support values than BI trees, especially at several nodes that involved the family-level relationships of Aderidae, Ciidae, Scraptiidae (Fig. 7). This is consistent with previous phylogenetic studies using mitogenomes [14, 31, 37].

Phylogenetic relationships within Tenebrionoidea are ambiguous. The inter-family relationships are also uncertain, especially for Aderidae, Ciidae, Oedemerdae, and Prostomedae, which are respectively represented by only one taxon. However, all tree topologies well recover the monophyly of Tenebrionidae, Meloidae, and Mordellidae (Fig. 7, Additional file 10: Figure S1). Tenebrioninae and Diaperinae are never recovered as monophyletic groups, as is reported by Gunter et al. [56]. To date, no phylogeny has successfully resolved the interfamilial relationships between tenebrionoids, either by using morphological or molecular characteristics. The comprehensive phylogeny of Coleoptera [57, 58] and the largest molecular phylogeny of Tenebrionoidea [56] were unable to recover strong support or definitively infer the phylogenetic relationships within the superfamilies. In contrast to these phylogenetic studies based on several genes, our phylogenetic inferences may be hindered by lack of across taxon sampling rather than dataset validity.

The genus Lytta is the sister group to Epicauta rather than grouped with Mylabris within Meloidae (Fig. 7). The placement of Lytta and whether Lytta is more closely related to Mylabris or Epicauta could not be inferred from previous molecular phylogeny of the family Meloidae based on 16S rRNA and ITS2 [32]. We confirmed this relationship by multiple inference approaches based on both nucleotide and amino acid sequences of 13 mitochondrial PCGs. Considering the high diversity and rapid radiations of insects [59, 60], mitochondrial genomes could be a better approach to resolve intractable phylogenetic relationships due to its relatively rapid rate of mutation and purely maternal inheritance [3, 61]. Consequently, we believe our conclusion that the genus Lytta is more closely related to Epicauta than Mylabris or Hycleus is reliable because it is based on complete mitogenomes and the preference of stop codon usage.

The data collected thus far regarding meloid mitogenomes could not resolve the phylogenetic relationships within Meloidae. In fact, no phylogeny of Meloidae based on either morphological or molecular characteristics has been able to successfully resolve the relationships at genera and species levels. Taxon sampling is known to be one of the most significant determinants of accurate phylogenetic inferences, particularly in species rich lineages [62, 63]. Considering the diversity of the family Meloidae and the limitation of the present molecular information, more conclusive phylogenetic results will be achievable as bio-information becomes increasingly available. This study will assist with these more conclusive phylogenetic results and future studies on taxonomy, phylogeny and systematics of Meloidae insects.

Conclusions

Our study presents the mitochondrial genomes of seven meloid beetles. All complete mitogenomes of meloids retain the typical gene content and organization of the ancestor. The evolutionary rates of all PCGs in the studied Meloidae indicate that their evolution is according to purifying selection. The comparison of tRNA secondary structures exhibit diverse substitution patterns in Meloidae. Two large intergenic spacers exist in the three studied Hycleus mitogenomes, and the sequence and structure of the two IGSs contributed to our conclusion regarding their possible evolutionary mechanisms. The phylogenetic results inferred from mitochondrial genomes support that the genus Lytta is more closely related to Epicauta than to Mylabris. Although data collected thus far could not resolve the phylogenetic relationships within Meloidae, this study will assist in future mapping of the Meloidae phylogeny.

Methods

Samples collection and DNA extraction

The specimens of the seven meloid species used for this study are listed in Table 1, with locality data. The fresh materials were immediately preserved in 100% ethanol and stored in a − 20 °C refrigerator. Total genomic DNA was extracted from a frozen adult using Tianamp Genomic DNA kit.

PCR amplification and sequencing

All mitochondrial genomes of these collected meloid species were generated by amplifications of nine overlapping PCR fragments. Eight fragments were amplified using common primers for all seven species designed from the aligned E. chinensis mitogenome (GenBank accession number KP692789) [37], only the primers of the control region were specifically designed for each species. Details of primers are given in Additional file 11: Table S10. The PCR was performed with Vazyme Taq DNA Polymerase (Mg²⁺ Plus Buffer) and carried out on a PTC-100 thermal cycler (BioRad, Hercules, CA). PCR conditions used were: 5 min denaturation at 94 °C; 35 cycles of 30s at 94 °C, 30s at 49–56 °C and 1–3 min (1 min/Kbp) at 72 °C; followed by 10 min extension hold at 72 °C. The PCR products were sequenced on an ABI PRISM 3730 DNA sequencer by Tsingke Biotechnology Company with primers walking on both strands.

Genome assembly and annotation

Sequences from overlapping fragments were assembled with the neighboring fragments using SeqMan program included in the Lasergene software package (DNAStar Inc., Madison, Wisc.). Protein-coding regions were identified by ORF Finder from the NCBI website with invertebrate mitochondrial genetic codes, and compared with published mitochondrial sequences by using MEGA 6.0 [64]. Most of the tRNA genes were identified using tRNAscan-SE 1.21 [65] with invertebrate genetic codon predictors; however, the trnS1 was predicted by alignments with other homologous genes because of its lack of dihydrouridine (DHU) arm. The rRNA gene boundaries were interpreted as the end of a bounding tRNA gene and alignment of sequences with homologous regions of known Coleopteran mitogenomes. The control regions were assumed between rrnS and trnI within all meloid mitochondrial genomes. The A + T contents, relative synonymous codon usage values, and evolutionary rates (number of variable sites, nucleotide diversity, and Ka/Ks ratios) for each PCG were calculated via MEGA 6.0 [64].

Phylogenetic analysis

Phylogenetic analyses were assessed using 29 Tenebrionoidea species representing 8 families, with the Chrysomelid Diabrotica barberi (GenBank accession number NC_022935.1) [8] employed as the outgroup. Species’ PCGs were extracted according to GenBank annotations by using GenScalpel [66]. All these nucleotide and amino acid sequences were aligned using MUSCLE [67] with the default setting. Gaps and ambiguous sites were removed from the protein alignment to generate a 10,356-bp nucleotide dataset and a corresponding amino acid dataset (3452 amino acids). The best partitioning schemes and corresponding evolutionary models were selected by PartitionFinder 1.1.1 [68] with 13 partitions defined by genes. We set the model selection as Bayesian information criterion (BIC), unlinked branch lengths, and greedy search algorithm to estimate the best fitting schemes. The best-fit partitioning schemes and corresponding models are shown in Additional file 12: Table S11.

Phylogenetic analyses with site-homogeneous model were conducted to reconstruct the phylogenetic trees of superfamily Tenebrionoidea by using Bayesian inference (BI) and maximum likelihood (ML) methods. Bayesian phylogenetic analysis was implemented using MrBayes 3.2.2 [69], and ran for 10,000,000 generations sampling per 1000 generations. Bayesian posterior probabilities were estimated using the Markov chain Monte Carlo (MCMC) sampling approach. ML analysis was carried out with RAxML-HPC2 on XSEDE 8.0.24 [70] using 1000 bootstraps to estimate the node support. Bayesian analyses with a site-heterogeneous model were performed using PhyloBayes MPI 1.5a [71] with two MCMC chains run after the removal of constant sites under the CAT-GTR model.