Novel Gene Rearrangement Pattern in Pachycrepoideus vindemmiae Mitochondrial Genome: New Gene Order in Pteromalidae (Hymenoptera: Chalcidoidea)

Simple Summary The mitochondrial genome is a reliable genetic marker for reconstructing phylogeny and Pteromalidae is a diverse and complex family of chalcid wasps, but its evolutionary history is still poorly understood. In this study, we sequenced the mitochondrial genomes of four species (Muscidifurax similadanacus, M. sinesensilla, Nasonia vitripennis, and Pachycrepoideus vindemmiae) of Pteromalidae. Additionally, a phylogenetic hypothesis was reconstructed for the subfamilies of Pteromalidae that includes newly acquired mitogenomes and those deposited in NCBI. We used pairwise breakpoint distances to infer this phylogeny. Our study enriches the overall knowledge on gene rearrangement in Pteromalidae, reveals the evolutionary relationships among several major groups of Pteromalidae, accumulates molecular data for a Pteromalidae phylogeny, and provides a genetic background basis for biological control in agriculture and forestry. Abstract The mitochondrial genomes of Muscidifurax similadanacus, M. sinesensilla, Nasonia vitripennis, and Pachycrepoideus vindemmiae were sequenced to better understand the structural evolution of Pteromalidae mitogenomes. These newly sequenced mitogenomes all contained 37 genes. Nucleotide composition was AT-biased and the majority of the protein-coding genes exhibited a negative AT skew. All 13 protein-coding genes (PCGs) initiated with the standard start codon of ATN, excepted for nad1 of N. vitripennis, which started with TTG, and terminated with a typical stop codon TAA/TAG or an incomplete stop codon T. All transfer RNA (tRNA) genes were predicted to fold into the typical clover-leaf secondary structures, except for trnS1, which lacks the DHU arm in all species. In P. vindemmiae, trnR and trnQ lack the DHU arm and TΨC arm, respectively. Although most genes evolved under a strong purifying selection, the Ka/Ks value of the atp8 gene of P. vindemmiae was greater than 1, indicating putative positive selection. A novel transposition of trnR in P. vindemmiae was revealed, which was the first of this kind to be reported in Pteromalidae. Two kinds of datasets (PCG12 and AA) and two inference methods (maximum likelihood and Bayesian inference) were used to reconstruct a phylogenetic hypothesis for the newly sequenced mitogenomes of Pteromalidae and those deposited in GenBank. The topologies obtained recovered the monophyly of the three subfamilies included. Pachyneurinae and Pteromalinae were recovered as sister families, and both appeared sister to Sycophaginae. The pairwise breakpoint distances of mitogenome rearrangements were estimated to infer phylogeny among pteromalid species. The topology obtained was not totally congruent with those reconstructed using the ML and BI methods.

