The complete mitochondrial genomes of four lagriine species (Coleoptera, Tenebrionidae) and phylogenetic relationships within Tenebrionidae

Sampling and DNA extraction

Cerogira janthinipennis and Crypsis chinensis specimens were collected in the Dayaoshan Mountains of the Guangxi Zhuang Autonomous Region, China. Luprops yunnanus and Spinolyprops cribricollis were collected using the sieve soil method in Jinghong City, Yunnan Province, China. A specimen of A. unidentasus was collected using the same method in Shiqian County, Guizhou Province, China. Gonocephalum kochi specimens were collected in Zunyi City, Guizhou Province, China, and Morphostenophanes yunnanus were collected in Yunlong County, Yunnan Province, China. The leg and thorax tissues were immediately preserved in 95% ethanol in the field and then were stored at −24 °C. The specimens were identified based on morphological characteristics. Total DNA was extracted from the muscle tissues using the Ezup Column Animal Genomic DNA Purification Kit (Shanghai, China) in accordance with the manufacturer’s instructions. This study was conducted under the approved guidelines of Animal Care and Use Committee of China West Normal University (approval number CWNU2022D021).

Sequencing, sequence assembly, annotation, and analysis

All of the mitogenomes in this study were sequenced with a whole genome shotgun strategy and Illumina sequencing technology using 150 bp paired-end reads. The adapter sequence and low-quality bases were deleted using Trimmomatic v. 0.36 (Bolger, Logse & Usadel, 2014). The target reads were assembled using IDBA UD v. 1.1.1 and Celera Assembler v. 8.3 (Crampton-Platt et al., 2015). The de novo assembly of the newly-generated sequences was conducted using Geneious 11.0.2 (Kearse et al., 2012). Protein-coding genes (PCGs) were identified as open reading frames. A total of 22 transfer RNA genes (tRNAs) were identified with the use of the MITOS Webserver, setting the default parameters with the Invertebrate Mito genetic code (Bernt et al., 2013). Their secondary structures were plotted manually from the MITOS predictions using Adobe Illustrator. Every sequence of tRNA genes was manually checked. The ribosomal RNA genes (rRNAs) and control region were identified by the boundaries of the tRNA genes. Mitogenome maps were produced using Organellar Genome DRAW COGDRAW (Greiner, Lehwark & Bock, 2019). The base composition and relative synonymous codon usage values of four Lagriinae mitogenomes were calculated using MEGA v 6.0 (Kumar, Stecher & Tamura, 2016). Strand asymmetry was calculated through the formula provided by Perna & Kocher (1995): AT-skew = (A – T)/(A + T), GC-skew = (G – C)/(G + C). The non-synonymous substitutions (Ka), synonymous substitutions (Ks), and Ka/Ks of PCG genes were calculated using DnaSP v. 5 (Rozas et al., 2003). Tandem repeat units in the control region of four Lagriinae species were identified using the Tandem Repeats Finder online tool (Benson, 1999).

Phylogenetic analyses

Phylogenetic trees for C. janthinipennis, L. yunnanus, A. unidentasus, S. cribricollis, and other tenebrionid beetles were reconstructed by aligning sequences of mitogenomes with those of 39 tenebrionid species (Table S1). Lymexylon navale and Hyloecetus dermestoides, belonging to Lymexyloidea, were chosen as outgroups, as they are phylogenetically closed to Tenebrionoidea in the Coleoptera (Timmermans et al., 2016; Cai et al., 2022). Phylogenetic trees were constructed using PhyloSuite v. 1.2.2 (Zhang et al., 2020) based on the maximum likelihood (ML) and Bayesian inference (BI) methods, respectively. The default parameters (ML: bootstrap: ultrafast; number of bootstrap: 5,000; maximum iterations: 1,000; minimum correlation coefficient: 0.90; replicates: 1,000; BI: generations: 2,000,000; sampling frequency: 100; number of runs: two; number of chains: four; burn-in fraction: 0.25) were used to construct phylogenetic trees. The default parameters (ML: bootstrap: ultrafast; number of bootstrap: 5,000; maximum iterations: 1,000; minimum correlation coefficient: 0.90; replicates: 1,000; BI: generations: 2,000,000; sampling frequency: 100; number of runs: two; number of chains: four; burn in fraction: 0.25) were used to construct phylogenetic trees. The resulting phylogenetic trees were edited and visualized using FigTree v. 1.43 (Horton et al., 2015).

