The Complete Mitogenome of Pyrrhocoris tibialis (Hemiptera: Pyrrhocoridae) and Phylogenetic Implications

We determined the complete mitogenome of Pyrrhocoris tibialis (Hemiptera: Heteroptera: Pyrrhocoridae) to better understand the diversity and phylogeny within Pentatomomorpha, which is the second largest infra-order of Heteroptera. Gene content, gene arrangement, nucleotide composition, codon usage, ribosomal RNA (rRNA) structures, and sequences of the mitochondrial transcription termination factor were well conserved in Pyrrhocoroidea. Different protein-coding genes have been subject to different evolutionary rates correlated with the G + C content. The size of control regions (CRs) was highly variable among mitogenomes of three sequenced Pyrrhocoroidea species, with the P. tibialis CR being the largest. All the transfer RNA genes found in Pyrrhocoroidea had the typical clover leaf secondary structure, except for trnS1 (AGN), which lacked the dihydrouridine arm and possessed an unusual anticodon stem (9 bp vs. the normal 5 bp). A total of three different phylogenetic relationships among the five super-families of Pentatomomorpha were obtained using three analytical methods (MrBayes and RAxML under site-homogeneous models and PhyloBayes under a site-heterogeneous CAT + GTR model) and two mitogenomic datasets (nucleotides and amino acids). The tree topology test using seven methods statistically supported a phylogeny of (Aradoidea + (Pentatomoidea + (Lygaeoidea + (Pyrrhocoroidea + Coreoidea)))) as the best topology, as recognized by both RAxML and MrBayes based on the two datasets.


Introduction
Pentatomomorpha (Insecta: Hemiptera) is the second largest among the seven infra-orders of Heteroptera, with >14,000 species in 40 approximate families [1]. Most of Pentatomomorpha insects are phytophagous and major pests in agriculture, forestry, and livestock, while a few species of them are bloodsucking or predatory [1]. Pentatomomorphahas been divided into five super-families: Aradoidea, Pentatomoidea, Coreoidea, Lygaeoidea, and Pyrrhocoroidea. The super-families, except the Aradoidea, are grouped as Trichophora [1,2]. During the past several decades, phylogenetic relationships among the five super-families within Pentatomomorpha have been extensively explored based on the morphological and molecular data. For example, using morphological evidence, Henrysuggested Failed and unsatisfactory sequencing fragments were cloned into the pEASY-T1 vector (TransGen Biotech, Beijing, China), and then conducted bidirectional sequencing.