The phylogeny of Pteromalidae has long been debated and is still controversial. The relationships among Pteromalidae and other chalcid families are still uncertain. Chen et al. considered that Agaonidae was closely related to Pteromalidae [15], whereas Wu et al. proposed that Pteromalidae and Eurytomidae were sister groups [16]. Pteromalidae has been considered as a polyphyletic group for decades [17]. This concept was supported by several studies [18][19][20], even though Pteromalidae was recovered as monophyletic in several other studies [16,21,22]. Rasplus et al. [23] suggested that Sycoecinae, Otitesellinae, and Sycoryctinae (used to be Agaonidae) should be transferred to Pteromalidae based on ribosomal and mitochondrial genes as well as morphology. Dzhanokmen [24] combined morphology and biological characters to reconstruct an unformal phylogeny and suggested Pteromalinae to be sister to all other Pteromalidae. Desjardins et al. [25] used four nuclear coding genes to reconstruct the phylogenetic relationship of Pteromalidae. Their results confirmed the monophyly of nearly all tribes represented by multiple specimens, and two subfamilies. Other subfamilies appeared para-or polyphyletic, and monophyly was significantly rejected for Miscogasterinae, Ormocerinae, and Colotrechninae. Munro et al. [18] conducted a molecular analysis based on 18S and 28S ribosomal genes, which supported the monophyly of seven subfamilies and the paraphyly or polyphyly of nine subfamilies. A similar result was recovered by Heraty et al. [19]. Pteromalidae, as historically understood, is both a highly diverse and a highly disparate group, these properties rendering it difficult to reconstruct a reliable phylogenetic hypothesis based on both morphological and molecular characters. Hereafter, we propose to investigate the relationships among Pteromalidae based on complete mitochondrial genomes.
In most insect groups, gene rearrangements in mitogenomes are relatively rare events and can be reliably used as phylogenetic markers [3]. To the opposite, gene rearrangements appear common in Chalcidoidea mitogenomes and may represent evolutionarily independent events. Gene rearrangements in Chalcidoidea have been hypothesized to correlate with the diversity of their life-history traits and with their phylogenetic positions [14]. In Pteromalidae, twenty nearly complete mitochondrial genomes belonging to four subfamilies have been reported so far, and a total of 14 gene rearrangements have been found in these subfamilies [14,16,26]. Oliveira et al. [26] first sequenced a partial mitochondrial genome of N. vitripennis. Compared with a hypothetical ancestral pattern, they identified six inverted protein-coding genes (PCGs). The incomplete mitochondrial genome of Nasonia vitripennis Walker, 1836 (Pteromalidae) was the first mitochondrial genome to be reported in Chalcidoidea; it showed novel gene rearrangements compared to other insects and was considered the typical model for Chalcidoidea [26]. However, nine genes were not sequenced and the gene arrangements of N. vitripennis were still partly undetermined. Subsequently, Xiao et al. [14] sequenced the mitochondrial genomes of two species of Sycoryctinae and showed that gene rearrangements mostly concerned tRNAs. Finally, Tang et al. and Wu et al. [16,27] performed mitochondrial sequencing of Pteromalidae species and summarized gene rearrangements.
The CREx program online is widely used to heuristically determine pairwise rearrangement events in mitochondrial genomes [28]. It considers transpositions, reverse transpositions, reversals, and tandem-duplication-random-loss (TDRL) events, which are based on common intervals that reflect genes that appear consecutively in several of the input gene orders. Wu et al. [16] proposed that "closely related groups in lower taxonomic levels tend to exhibit more similar gene orders, which indicates that gene rearrangements may still be useful for phylogenetic analysis at higher taxonomic levels. Pairwise breakpoint distances can be used to analyze the rates of mitochondrial gene rearrangement for Chalcidoidea using the CREx web server".
We measured complete or nearly complete mitochondrial genomes of four Pteromalidae, Muscidifurax similadanacus Xiao and Zhou, 2018 [29], M. sinesensilla Xiao & Zhou, 2018 [29], and N. vitripennis as primary parasitoids and Pachycrepoideus vindemmiae Rondani, 1875 which is both a primary parasitoid and hyperparasitoid [29][30][31][32][33][34][35][36]. Using these newly acquired mitogenomes of Pteromalidae as well as those deposited in GenBank, we propose a phylogenetic hypothesis based on mitogenomes for Pteromalidae. A second aim of this study is to analyze the characteristics of the Pteromalidae mitogenomes and to test the ability of pairwise breakpoint distances to infer reliable phylogenetic relationships.

Sample Preparation and DNA Extraction
All specimens were collected from the wild, including the three primary parasitoids and Pachycrepoideus vindemmiae as primary parasitoids and hyperparasitoids [37]. M. similadanacus and M. sinesensilla were collected in Xinjiang province of China, N. vitripennis was collected in Shandong province, and P. vindemmiae in Anhui province (Table 1), and then they were reared in the laboratory for more than three years. All the specimens were stored at −20 • C in absolute ethanol prior to DNA extraction. Total genomic DNA was extracted using the cetyltrimethyl ammonium bromide (CTAB) method [38].

