Mitogenomes of Nine Asian Skipper Genera and Their Phylogenetic Position (Lepidoptera: Hesperiidae: Pyrginae)

Simple Summary This is the first attempt to test the validity of morphological characters to diagnose the tribes used in the conventional taxonomy of the pyrginae skippers. Current studies on the group deal only with molecular phylogeny and its resulting classification without considering morphology. The diagnostic characters currently used in distinguishing the two tribes cannot be adopted. When a new taxonomic framework is proposed based on molecular data, reevaluation of morphological characters is suitable. Abstract In this study, complete mitochondrial genomes of nine species representing three tribes in the subfamily Pyrginae sensu lato were newly sequenced. The mitogenomes are closed double-stranded circular molecules, with the length ranging from 15,232 bp to 15,559 bp, which all encode 13 protein-coding genes (PCGs), two ribosomal RNA (rRNA) genes, 22 transfer RNA (tRNA) genes, and a control region. The orientation and gene order of these nine mitogenomes are identical to the inferred ancestral arrangement of insects. All PCGs exhibit the typical start codon ATN except for cox1 (using CGA) and cox2 (using TTG) in Mooreana trichoneura. Most of the PCGs terminate with a TAA stop codon, while cox1, cox2, nad4, and nad5 end with the incomplete codon single T. For the different datasets, we found that the one comprising all 37 genes of the mitogenome produced the highest nodal support, indicating that the inclusion of RNAs improves the phylogenetic signal. This study re-confirmed the status of Capila, Pseudocoladenia, and Sarangesa; namely, Capila belongs to the tribe Tagiadini, and Pseudocoladenia and Sarangesa to the tribe Celaenorrhini. Diagnostic characters distinguishing the two tribes, the length of the forewing cell and labial palpi, are no longer significant. Two populations of Pseudocoladenia dan fabia from China and Myanmar and P. dan dhyana from Thailand are confirmed as conspecific.


