M. pulchricornis, M. sp. 1 and M. sp. 2 were all collected in the Chinese provinces of Zhejiang (Ningbo city), Guizhou (Guiyang city), and Hebei (Shijiazhuang city), respectively. The specimens were morphologically identified by Prof. Cornelis van Achterberg (Zhejiang University, China). All specimens were initially preserved in 100% ethanol and then stored at 4 °C before DNA extraction. Whole genomic DNA was extracted from every sample using the DNeasy tissue kit (Qiagen, Hilden, Germany).

The libraries were prepared for each DNA sample using the VAHTS^® Universal DNA Library Prep Kit. All constructed libraries were then sequenced as 150 bp paired-end on a full run (2 × 150 PE) using MGISEQ2000 platform. FastQC v0.11.9 (Andrews, 2015) was used to check the data quality, and fastp v0.23.1 (Chen et al., 2018) was used to trim adaptors and remove low quality reads with default parameters. More than 5 GB of clean data for each sample was used in de novo assembly. The mitogenomes were assembled using MitoZ v2.3 (Meng et al., 2019), IDBA v1.1.3 (Peng et al., 2012) and SPAdes v3.13.0 (Bankevich et al., 2012) with default parameters, respectively. Subsequently, to verify the accuracy of de novo assembly, one fragment in the cox1-cox2 junction was amplified with primers (C1-J-2195: 5’-TGA​TTT​TTT​GGG​CAT​CCT​GAA​GT-3’, C2-N-3665: 5’-CCA​CAA​ATT​TCA​GAA​CAT​TGA​CC-3’) for three species by polymerase chain reaction (PCR) and was then Sanger-sequenced. Finally, all assemblies were then integrated with GENEIOUS v2020.0.5 (Biomatters Ltd. San Diego, CA, United States).