Sequencing and Assembly
Sequencing was performed using a whole-genome shotgun (WGS) strategy on the Illumina Novaseq platform. The PCR library was constructed with the use of a TruSeqTM DNA Sample Prep Kit for each species using genomic DNA with an insert size of 350 bp and was sequenced on the Illumina platform at Berry Genomics, Beijing; a total of 10 Gb clean data (150 bp pair-end reads) were obtained for each species. The fragments were repaired by the combined action of 3 -5 exonuclease and polymerase in an End Repair Mix. Thereafter, magnetic beads purified the connection products, which were incubated with the DNA fragment to remove free and self-connecting joint sequences. Finally, the homogenized and mixed libraries were gradually diluted and quantified to 4-5 pM for sequencing. The quality of data was checked with FastQC [6,41]. The original data adapter had been removed by AdapterRemoval version 2 [42]. SOAPec version 2.01 was used for quality correction with K-mer set to 17. Reads with a length of less than 50 bp were excluded. Geneious v 11.0.2 was used to assemble and annotate the mitochondrial (mt) genomes [6,43]. The COX1 gene of N. vitripennis (MT635402) was used as a reference to map and to identify the COX1 gene of these four pteromalid species in our mitogenome sequences [44]. Subsequently, the sequence was extended on each side using a de novo method based on the local alignment method.

Mitochondrial Genome Annotation
The tRNA genes and protein-coding genes were identified using the MITOS Webserver. The secondary structure was also predicted using the MITOS WebSever, setting the parameters with the Invertebrate Mito genetic code [6,45]. Protein-coding genes (PCGs) were identified as open reading frames corresponding to the 13 PCGs in the metazoan mitogenomes. Mitogenome maps were produced using Organellar Genome DRAW (OG-DRAW) [46].

Comparative Analysis
Base composition and relative synonymous codon usage (RSCU) were analyzed using MEGA X [47]. Geneious v 11.0.2 [43] was used to check all the genes in the mitochondrial genome. Comparative analyses of codon usage for these four mitogenomes were calculated using PhyloSuite [22,45]. The predicted secondary structures of all tRNAs were drawn using Adobe Illustrator CC 2018. Multiple-substitution correction of Non-synonymous (Ka)/synonymous (Ks) mutation rate ratios among the 13 PCGs was calculated with DnaSP v5, using that of Pachyneuron aphidis (Pteromalidae: Pachyneurinae) as a reference sequence [48,49]. AT/GC skewness was calculated as AT skew= (A − T)/(A + T) and GC skew = (G − C)/(G + C) [50].

Phylogenetic Analysis
To reconstruct the phylogeny of Pteromalidae, 4 newly generated mitogenomes and 18 others downloaded from GenBank were analyzed. Species of Eupelmus (Eupelmidae) and Torymus (Torymidae) were used as outgroups ( Table 1). All protein-coding genes were aligned individually based on codon optimized multiple alignments using the MAFFT 7.3.1 and G-INS-I algorithm [51]. Aligned sequences were then concatenated into two different datasets (PCG12, including the first and second codon of 13 protein-coding genes, and AA: 13 PCGs translated into amino acids). Maximum-likelihood (ML) analysis was conducted in W-IQtree [52] using the best-fit substitution model. An ultrafast bootstrap (UFB) [53] of 1000 replications and the SH-aLRT test [54] was used to assess branch supports. Bayesian inference was implemented in PhyloBayes MPI 1.5a [55] and we used a site-heterogeneous mixture model (CAT+GTR). The trees were sampled every 1000 generations, and the burning rate was set as 0.25 of the sampled value. FigTree v.1.3.1 [56] was used to view topologies.

Ancestral Character Reconstruction and Gene Rearrangement in Pteromalidae
To explore the effect of parasitic lifestyles (primary parasitoids or/and hyperparasitoids) and of the phylogenetic position in gene rearrangement, we performed an ancestral character reconstruction using Mesquite 3.04 [57]. Because some mitogenome sequences were incomplete and pairwise breakpoint distances could not be calculated, we used mitogenomes with 37 genes (13 PCGs, 2 rRNAs, 22 tRNAs) and the complete/incomplete control region to performed our test.
Four newly generated mitogenomes and three obtained from GenBank were analyzed to reconstruct the ancestral character of Pteromalidae. The parasitic lifestyle of the species included in our dataset was coded as follows: 1: primary parasitoids; 2: both primary parasitoids and hyperparasitoids; and 3: undefined. The results obtained were organized in figures using Adobe Illustrator CC 2018. A map of the mitochondrial gene rearrangement was depicted with the Illustrator of biological sequences (IBS) [58]. Pairwise breakpoint distances (PBD) were used to analyze the rate of mitochondrial gene rearrangement for Pteromalidae (complete mitogenomes only) using the CREx web server [28]. A heatmap of the pairwise breakpoint distance was drawn using TBtools [59].