Introduction
The family Hesperiidae is one of species-richest groups of Lepidoptera, which account for one-fifth of the world's butterfly species, though the number is significantly underestimated [1]. Traditionally, they had been classified into four or five subfamilies, and genus group was adopted instead of tribe, until the groups were divided [2] and members were re-classified based on molecular phylogeny [1,3]. Currently, as many as 10 subfamilies with numerous tribes are proposed [4,5].
The conventional concept of the subfamily Pyrginae is considered polyphyletic because members of the subfamily Eudaminae, which are morphologically significant, are often embedded within Pyrginae. As a result, Pyrginae is divided into three subfamilies: Tagiadinae (with the tribes Tagiadini, Netrocorynini, and Celaenorrhini), Pyrrhopygini (with the tribes Pyrrhopygini and Jerini), and Pyrginae (with the tribes Pyrgini, Carcharodini, Achlyodini, and Erynnini). These taxonomic arrangements, however, are purely or almost exclusively based on molecular data, and subfamily-and tribe-level assignments are arbitrary. On the one hand, biological taxonomy must reflect the phylogenetic relationship. On the other hand, however, it is a general reference system, and such a splitting of groups may result in taxonomic confusion. For this reason, herein we adopt conventional Pyrginae sensu lato, the subfamily currently composed of 646 species in 86 genera worldwide [6][7][8][9][10].
Most of the phylogenetic analyses of the family Hesperiidae rely on single-gene loci analyses, both mitochondrial (such as cox1, cox2, and 16S rRNA) and nuclear (wingless and EF-1α). Although the mitogenome, in whole or part, has been widely used as a molecular marker in population genetics as well as evolutionary and phylogenetic studies, relatively few studies have attempted to study the phylogeny of Hesperiidae using the mitogenomes [12][13][14]. Data from a new perspective, such as complete mitogenomes, may help improve phylogenetic resolution for these groups. In this study, we sequenced nine additional mitogenomes of species belonging to nine genera in three tribes of the subfamily Pyrginae. In addition, we analyzed the characteristics of the nine mitogenomes in detail. Phylogenetic relationships of the family Hesperiidae were explored in combination with all the available 90 complete Pyrginae mitogenomic sequences available in GenBank. This enabled us to test the monophyly of Pyrginae, explore relationships within Pyrginae, and examine the taxonomic status of some ambiguous genera in Pyrginae.
Due to the compositional heterogeneity of the mitochondrial genome of Lepidopteran insects, most previous studies have attempted to address it by removing the third codon position of protein-coding genes, as well as excluding rRNAs and tRNAs [12][13][14][15][16][17]. The latest phylogenetic analysis by Ma et al. [18] established a dataset that contained RNAs and produced the most consistent topologies and higher node support values [19,20]. Therefore, the RNAs increased the phylogenetic resolution. Among the RNAs, tRNA is quite conservative, especially in the sequence of the stem area; its evolutionary rate is slower than the rate of other components of mitochondrial genomes. rRNAs also have higher conservatism and slower evolutionary rate. Thus, tRNA and rRNA are often used in phylogenetic analysis of various taxa, which has a great influence on the phylogenetic results [18,[21][22][23][24]. Thus, we performed phylogenetic analysis of the family Hesperiidae using several different datasets to explore the impact of the inclusion or exclusion of tRNA on phylogenetic resolution.
Current studies deal only with molecular phylogeny and its resultant classification without reevaluating diagnostic morphological characters used in conventional taxonomy. We examined the wing venation and labial palpi of five genera to test the validity of those morphological characters as diagnostic features under the current phylogenetic framework.

Sample Collection and Genomic DNA Extraction
Mitochondrial genomes of nine species belonging to nine genera of Pyrginae were sequenced. Of those nine genera, three (Satarupa, Mooreana, and Abraximorpha) were newly sequenced genera and the rest were genera sequenced previously; therefore, we used different species in those genera, except for Pseudocoladenia, to confirm the monophyly of the genera. All the species/specimen used in this study are listed in (Table 1) and were collected and stored in 100% ethanol at −20 • C in the Entomological Museum, NWAFU. The specimens were initially examined using morphological characteristics, particularly the genitalia, and confirmed via cox1 barcoding using the BOLD database [15,25,26]. The extraction of the Genomic DNA was done from the thoracic muscle (mitogenomes) using the Biospin Insect Genomic DNA Extraction Kit (Qiagen, Hilden, Germany). The NGS (Illumina HiSeq X; Biomarker Tech, Beijing, China) was employed to determine the nine mitogenomes of Pyrginae.

Bioinformatics Analyses
The extraction of the complete mitogenome sequences of the nine Pyrginae species was done using the Illumina HiSeq 2000 system by Genesky Biotechnologies Inc. (Shanghai, China). The correct identification rate of bases was very high, reaching 99.9%. First, the raw paired reads were retrieved and quality-trimmed using CLC Genomic Work bench v10.0 (CLC Bio, Aarhus, Denmark) with default parameters. The basic statistics of sequencing for each mitochondrial genome are presented in the Supplementary Materials (Table S1). The format of sequencing was Illumina, and the length of reads was 150 bp paired-ends. Then, the clean paired reads were used for mitogenome reconstruction using MITObim v1.7 software with default parameters and the mitogenome of Tagiades vajuna (KX865091) as the reference [16,27]. We selected a mitogenome, T. vajuna, as the reference and compared it with the nine mitogenomic sequences using MAFFT integrated into Geneious [28,29]. We conducted the annotation of the mitogenomes and comparative analyses following the methodology outlined above. We selected the complete mitogenome sequence of T. vajuna as a reference and used the Geneious 8.1.3 software to annotate all the various genomic features. Protein-coding genes (PCGs) were found by searching for ORFs (employing the invertebrate mitochondrial genetic code translation table5) and checking nucleotide alignments against the reference genome in Geneious 8.1.3. All RNAs (rRNAs and tRNAs) were identified using the MITOS Web Server (http://mitos.bioinf.uni-leipzig.de/index.py (accessed on 1 July 2021)), and tRNA secondary structures were visualized according to these results [27]. Finally, we used the Geneious software to visualize all gens by inspecting against the reference mitogenome. Nucleotide composition, codon usage, comparative mitogenomic architecture tables for the nine mitogenomes, and data used to plot RSCU (relative synonymous codon usage) figures were all calculated using PhyloSuite [28]. The nine newly sequenced mitogenome sequences were uploaded onto GenBank with the accession numbers as specified in Table 1.

Sequence Read Archive (SRA) Data Extraction
We downloaded and extracted the raw data of mitochondrial genomes of 86 hesperiid species from GenBank to reconstruct the phylogenetic relationships. The raw data were assembled and annotated by Geneious 8.1.3 [26]. The SRA data used in this study were obtained from GenBank and are listed in Table 2.

Phylogenetic Analysis
Phylogenetic analyses were conducted based on three datasets:  (Tables S2 and S3). Nucleotide sequences were aligned using the G-INS-i (accurate) strategy and codon alignment mode. All rRNAs were aligned in the MAFFT with the Q-INS-i strategy [29]. Poorly matched sites in the alignments were removed using Gblocks v0.91b [30]. Individual genes were also concatenated using PhyloSuite v1.2.2.

Morphological Comparison
For the morphological comparison, we chose two types of the diagnostic characters one of which is the wing venation of five species representing five tribes. Another diagnostic characters is the lateral view of the head, showing the labial palpus of the five species representing five tribes. First, we rinsed the scales of Celanorrhinus maculosus and Saranges dasahara with 10% sodium hypochlorite and then traced the structure of the wing venation with digital boards, finally obtaining pictures with a simplified drawing. We examined the labial palpus of the five species using an electron microscope. Then, we constructed all the lateral views of head showing the labial palpus through freehand drawings.

Genome Organization and Base Composition
The mitogenome of nine species are a single, covalently closed circular doublestranded DNA molecule ( Figure 1) composed of 37 coding genes [48][49][50]. The mitogenome sizes are shown in Table 3. Including the newly sequenced mitogenomes of nine species in the present study, 49 species of Hesperiidae have mitogenome data public available, with the length ranging from 15,232 bp (M. trichoneura) to 15,559 bp (Erynnis popoviana).

Genome Organization and Base Composition
The mitogenome of nine species are a single, covalently closed circular doublestranded DNA molecule ( Figure 1) composed of 37 coding genes [48][49][50]. The mitogenome sizes are shown in Table 3. Including the newly sequenced mitogenomes of nine species in the present study, 49 species of Hesperiidae have mitogenome data public available, with the length ranging from 15,232 bp (M. trichoneura) to 15,559 bp (Erynnis popoviana).  All the nine mitogenomes contain 13 protein coding genes (PCGs), 22 transfer RNA genes (tRNAs), 2 ribosomal RNA genes (rRNAs), and an AT-rich region. Like many other insect mitogenomes, its major strand codes for 23 genes (9 PCGs and 14 tRNAs), while the minor strand codes for the remaining 14 genes (4 PCGs, 8 tRNAs, and 2 rRNA genes). The orientation and gene order of these mitogenomes of nine species are identical to the hypothesized ancestral arthropod arrangement found in the insect Drosophila yakuba [51].
The nucleotide composition of the nine mitogenomes is significantly biased towards A and T, with a relative A + T content of 79.5% to 82.4% in the whole genome, 77.5% to 80.9% in the PCGs, 81.2% to 83% in the transfer RNAs, and 84.3% to 86% in the ribosomal RNAs. Nearly all AT-skew and GC-skew are negative, except for the AT-skew in Ce. aspersus (0.003) and Sat. nymphalis (0.026), showing that there are more TC than AG across the whole mitogenome (Table 3).

Protein-Coding Genes
Of all the 13 PCGs in the nine mitogenomes of the subfamily Pyrginae, nine were located on the majority strand (J-strand), and the other four PCGs were located on the minority strand (N-strand). The total length of the 13 PCGs of nine species is 11,178 bp to 11,199 bp. The A + T content of the third codon positions of the PCGs (88-96.2%) was much higher than either the first (72.8-75.4%) or second codon positions (69.6-71.3%) ( Table 4). Most PCGs are initiated by a typical ATN as the start codon; the cox1 of nine species have CGA as the start codon, and cox2 in M. trichoneura is started with TTG. The canonical TAN stop codon occurs in most PCGs. In most of the nine species, cox1, cox2, nad4, and nad5 genes use T as a truncated stop codon, except for the nad4 gene of P. dan fabia, and nad5 gene of E. popoviana, M. trichoneura, and P. dan fabia-all use TAA as a stop codon. Truncated stop codons are common in insect mitogenomes and might be completed by post-transcriptional polyadenylation [52].
Relative synonymous codon usage (RSCU) values for the nine mitogenomes are calculated and summarized in Figure 2. A + T bias is also reflected in the relative codon usage by the PCGs. The amino acid frequencies, excluding stop codons, are similar amongst the different skipper mitogenomes. The most frequently used codons across all the species are UUU (Phe), UUA (Leu), AUU (Ile), AUA (Met), and AAU (Asn), all of which are composed wholly of A or T. The results indicate a preference for NNU and NNA codons in skipper mitogenomes, which has been observed before [18,24,32].

Transfer and Ribosomal RNA Genes
All 22 standard tRNAs of the mitogenomes of the arthropod were found in the nine mitogenomes. The total length of tRNAs in the nine mitogenomes was 2121 bp (A. davidii) to 2165 bp (Ce. aspersus) ( Table 3). Most tRNAs could be folded into the canonical cloverleaf structure, except for trnS1 (AGN), with its dihydrouracil (DHU) arm forming a simple loop (all the secondary structure of tRNAs are shown in Supplementary Material Figures S1-S9), which was considered a typical feature in metazoan mitogenomes (Wolstenholme 1992). Fourteen tRNA genes were encoded by J-strand and the remaining eight were encoded by N-strand. There were base pair mismatches in the receptor arm, DHU arm, anticodon arm, and TΨC arm on the tRNAs of the nine mitogenomes, most of which were G-U mismatches, followed by U-U mismatches, and U-C, A-C, A-A, and A-G pairs (Table 5). Table 5. Mismatch among the tRNAs of the nine mitogenomes.

Species
Acceptor Arm DHU Arm Anticodon Arm TΨC Arm The two rRNA genes, rrnL and rrnS, were located between trnL (CUN) and trnV, and trnV and the A + T-rich region, respectively (Table 6). In the nine newly sequenced mitogenomes, the lengths ranged from 1347 bp (A. davidii) to 1390 bp (M. trichoneura) for the rrnL, and from 772 bp (T. menaka) to 788 bp (M. trichoneura) for the rrnS ( Table 6). The A + T content of total rRNA genes was 84.3% (Sat. nymphalis) to 86% (E. popoviana), which was higher than that in the whole genome, indicating a moderate A + T preference in the total rRNA genes.
The control region (A + T-rich region) was located between rrnS and trnM, which ranged from 278 bp to 466 bp. This region contained the highest proportion of A and T, ranging from 92.7% to 96.3%. The high A+T could be involved in the regulation of transcription and replication of the mitogenome [52,53]. Both the AT-skew and GC-skew were negative in the control region of most species (except for the AT-skew of Ce. aspersus, Sat. nymphalis), indicating a clear bias towards the utilization of T and C.

Phylogenetic Relationships
Only one phylogenetic hypothesis based on PRT datasets and using the BI method was proposed here (Figure 3); the rest are shown in the Supplementary Materials (Figures S10-S14). The resulting tree topologies are all congruent at the subfamily level, and nodal support values vary slightly for different analyses ( Figure S15). The phylogenetic relationship among subfamilies is (Coeliadinae + (Euschemoninae + (Eudaminae + Pyrginae sensu lato) + (Heteropterinae + (Barcinae + Hesperiinae)))).
The monophyly of the subfamily Pyrginae is strongly supported (nodal support value is 0.867) as is the monophyly of Eudaminae + Pyrginae (nodal support value is 0.914), which means that these two subfamilies are not necessarily divided if based on this phylogeny. There are some differences in the phylogenetic topological structures obtained based on different datasets and different methods in this study, which are mainly reflected in the positional relationships of the subfamily Eudaminae, and the tribe Tagiadini and the phylogenetic relationships among the tribes Erynnini, Achlyodidini, Carcharodini, and Pyrgini. Since the BI tree based on PRT datasets has the high nodal support, only BI trees of PRT datasets are shown in this paper.
All the nine samples are placed within the subfamily Pyrginae sensu lato (Figures 3 and 4). As in the traditional taxonomy of the group, Abraximorpha, Gerosis, Mooreana, and Tagiades are placed in Tagiadini, Celaenorrhinus in Celaenorrhini, and Erynnis in Erynnini.

Phylogenetic Relationships
Only one phylogenetic hypothesis based on PRT datasets and using the BI method was proposed here (Figure 3); the rest are shown in the Supplementary Materials ( Figures  S10-S14). The resulting tree topologies are all congruent at the subfamily level, and nodal support values vary slightly for different analyses ( Figure S15). The phylogenetic relationship among subfamilies is (Coeliadinae + (Euschemoninae + (Eudaminae + Pyrginae sensu lato) + (Heteropterinae + (Barcinae + Hesperiinae)))). The monophyly of the subfamily Pyrginae is strongly supported (nodal support value is 0.867) as is the monophyly of Eudaminae + Pyrginae (nodal support value is 0.914), which means that these two subfamilies are not necessarily divided if based on this phylogeny. There are some differences in the phylogenetic topological structures obtained based on different datasets and different methods in this study, which are mainly reflected In the tribe Tagiadini, Tagiades menaka forms a monophyletic clade with T. vajuna. Our results support the view that the genus Daimio Murray, 1587, is a synonym of the genus Tagiades Hübner, 1819 [54]. The genus Abraximorpha is a sister group of the genus Odontoptilum. Gerosis phisara forms a monophyletic clade with G. bhagava. The genus Ctenoptilum and the genus Gerosis are sister groups with high nodal support and a stable relationship. Abraximorpha, Odontoptilum, and Ctenoptilum share a significant character, asymmetric genitalia. Based on our phylogeny, the character was lost secondarily in Gerosis. Another character also lost in Gerosis is a hair tuft at the tip of the female abdomen. The genus Satarupa forms a monophyletic clade with the genus Mooreana. In the tribe Celaenorrhini, Celaenorrhinus aspersus forms a monophyletic clade with C. maculosus and C. syllius with strong support. In the tribe Erynnini, Erynnis popovina forms a monophyletic clade with E. montanus and E. brizo brizo.
This study re-confirmed the status of Capila, Pseudocoladenia, and Sarangesa. Capila belongs to the tribe Tagiadini, and Pseudocoladenia and Sarangesa belong to the tribe Celaenorrhinini [5].
in the positional relationships of the subfamily Eudaminae, and the tribe Tagiadini and the phylogenetic relationships among the tribes Erynnini, Achlyodidini, Carcharodini, and Pyrgini. Since the BI tree based on PRT datasets has the high nodal support, only BI trees of PRT datasets are shown in this paper.
All the nine samples are placed within the subfamily Pyrginae sensu lato (Figures 3  and 4). As in the traditional taxonomy of the group, Abraximorpha, Gerosis, Mooreana, and Tagiades are placed in Tagiadini, Celaenorrhinus in Celaenorrhini, and Erynnis in Erynnini. In the tribe Tagiadini, Tagiades menaka forms a monophyletic clade with T. vajuna. Our results support the view that the genus Daimio Murray, 1587, is a synonym of the genus Tagiades Hübner, 1819 [54]. The genus Abraximorpha is a sister group of the genus Odontoptilum. Gerosis phisara forms a monophyletic clade with G. bhagava. The genus Ctenoptilum and the genus Gerosis are sister groups with high nodal support and a stable relationship. Abraximorpha, Odontoptilum, and Ctenoptilum share a significant character, asymmetric genitalia. Based on our phylogeny, the character was lost secondarily in Gerosis. Another character also lost in Gerosis is a hair tuft at the tip of the female abdomen. The genus Satarupa forms a monophyletic clade with the genus Mooreana. In the tribe Celaenorrhini, Celaenorrhinus aspersus forms a monophyletic clade with C. maculosus and C. syllius with strong support. In the tribe Erynnini, Erynnis popovina forms a monophyletic clade with E. montanus and E. brizo brizo.
This study re-confirmed the status of Capila, Pseudocoladenia, and Sarangesa. Capila belongs to the tribe Tagiadini, and Pseudocoladenia and Sarangesa belong to the tribe Celaenorrhinini [5].
Capila translucida forms a monophyletic clade with C. zennara in the tribe Tagiadini. In Celaenorrhini, Sarangesa forms a monophyletic clade with Eretis. The intergeneric relationship within the tribe Celaenorrhinini is (Celaenorrhinus + ((Sarangesa + Eretis) + Pseudocoladenia)). Evans (1949) adopted the length of the forewing cell and labial palpi to distinguish his Celaenorrohinus-group and Tagiades-group. Then, Evans (1952) more precisely stated Evans (1949) adopted the length of the forewing cell and labial palpi to distinguish his Celaenorrohinus-group and Tagiades-group. Then, Evans (1952) more precisely stated that the forewing cell is equal to two-thirds of the length of the costa and equal to or longer than the dorsum in Pyrginae Section I, including current Celaenorrhini. In contrast, the forewing cell of Celaenorrhinus maculosus (Figure 6c) is equal to the dorsum but much shorter than two-thirds of the costa. On the other hand, in Tagiades vajuna (Figure 6a), the forewing cell is shorter than two-thirds of the length of the costa but equal to the dorsum. In Capila omeia (Figure 6b), which is shifted from Celaenorrhini to Tagiadini, the forewing cell is not equal to or longer than two-thirds the length of the costa or dorsum. This is also the case in Sarangesa dasahara (Figure 6d) and Pseudocoladenia dan (Figure 6e), both of which are shifted from Tagiadini to Celaenorrhinni. To summarize, the length of the forewing cell cannot be adopted as a diagnostic character to distinguish the tribes Celaenorrhini and Tagiadini.
The other diagnostic character that Evans (1949Evans ( , 1952 used is the labial palpi. In Section I in America or the Celaenorrhinus-group in Asia, the labial palpi are erect. In other words, the second segment is appressed to the face, and the third segment is not protruding in front of the second segment. However, in Section II in America or the Tagiades-group in Asia, the labial palpi aren't erect, and the third segment is protruding. The character is obviously applicable for Celaenorrhinus ( Figure 5c) and Tagiades (Figure 5a), respectively. In Capila (Figure 5b), shifted from Celaenorrhini to Tagiadini, the third segment protrudes in front of the second segment, but this is also the case in Sarangesa (Figure 5d) and Pseudocolodenia (Figure 5e), shifted the other way. Therefore, in conclusion, the character cannot be adopted for either. Thus, we were not successful in finding a good diagnostic character for distinguishing between these two tribes. Evans (1949) classified Coladenia (now Pseudocoladenia) dan into 11 subspecies in four groups. Subsequently, because of their sympatric distribution and the differences in male and female genitalia, Huang and Xue (2004) raised the taxonomic status of two taxa from China to distinct species. In order to clarify the species limitation of some taxa, we calculated p-distance between our sample, P. dan fabia from Yunnan Province, China, and that from Myanmar (GenBank #SRR7174480) and P. dan dhyana from southern Thailand (#KY019868). The distance between P. dan fabia from the two different localities was 0.022, whereas the distance between P. dan fabia from China and P. dan dhyana was 0.027. Though the latter number reaches around the border of species and subspecies [55], we retained the taxonomic status, pending inclusion of the nominate subspecies dan from South India in our analysis.     Chou (1998) [56].

Supplementary Materials:
The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/insects13010068/s1, Figure S1: Predicted secondary cloverleaf structure for the tRNAs of A. davidi. Lines (-) indicate Watson-Crick base pairings, whereas dots (·) indicate unmatched base pairings. Figure S2: Predicted secondary cloverleaf structure for the tRNAs of Ca. translucida. Lines (-) indicate Watson-Crick base pairings, whereas dots (·) indicate unmatched base pairings. Figure S3: Predicted secondary cloverleaf structure for the tRNAs of Ce. aspersus. Lines (-) indicate Watson-Crick base pairings, whereas dots (·) indicate unmatched base pairings. Figure S4: Predicted secondary cloverleaf structure for the tRNAs of E. popoviana. Lines (-) indicate Watson-Crick base pairings, whereas dots (·) indicate unmatched base pairings. Figure S5: Predicted secondary cloverleaf structure for the tRNAs of G. phisara. Lines (-) indicate Watson-Crick base pairings, whereas dots (·) indicate unmatched base pairings. Figure S6: Predicted secondary cloverleaf structure for the tRNAs of M. trichoneura. Lines (-) indicate Watson-Crick base pairings, whereas dots (·) indicate unmatched base pairings. Figure S7: Predicted secondary cloverleaf structure for the tRNAs of P. fabia. Lines (-) indicate Watson-Crick base pairings, whereas dots (·) indicate unmatched base pairings. Figure S8: Predicted secondary cloverleaf structure for the tRNAs of Sat. nymphalis. Lines (-) indicate Watson-Crick base pairings, whereas dots (·) indicate unmatched base pairings. Figure S9: Predicted secondary cloverleaf structure for the tRNAs of T. menaka. Lines (-) indicate Watson-Crick base pairings, whereas dots (·) indicate unmatched base pairings. Figure S10: Phylogenetic tree produced by Bayesian inference analysis of the PCG dataset. Bayesian posterior probability (BPP) support values are indicated above the branches. Figure S11: Phylogenetic tree produced by maximum likelihood analyses of PCG dataset. Bootstrap support values (BS) are indicated above the branches. Figure S12: Phylogenetic tree produced by maximum likelihood analysis of the PRT dataset. Bootstrap support values (BS) are indicated above the branches. Figure S13: Phylogenetic tree produced by Bayesian inference analysis of the PCG12RT dataset. Bayesian posterior probability (BPP) support values are indicated above the branches. Figure S14: Phylogenetic tree produced by maximum likelihood analyses of PCG12RT dataset. Bootstrap support values (BS) are indicated above the branches. Figure S15: Phylogenetic tree produced by Bayesian inference analysis of the PRT dataset. Bayesian posterior probability (BPP) support values are indicated above the branches. Table S1: The basic statistics of sequencing for nine mitochondrial genomes. Table S2: The best partitioning schemes and models for Bayesian inference (BI) method based on the three datasets selected by PartitionFinder. Table S3: The best partitioning schemes and models for Maximun likelihood (ML) method based on three datasets selected by PartitionFinder.