Specimens of Ooencyrtus plautus were collected from Fuzhou city, Fujian province, in September 2020. They were reared in the laboratory, then processed for DNA extraction. Total genomic DNA was extracted using the cetyltrimethyl ammonium bromide (CTAB) method (Shahjahan et al. 1995).

Sequencing was performed using a whole genome shotgun (WGS) strategy on the Illumina Novaseq platform. The quality of data was checked using FastQC (Andrews 2013). Raw data were filtered into high-quality data after processing with AdapterRemoval v. 2. Then, Single Nucleotide Polymorphism (SNP), Insertion and Deletion (InDel), Copy Number Variation (CNV), and Structural Variation (SV) analysis were conducted to ensure the reliability of nucleotides. The assembly of the mitochondrial genome was accomplished with Novoplasy v. 2.7 (Dierckxsens et al. 2017).

Gene annotation of mitochondrial genome was performed using MitoZ and 13 protein-coding genes as well as 22 tRNA genes are annotated (Meng et al. 2019). Twenty-two tRNA genes were verified with the use of MITOS WebServer, setting the parameters with the Invertebrate Mito genetic code (Minh et al. 2013). Every sequence of tRNA genes was manually checked separately. According to the secondary structure of tRNA predicted by MITOS, it was drawn manually using Adobe Illustrator. Annotation of the 13 protein-coding genes were refined manually by identifying the corresponding open reading frames using the invertebrate mitochondrial code. The mitogenome maps were produced using Organellar Genome DRAW (OGDRAW) (Marc et al. 2013).

Geneious v. 11.0.2 was used to examine all genes in the mitochondrial genome (Kearse et al. 2012). MEGA X (Sudhir et al. 2018) was used to analyze base composition and relative synonymous codon usage (RSCU). To calculate the number of synonymous substitutions for each synonymous site (Ks) and the number of non-synonymous substitutions for each non-synonymous site (Ka) through DnaSP 5. AT/GC skewness is calculated as AT skewness = (A − T) / (A + T) and GC-skew = (G − C) / (G + C) (Rozas et al. 2003). To measure the relative composition of different bases was measured based on GC and AT skew.

In this study, a total of 10 species of mitogenomes were analyzed, of which Encarsia formosa and Encarsia obtusiclava of family Aphelinidae were selected as outgroups (Table 1). MAFFT v. 7.3.1 was used to align the PCGs (Kazutaka and Standley 2016). Alignments of individual genes were concatenated to generate two kind of data sets: 1) the PCG matrix, including all three codon positions of protein-coding genes; 2) the AA matrix, translated 13 PCG into amino acids.

Both ML (maximum likelihood) and BI (Bayesian inference) analyses were performed on the concatenated data set used for phylogenetic reconstruction. In W-IQtree (Trifinopoulos et al. 2016), the best-fit model was used for maximum likelihood analysis. The site-heterogeneous mixture model (CAR+GTR) was used in PhyloBayes analysis. The trees were sampled every 1000 generations (Huelsenbeck 2012). FigTree v. 1.3.1 was used to view the generated tree.

The nearly complete mitogenome of O. plautus was 15,730 bp in size, including 13 PCGs, 22 tRNAs, 2 rRNAs, and a nearly complete control region. In Ooencyrtus plautus, 27 genes (15 tRNAs, 2 rRNAs, and 10 PCGs) were encoded by the majority-strand (J-strand), and 10 genes (3 PCGs and 7 tRNAs) were encoded by the minority-strand (N-strand). The values of AT skew and GC skew are often used to reveal the nucleotide composition of the mitochondrial genome (Alexandre et al. 2005). In this study, the nucleotide compositions of eight complete or nearly complete mitogenomes in Encyrtidae were investigated by calculating the percentages of AT-skew and GC-skew (Fig. 3). The results of the nucleotide skew statistics showed that the AT-skews in PCGs, tRNAs, and rRNAs of encyrtid mitogenomes were almost all positive, while the GC-skews were almost all obviously negative.

Figure 3. 

Nucleotide composition of various data sets of mitogenomes. Hierarchical clustering of Encyrtidae species based on the AT-skew and GC-skew.

The nucleotide composition O. plautus was significantly biased toward adenine and thymine, with an A + T content of 84.6% (Table 2). The PCGs on the N strand had a T skew and slight G skew. However, the PCGs on the J strand exhibited a C skew. The tRNAs was negative for AT skews, besides the tRNAs on the N strand, which had an equal G and C distribution. Two rRNA genes (lrRNA and srRNA) were encoded on the J strand and exhibit an A skew and slight C skew.

Intergenic spacers and overlapping genes are very common in arthropod mitochondrial genomes (Ma et al. 2015; Chen et al. 2018). In the mitochondrial genome of O. plautus, a total of 166 bp of intergenic spacers ranging from 1 to 47 bp were found in eight locations. The minimum intergenic spacers (1 bp) located at trnT-trnP, trnQ-trnW, nad6-cob, and trnS2-nd1. The maximum (47 bp) gene spacer was located at trnW-trnC. There were 18 overlapping gene regions, ranging from 2 bp to 7 bp in length in the nine mitogenomes. The longest overlapping sequence (7 bp) was between atp6-atp8 and nad4-nad4l.