Characteristics of Base Composition
The nucleotide composition of the mitogenomes exhibits a high A + T content, with an average of 83.98%, which is usual in mtDNA in arthropods and Hymenoptera [60][61][62]. Among the newly generated mitogenomes, P. vindemmiae had the highest AT content (85.3%) and N. vitripennis had the lowest (83.0%) (Tables 1 and 2). Based on the analysis of base content, the range of the AT (GC) skew value varies from 0.018 (−0.112) to 0.047

Characteristics of Base Composition
The nucleotide composition of the mitogenomes exhibits a high A + T content, with an average of 83.98%, which is usual in mtDNA in arthropods and Hymenoptera [60][61][62]. Among the newly generated mitogenomes, P. vindemmiae had the highest AT content (85.3%) and N. vitripennis had the lowest (83.0%) (Tables 1 and 2). Based on the analysis of base content, the range of the AT (GC) skew value varies from 0.018 (−0.112) to 0.047 (−0.091); for P. vindemmiae it was 0.036 (−0.108), indicating a high A-T bias in Pteromalidae mitogenomes.
Among the 13 PCGs, the highest AT content was 84.2% for P. vindemmiae, and the lowest was 81.7% for N. vitripennis. Among the newly generated mitogenomes, P. vindemmiae had the highest AT contents for the full genome, PCGs, tRNAs, and rRNAs, which were 85.3%, 84.2%, 88.9%, and 88.0%, respectively (Table 4). There are three hydrogen bonds between G and C and two hydrogen bonds between A and T, so GC was more stable than AT. High AT content is likely to be one of the most important factors explaining rearrangement in mitochondrial genomes.

Overlap and Gap
A total length of 370 bp of gaps were found in the newly sequenced mitogenomes, of which M. similadanacus had 58 bp of gaps, M. sinesensilla 121 bp, N. vitripennis 94 bp, and P. vindemmiae 97 bp. Among these species, the longest gap was 35 bp and was located between trnI and nad2 of P. vindemmiae. The shortest gap was 1 bp and was located between trnQ and nad1 of M. similadanacus, trnN and trnR and trnH and nad4 of M. sinesensilla, trnD and atp8 and cytb and nad6 with N. vitripennis, trnS1 and trnY and trnC and trnN, cox2, and trnK, trnQ, and nad1 of P. vindemmiae.
Overlapping genes are common in mitochondrial genomes of arthropods. A total of 98 bp of overlaps were observed in our sequenced mitogenomes. Overlapping genes might reflect the selection for a short and economic mitogenome, and they usually involve trn genes, because their sequences are constrained by fewer mutations [63,64]. Most of the discovered overlaps appeared in tRNA. M. similadanacus had 31 bp of overlaps, M. sinesensilla 23 bp, N. vitripennis 21 bp, and P. vindemmiae 23 bp. The length of overlap ranged from 1 to 10 bp, and the longest was located between atp8 and atp6 of M. similadanacus and M. sinesensilla. The shortest was 1 bp and were located between atp6 and cox3, cytb and nad6, and trnS2 and cytb of M. similadanacus; nad3 and trnC and atp6 and cox3 of M. sinesensilla; nad3 and trnC and nad5 and trnF of N. vitripennis; and atp6 and cox3, trnT and nad4l, and cytb and nad6 of P. vindemmiae.

