A Molecular Phylogeny of Hemiptera Inferred from Mitochondrial Genome Sequences

Classically, Hemiptera is comprised of two suborders: Homoptera and Heteroptera. Homoptera includes Cicadomorpha, Fulgoromorpha and Sternorrhyncha. However, according to previous molecular phylogenetic studies based on 18S rDNA, Fulgoromorpha has a closer relationship to Heteroptera than to other hemipterans, leaving Homoptera as paraphyletic. Therefore, the position of Fulgoromorpha is important for studying phylogenetic structure of Hemiptera. We inferred the evolutionary affiliations of twenty-five superfamilies of Hemiptera using mitochondrial protein-coding genes and rRNAs. We sequenced three mitogenomes, from Pyrops candelaria, Lycorma delicatula and Ricania marginalis, representing two additional families in Fulgoromorpha. Pyrops and Lycorma are representatives of an additional major family Fulgoridae in Fulgoromorpha, whereas Ricania is a second representative of the highly derived clade Ricaniidae. The organization and size of these mitogenomes are similar to those of the sequenced fulgoroid species. Our consensus phylogeny of Hemiptera largely supported the relationships (((Fulgoromorpha,Sternorrhyncha),Cicadomorpha),Heteroptera), and thus supported the classic phylogeny of Hemiptera. Selection of optimal evolutionary models (exclusion and inclusion of two rRNA genes or of third codon positions of protein-coding genes) demonstrated that rapidly evolving and saturated sites should be removed from the analyses.