Genome organization and base composition

We sequenced and annotated the complete mitogenomes of four lagriine species Cerogira janthinipennis (ON303727, length: 15,328 bp), Luprops yunnanus (ON303728, length: 16,437 bp), Anaedus unidentasus (ON303730, length: 15,387 bp), and Spinolyprops cribricollis (ON303729, length: 15,454 bp), as well as three other tenebrionid species: Gonocephalum kochi (ON303729, length: 15,825 bp), Crypsis chinensis (ON303729, length: 16,724 bp), and Morphostenophanes yunnanus (MW822745, length: 15,616 bp). The complete mitogenomes of these seven species contained 13 protein-coding genes (PCGs), two ribosomal RNA genes (rRNAs), 22 transfer RNA genes (tRNAs), and a control region (CR). Four of the 13 PCGs (nad1, nad4L, nad4, nad5), two rRNAs, and eight of the 22 tRNAs (trnY, trnC, trnQ, trnV, trnL1, trnP, trnH, trnF) are encoded on the N-strand, while the other 23 genes (nine PCGs and 14 tRNAs) are encoded on the J-strand (Table S1; Figs. S1, S2). The mitogenome sequences of these four lagriine species were determined to be medium- to maximum-sized compared to those of the Tenebrionidae species as a whole, which ranged from 15,328 bp (C. janthinipennis) to 16,437 bp (L. yunnanus).

The AT nucleotide contents of the four lagriine mitogenomes were similar (79.81% in C. janthinipennis, 72.12% in L. yunnanus, 74.90% in A. unidentasus, and 75.43% in S. cribricollis). The entire mitogenome had high A + T contents (ranging from 72.12–79.81%; 69.86–78.54% for PCGs, 77.34–81.16% for tRNAs, 76.79–82.56% for rRNAs, and 77.62–87.33% for the CR) (Table 1). The AT-skews in the four mitogenome sequences were −0.048, 0.127, 0.091, and 0.044, respectively, of which the C. janthinipennis had a negative value. All of the GC skews were negative (Table 1), which indicated that the content of G is lower than C.

Protein-coding regions and codon usage

In these four lagriine complete mitogenomes, the lengths of PCGs ranged from 11,030 to 11,087 bp, accounting for 67.36–72.22% of the full mitogenomes. In 13 PCGs, the nad5 (1,699–1,707 bp) and atp8 (147–192 bp) were the largest and smallest genes, respectively. All the PCGs had a typical ATN start codon, except nad1, which began with TTG. For stop codons, most PCGs terminated with TAR (TAA/TAG) and few had an incomplete stop codon T- (the TAA stop codon is completed by the addition of 3′ A residues to the mRNA) (Table S1).

In PCGs, the numbers of amino acids and relative synonymous codon usage (RSCU) were calculated. The 13 PCGs of A. unidentasus, C. janthinipennis, L. yunnanus, and S. cribricollis contained 3,684, 3,682, 3,681, and 3,667 codons, respectively, excluding stop codons. In these four species, the most frequently used codons were UUU (276, 379, 224, 280), UUA (299, 427, 274, 348), AUU (334, 401, 266, 350), and AUA (223, 243, 200, 231) (Fig. 1). Accordingly, F, L2, I, and N are the most frequently used amino acids in these four species.

The ratio of nonsynonymous/synonymous (Ka/Ks) of 13 PCGs was also calculated; the result is representative of mutations and the evolutionary rate (Hurst, 2002). The Ka/Ks of nad1 (3.157) are distinctly higher than 1 (Fig. 2), indicating that it is under strong positive selection (Mori & Matsunami, 2018), whereas cox1 (0.088), cox2 (0.114), and cytb (0.129) had values lower than the other genes (Fig. 2). These ratios suggest that the cox1, cox2, and cytb genes can be used as barcodes for deducing the phylogenetic relationships of Tenebrionidae.