Transfer RNA and Ribosomal RNA Genes
In total, 22 tRNA genes were interspersed throughout the mitochondrial genomes of Pteromalidae. Within the tRNA secondary structures of pteromalid mitogenomes, the DHU arm of trnS1 was missing, while in P. vindemmiae, the DHU arm and TΨC arm were missing in trnR and trnQ, respectively. Many insects lack DHU arms in the trnS1 secondary structures, while other secondary structures are typical clover structures. This folding probably happened early in the evolution of Metazoa [65]. The TΨC arm of some insects is also missing [66]. Among the 22 tRNAs, 16 tRNAs were in the J chain and the six others were in the N chain. Additionally, the size of tRNAs ranged from 55 to 74 bp (Figure 2). The average nucleotide composition of these tRNAs was A: 44.5%, T: 43.7%, C: 6.9%, and G: 4.9%, with a total average A + T content of 88.2%. Most AT skews were positive, except in P. vindemmiae (−0.006). In the same manner, all GC skews were negative, which indicates a slight bias towards the use of A and T in tRNAs.
The positions of large and small ribosomal RNA genes (rrnL and rrnS) were consistent with those observed in most other insects (between trnL1-trnA and trnA-trnV). The lengths of rrnS and rrnL for N. vitripennis were the longest among the four species, and the lengths of rrnS and of rrnL for P. vindemmiae were the shortest. Among these four species, the length of rrnS ranged from 763 to 787 bp, and the average AT content was 87.6% (P. vindemmiae had the highest AT content); the length of rrnL ranged from 1300 to 1325 bp, and the average AT content was 87.0%. In P. vindemmiae, the length of the rrnS gene was 763 bp with an AT content of 88.9%, while the rrnL gene was 1300 bp with an AT content of 87.5%.

Protein-Coding Genes
In our four new mitogenomes, 10 of 13 protein-coding genes were encoded by the majority strand, while 3 genes (nad2, nad6 and cytb) were encoded by the minority strand; the length of the 13 protein-coding genes ranged from 161 bp to 1675 bp. Among them, the length of the 13 protein-coding genes of P. vindemmiae was 159-1672 bp. Among the 13 protein-coding genes, the shortest and longest encoding genes were atp8 and nad5. In M. similadanacus and M. sinesensilla, the protein-encoding genes had the same length excepted for nad1, nad5, and nad6. The length of cox3, atp6, and nad4l were the same in the four newly sequenced species, being, respectively, 786 bp, 675 bp, and 288 bp.
The total length of the PCGs of the four new mitogenomes ranged from 11092 bp to 11117 bp. The shortest was 11092 bp in P. vindemmiae and accounted for 74.69% of the entire genome. The highest average AT content of the 13 protein-coding genes was 84.20% in P. vindemmiae and the lowest was 81.70% in N. vitripennis. In all new mitogenomes, cox1 had the lowest AT content in the 13 PCGs, ranging from 76.9% (P. vindemmiae) to 74.5% (N. vitripennis), and atp8 had the highest AT content in these 13 PCGs, ranging from 92.50% (P. vindemmiae) to 88.10% (N. vitripennis).
The preferred initiation codon of the pteromalid mitogenomes was ATN, as observed in most other insect mitogenomes [47], except for nad1 of N. vitripennis, which starts with TTG. In these four pteromalid species, the start codon of six genes (cox3, atp6, cox1, nad4, nad6, and cytb) was ATG and it was ATT for the other three genes (cox2, nad5, nad4l). The stop codon of eight genes (nad2, cox3, atp6, atp8, cox2, nad4l, nad6, cytb) was TAA. TAG was the stop codon of nad5 of N. vitripennis. Four genes of P. vindemmiae (nad3, cox1, nad4, and nad5), two genes of N. vitripennis (cox1 and nad4), two genes of M. sinesensilla (nad5 and nad1), and one gene of M. similadanacus (nad5) use the incomplete T-as stop codons (Tables 2 and 3).  Codons with high AT content were preferred, which was consistent with most other insect mitogenomes [67]. The relative synonymous codon usages (RSCU) of these four species are shown in Figure 3. Taken together, the most frequently used codons were UUA (Leu2), CGA (Arg), GGA (Gly), UCA (Ser2), CCU (Pro), GUU (Val), and GCU (Ala), whereas codons ending with G or C, CUG, CUC, CAG, and GGC, were the less frequently used codons. The third codon position of A/T occurred more frequently than that of G/C, reflecting AT nucleotide bias in the mitochondrial PCGs among Pteromalidae.