In the mitochondrial genome of O. plautus, the total length of protein-coding genes was 11,072 bp, accounting for 70.39% of the entire genome. Most of them were encoded on the J strand, and only nad2, nad6, and cob were encoded on the N strand. The average A + T content of the 13 protein-coding genes was 82%, ranging from 73.7% (cox1) to 93.80% (atp8) for individual genes.

To investigate further this high A and T content, and the frequency of synonymous codon usage, we calculated relative synonymous codon usage (RSCU) values. The relative synonymous codon usages (RSCU) of the O. plautus are shown in Fig. 4. Taken together, the most frequently used codons are UUA (Leu2), AGA (Ser1), and CGA (Arg), whereas those ending in G or C, CUG, CUC, CAG, and GGC were the less frequently used codons. The codons ending with A or T are predominant, which at least partly leads to the bias towards A and T.

Figure 4. 

Relative synonymous codon usage (RSCU) of the mitochondrial genomes of Ooencyrtus plautus.

The predicted initiation codons are ATN as in most other insect mitochondrial genomes (Dowton and Austin 1999). There were six genes (cox1, cob, nad4, nad6, atp6, and cox3) starting with ATG and seven genes (nad1, nad3, nad6, nad4L, nad5, cox2, and atp8) starting with ATT. All protein-coding genes terminated at the most common stop codon TAA, except for nad1 and nad5, which stopped with single T (Table 3).

Previous work (Zhang et al. 2020) showed that Eulophidae and Encyrtidae have close relationships, so we chose Eulophidae species as the reference sequence. Using Chouioia cunea as a reference sequence, we calculated the non-synonymous substitution rate (Ka), synonymous substitution (Ks), and Ka/Ks ratio of each PCG for O. plautus (Fig. 5). Ka and Ks are non-synonymous and synonymous substitution rates, respectively. They are controlled by functionally related sequence contexts, such as encoding amino acids and participating in exon splicing (Parmley and Hurst 2007). The ratio of the two parameters Ka/Ks (a measure of the strength of selection) is defined as the degree of evolutionary change (Wang et al. 2011). A value of Ka/Ks greater than 1 means positive selection exists, indicating that non-synonymous mutations are more favored by Darwinian selection, and they will be retained at a rate greater than synonymous mutations. In O. plautus, the Ka/Ks ratio of nad2 was greater than 1, indicating that the nad2 gene had a positive selection. This phenomenon was first discovered in Encyrtidae. In addition, the Ka/Ks of atp8 is as high as 0.9389, which is the highest except for nad2. The high Ka/Ks phenomenon of nad2 and atp8 is also found in other species (Jia et al. 2020; Guo et al. 2021; Xu et al. 2021). The reason for this phenomenon may be that the evolution speed of a gene is related to its function (Wang et al. 2011).

Figure 5. 

Evolutionary rates of mitochondrial genomes in Ooencyrtus plautus. The numbers of nonsynonymous substitutions per nonsynonymous site (Ka), the number of substitutions per synonymous site (Ks), and the ratio of Ka/Ks for each gene are given, using Chouioia cunea as the reference sequence.

In total, 22 transfer RNA genes were found, ranging in size from 59 bp (trnS1) to 71 bp (trnT). Most of the tRNAs were encoded on the J strand and only 7 tRNAs (trnT, trnE, trnK, trnC, trnQ, trnR, and trnS2) were encoded on the N strand. The average nucleotide composition of these tRNAs was A: 43.7%, T: 45.5%, C: 6.6%, and G: 4.2%, with a total average A + T content of 89.2%. All tRNA sequences can be folded into the canonical cloverleaf secondary structure, except for trnS1 which lacked the dihydrouridine (DHU) arm. A lack of the DHU arm in trnS1 was found in the mitochondrial genomes of most insects (Dowton et al. 2002b). Changes in the length of the DHU and TΨC arms lead to differences in the size of the tRNA sequence (Shao et al. 2001). In addition, some mismatches (G-U in trnD, trnA, and trnV, U-U in trnS1, and two U-U in trnA) were found in O. plautus.

The two rRNA genes, the larger ribosomal gene (rrnL) and the smaller ribosomal gene (rrnS), were located between trnA and trnL1, and trnV and trnA, respectively. The average of the total size of two rRNAs was 2,103 bp and the average A + T content was 89%.

Figure 6. 

Predicted secondary cloverleaf structure for the tRNAs of Ooencyrtus plautus.

Gene rearrangement in the mitochondrial genome is a relatively rare event in the evolutionary history of insects (Cameron 2014). However, in the Hymenoptera lineage, more and more species have mitochondrial gene rearrangements (Shao et al. 2003). Gene rearrangements are a common phenomenon found in almost all the sequenced Chalcidoidea mt genomes (Fig. 2). In the Chalcidoidea mitochondrial genomes studies, most of the tRNA genes are rearranged (Tang et al. 2018). In the Hymenoptera, numerous rearrangements of protein-coding genes have been identified in several groups (Xiao et al. 2011; Wei et al. 2014). Compared with the putative ancestral pattern of the insect mitochondrial genome, dramatic gene rearrangements, not only in tRNA genes but also in protein-coding genes, were found in O. plautus mitochondrial genomes. In the mitogenome of O. plautus, a total of 17 genes, including six PCGs, 10 tRNAs, and one rRNA gene are rearranged (Fig. 7).