Hemiptera is one of the largest insect orders, comprising more than 50,000 described species. The extraordinary diversity in terms of morphology and lifestyle adaptions has long attracted the attention of evolutionary biologists and systematists. For example, some hemipteran species display bizarre morphology, some are brilliantly colored, and some produce cuticular waxes (e.g., the strangely protruding head of Fulgora laternaria and the white wax of Geisha distinctissima) [23]. Triatoma dimidiata is the vector of Chagas disease, which is a predominantly chronic disease affecting millions of people [24]. Due to high reproductive potentials, capabilities of dispersal, and transmission of plant viral diseases, some delphacids have caused considerable damage to grain production and are identified as the causes of rice famines in several Asian countries [25].
The hemipteran infraorder Fulgoromorpha includes more than 9000 described species in ,20 families, and all taxa in Fulgoromorpha are terraneous and plant-feeding. The monophyly of Fulgoromorpha is well supported by morphological and molecular data [36,37]. However, the taxonomic position of Fulgoromorpha within Hemiptera is still controversial. Previous morphological studies suggested that Fulgoromorpha and Cicadomorpha formed Auchenorrhyncha, and that Auchenorrhyncha is more closely related to Coleorrhyncha and Sternorrhyncha than to Heteroptera [38]. Cobben [39] suggested that both Heteroptera and Fulgoromorpha form the sister clade to (Sternorrhyncha,Cicadomorpha) according to a cladistic study of morphological traits. Hamilton [33] examined the phylogenetic affiliations using mouthparts and features of the head and proposed (Fulgoromorpha,(Sternorrhyncha,Cicadomorpha). Emel'yanov [40] placed the clade (Cicadelloidea + Fulgoroidea) in a sister position to the clade (Cercopoidea + Cicadoidea) and suggested that Fulgoromorpha is closer to Cicadomorpha than to other hemipterans. Although the position of Fulgoromorpha is crucial in determining the phylogenetic framework for Hemiptera, few mitogenome studies have addressed this. Fragments of the mitochondrial genes encoding 16S rRNA, 12S rRNA, cytb, and cox1 of some planthopper species [37,41,42,43,44,45], and mitogenomes from three fulgoroid species [46,47,48] have already been sequenced and utilized for phylogenetic studies. However, these studies usually have neither sufficient genetic information nor broad taxonomic samples.
Here, we describe three new complete mitogenome sequences, which add a major group (Fulgoridae) and a second representative of the highly derived clade Ricaniidae. In addition, we analyzed all 49 available complete or nearly complete mitogenomes from Hemiptera, with the aim of estimating a phylogeny of the order.

Ethics Statement
No specific permits were required for the insect specimens collected for this study in China. These specimens were collected on the roadside. The field studies did not involve endangered or protected species. Pyrops candelaria, Lycorma delicatula, and Ricania marginalis are all common fulgoroid species in China and are not included in the ''List of Protected Animals in China''.

Insects
Adult specimens of Pyrops, Lycorma, and Ricania were collected in Fujian, Henan, and Zhejiang provinces, China, respectively. They were preserved in 100% ethanol and stored at 280uC in the Key Laboratory of Zoological Systematics and Evolution, Institute of Zoology, Chinese Academy of Sciences. Flies were identified by Nan Song with reference to Chou et al. [49].

DNA Extraction, PCR, Cloning, and Sequencing
After an examination of external morphology for identification, the muscle tissue under the pronotum of each specimen was used for DNA extraction. A modified method of salt-extraction protocol [50] was employed to extract DNA.
The whole mitogenome was amplified in overlapping PCR fragments. Initial PCR primers are based on Simon et al. (2006) [51]. Short regions within individual genes were amplified using QIAGEN Taq DNA polymerase (QIAGEN, China) in PCR reaction under the following conditions: 5 min at 94uC, followed by 30 cycles of 50 s at 94uC, 50 s at 50uC, and 1-3 min at 72uC. The final elongation step was continued for 10 min at 72uC. The sequences obtained from these regions were then used to design specific primers for long PCRs that allowed us to link all of the shorter regions. The large fragments (.2000 bp) were obtained using the QIAGEN Long Taq DNA polymerase (QIAGEN, China) under the following conditions: 2 min at 96uC, followed by 30 cycles of 10 s at 98uC, and 3 min at 68uC. The final elongation was continued for 10 min at 72uC. These PCR products were analyzed by 1.0% agarose gel electrophoresis.
PCR products of ,1200 bp were directly sequenced after purification. Whereas products of 1.2-3 kb were cloned into pBS-T Easy vector (QIAGEN, China), and the resultant plasmid DNA was isolated using the TIANprp Midi Plasmid Kit Purification System (QIAGEN, China) and sequenced by means of primer walking. DNA sequencing was performed using BigDye terminator chemistry and ABI 3730xl Genetic Analyzer (PE Applied Biosystems, USA). After the entire mitogenome for Pyrops had been sequenced, the primers were reused to amplify other species; however, some primers were redesigned for efficient sequencing because of minor sequence differences between species.

Sequence Assembly, Annotation, and Analysis
Raw sequence files were proofread and aligned into contigs in BioEdit 7.0.5.3 [52]. Contig sequences were checked for ambiguous base calls, and only non-ambiguous regions were used for annotation. Sequence alignment, genome assembly, and calculations of nucleotide composition were all conducted with MEGA 5 [53]. The 22 tRNA genes were identified using the tRNAScan-SE server [54]. The tRNAs not found by tRNAScan-SE were identified through comparison with the regions coding these tRNAs in other insects. Protein-coding genes and rRNA genes were determined by comparison with those of published insect mitochondrial sequences. Potential secondary structure folds in the A+T-rich region of the genome were predicted using Mfold 3.2 [55]. New mitogenome sequences obtained in this study were deposited in GenBank under accession numbers FJ006724, EU909203, and JN242415 for Pyrops, Lycorma, and Ricania, respectively (Table S1).

Phylogenetic Analyses
We conducted phylogenetic analyses using all of the currently available mitogenomes of Hemiptera along with the three newly sequenced ones (Table S1).Three orthopterans and a psocid served as outgroups. The orthopterans form a basal lineage and the psocid is in a sister clade to Hemiptera.
The nucleotide sequences of 13 protein-coding genes and 2 rRNA genes were used to reconstruct the phylogenetic relationships in Hemiptera. All protein-coding genes were aligned at the amino acid level using the default settings in ClustalW (as implemented in MEGA 5). The alignments were back-translated into the corresponding nucleotide sequences. The stop codons of protein-coding genes were excluded when aligned. Two rRNAs were respectively aligned in MEGA 5 as nucleotides. Ambiguously aligned regions of protein-coding genes and rRNA genes were checked by eye. Potential saturation in each protein-coding gene and rRNA was individually assessed, with transitions and transversions plotted against corrected genetic distance using DAMBE 5.1.2 [56]. Both average p-distances and pairwise relative -rates were calculated by PHYLTEST [57].
Phylogenetic trees were estimated from each data set using maximum parsimony (MP), maximum likelihood (ML). and Bayesian inference (BI). MP analyses were performed using PAUP* 4.0b10 [58], with gaps treated as missing data. A total of 1,000 random-addition searches using tree-bisection-reconnection were performed for each MP analysis. Bootstrap support was calculated from 1000 bootstrap replicates with 100 random additions per replicate in PAUP*. Tree statistics were also calculated in PAUP* ( Table 1). All four outgroup taxa, partial outgroup (psocid only), or long-branch taxa were respectively removed from MP analyses to search for potential long-branch attraction artifacts.
ML analyses were conducted using the program TreeFinder [59]. For each dataset, the GTR+I+G model was identified as the best-fit one under the Akaike information criterion (AIC) by the ''model proposer'' in TreeFinder. For the datasets with 3rd codon positions, each codon position was treated as a separate partition. A search-depth level of 2 was selected and bootstrap analysis was performed with 1000 replicates.
Bayesian analyses were conducted with MrBayes 3.2 [60]. Prior to the analyses we tested for the appropriate nucleotide substitution model via AIC with MrModeltest v2.3 [61] for each gene, for each codon and for all concatenated datasets (Table S2). GTR+I+G were estimated as the best-fit substitution model for all partitions. Each analysis comprised four independent Markov chains of 3,20 million generations each, sampling every 100 generations, and the first 25% discarded as burn-in. The datasets ALL_123 and PCG_123 were partitioned first by gene, then by codon position. All model parameters were unlinked across partitions. Markov chain stationarity was considered to be reached when the average standard deviation of split frequencies fell below 0.01 and potential scale reduction factor values approached 1.0 [62].

Genome Features
The complete mtDNA sequences of Pyrops, Lycorma, and Ricania were determined to be 16021 bp, 15946 bp, and 15698 bp in size, respectively. The three genome sizes are well within the observed range of insect mitogenomes (14-19 kb) and the length variation occurs in the A+T-rich regions, which range from 1324 to 1642 nucleotides. Like other insect mitogenomes, three fulgoroid mitogenomes contain the typical 13 protein-coding genes (PCGs), 22 tRNAs, two rRNAs, and A+T-rich region (Fig. 1). The position and orientation of the mitochondrial genes are the same as those found in the putative ancestral insect mitogenome [62]. A summary of the mitochondrial genes of Pyrops, Lycorma, and Ricania is given in Table S3.
The A+T content at the third codon position is higher than that at the first or second codon position. All of the majority-strand genes favor A and C (A-skew and C-skew are 0.24-0.28 and 0.20-0.28, respectively).
We determined the codon usage of 13 protein-coding genes ( Table 3) and found that nine codons (AAT-Asn, TTT-Phe, ATT-Ile, TTA-Leu, AAA-Lys, ATA-Met, TCA-Ser, ACA-Thr, TAT-Tyr) are most frequently used in the three fulgoroids. All of these codons have an A and/or T at the third codon position. This codon bias is also caused by the A+T-rich composition of the mitogenome.
Transfer RNA and Ribosomal RNA Genes All 22 tRNA coding genes usually found in mitogenome of metazoans are present in these three fulgoroids. The anticodon nucleotides for the corresponding tRNA genes are identical to those of other available hemipteran mitogenomes [46,47,48,63,64]. All tRNA genes have the typical clover-leaf structure with one exception: tRNA-ser(AGN), in which the dihydrouridine arm formed a simple loop as in some other metazoan species, including most insects [1,23,65,66]. The tRNAs are all found to be between 60 bp and 75 bp in length.
The arrangements of both 16S rRNA and 12S rRNA in the fulgoroid mitogenomes are conserved. They are located between tRNA-Leu(CUN) and tRNA-Val and between tRNA-Val and the A+T-rich region, respectively (Fig. 1). The lengths of 16S rRNA and 12S rRNA are determined to be 1214-1216 bp and 717-736 bp, respectively. These lengths are similar to those of other sequenced fulgorids (1192-1219 bp for 16S rRNA and 711-747 bp for 12S rRNA in Sivaloka damnosus, Geisha, and Laodelphax striatellus) [46,47,48].

A+T-rich Region
The A+T-rich regions can be roughly divided into five sections, which is similar to other fulgoroid species [46,47,48]. In the A+T- rich region, the poly-T stretch may be important to the initiation of mtDNA replication. The stem-loop secondary structure was suggested to be the site of initiation of the secondary strand synthesis in Drosophila [67]. Such poly-T and stem-loop structures also exist in the fulgoroids' A+T-rich region. Repetitive sequences have been commonly found in insect mitogenomes [1,68], and the length variations are due to the variable number of repeat copies [69,70,71]. In the six sequenced fulgoroids, all A+T-rich regions contain tandem repeat units. But there are no any similarities in the size and nucleotide composition of repeat units among them.

Saturation Test
Because saturation in substitutions can have negative effects on phylogenetic inference [72], the levels of saturation in protein-coding genes and mitochondrial rRNA genes were separately explored (Table S4). Xia's saturation idex (Iss) was estimated and compared to the critical values assuming symmetric (Iss.cSym) and asymmetric (Iss.cAsym) topologies, and a P-value was obtained to assess statistical significance [73,74]. No saturation was detected in the 1st and 2nd codons of protein-coding genes [P(Iss , Iss.c) ,0.05]. For the 3rd codon positions, saturation was detected, with the resulting Iss (0.861) higher than Iss.cSym (0.809). We found similar results in comparisons with Iss.cAsym values. With respect to rRNAs, no saturation was detected in the 16S rRNA alignment, but Iss was greater than both Iss.cSym and Iss.cAsym (NumOUT = 16 or 32) in the 12S rRNA alignment. Therefore, there is some saturation in the third codon positions and rRNAs.

Sequence Diversity and Relative-rates Tests
We calculated average p-distances from PCG_123 between major groups using PHYLTEST (Table S5). The average genetic distances between six whiteflies and other hemipteran groups are higher than those without whiteflies involved. The distances between whiteflies and the psocopteran are close to that between whiteflies and aphids. Furthermore, we carried out pairwise relative-rates tests on PCG_123 using PHYLTEST (Table S6). The family Aleyrodoidea displays higher rates of nucleotide change in the mitochondrial protein-coding genes. The hierarchy of the rate of nucleotide substitution was whiteflies .psyllids . fulgoroids . aphids . cicadas . true bugs. This hierarchy did not change whether Orthoptera or Lepidoptera was used as outgroup.

Phylogenetic Analyses
In this study, we included 49 taxa of Hemiptera representing four higher groups. After alignment and concatenation, the protein-coding genes totalled 11205 bp and the rRNAs (16S rRNA +12S rRNA) 2173 bp (Table 1).
Five different datasets with three inference methods produced 15 phylogenetic trees. These tree topologies were highly compatible with each other. Fig. 2, 3, 4 illustrates the results of analyses from the ALL_12 dataset.
In most trees provided in this study, Hemiptera was divided into two groups (Cicadomorpha,(Fulgoromorpha,Sternorrhyncha)) and Heteroptera, rendering Auchenorrhyncha paraphyletic. Heteroptera and Sternorrhyncha were well supported (both MP and ML bootstrap support values were 100, and Bayesian posterior probability was 1.00).Within Hemiptera, Fulgoromorpha was closer to Sternorrhyncha than to Heteroptera in 11 trees. In the remaining four trees (MP_ALL_123, MP_PCG_123, MP_rRNAs and BI_rRNAs), the group (Cicadomorpha,Fulgoromorpha) was sister to Heteroptera, and Sternorrhyncha was sister clade to the remaining taxa.
Monophyly at the superfamily level within Sternorrhyncha was well supported in all of the topologies. In the MP trees, Aphidoidea was placed as the sister clade to all other taxa within Sternorrhyncha, whereas Psylloidea was found to be the sister group to the remaining sternorrhynchan lineages. But in all model-based analyses, Psylloidea was recovered as the sister group to all other taxa within Sternorrhyncha, followed by Aphidoidea.
The interfamilial relationships within Fulgoromorpha were consistently supported. The relationships (Delphacidae,(Fulgoridae,(Flatidae,(Issidae,Ricaniidae))) were resolved in all analyses except for the datasets using rRNAs. The Laodelphax was placed as the sister to all other taxa in Fulgoromorpha on all trees, and the monophyly of Pyrops + Lycorma was also supported irrespective of the analytical method used.

Effects of Long-branch Attraction
When parsimony and model-based estimates of phylogeny infer different trees, it is desirable to explore the causes of the discrepancy [75]. Long-branch attraction (LBA) is commonly cited as the reason for incongruence between different treereconstruction methods [76,77]. Maximum parsimony is considered susceptible to attraction between long branches [78]. In the present study, six representatives of Aleyrodoidea always display branches significantly longer than other hemipteran species (Fig. 2, 3, 4 and Figure S1). This raised the suspicion that the resulting topologies may be compromised by a long-branch attraction artifact.
To test whether LBA affected our phylogenetic analyses, we followed several strategies suggested by Bergsten [79]:  long branches yielded almost identical internal structure of Hemiptera: Fulgoromorpha formed a sister group to Sternorrhyncha with high support values (bootstrap value = 100), and Heteroptera alone constituted a major clade within Hemiptera. The topology was very similar to that obtained with the datasets excluding 3rd codon positions, but most higher taxa were recovered with higher support (Fig. 5). (iii) Removal of outgroups. MP analyses of datasets ALL_123 and PCG_123, excluding four outgroup taxa, produced an identical tree topology for Hemiptera (Fig. 6). The only difference was that the MP tree on PCG_123 had an increase in support for a few nodes. In our analyses of rRNAs, we found that the outgroup psocid was always nested within the ingroup. Accordingly, we also estimated trees using PCG_123 with the psocid removed. We estimated a congruent ingroup topology with those removing all outgroups.
The analyses performed above indicated that there is a strong attraction of the clade Sternorrhyncha to the outgroups, and this probably led to LBA artifacts.

rRNA Dataset
The rRNA dataset produced less-resolved topologies. The outgroup Lepidopsocidae sp. tended to group within the ingroup close to Sternorrhyncha in ML and BI analyses. In terms of the family-level relationships in Fulgoromorpha, Flatidae+Ricaniidae formed a moderately supported clade; the position of Issidae varied across different analytical methods. Coreoidea was paraphyletic in the three rRNA analyses and several other datasets that included rRNA (MP_ALL_12, ML_ALL_12, and ML_ALL_123). It is possible that the rRNA data support paraphyly of Coreoidea. Similarly, the monophyly of Cicadomorpha could not be confirmed by the rRNA dataset in MP and ML analyses. These problematic topologies might result from alignment errors or from long-branch attraction.

Mitochondrial Genome Structure
The start codon for cox1 is highly variable across insects, and frequently uses noncanonical start codons that code for amino acids other than methionine [80,81,82]. However, fulgoroid species do not share this feature. Pyrops uses ATC and Ricania uses ATA as start codon for cox1, whereas conventional methionine (start codon with ATG) is used in Lycorma. Although ATG, ATT, ATA, or ATC is universally used as a start codon in mitochondrial protein-coding genes in vertebrates and insects [83,84], there are apparent exceptions to the ATN rule [1]. For example, GTG has been found as the start codon in some insect mitogenomes: for atp8 in Bactrocera dorsalis [66] and for nad5 in Triatoma [24] and Pteronarcys princes [84]. Similarly in the mitogenome of Pyrops, we found GTG as start codon for nad1.
In accordance with other fulgoroid species, overlapping proteincoding genes are present in Pyrops, Lycorma and Ricania; a 7-bp overlap exists not only between atp8 and atp6 but also between nad4l and nad4. In this case, hairpin structures forming at the 39 end of the upstream protein's mRNA may act as a signal for the cleavage of the polycistronic primary transcript [62].

Phylogenetic Relationships within Hemiptera
Maximum-parsimony analysis, in which all characters were included, produced poorly resolved and poorly supported trees.
Excluding third codon positions tended to increase both resolution and support for nodes. This may be attributed to the generally high rate of evolution at third codon positions for Hemiptera. Maximum-likelihood and Baysian analyses produced broadly similar topologies to the MP analyses without 3rd codon positions. Concordant results between the different analytical approaches provide some confidence in the subordinal structure obtained for   Hemiptera: (((Fulgoromorpha,Sternorrhyncha),Cicadomorpha),-Heteroptera).
The monophyly of Hemiptera is well supported by nucleotide sequences of protein-coding genes, which shows that mitogenome data are very effective in resolving relationships within this group. The relationships (Cicadomorpha,(Fulgoromorpha,Sternorrhyncha)) were consistently recovered by the analyses using datasets ALL_12 and PCG_12. This suggested a paraphyletic Auchenorrhyncha, which has been proposed by previous studies [33,85,86].
The monophyletic Fulgoromorpha was well recovered in the analyses of the full dataset (PCGs + rRNAs). The phylogenetic hypothesis of (Laodelphax,((Pyrops,Lycorma),(Geisha,(Sivaloka,Ricania))) was strongly supported in the analyses without rRNAs. This result was consistent with those of other molecular studies [23]. The rRNA dataset produced slightly differing topologies ( Figure S1). Common to the trees was the placement of Laodelphax as sister taxon to other fulgoroids, which was in line with previous research [25,37,41,87,88]. On the other hand, a sister relationship between Pyrops and Lycorma was consistently resolved in all analyses. This was in accordance with the findings of Yeh et al. [37].

Long-branch Attraction Artifact
Long-branch attraction refers to the erroneous grouping of two or more unrelated branches as sister groups due to undetected parallel evolution (homoplasy) [90]. These homoplasies are more likely to occur along long branches. Parsimony analysis has been demonstrated to be vulnerable to the high level of homoplasy, and to be particularly sensitive to LBA [79]. In contrast, by correcting for saturation of substitutions, model-based methods may be relatively insensitive to LBA and infer topologies more accurately [79].
In this work, the conflict between parsimony and model-based methods was detected. We followed a similar protocol to that of Bergsten [79] to investigate possible LBA artifacts. Our results showed that the topological affinities (Sternorrhyncha,(Heteroptera,(Cicadomorpha,Fulgoromorpha))) is most likely to be the result of the attraction between the whiteflies and the distant psocid outgroup. This is supported by the fact that these two groups share similar rates of evolution and have less sequence divergence. By contrast, the topologies from MP analyses on datasets without 3rd codon positions (removal of fast evolving sites) and those from the model-based methods (which are less commonly affected by LBA) are correct.

Mitochondrial rRNA Genes
In addition to the datasets comprising protein-coding genes, we attempted to gain a better resolution of hemipteran relationships using different genes evolving at different rates, 16S rRNA and 12S rRNA. However, these data yielded the worst topologies and the weakest nodal support, indicating that they might be unsuitable for reconstructing the evolutionary relationships of higher taxa in Hemiptera. This may be due to the substitutional saturation observed in mitochondrial rRNA genes. In addition, it was previously seen that mitochondrial rRNA genes displayed significant differences among lineages in their evolutionary rates [90]. Thus, the phylogenetic inferences based on mitochondrial rRNAs were likely to be subject to a LBA artifact, and the lineage leading to Sternorrhyncha might have been incorrectly pulled towards the root of the tree.

Comparison to Previous Phylogenetic Results
Complete mtDNA sequences can be informative at deep phylogenetic levels [91], and their phylogenetic utility has been demonstrated in several insect orders (see Introduction). Our results show the paraphyly of Auchenorrhyncha, which is consistent with analyses of full-length or partial sequences of nuclear 18S rDNA [87,92,93]. The discordance between the trees we present here and those reported in previous molecular phylogenetic studies [27,92,93] relates to the position of Fulgoromorpha in Hemiptera. This may be attributed to different methods, molecular markers, or outgroups used among the phylogenetic analyses. In addition, previous molecular studies were mainly based on single 18S rDNA fragments (,1000 bp, which seems to contain insufficient information) and parsimony analysis. Both the short sequences and methods sensitive to LBA may cause problems in resolving phylogenetic relationships at the intraordinal level of insects.
Although it has been used extensively for studies inferring phylogenies, mtDNA has shortcomings that can limit its potential in recovering the phylogenetic signal [94]. In contrast, nuclear genes evolve more slowly, making them effective sources of information for the analysis of deep phylogenetic relationships. Therefore, a combination of mitochondrial genomes and nuclear genes is expected to provide more precise estimates of phylogenetic trees.

Conclusions
Hemiptera is the largest nonholometabolan insect assemblage and still contains many species that are underrepresented by complete mitochondrial sequences. Although our limited taxon sampling only can provide a preliminary phylogenetic picture of Hemiptera, our result confirms that mitogenome data are very effective in resolving deeper relationships within this order, and the inclusion or exclusion of third codon positions has a strong influence on phylogenetic reconstruction. Moreover, our results are compatible with our previous findings based on more limited taxon sampling [47,48]. Most topologies inferred in the current study appear to be more consistent with the classical hypothesis by Hamilton [33]. This work has added to current knowledge on the hemipteran phylogeny inferred from mitogenomes. Future research that integrates more taxon sampling, more mitogenome sequences, and data from other molecular markers will provide greater insight into the evolution of Hemiptera. Figure S1 Tree topologies of MP, ML, and BI inferred from different data partitions except ALL-12.

(RAR)
Table S1 List of taxa used in the phylogenetic analysis.