Evolutionary Rate Analysis
The non-synonymous/synonymous substitution ratio (Ka/Ks) can be used to estimate whether a sequence is undergoing purifying, neutral, or positive selection. Among the four new mitogenomes, the Ka/Ks value of cox1 was the lowest, at 0.12, 0.13, 0.14, and 0.14, respectively; among M. sinesensilla and P. vindemmiae, the Ka/Ks value of atp8 was the highest at 0.97 and 1.59 respectively, but the highest Ka/Ks value of M. similadanacus was that of nad2 (0.76), and the highest of N. vitripennis was that of nad4l (0.96).
In total, only the Ka/Ks value of atp8 of P. vindemmiae was greater than one. A Ka/Ks greater than one indicates that this protein-coding gene has been subjected to positive selection effects during evolution. This may suggest that non-synonymous mutations of atp8 were selectively retained, which is rarely observed in insects. Following non-synonymous mutations, the mutated protein-coding gene was retained. Additionally, most of the 13 protein-coding genes were under purifying selection (Figure 4).

Evolutionary Rate Analysis
The non-synonymous/synonymous substitution ratio (Ka/Ks) can be used to estimate whether a sequence is undergoing purifying, neutral, or positive selection. Among the four new mitogenomes, the Ka/Ks value of cox1 was the lowest, at 0.12, 0.13, 0.14, and 0.14, respectively; among M. sinesensilla and P. vindemmiae, the Ka/Ks value of atp8 was the highest at 0.97 and 1.59 respectively, but the highest Ka/Ks value of M. similadanacus was that of nad2 (0.76), and the highest of N. vitripennis was that of nad4l (0.96).
In total, only the Ka/Ks value of atp8 of P. vindemmiae was greater than one. A Ka/Ks greater than one indicates that this protein-coding gene has been subjected to positive selection effects during evolution. This may suggest that non-synonymous mutations of atp8 were selectively retained, which is rarely observed in insects. Following nonsynonymous mutations, the mutated protein-coding gene was retained. Additionally, most of the 13 protein-coding genes were under purifying selection (Figure 4).

Phylogenetic Analysis
Two datasets (PCG12 and AA) and two inference methods, namely maximum-likelihood (ML) and Bayesian inference (BI), were used to reconstruct the phylogeny of pteromalid mitogenomes, and four topologies of ML-AA, ML-PCG12, BI-AA, and BI-PCG12 were reconstructed. All topologies were highly supported and both methods yield similar topologies. All pteromalid subfamilies were recovered as monophyletic. In all topologies, Pachyneurinae was recovered as sister to Pteromalinae, and this clade was then sister to Sycophaginae. Pteromalinae were subdivided into two monophyletic tribes, Otitesellini and Pteromalini ( Figure 5).

Phylogenetic Analysis
Two datasets (PCG12 and AA) and two inference methods, namely maximum-likelihood (ML) and Bayesian inference (BI), were used to reconstruct the phylogeny of pteromalid mitogenomes, and four topologies of ML-AA, ML-PCG12, BI-AA, and BI-PCG12 were reconstructed. All topologies were highly supported and both methods yield similar topologies. All pteromalid subfamilies were recovered as monophyletic. In all topologies, Pachyneurinae was recovered as sister to Pteromalinae, and this clade was then sister to Sycophaginae. Pteromalinae were subdivided into two monophyletic tribes, Otitesellini and Pteromalini ( Figure 5).

Ancestral Character Reconstruction and Gene Rearrangement
Based on the BI topology, the ancestral biology of Pteromalidae included in the analysis was inferred to be both primary parasitoid and hyperparasitoid ( Figure 6).