Figure 7. 

Mitochondrial genome organization in Encyrtidae. Colored boxes indicate rearranged gene clusters, and gray boxes indicate conserved gene blocks.

In Encyrtidae, the gene order of all species of Encytus are identical. We used their gene order as a template to analyze the gene rearrangement of Encyrtidae. It can be summarized into four obvious rearrangements. They are the rearrangement of gene clusters trnI-nad2-trnW, trnY-trnS1-trnN-trnC-trnR, and trnA-trnQ-rrns-trnV-trnM. In addition, there are inversions of trnE-trnF and trnP-trnT.

Firstly, the trnW-trnI-nad2 gene cluster in these four species (Aenasius arizonensis, Diaphorencyrtus aligarhensis, Metaphycus eriococci, and Platencytus parkeri) was rearranged as trnI-nad2-trnW, most species in Chalcidoidea of which were in this order. However, in O. plautus, trnW was translocated, causing the gene cluster trnI-nad2-trnW to be divided into trnI-nad2 and trnW, which is the first of its kind to be reported in Chalcidoidea. Besides, the trnI-nad2 gene cluster was translocated to upstream. Secondly, the trnN-trnS1-trnY-trnC-trnR gene clusters in four species (A. arizonensis, D. aligarhensis, M. eriococci, and O. plautus) have undergone disorderly rearrangements. Thirdly, trnE-trnF inversion occurred in both A. arizonensis and D. aligarhensis. In addition, trnP-trnT inversion also occurred in D. aligarhensis. Finally, in D. aligarhensis and O. plautus, trnQ in the gene cluster trnA-trnQ-rrnS-trnV-trnM has a long-distance transposition. In M. eriococci, trnM and trnQ were transposed in the gene cluster trnA-trnQ-rrnS-trnV-trnM.

Phylogenetic relationships were analyzed using the concatenated nucleotides and amino acids sequences of 13 PCGs from eight encyrtid species and two outgroups. Four topologies were constructed using both ML and BI methods and two different data sets.

The topological structures of these two phylogenetic trees reconstructed by ML analysis were identical with the phylogenetic tree reconstructed by BI based on PCG (Fig. 8). While Bayesian tree reconstructed on the AA data set showed different phylogenetic relationships in the clade consist of Encyrtus rhodococcusiae, E. eulecaniumiae and E. sasakii which was E. eulecaniumiae + (E. rhodococcusiae + E. sasakii) rather than E. rhodococcusiae + (E. eulecaniumiae + E. sasakii). The deviation and rate heterogeneity of the nucleotide composition of the mitogenomes are the basis of the rapid evolution of Chalcidoidea, which leads to inconsistent topologies based on different data types and analysis strategies (Wu et al. 2020). Currently, the CAT + GTR model implemented in PhyloBayes software was found to be the best fitting model for all data sets. In addition, the AA data set is considered as the best fit matrix for reconstructing phylogenetic trees (Li et al. 2015, 2017).

Figure 8. 

Phylogenetic tree produced by maximum likelihood and Bayesian inference analyses based on the 13 PCG and the AA data set A phylogenetic tree produced by Bayesian inference analyses based on 13 PCG data set B phylogenetic tree produced by maximum likelihood analyses based on the AA data set C phylogenetic tree produced by maximum likelihood analyses based on the 13 PCG data set.

In the BI topology based on the AA data set, A. arizonensis and M. eriococci formed a clade representing Tetracneminae. The remaining six species formed a monophyletic clade representing Encyrtinae. In Encyrtinae, Encyrtus formed a monophyletic clade as a sister group to the clade formed by O. plautus and D. aligarhensis. Encyrtus sasakii and E. rhodooccisiae are most closely related in this monophyletic clade. In the BI topology based on the PCG data set, E. eulecaniumiae and E. sasakii form a sister group then to E. rhodooccisiae. Encyrtus infelix was the first diverged in Encyrtus.

Pairwise breakpoint distances (PBD) between the mitochondrial genomes of each species in Encyrtidae were calculated using the web server CREx, and heatmaps were constructed (Fig. 9). In the BI topology based on the AA data set, E. infelix was the first divergence branch in Encytus, which is consistent with the structure of the rest phylogenetic trees. The value of pairwise breakpoint distances between E. infelix and D. aligarhensis was 14, and the value of E. infelix and O. plautus was 10. Lower values of PBD indicated closer relationships which was consistent with the topology among E. infelix, O. plautus, and D. aligarhensis on the phylogenetic tree.

Figure 9. 

Mitogenomic phylogeny of eight Encyrtidae species and two outgroups based on the AA data set using the Bayesian inference. Heatmap of PBD values among encyrtid species is nested within the phylogeny, and the PBD values are listed.