Three assembled mitogenomes were annotated using MITOS web server (http://mitos.bioinf.uni-leipzig.de/index.py; accessed on 15 November 2022) (Bernt et al., 2013). The start and stop positions of 13 PCGs were manually adjusted and corrected by aligning published data of the Braconidae species in GenBank. tRNAScan-SE (Chan et al., 2021) was used to verify the results of putative tRNA genes. Nucleotide composition, codon usage, and relative synonymous codon usage (RSCU) values were estimated using MEGA 11 (Tamura et al., 2021). The bias of nucleotide composition was measured as AT-skew = (A − T)/(A + T) and GC-skew = (G − C)/(G + C) (Perna and Kocher, 1995). The gene rearrangements in Meteorus mitogenomes were analyzed by comparison with the ancestral insect mitogenomes (Boore, 1999).

A total of 21 Braconidae mitogenomes were used for phylogenetic analyses, including three newly obtained mitogenomes. Two Ichneumonidae species, Euceros kiushuensis and Diadegma semiclausum, were used as outgroups (Supplementary Table S1). The PCGs were aligned using MAFFT v7.464 (Katoh and Standley, 2013) with the default algorithm, and the best partition schemes and substitution models (Supplementary Table S2) for the datasets were analyzed by PartitionFinder v1.1.1 (Lanfear et al., 2012). Bayesian inference (BI) and maximum likelihood (ML) were selected to reconstruct the phylogenetic trees. For BI analysis, MrBayes v3.2.7a (Ronquist and Huelsenbeck, 2003) was used to run four independent Markov chains for 100 million generations, with tree sampling occurring every 1,000 generations and a burn-in of 25% of the trees. The stationarity of the run was assessed by Tracer v1.7. (ESS values >200) (Rambaut et al., 2018). Maximum likelihood (ML) analysis was performed with RAxML-HPC2 v8.2.12 (Stamatakis, 2014) under the GTRGAMMA model. A total of 200 runs for different individual partitions were conducted with 1,000 bootstrap replicates.

The mitogenomes of three Meteorus species were successfully obtained, and the newly sequenced mitogenomes were submitted to GenBank (accession numbers: OP832526 - OP832528) (Supplementary Table S1). All sequences were nearly complete mitogenomes, measuring 15,590 bp for M. pulchricornis, 15,799 bp for M. sp. 1, and 16,639 bp for M. sp. 2. In all assembled mitogenomes, all 13 protein-coding genes (PCGs), 22 tRNA genes, and two rRNA genes were found (Figure 1). Non-etheless, the A-T control region was incomplete in all genomes due to the high A+ T content in the Hymenoptera mitogenomes (Zheng et al., 2021; Jasso-Martinez et al., 2022a).

The mitogenome of M. pulchricornis contained 19 non-coding regions ranging in size from 1 to 67 bp, with a total length of 407 bp. The nucleotides from 8 overlap regions were up to 22 bp in total. The maximum overlap length was 7 bp, located at two junctions (atp8-atp6, nad4-nad4l), while the minimum overlap length was 1 bp and occurred at four junctions (atp6-cox3, trnA-trnR-trnN, trnT-trnP). M. sp. 1 had 18 non-coding regions with lengths ranging from 1 to 83 bp, for a total of 303 bp. The nucleotides from 5 overlap regions amounted to 18 bp in total. Two junctions (atp8-atp6, nad4-nad4l) had overlaps of 7 bp in length. M. sp. 2 possessed 23 non-coding regions ranging from 1 to 465 bp in length and a total length of 1,138 bp, which had the longest intergenic nucleotides among the three Meteorus mitogenomes. It had 3 overlap regions with a total of 12 bp in length. Three junctions (nad4-nad4l, atp8-atp6, and trnE-trnF) had overlaps of 7 bp, 4 bp and 1 bp in length, respectively. The length of overlap (4 bp) at the junction atp8-atp6 in the mitogenomes of M. sp. 2 was unusually short compared to that of other wasps, which was typically 7 bp (Zhu et al., 2018; Tang et al., 2019).

The A + T content for the sequenced region of the mitogenomes was 84.41% (M. pulchricornis), 82.70% (M. sp. 1) and 84.31% (M. sp.2) (Supplementary Table S3), which was similar to other wasp species (Tang et al., 2019). All 13 PCGs were detected in the newly generated mitogenomes, with sizes ranging from 11,063 bp (M. sp. 1) to 11,120 bp (M. sp. 2). The entire A + T content of all the PCGs was from 80.48% (M. sp. 1) to 83.37% (M. pulchricornis) (Supplementary Table S3). The size of the PCGs in the three mitochondrial genomes was similar to other wasps (Zheng et al., 2021; Shu et al., 2022; Zheng et al., 2022). The AT-skew in three Meteorus mitogenomes was negative (−0.1574 to −0.1479), and the GC-skew was positive (0.0662–0.1088) (Supplementary Table S3). All 22 typical tRNA genes were found in three mitogenomes. The size of all tRNA genes identified ranged from 57 bp to 72 bp (Supplementary Table S4). Two rRNA genes (rrnS and rrnL) were identified in all mitogenomes. The length of rrnS was 764 bp (M. pulchricornis), 735 bp (M. sp. 1) and 790 bp (M. sp. 2), and the size of rrnL was 1,253 bp, 1,325 bp and 1,347 bp, respectively (Supplementary Table S3).

The majority of the protein-coding genes used ATN (ATG, ATT, or ATA), except ATC, as an initiation codon. Although most PCGs terminated with the conventional TAA and TAG as stop codons, some of them have incomplete stop codons (T or TA) in some PCGs (Supplementary Table S5). The start and stop codons were typical of insect mitogenomes (Li et al., 2021; Ge et al., 2022). Relative synonymous codon usage (RSCU) values of three Meteorus species were analyzed (Figure 2; Supplementary Table S6). The total codon number of these three species was 3,689, 3,677 and 3,695, respectively. The codon GCG was not detected in all species, while codons UGC, CGC, AGC and GGC appeared in some species. Meanwhile, the five most frequently used codons UUA, AUU, UUU, AUA and AAU were observed in the mitogenomes due to the high A + T content in the mitogenomes. These results are consistent with other published wasp mitogenomes (Zhu et al., 2018; Zheng et al., 2022).

The order of PCGs was relatively conservative in the mitogenomes of Braconidae. So far, only Stenocorse bruchivora (Doryctinae), two Chelonus spp. (Cheloninae), and two Cotesia spp. (Microgastrinae) were found PCGs rearrangement in mitogenomes (Wei et al., 2010; Jasso-Martinez et al., 2022a; Yuan et al., 2022). In contrast, tRNA rearrangements occurred in all known Braconidae mitogenomes (Li et al., 2016; Jasso-Martinez et al., 2022a). PCG rearrangements did not occur in any of the four Meteorus mitogenomes studied, but tRNA rearrangements varied (Figure 3). Only seven tRNAs (trnW, trnY, trnL2, trnH, trnT, trnP and trnV) were conserved and trnG had different locations in each of the four mitogenomes. trnL1 and trnS2 were inverted and exchanged. Furthermore, trnC was translocated from nad2-cox1 junction to cox3-nad3 junction. The tRNA cluster (trnA-trnR-trnN-trnS1-trnE-trnF) between nad3 to nad5 exhibited two distinct patterns. In M. sp. USNM and M. pulchricornis, trnE was translocated upstream of trnA, and trnF was rearranged to cox3-nad3 junction, which reduced one tRNA in the tRNA cluster. In M. sp. 1 and M. sp. 2, the number of tRNAs between nad3 to nad5 remained at six, and trnN was translocated downstream of trnF. In addition, the positions of trnA and trnR were interchanged in M. sp. 2. Interestingly, trnI and trnM, which were usually conservative in Braconidae mitogenomes, were translocated from upstream of nad2 to the nad1-rrnL junction, and trnL2 was duplicated into the nad1-rrnL junction as well in M. sp. USNM. In summary, the pattern of tRNAs suggested that M. sp. USNM and M. pulchricornis were more closely related than the other two species.

The rich and diverse rearrangements presented by the tRNAs in the genus Meteorus were the first to be found in the Braconidae. tRNA rearrangements are usually conserved in the same genus in Hymenoptera (Jasso-Martinez et al., 2022a). Although different rearrangement of tRNAs within the same genus in Hymenoptera has been reported in Aphelinidae, Chrysididae and Ichneumonidae (Zhu et al., 2018; Zheng et al., 2021; Zheng et al., 2022), the tRNA rearrangements within Meteorus in Braconidae were more dramatic than in the aforementioned families, which exhibited only one or two different tRNA rearrangements in the same genus. To our knowledge, the various tRNA rearrangement patterns within the same genus had not been found in the other insects.

In this study, a phylogenetic analysis based on 13 PCGs was constructed using maximum likelihood and Bayesian methods, in order to validate the phylogenetic position of Meteorus within Braconidae (Figure 4). The consensus tree derived from the above two inference methods generated the congruent results. Most nodes had significantly supported bootstrap and posterior probability values. The tree recovered the well-accepted major lineages within the family and supported the division of Braconidae into cyclostomes and non-cyclostomes as generally accepted (Jasso-Martinez et al., 2022a; Jasso-Martinez et al., 2022b). Meanwhile, consistent with previous studies, within the non-cyclostomes, Euphorinae was recovered as a sister clade to the remaining non-cyclostomes (Jasso-Martinez et al., 2022b). The four Meteorus species formed a clade within Euphorinae, and were close to Zele chlorophthalmus rather than Dinocampus coccinellae. In the Meteorus, the species were clustered in two clades: one comprised of M. sp. USNM and M. pulchricornis while M. sp. 1 and M. sp. 2 grouped together in the second clade. This structure matched the tRNA rearrangement result presented above.