Ancestral Character Reconstruction and Gene Rearrangement
Based on the BI topology, the ancestral biology of Pteromalidae included in the analysis was inferred to be both primary parasitoid and hyperparasitoid ( Figure 6).
Gene rearrangement patterns in Hymenoptera are usually more complex and variable than in other insect orders. Compared to the putative ancestral pattern [69,70] of the insect mitogenome, dramatic gene rearrangements were observed in Pteromalidae mitogenomes. These rearrangements occurred in tRNA genes but also in protein-coding genes. Five gene blocks were found in Pteromalidae (only considering complete mitogenomes). The first gene block was cox3-atp6-atp8, which existed in all sequenced species. With the exception of Philotrypesis tridentata, all other species shared the gene block cox2-trnL2-cox1. All species also shared the gene block of nad5-trnH-nad4-nad4l-trnT-trnP-nad6-cytb, but P. Gene rearrangement patterns in Hymenoptera are usually more complex and variable than in other insect orders. Compared to the putative ancestral pattern [69,70] of the insect mitogenome, dramatic gene rearrangements were observed in Pteromalidae mitogenomes. These rearrangements occurred in tRNA genes but also in protein-coding genes. Five gene blocks were found in Pteromalidae (only considering complete mitogenomes). The first gene block was cox3-atp6-atp8, which existed in all sequenced species. With the exception of Philotrypesis tridentata, all other species shared the gene block cox2-trnL2-cox1. All species also shared the gene block of nad5-trnH-nad4-nad4l-trnT-trnP-nad6-cytb, but P. tridentata showed an inversion segment (cytb-nad6-trnP-trnT-nad4l-nad4-trnH-nad5). Finally, most species exhibited a gene block of trnS2-trnQ-nad1-trnL1-rrnL-trnA. In addition, gene block trnI-nad2-trnW-trnY-trnS1 was shared by all sequenced genomes.

Figure 6.
Left side is the ancestral state reconstruction based on the BI tree. Red pies represent primary parasitoids, half-green and half-red pies represent primary parasitoids and hyperparasitoids, and white pies represent undefined biology. Gene orders were arranged on the right side of pteromalids. Different gene rearrangements are highlighted in different colors. The heat map of the pairwise breakpoint distance of Pteromalidae is on the right side.

Discussion
Gene rearrangement information is highly valuable for phylogenetic reconstruction of specific lineages [8,64], especially in the classification of low-level elements [27,74]. Gene rearrangement has been reported in most insect orders and an increased rate of gene rearrangements has been observed in Hemiptera [10,73], Protura [75], and Hymenoptera [4,14,62]. Since the pioneering work of Boore et al. [2], it has been acknowledged that the gene structure of mitochondrial genomes contains phylogenetic signal [3]. Detection of gene rearrangements in lower-level elements of insect mitogenomes is expected to shed light on the evolution of these taxa [62,74]. The mitochondria of Hymenoptera generally exhibit extremely high molecular evolution rates and extensive gene rearrangements due to their parasitic lifestyles. However, our understanding of the mitochondrial genome of pteromalids is still limited, and the mechanism of gene rearrangement is still unclear.
The site-heterogeneous mixture model (CAT+GTR) implemented in PhyloBayes has been shown to correct the wrong grouping of unrelated taxa that share a similar base composition and an accelerated evolutionary rate [77]. ML and BI analyses yield similar Figure 6. Left side is the ancestral state reconstruction based on the BI tree. Red pies represent primary parasitoids, half-green and half-red pies represent primary parasitoids and hyperparasitoids, and white pies represent undefined biology. Gene orders were arranged on the right side of pteromalids. Different gene rearrangements are highlighted in different colors. The heat map of the pairwise breakpoint distance of Pteromalidae is on the right side.
Comparison with other pteromalid mitogenomes revealed a novel inversion in the trnE-trnF gene cluster of Pachyneurinae and a novel transposition of trnR in P. vindemmiae, which was the first of this kind to be reported in Pteromalidae ( Figure 6). The accelerated rate of gene rearrangement may be the result of the fast evolution of this group [71][72][73].
In Pteromalidae, most species have a high frequency of mitochondrial gene rearrangements, and the breakpoint distances range from 0 to 21 ( Figure 6). According to our phylogeny, M. similadanacus and M. sinesensilla were the most closely related species, and the pairwise breakpoint distance between them was 0, which is in agreement with the assumption that a lower value of the pairwise breakpoint distance indicates close relationship in the topology ( Figure 6). P. vindemmiae was sister to P. aphidis and the pairwise breakpoint distance between them was 10. P. vindemmiae and M. similadanacus were more distantly related in the topology; however, the pairwise breakpoint distance observed between them was only 5, which is a lower value than that observed between P. vindemmiae and P. aphidis. N. vitripennis and A. calandrae were closely related on the phylogenetic tree, and the pairwise breakpoint distance value was 7, whereas the value of N. vitripennis and P. vindemmiae was only 5. Globally, the pairwise breakpoint distance of mitogenome rearrangements appeared not to be consistent with the relationships and the proximity between Pteromalidae species observed in our topology.