In PCGs, the nucleotide diversity (Pi) with different data ranged from 0.211 (cox1) to 0.978 (atp8) (Fig. 3). Among the 13 PCGs, the gene atp8 (Pi = 0.978) is the most diverse nucleotide, followed by nad2 (Pi = 0.407) and nad6 (Pi = 0.375), however, cox1 is the most conserved gene with the lowest value (Pi = 0.211).

Transfer and ribosomal RNA genes

The total tRNAs lengths of A. unidentasus, C. janthinipennis, L. yunnanus, and S. cribricollis were 1,403, 1,391, 1,421, and 1,395 bp, respectively. Eight of the 22 tRNAs (trnY, trnC, trnQ, trnV, trnL1, trnP, trnH, trnP) were encoded on the N-strand, while the other 14 tRNAs were on the J-strand (Table 1). The tRNAs of four species were between 52 and 70 bp in length. The trnK is the longest of four species, while trnS1, trnA, trnS1, and trnS1 are the shortest of A. unidentasus, C. janthinipennis, L. yunnanus, and S. cribricollis, respectively. Among the 22 tRNAs, 21 tRNAs showed the typical clover-leaf secondary structure (Figs. S2–S5), however, trnS1 had a simple loop (Fig. 4). In L. yunnanus, a 526 bp insertion was present at trnD and atp8 junction.

The total length of rRNAs ranged from 1,961 bp (C. janthinipennis) to 2,011 bp (A. unidentasus). The AT content of four species ranged from 76.79% (L. yunnanus) to 82.56% (C. janthinipennis) (Table 2). The two rRNAs (rrnL and rrnS) of four species are encoded on the N-strand. The rRNAs of four species were located between the control region and trnL1, and separated by trnV.

Control region

The control region (CR) was the longest non-coding region. It was located between trnI and rrnS. The length of CR ranged from 831 bp (A. unidentasus) to 1,269 bp (L. yunnanus). The AT content of the four species ranged from 77.62% (L. yunnanus) to 87.33% (C. janthinipennis), which is distinctly higher than that of the full mitogenome, PCGS, rRNAs, and tRNAs in the same species. These four sequences had a positive AT skew. The tandem repeat element numbers of A. unidentasus, C. janthinipennis, L. yunnanus, and S. cribricollis are 1, 5, 3, and 3, respectively.

Phylogenetic analyses

The phylogenetic relationships among 39 Tenebrionidae species (Table S1) were reconstructed based on the nucleotide sequences of mitogenomes using the maximum likelihood and Bayesian inference methods. The outgroup species, Mylabris sp. (Meloidae) and Oedemera virescens (Oedemeridae), were the first to be separated from the tenebrionid clade. Phylogenetic trees based on the ML and BI methods showed approximately identical topologies (Figs. 5, S10). All the Tenebrionidae species were clustered together, which suggested that Tenebrionidae is monophyletic. The target species (A. unidentasus, C. janthinipennis, L. yunnanus, and S. cribricollis) and other lagriine species were clustered into single branch with high support (Fig. 5, BI = 0.87) suggesting that Lagriinae is monophyletic. These results are consistent with previous studies (Doyen & Tschinkel, 1982; Doyen, Matthews & Lawrence, 1989; Kergoat et al., 2014; Gunter et al., 2014; Timmermans et al., 2016; Aalbu, Kanda & Smith, 2017; Wu et al., 2022). The three other target species, Crypsis chinensis, Gonocephalum kochi, and Morphostenophanes yunnanus, were clustered together with Diaperinae, Blaptinae, and Stenochiinae species, respectively. In this study, the monophyly of the subfamilies Pimelinae (BI = 1), Blaptinae (BI = 0.89), Stenochiinae (BI = 1), Alleculinae (BI = 1) was also supported (Fig. 5); the subfamily Diaperinae was paraphyletic, and Tenebrioninae appears polyphyletic. The subfamily Pimelinae was the sister group to the remaining subfamilies, which was consistent with the results of Kergoat et al. (2014) but inconsistent with the results of other studies (Doyen & Tschinkel, 1982; Gunter et al., 2014; Timmermans et al., 2016; Wu et al., 2022). The results also suggested that the classification of Tenebrioninae and Diaperinae needs to be further studied.