Annotation and Sequence Analysis
To ensure the accurate reading of the sequence, the sequencing results were manually corrected, which ensures the consistency of base recognition and an original peak map. After removing the vector and primer sequences, the bidirectional sequencing results of each PCR fragment were assembled to obtain the correct sequence of the target fragment using BioEdit 7.2 (https://bioedit.software.informer. com/). Mitochondrial protein-coding genes (PCGs), ribosomal RNA genes (rRNAs), and transfer RNA (tRNA) genes of P. tibialis were annotated by using the methods described in our previous studies [7,22]. MEGA 6.06 [23] was used to analyze base composition, codon usage, synonymous mutations (Ks), and non-synonymous mutations (Ka). Asymmetry of base composition between two chains was calculated as the following formula [24]:

Phylogeneticsand Analyses of Sequence Heterogeneity
Twenty-eight Pentatomomorpha species with complete or nearly complete mitogenomes were used in phylogenetic analyses, which represented five super-families and 17 families. Two species of Cimicomorpha, Lygus lineolaris and Apolygus lucorum, were used as outgroups. Details of the species used in this study are listed in Table S2. The complete sequences of each gene were used for phylogenetic analysis (excluding stop codons of the PCGs). All PCGs were aligned based on amino acid sequence alignments in MAFFT (http://mafft.cbrc.jp/alignment/server/). The aligned sequences were concatenated as two matrices used in phylogenetic analyses: (1) The P123 matrix, including all the three codon positions of PCGs, and (2) the amino acid (AA) matrix, including the amino acid sequence of PCGs. Next, GBlocks [25] was used to remove vacancies and blur sites for two aligned sequence matrices. The aligned single PCGs were combined to obtain the final mitochondrial gene set. DAMBE analysis [26] results showed that the three codon sites of 13 PCGs were not significantly saturated (Table S3). All codons, thus, were used in phylogenetic analysis. PartitionFinder 1.1.1 [27] was adopted to test the optimal partition and evolutionary model for sequence data, and the results were used for downstream phylogenetic analyses (Table S4).
Subsequently, for the nucleotide (P123) and amino acid (AA) datasets of 13 PCGs, three methods (RAxML, MrBayes, and PhyloBayes) were used to construct phylogenetic relationships within Pentatomomorpha. All phylogenetic analyses were performed on the CIPRES Science Gateway 3.3 [28] online platform. Machine learning analysis was performed using the GTRGAMMAI model in RAxML-HPC2 on XSEDE 8.0.24 [29], and the branch reliability was evaluated with 1000 rapid bootstraps. Bayesian analysis was implemented using MrBayes 3.2.2 [30], with four independent Markov chains, including three hot chains and one cold chain, running a 1 × 10 8 generation simultaneously. For every 1000 runs of collection sampling, when estimated sample size is greater than 100 and the potential scale reduction factor is close to 1.0. The two analytical processes are considered to be a stable state. Removing 25% aged samples, the remaining samples were used to construct the 50% consensus tree, and the Bayesian posterior probability (BPP) was calculated. Another Bayesian analysis under a site-heterogeneous model was implemented using PhyloBayes MPI 1.5a on the CIPRES webserver [28]. After removing constant sites from the alignment, two independent chains starting from a random tree were run under the CAT + GTR model.
The heterogeneity of sequence divergence within the two datasets (AA and P123) was analyzed using AliGROOVE [31] with the default sliding window size. Indels in the nucleotide datasets were treated as ambiguities and a BLOSUM62 matrix was used as the default amino acid substitution matrix. The metric establishes pairwise sequence distances between the individual terminal branch or subclades with the terminal branch outside of the focal group. The distances were then compared to distances over the entire data matrix. The metric values can vary between -1, if distances are very different from the average for the entire data matrix to +1 for distances that match the average for the entire matrix. This provides an indirect measure of heterogeneity of a given sequence or clade with respect to the full dataset.

Genome Organization
The P. tibialis mitogenome was a typical closed-circular DNA molecule with 16,577 bp in size, and contained 37 typical mitochondrial genes, i.e., 13 PCGs, 22 tRNAs, and two rRNAs ( Figure 1 and Table S5). The size of the P. tibialis mitogenome was similar that of the same-family species Dysdercus cingulatus (16,249 bp), but much larger than Physopelta gutta (Largidae, 14,935 bp). The size difference among the three pyrrhocorid mitogenomes was primarily due to variable length of the putative control region (CR), which ranges from 224 bp in P. gutta to 1620 bp in P. tibialis. The P. tibialis mitogenome is highly compact in genome size as that in other animals, with seven gene overlapping regions involved in a total of 33 nucleotides. The largest gene overlapping region (8 bp) is located between trnW and trnC, as observed in the other two Pyrrhocoroidea species [2]. Except for a large non-coding region (mitochondrial CR), there were 14 small non-coding intergenic regions. The intergenic spacer between trnP and trnT was 345 bp. The remaining 13 regions ranged from 1 bp to 9 bp, with a total of 63 bp. Because of the inconsistency of phylogenetic relationships obtained by different data sets and 143 phylogenetic analytical methods, in order to determine which phylogenetic tree is trustworthy, the 144 IQ-TREE web server (http://iqtree.cibiv.univie.ac.at/) was used to test the tree topological structure 145 [32]. The data sets of P123 and AA were analyzed by KH (Kishino-Hasegawa test) [33], SH

148
Furthermore, 1000 replicates were set.    157 tibialis mitogenome is highly compact in genome size as that in other animals, with seven gene 158 overlapping regions involved in a total of 33 nucleotides. The largest gene overlapping region (8 bp) 159 is located between trnW and trnC, as observed in the other two Pyrrhocoroidea species [2]. Except

170
With some notable exceptions, within the class Insecta, the order of the mitochondrial genes is 171 highly conserved and has led to the proposal of an ancestral gene order. A translocation that trnT 172 was followed by trnP was observed in the mitogenome of three species of Pyrrhocoroidea.

173
Comparing with the original gene order in mitogenomes of arthropods, the original trnT-trnP is 174 rearranged as trnP-trnT in the P. tibialis mitogenome, and there was a large non-coding sequence 175 between trnPand trnT (Figure 1 and 2). This rearrangement was also observed in the other two 176 pyrrhocorid species, which indicates that this kind of T-P rearrangement is specific to the With some notable exceptions, within the class Insecta, the order of the mitochondrial genes is highly conserved and has led to the proposal of an ancestral gene order. A translocation that trnT was followed by trnP was observed in the mitogenome of three species of Pyrrhocoroidea. Comparing with the original gene order in mitogenomes of arthropods, the original trnT-trnP is rearranged as trnP-trnT in the P. tibialis mitogenome, and there was a large non-coding sequence between trnP and trnT (Figures 1 and 2). This rearrangement was also observed in the other two pyrrhocorid species, which indicates that this kind of T-P rearrangement is specific to the Pyrrhocoroidea. Due to the unique sort of model of mitochondrial genes in the Pyrrhocoroidea, we speculate that this phenomenon in Pyrrhocoroidea likely occurred after the Pyrrhocoroidea differentiated from ancestors. The gene rearrangement was seldom discovered in Heteropteraso far. It was limited to Aradoidea, Reduvioidea, and Pyrrhocoroidea. On the contrary, gene rearrangements are common in Auchenorrhyncha.

187
The unexpected non-coding region between trnP and trnT was further verified by PCR 188 amplification and sequencing using species-specific primers. Non-coding regions were also found 189 in the other two pyrrhocorid species, despite their relatively small size ranging from 1 bp to 60 bp 190 in P. gutta and D. cingulatus, respectively. This rearrangement can be explained by the 191 duplication-random deletion model that is frequently used to explain the diversity of gene 192 rearrangements in metazoan mitogenomes [37,38]. This model considers that some genes first use

198
The nucleotide composition of the P. tibialis mitogenome was significantly biased toward A 199 and T. The total A+T content of the J (H)-strand was 75.9%, which is similar to that of the other two 200 pyrrhocorid species (Table 1) The unexpected non-coding region between trnP and trnT was further verified by PCR amplification and sequencing using species-specific primers. Non-coding regions were also found in the other two pyrrhocorid species, despite their relatively small size ranging from 1 bp to 60 bp in P. gutta and D. cingulatus, respectively. This rearrangement can be explained by the duplication-random deletion model that is frequently used to explain the diversity of gene rearrangements in metazoan mitogenomes [37,38]. This model considers that some genes first use slipped-strand mispairing to produce a repetitive gene, and then delete a redundant gene or transfer to a pseudogene randomly during the process of duplication. Thus, the model can explain gene translocation, the small non-coding region generated around the start point, length heterogeneity, the extreme variation of the replacing rate, and a repetitive copy of the tRNA gene.

Nucleotide Composition and Codon Usage
The nucleotide composition of the P. tibialis mitogenome was significantly biased toward A and T. The total A + T content of the J (H)-strand was 75.9%, which is similar to that of the other two pyrrhocorid species (Table 1). Among 13 PCGs, cox1 had the lowest A + T content (69.02%), while the highest was 86.27% in atp8. The analysis of the nucleotide composition at each codon position of the concatenated 13 PCGs of C. tetraspilus demonstrated that the third codon position (86.19%) had an A + T content higher than that of the first (70.34%) and second (67.46%) positions. The similar nucleotide composition patterns were also observed in other Coreoidea species [7]. Similar to other sequenced Pentatomomorpha insects, the base composition of the P. tibialis mitogenome was clearly biased toward A + T ( Table 1). The A + T content of the J (H)-strand in the P. tibialis mitogenome was 75.91%, with 75.15% in 13PCGs, 79.94% in rrnL, 77.91% in rrnS, and 73.46% in the CR. The A + T content of three codon sites in 13 PCGs was significantly different, with the third codon site showing much higher A + T content than that of the first and the second sites. AT-skew of the J (H)-strand in the P. tibialis mitogenome was positive (0.101), while the GC-skew was negative (−0.181).
The high A + T content and nucleotides bias in the P. tibialis mitogenome were also reflected in codon usage of PCGs ( Figure 3). Relative synonymous codon usage (RSCU) analysis indicated that all the 62 mitochondrial codons of invertebrates were used in the P. tibialis mitogenome. However, the appearance frequency of codons that ended with A or T was much higher than that of the other synonymous codons, which indicates that codons with rich AT content were frequently used ( Figure 3). Four AT-rich codons (TTT, TTA, ATT, and ATA) were the most frequently used codons in the P. tibialis mitogenome, which accounts for 31.8%. This pattern of codon usage in the P. tibialis mitogenome was highly similar to that of the other two species (D. cingulatus and P. gutta) within the super-family Pyrrhocoroidea.
Genes 2019, 10, x FOR PEER REVIEW 6 of 13 in the P. tibialis mitogenome, which accounts for 31.8%. This pattern of codon usage in the P. tibialis 219 mitogenome was highly similar to that of the other two species (D. cingulatus and P. gutta) within 220 the super-family Pyrrhocoroidea.   Except for cox1 that started with 'TTG', the remaining 12 PCGs began with 'ATN' in the P.
To analyze the evolutionary patterns of 13 PCGs in three species of Pyrrhocoridae, Ka, Ks, Ka/Ks and GC percentage of each PCG were calculated, respectively ( Figure 4). The average Ks values of the three species were similar among 13 PCGs, and all values were around one. Conversely, Ka showed a great variation, of which atp8 was the largest, which indicates that this gene had the fastest evolutionary rate. Moreover, nad2 and nad6 also presented a faster evolutionary rate, while cox1, cox3, and cob showed the slowest evolutionary rate. Ka/Ks values of all PCGs were less than 1 (<0.53), which suggested that these genes evolved under a strong purifying selection due to a functional constraint. Ka/Ks values of 13 PCGs were significantly negatively correlated with the GC percentage (R 2 = 0.72, p < 0.01), which indicates that the variation of GC content may have caused different evolution patterns of mitochondrial PCGs. Alternatively, various evolutionary patterns promoted changes of GC content.   arm, the amino acid acceptor arm, the anticodon arm, and the dihydrouridine arm ( Figure S1).

Ribosomal and Transfer RNAs
The size of rrnL and rrnS genes in the P. tibialis mitogenome were 1266 bp and 824 bp, respectively, and the content of AT was 79.94% and 77.91%, which were very similar to that of the other two pyrrhocorid species (Table 1). The total length of the 22 tRNA genes in the P. tibialis mitogenome was 1449 bp, among which the longest was 72 bp and the shortest is only 63 bp, with an average length of 66 bp. All 22 tRNAs can form a classical clover structure, including the TψC arm, the amino acid acceptor arm, the anticodon arm, and the dihydrouridine arm ( Figure S1). Some of tRNA's dihydrouridine (DHU) arm (trnQ, trnY, trnC, trnG, trnF, trnH, and trnP), amino acid acceptor arm (trnQ, trnC, trnG, trnA, and trnF), and anticodon arm (trnH and trnT) showed individual base mismatches, which was common in insect mitogenomes. Nevertheless, it was worth noting that the DHU arm of trnS1 not only existed mismatches of G-U and A-C, but also the length of the stem was only 3 bp. The loop was only composed of three nucleotides, which suggested that the arm may not exist. In fact, the lack of a DHU arm in trnS1 is common in sequenced mitogenomes of Metazoans [11]. DHU arms of trnS1 seemed to be missed in many sequenced insects including the other two pyrrhocorid species [2].

Non-Coding Regions
The CR of the P. tibialis mitogenome was located between rrnS and trnI, with a length of 1620 bp and an AT content of 73.46% (Table S5). P. tibialis and D. cingulatus had similar CR length (1620 bp and 1617 bp), but CR of P. gutta was only 224 bp ( Figure S2). The CRs of these three mitogenomes had tandem repeats, and there were two types of tandem repeats in P. tibialis and D. cingulatus, while there was only one type of tandem repeat sequence in P. gutta. There was a large non-coding region close to rrnS in P. tibialis and D. cingulatus, but it was not found in P. gutta. The absence of both non-coding regions close to rrnS and a small number of tandem repeats may be the reason that length of CR in P. gutta was significantly shorter than that of the other two species of the same super-family Pyrrhocoroidea.

Phylogenetic Analysis
Based on two datasets (P123 and AA) and three analytical methods (RAxML, MrBayes, PhyloBayes), a total of six phylogenetic trees were obtained ( Figure 5, Figures S4-S9). These phylogenetic trees highly supported the monophyly of each of Pentatomoidea, Aradoidea, Lygaeoidea, Pyrrhocoridae, and Coreoidea, which has been revealed by previous studies [7,41]. A phylogeny of (Aradoidea+ (Pentatomoidea+ (Lygaeoidea+ (Pyrrhocoroidea + Coreoidea)))) was consistently supported by four phylogenetic analyses based on two methods (RAxML and MrBayes) and two datasets (Phylogeny 1 in Figure 5). However, phylogenetic trees constructed by the PhyloBayes method were different from that of both RAxML and MrBayes. The AA-tree was (Aradoidea + (Pentatomoidea + (Coreoidea + (Pyrrhocoroidea + Lygaeoidea)))), and the P123-tree was ((Aradoidea + Pentatomoidea) + (Coreoidea + (Pyrrhocoroidea + Lygaeoidea))). Our results indicated that the degrees of heterogeneity of the amino acid dataset were lower than those of the nucleotide datasets ( Figure S4), as similarly reported in previous studies [19]. However, both AA and P123 datasets generally showed a high site-homogeneous pattern (or no significant site heterogeneity) within datasets used in this study (all AliGROOVE scores between two species > 0.0), despite a lower heterogeneity of the AA dataset than that of the P123 dataset. Furthermore, all of topology test results presented the greatest test values in phylogeny 1 based on both P123 and AA datasets (Table 2), which supports a phylogeny of (Aradoidea+ (Pentatomoidea+ (Lygaeoidea+ (Pyrrhocoroidea + Coreoidea)))). In general, the use of the CAT + GTR model, the site-heterogeneous mixture model (CAT-based model), implemented in PhyloBayes tends to reduce tree reconstruction artifacts [42][43][44], and then shows significant improvement over site-homogenous models in the reconstruction of the phylogeny of Heteroptera [14,19]. Nevertheless, in this study, due to a low heterogeneous sequence divergence between heteropteran mitochondrial AA and P123 datasets, low statistical values of topology tests were consistently presented in seven test methods using PhyloBayes. This observation suggested that robustness of tree topology was not improved by the CAT + GTR model of PhyloBayes. RAxML and MrBayes trees, which are phylogenetic analyses under a site-homogenous model, were presented to be stable. These trees consistently presented a phylogenetic relationship at the super-family level within Pentatomomorpha, which is similar to previous studies [7,14,19,45]. The trees were largely congruent with previous results based on the morphological data [6]. Therefore, compared to PhyloBayes, RAxML and MrBayes may be more plausible and phylogenetic within the suborder Pentatomomorpha, i.e., the phylogeny of (Aradoidea + (Pentatomoidea + (Lygaeoidea + (Pyrrhocoroidea + Coreoidea)))) was more reliable than the remaining two phylogenies obtained in this study. phylogenetic relationships within Lygaeoidea, only one family, Lygaeidae [19] or Geocoridae [14], 337 was included in their analyses. For the two families Berytidae and Colobathristidae, which family 338 was closer to the internal-branch families was still controversial, as reported in previous studies 339 [14,19]. Herein, excluding the analysis based on P123 and PhyloBayes, the remaining five analyses found that Geocoridae did not cluster together with Berytidae, but it was a sister lineage with the 348 family Malcidae [19]. Therefore, the present study based on the limited taxa was difficult to well   [34]. SH, Shimodaira-Hasegawa test [35]. WKH, weighted KH test [34]. WSH, weighted SH test [35]. ELW, Expected Likelihood Weight [36]. AU, approximately unbiased test [37].
Internal phylogenetic relationships within each of three super-families (Pentatomoidea, Aradoidea, and Pyrrhocoroidea) were stable, but internal phylogenetic relationships within Lygaeoidea and Pentatomoideawereunstable (Table S7). Within the super-family Lygaeoidea, phylogenetic relationships among five families (Malcidae, Lygaeidae, Geocoridae, Berytidae, and Colobathristidae) were variable depending on different analytical methods and datasets, which indicates the necessity of further study for phylogenetic relationships within Lygaeoidea. Three of six phylogenetic analyses (both P123 and AA trees constructed by MrBayes, P123 tree constructed by RAxML) consistently supported a phylogeny of ((Malcidae + Lygaeidae) + the other families) Figures S3-S8), while two AA-trees based on RAxML and PhyloBayes consistently presented a closer relationship between Lygaeidae and Geocoridae. Despite recent studies explored phylogenetic relationships within Lygaeoidea, only one family, Lygaeidae [19] or Geocoridae [14], was included in their analyses. For the two families Berytidae and Colobathristidae, which family was closer to the internal-branch families was still controversial, as reported in previous studies [14,19]. Herein, excluding the analysis based on P123 and PhyloBayes, the remaining five analyses consistently supported a closer phylogenetic relationship of the family Berytidae with internal branches within Lygaeoidea. Notably, it was difficult to demonstrate the phylogenetic location of Geocoridae in the super-family Lygaeoidea in this study, due to a severely confused topology between Geocoridae and the other families. Within the super-family Lygaeoidea, the phylogenetic location of Geocoridae was unstable. In the current study, only one species was included for each family (Geocoridae or Berytidae). However, Liu et al. added one species (Metatropis longirostris) belonging to the family Berytidae of phylogenetic trees in the super-family Lygaeoidea, which found that Geocoridae did not cluster together with Berytidae, but it was a sister lineage with the family Malcidae [19]. Therefore, the present study based on the limited taxa was difficult to well systematically infer the phylogenetic relationships at the family level, but several interesting issues involved in phylogenetics mentioned above within the super-family Lygaeoidea deserve to be further clarified in the future.
For the super-family Pentatomoidea, incongruent phylogenetic relationships were observed between topologies constructed based on the P123 and AA datasets. The former consistently presented phylogenetic relationships among two analytical methods (RAxML and MrBayes): ((Plataspidae, + (((Dinidoridae + Tessaratomidae) + Cydnidae) + Pentatomidae)) + Urostylididae), while the latter consistently exhibiteda phylogeny of (((Plataspidae + ((Dinidoridae + Tessaratomidae), Cydnidae)) + Pentatomidae) + Urostylididae). On account of lower degrees of heterogeneity of the AA dataset than that of the P123 dataset, it is, thus, speculated that phylogenetic relationships among families based on the nucleic dataset may be unfavorable in the current study. In addition, the PhyloBayes analysis based on the AA dataset consistently presented the same phylogenetic relationship with the AA-tree based on RAxML and MrBayes, which further supports (((Plataspidae + ((Dinidoridae + Tessaratomidae), Cydnidae)) + Pentatomidae) + Urostylididae) within Pentatomoidea in this study. Notably, the confusion of taxonomic relationships at the family level within Pentatomoidea needs to be further resolved with a denser taxa sampling by including molecular and morphological data in the future.

Conclusions
In conclusion, the complete mitogenome of P. tibialis was determined in this study, which further enriched the number of mitogenomes of Pyrrhocoroidea. We analyzed the main features of the P. tibialis mitogenome, and provided a comparative analysis with two other pyrrhocorid species. Furthermore, phylogenetic analyses were conducted for the five super-families within Pentatomomorpha based on nucleic and amino acid datasets. In this study, a low heterogeneity in base composition and contrasting evolutionary rates were detected among the five super-families within Pentatomomorpha, which results in no improvement of phylogenetic inference under a site-heterogeneous mixture model. Phylogenetic analyses based on mitogenomic data supported the monophyly of each super-family within Pentatomomorpha. Aphylogenetic relationship within this suborder was preliminarily proposed. This study is valuable for further understanding phylogenetic relationships among super-families within Pentatomomorpha. Notably, due to single species in several families and a single sequenced mitogenome per species, it likely limited the achievement of more phylogenetic information. Thus, to obtain more clarity, it is important to sequence more species per family and multiple mitogenomes per species in the future.
Supplementary Materials: The following are available online at http://www.mdpi.com/2073-4425/10/10/820/s1. Table S1: PCR primers used in this study. Table S2: List of the species included in the present study. Table S3: Saturation test implemented in DAMBE. Table S4: The best partitioning schemes and substitution models selected by the PartitionFinder for the P123 and AA datasets. Table S5: Annotation and organization of the complete mitochondrial genome of Pyrrhocoris tibialis. Table S6: Start and stop codons and lengths of protein-coding genes (PCGs), and two rRNA genes from three Pyrrhocoroidea species. Table S7: Phylogenetic relationships within Lygaeoidea and Pentatomoidea based on different datasets and analytical methods. Figure S1: Putative secondary structures of the 22 tRNA genes identified in the mitochondrial genome of Pyrrhocoris tibialis. Figure  S2: Organization of the A + T-rich region in Pyrrhocoroidea mitochondrial genomes. The location and copy number of tandem repeats are shown by the colored oval with Arabic numerals inside. Non-repeat regions are indicated by the colored box with the sequence size inside. Figure S3: Heterogeneous amino acid (A) and nucleotide (B) sequence divergence within heteropteran mitochondrial genomes. The mean similarity score between sequences is represented by a colored square, which ranges from −1 to +1. A lower value indicates a greater difference in rates from the remainder of the data set, i.e., heterogeneity (red coloring), while a greater value indicates rates that match all other comparisons (blue coloring). Figure S4: MrBayes phylogenetic relationships among five Pentatomomorpha super-families. Phylogenetic analysis is based on the concatenated nucleotide sequences of 13 mitochondrial PCGs. Numbers on branches are Bayesian posterior probabilities. Figure S5: RAxML phylogenetic relationships among five Pentatomomorpha super-families. Phylogenetic analysis is based on the concatenated nucleotide sequences of 13 mitochondrial PCGs. Numbers on branches are bootstrap values. Figure S6: PhyloBayes phylogenetic relationships among five Pentatomomorpha super-families. Phylogenetic analysis is based on the concatenated nucleotide sequences of 13 mitochondrial PCGs. Numbers on branches are Bayesian posterior probabilities. Figure S7: MrBayes phylogenetic relationships among five Pentatomomorpha super-families. Phylogenetic analysis is based on the concatenated amino acid sequences of 13 mitochondrial PCGs. Figure S8: RAxML phylogenetic relationships among five Pentatomomorpha super-families. Phylogenetic analysis is based on the concatenated amino acid sequences of 13 mitochondrial PCGs. Figure S9: PhyloBayes phylogenetic relationships among five Pentatomomorpha super-families. Phylogenetic analysis is based on the concatenated amino acid sequences of 13 mitochondrial PCGs.