Discussion
Gene rearrangement information is highly valuable for phylogenetic reconstruction of specific lineages [8,64], especially in the classification of low-level elements [27,74]. Gene rearrangement has been reported in most insect orders and an increased rate of gene rearrangements has been observed in Hemiptera [10,73], Protura [75], and Hymenoptera [4,14,62]. Since the pioneering work of Boore et al. [2], it has been acknowledged that the gene structure of mitochondrial genomes contains phylogenetic signal [3]. Detection of gene rearrangements in lower-level elements of insect mitogenomes is expected to shed light on the evolution of these taxa [62,74]. The mitochondria of Hymenoptera generally exhibit extremely high molecular evolution rates and extensive gene rearrangements due to their parasitic lifestyles. However, our understanding of the mitochondrial genome of pteromalids is still limited, and the mechanism of gene rearrangement is still unclear.
The site-heterogeneous mixture model (CAT+GTR) implemented in PhyloBayes has been shown to correct the wrong grouping of unrelated taxa that share a similar base composition and an accelerated evolutionary rate [77]. ML and BI analyses yield similar topologies ( Figure 5). Because Pteromalidae mitogenomes show a high evolutionary rate, we selected the BI topology to reconstruct the ancestral states. In previous classification, Walkerella microcarpae and Micranisa ralianga belonged to Otitesellinae and Apocrypta bakeri, Philotrypesis tridentata, Philotrypesis pilosa, and Philotrypesis sp. to Sycoryctinae. Based our results, Otitesellinae were recovered as nested within Sycoryctinae, which made Sycoryctinae paraphyletic, a result consistent with Zhao et al. [40] and Cruaud et al. [78]. The clade comprising Otitesellinae and Sycoryctinae appeared sister to other Pteromalinae included in our analysis, a result consistent with Rasplus et al. [23]. In addition, Sycoryctinae were not recovered as monophyletic, a result consistent with Munro et al. [18]. However, according to the classification system recently proposed by Burks et al. [68] and based on the analyses by Cruaud et al. [78], our results also corroborated that members of the previous Otitesellinae and Sycoryctinae belong in fact to Pteromalinae and represent at most a tribe of this subfamily. All other pteromalid subfamilies were recovered as monophyletic. Therefore, our study supports the classification proposed by Burks et al. [68] based on nuclear Ultra-Conserved Elements and exons [78]. Our result shows the power of mitogenomes to reconstruct family-level phylogenies in Chalcidoidea, only for three of the eight pteromalid subfamilies and 11 of the 500+ recognized genera. The monophyly of Pteromalidae was not directly tested in the present study and this will be an important issue for future studies.
The pairwise breakpoint distances of mitogenome rearrangements appeared not to be completely consistent with the relationships observed among our Pteromalidae species; see Figure 6. However, the number of mitogenomes used in this study is still limited and increasing sampling is necessary to confirm this result.
Our studies highlighted new Pteromalidae gene rearrangements, revealed the evolutionary relationship between the main groups of Pteromalidae, accumulated molecular data for the study of Pteromalidae phylogeny, and provided a genetic background for biological control in agriculture and forestry.