Adult specimens of Menida musiva (Jakovlev, 1876) were collected from Gaoleshan National Nature Reserve (32°39.90'N, 113°37.37'E), Tongbai County, Nanyang City, Henan Province, China, in August 2019. Adult specimens of M. lata Yang, 1934 were collected from Buddhist College of Tongbo County (32°21'N, 113°23'E), Nanyang City, Henan Province, China, in August 2019. Adult specimens of M. metallica Hsiao & Cheng, 1977 were collected from Wuli Village (30°52'N, 103°35'E), Qingchengshan Town, Dujiangyan City, Sichuan Province, China, in September 2020. All samples were immediately placed in absolute ethanol and stored in a freezer at −20 °C until DNA extraction. The total DNA was extracted from thoracic tissue using the HiPure Universal DNA Kit (Jisi Huiyuan biotechnology, Nanjing, China).

The complete mitochondrial genomes of the three species were sequenced on Illumina Novaseq 6000 Sequencing System with a read length of PE150. Fastp (Chen et al. 2018) software was used to filter the original data and remove the joint sequences and low-quality reads to obtain high-quality, clean data. Three mitochondrial genomes were assembled using SPAdes v. 3.10.1 (Bankevich et al. 2012), and the assembly of the genomes did not depend on the reference genome. After assembly, the complete mitogenomes were manually annotated using Geneious v. 11.0 (Kearse et al. 2012) software. A reference sequence of M. violacea for annotation was obtained from the basic local alignment search tool (BLAST) in the NCBI database. PCGs were identified by open reading frame (ORF) Finder (http://www.ncbi.nlm.nih.gov/gorf/gorf.html) implemented through the NCBI website using invertebrate mitochondrial genetic codes. The position and structure of 22 tRNAs were predicted using the MITOS Web Server (http://mitos.bioinf.uni-leipzig.de/index.py/) (Bernt et al. 2013). The exact locations of rRNA adjacent genes and the control regions were determined by confirming the boundary between them. In addition, tandem repeats of the control region were identified with the Tandem Repeats Finder server (http://tandem.bu.edu/trf/trf.html) (Benson 1999).

The circular maps of mitogenomes were produced by the CGView Server (Grant and Stothard 2008). Nucleotide composition and codon usage were analyzed with MEGA v. 11 (Tamura et al. 2021). To investigate the evolutionary patterns among the mitochondrial PCGs in Pentatominae species, DnaSP5 software (Librado and Rozas 2009) was used to count the non-synonymous substitutions (Ka) and synonymous substitutions (Ks) of 13 PCGs of Pentatominae, and to calculate the Ka/Ks values. The skew of the nucleotide composition was calculated with the formulas: AT-skew = (A − T) / (A + T) and GC-skew = (G − C) / (G + C) (Perna and Kocher 1995).

We selected three newly sequenced species of Menida and 37 available mitogenomes of related taxa (including all available Pentatominae sequences and two Acanthosomatidae sequences as outgroups) from GenBank to determine the phylogenetic status of Menida and to discuss the phylogenetic relationships of tribes within the subfamily Pentatominae (Table 1). The phylogenetic relationships were reconstructed based on two datasets: (1) 13 PCGs + 2 rRNAs (PR) and (2) 13 PCGs + 2 rRNAs + 22 tRNAs (PRT). The two data sets represent relatively complete genetic evolution information of mitochondrial genomes.

The nucleic acid sequences of the PCGs and RNA genes were extracted using Geneious v. 11.0 and aligned using the MUSCLE strategy in MEGA v. 11. Multiple sequences for each species were then connected using SequenceMatrix v. 1.7.8 (Vaidya et al. 2011), protein-coding genes were optimized using MACSE (Ranwez et al. 2011), ambiguous loci were deleted using Gblocks (Talavera and Castresana 2007), and converted into Nexus and Phylip formats in Mesquite v. 3.7 (Maddison 2008). To determine the best model for partitioning, four datasets were analyzed using PartitionFinder v. 2.1.1 (Lanfear et al. 2017). The maximum likelihood (ML) and Bayesian inference (BI) methods were used for phylogenetic analysis based on two datasets. The ML trees were constructed by IQ-TREE v. 2.2.0 (Minh et al. 2020), and the support value for each node was evaluated by the standard bootstrap (BS) algorithm, which was tested 500,000 times. The Bayesian inference (BI) method was used for phylogenetic analysis based on four datasets. The BI tree was constructed by MrBayes v. 3.2.7 (Ronquist et al. 2012). Two independent runs were run for 10 million generations, and samples were taken every 1000 generations. Four independent Markov Chains (including three heated chains and a cold chain) were run. A consensus tree was obtained from all the trees after the initial 25% of trees from each MCMC run were discarded as burn-in, with the chain convergence assumed after the average standard deviation of split frequencies fell below 0.01.

We studied the relationship among four species of Menida (three newly sequenced species and one species downloaded from NCBI). All four mitogenomes are double-strand circular DNA molecules. The total lengths of the mitogenomes of M. musiva, M. lata, M. metallica, and M. violacea are 16,663bp, 16,463bp, 16,418bp, and 15,379bp, respectively. The mitogenomes of the four species each contain 37 genes (13 protein-coding genes (PCGs), 22 tRNA genes, and 2 rRNA genes) and a control region (Fig. 1), with 23 genes located on the J-strand and 14 genes on the N-strand. The sequence of genes was consistent with the original gene arrangement of Drosophila yakuba Burla, 1954 (Clary and Wolstenholme 1985; Hua et al. 2008) without rearrangement. Nucleotide composition of the complete mitogenome of M. musiva: A 42.51%, T 33.70%, C 14.18%, G 9.60%; nucleotide composition of the complete mitogenome of M. lata: A 41.95%, T 32.92%, C 15.08%, G 10.05%; nucleotide composition of the complete mitogenome of M. metallica: A 41.39%, T 33.51%, C 13.77%, G 11.33%; nucleotide composition of the complete mitogenome of M. violacea: A: 42.19%, T: 33.32%, C: 13.86%, G: 11.33%. The four species show similar nucleotide composition (Suppl. material 1: table S1). All the mitogenomes exhibit a strong base composition bias toward AT, ranging from 74.86% to 76.22% in the four species (mean value = 75.37%). Moreover, all mitogenomes have a slightly positive AT-skew (ranging from 0.11 to 0.12, mean = 0.11) and a negative GC-skew (ranging from 0.20 to −0.10, mean = −0.16) (Suppl. material 1: table S1).

Figure 1. 

Gene arrangements of the four complete mitochondrial genomes.

The four mitogenomes have similar overlapping regions and gene spacers. The longest intergenic region (31–34bp) of the four species of the genus Menida appeared between trnS2 and nad1, and there were mainly three conserved overlaps, with a 8 bp overlap between trnC/trnW (AAGCTTTA) and a 7 bp overlap between atp8/atp6 and nad4/nad4L (ATGATAA) (Suppl. material 1: table S2).

For the four studied species, nine PCGs (nad2, cox1, cox2, atp8, atp6, cox3, nad3, nad6, and cytb) were found to be coded on the majority strand (J-strand) and four PCGs (nad5, nad4, nad4L, and nad1), on the minority strand (N-strand). The longest PCG is nad5 (1707–1710 bp), while the shortest is atp8 (159 bp). Five PCGs (cox1, cox2, atp8, atp6, and nad3) did not vary in length among the four species. Most of the PCGs use ATN (ATT/ATA/ATG/ATC) as initiation codon. TTG was the second most used initiation codon, and was found in cox1, atp8, nad1, and nad6 (except in M. musiva). The coding region of most PCGs ends with the complete termination codon TAA, except cox1, cox2, and nad3, which ended with the incomplete stop codon T (Suppl. material 1: table S2).

Statistics on the relative synonymous codon usage (RSCU) of the four species showed distinct bias and similar codon usage patterns. The most frequently used codons are UUA (Leu2), while the least commonly used codons are AAC (Asn), GAC (Asp), UGC (Cys), CAC (His), AUC (Ile), UUC (Phe), and UAC (Tyr) (Fig. 2). These results indicate that the codons of the mitochondrial protein-coding genes of Menida prefer the codon ending with A/T.

Figure 2. 

Relative synonymous codon usage (RSCU) in the mitogenomes of four Menida species.

To further investigate the codon usage bias among Pentatominae species, we analyzed the correlations between ENC (effective number of codons), the GC content of all codons, and the GC content of the third codon positions. We found a positive correlation between ENC and GC content for all codons (R² = 0.9199) and the third codon positions (R² = 0.959) (Fig. 3). These results are consistent with prevailing neutral mutational theories, in which genomic GC content is the most significant factor in determining codon bias among organisms.

Figure 3. 

Evaluation of codon bias in the mitochondrial genomes of 40 Pentatominae species.

The values of Ka (the number of non-synonymous substitutions per nonsynonymous site), Ks (the number of synonymous substitutions per synonymous site), and Ka/Ks were calculated for each PCG, respectively (Fig. 3). The Ka/Ks ratio for all 13 PCGs were below 1.0, indicating evolution under purifying selection. The Ka/Ks ratio of atp8 was the highest, while that of cox1 was the lowest. We also observed lower Ka/Ks ratios in the genes that are usually used as a barcode, such as cox2, cox3, and cytb; it is showed that at the nucleotide and amino acid levels, these four genes had the lowest evolutionary rates (Fig. 4).

Figure 4. 

The Ka, Ks, and Ka/Ks values of protein-coding genes within Pentatominae.

The total lengths of the 22 tRNAs of the four species range between 1464 bp (M. musiva) and 1484 bp (M. metallica), and the length of 22 tRNA genes ranged from 63 to 72 bp. Fourteen tRNA genes (trnI, trnM, trnW, trnL2, trnK, trnD, trnG, trnA, trnR, trnN, trnS1, trnE, trnT, trnS2) are coded on the J-strand and eight (trnQ, trnC, trnY, trnF, trnH, trnP, trnL1, trnV) on the N-strand. We found that only trn S1 lacked the dihydrouridine (DHU) arm, and the remaining 21 tRNA genes can form a typical cloverleaf structure in the four species. All tRNAs in the four mitogenomes use the standard anticodon. Among all the tRNAs of the four species in Menida, trnH has the weakest conservatism compared with other genes. In addition, 16 wobble G-U pairs were found in 22 tRNAs of Menida (Fig. 5), which usually need three-dimensional structure to stabilize.

Figure 5. 

Potential secondary structure of tRNA in Menida musiva. The conserved sites within Menida were marked in orange.

The two rRNA genes (12S rRNA and 16S rRNA) are encoded on the N-strand in these species. The 16S rRNA gene, ranging from 1277 to 1285 bp in size, is located at a conserved position between trnL1 and trnV. The 12S rRNA (795–804 bp) was found between trnV and the control region. The complete secondary structures of the 12S rRNA and 16S rRNA genes are shown in Figs 6, 7, respectively. In Menida, 16S rRNA contained 78.49% conserved sites and 12S rRNA contained 78.17% conserved sites.

Figure 6. 

Potential secondary structure of 16S rRNA in Menida musiva. The conserved sites within Menida were marked in orange.

Figure 7. 

Potential secondary structure of 12S rRNA in Menida musiva. The conserved sites within Menida were marked in orange.

The control regions located between 12S rRNA and trnI of the four species showed more variation in length, and the length ranged from 686 to 2,002 bp. This variation leads to the difference in the total length of its mitochondrial genome. The AT content in the control area of M. musiva (82.82%) was significantly higher than that of the other three species. The longest repeating unit length (284 bp) was found in M. metallica. However, no tandem repeats were detected in M. violacea (Fig. 8).

Figure 8. 

Organization of the control region in the four mitochondrial genomes. The tandem repeats are showed by the magenta circle with repeat length inside. The orange boxes indicate the length of the sequence of the control region.

Before constructing the phylogenetic tree, we performed saturation and heterogeneity analysis on two data sets. Saturation analysis showed that the sequence was not saturated (Iss < Iss. c, and p < 0.05) (Suppl. material 1: fig. S1). Heterogeneity analysis of the two data sets shows that the composition of the sequence has low heterogeneity (Fig. 9). These two data sets are suitable for further phylogenetic studies.

Figure 9. 

Analysis of heterogeneity of sequence divergence for two datasets (PRT and PR). The heterogeneity of the corresponding sequence relative to other sequences increases as the indicated color becomes lighter. The species with relatively higher sequence heterogeneity are shown.

We constructed phylogenetic trees of Pentatominae based on the two data sets using BI and ML (Fig. 10). The topological structure of the four trees was highly consistent, and most clades had high posterior probabilities. The phylogenetic positions of the Pentatominae are as follows: (Neojurtina + ((Caystrini + Halyini) + (Eysarcorini + (Carpocorini + ((Palomena + Nezara) + (Anaxilaus + (Glaucias + Plautia))))) + ((Placosternum + Cappaeini) + (Sephelini + ((Catacanthini + Strachiini) + (Pentatoma + (Hoplistoderini + Menidini))))))). The species Neojurtina typica Distant, 1921 was the earliest diverged lineage within Pentatominae. Other species of Pentatominae were scattered in the phylogenetic tree. Placosternum and Cappaeini form a sister-group relationship, and the phylogenetic tree also strongly supports the monophyly of Pentatoma. Caystrini and Halyini form a sister group relationship. At the same time, our phylogenetic relationship also shows that the genus Menida and Hoplistoderini are closely related within Pentatominae. The four Menida species are well grouped; M. metallica and M. violacea are closely related, and M. lata has the longest differentiation time compared to the other species.

Figure 10. 

Phylogenetic relationships inferred by the BI and ML method based on the PRT and PR datasets. Numbers on nodes are the posterior probabilities (PP).