Comparative Mitogenomes of Two Coreamachilis Species (Microcoryphia: Machilidae) along with Phylogenetic Analyses of Microcoryphia

Simple Summary Bristletails (Insecta: Microcoryphia) are primarily wingless insects, some of which have been found to exhibit parthenogenesis. In the genus Coreamachilis, parthenogenesis occurs in C. coreanus, whereas sexual reproduction is found in C. songi. Therefore, after obtaining mitochondrial genome sequences of these two species, we analyzed their selection pressure, based on phylogenetic trees of Microcoryphia. However, no positive selection was found in the mitochondrial protein coding genes of either C. coreanus or C. songi. In addition, a long hairpin structure was found between ND1 and 16S rRNA genes in Machilinae and Petrobiinae, which was highly consistent with the phylogenetic results. Abstract The order Microcoryphia, commonly known as bristletails, is considered as the most primitive one among living insects. Within this order, two species, Coreamachilis coreanus and C. songi (Machilidae: Machilinae), display the following contrasting reproductive strategies: parthenogenesis occurs in C. coreanus, whereas sexual reproduction is found in C. songi. In the present study, the complete mitogenomes of C. coreanus and C. songi were sequenced to compare their mitogenome structure, analyze relationships within the Microcoryphia, and assess adaptive evolution. The length of the mitogenomes of C. coreanus and C. songi were 15,578 bp and 15,570 bp, respectively, and the gene orders were those of typical insects. A long hairpin structure was found between the ND1 and 16S rRNA genes of both species that seem to be characteristic of Machilinae and Petrobiinae species. Phylogenetic assessment of Coreamachilis was conducted using BI and ML analyses with concatenated nucleotide sequences of the 13 protein-coding genes. The results showed that the monophyly of Machilidae, Machilinae, and Petrobiinae was not supported. The genus Coreamachilis (C. coreanus and C. songi) was a sister clade to Allopsontus helanensis, and then the clade of ((C. coreanus + C. songi) + A. helanensis) was a sister clade to A. baii, which suggests that the monophyly of Allopsontus was not supported. Positive selection analysis of the 13 protein-coding genes failed to reveal any positive selection in C. coreanus or C. songi. The long hairpin structures found in Machilinae and Petrobiinae were highly consistent with the phylogenetic results and could potentially be used as an additional molecular characteristic to further discuss relationships within the Microcoryphia.


Introduction
Bristletails (Insecta: Microcoryphia) are primarily wingless insects, living in places such as under fallen leaves, the bark of trees, moss on wet logs, etc. They can be found at the Insects 2021, 12, 795 2 of 17 seashore or inland all over the world. There are over 500 species in two families (Machilidae and Meinertellidae), and three subfamilies of Machilidae (Petrobiellinae, Petrobiinae, and Machilinae) have been reported [1,2]. The genus Coreamachilis belongs to the subfamily Machilinae. At present, there are two reported species of Coreamachilis; C. coreanus is distributed in Northern Korea and China [3], along with C. songi in China [4]. When Mendes described C. coreanus from North Korea, he did not find any male individuals, suggesting that parthenogenesis occurred in C. coreanus [5]. By contrast, C. songi was reported to show sexual reproduction [4].
Many studies have reported parthenogenetic lineages, particularly among species living at higher altitudes or in extreme habitats [6], indicating an association between parthenogenesis and environmental selection pressures that lead to specialized life cycles [7]. It was suggested that the loss of sexual reproduction was driven by a very general selective force [8], and selective factors promoted by variations in environmental conditions may be advantageous to parthenogenic reproduction [9]. By contrast, sexual reproduction is common among animals, and there have been many studies suggesting that natural selection can affect the characteristics of sexual reproduction [10,11]. In addition, parthenogenesis is rare in Microcoryphia [5]. Hence, we wondered if the mechanism of parthenogenesis in C. coreanus is under selective pressure. Also, we aimed to determine if there is a relationship between sexual reproduction and parthenogenesis within Coremachilis, related to selective pressure.
Insect mitogenomes are usually double-stranded circular molecules with a length of 14-20 kb, encoding 13 protein-coding genes (PCGs), 22 transfer RNAs (tRNAs), two ribosomal RNAs (rRNAs), and a control region (CR; or the AT-rich region) [12,13]. Because of their characteristics of fast evolution rates, small genome sizes, and low sequence recombination, mitogenomes have been used widely as molecular markers for phylogenetic analyses [13][14][15]. Many studies have now shown that specific gene rearrangements can occur in the mitogenomes of different families or genera [16][17][18]. Intergenic regions can also be used as a synapomorphy for a genus [18,19], and these regions can sometimes fold into stem-loop secondary structures or hairpin structures. Hairpin structures have been observed at PCG junctions in the mitochondrial genomes of various metazoans, including insects [20][21][22][23][24][25]. Kim et al. [25] suggested that hairpin structures are important for precise cleavage of the mature protein-coding genes. In addition, potential hairpin structures were also found in the A + T-rich region of insect mitogenomes, which was considered to be the initiation site of secondary strand synthesis [21,24,26]. However, the hairpin structures found in insects were relatively short (about 10 bp), and only a few studies have inferred an association between phylogenetic relationships and hairpin structures. We were interested in whether the hairpin structures were associated with synapomorphy at genus or family levels within Microcoryphia. It is well known that the mitogenome is often assumed to be an important neutral marker [27]. However, recent research has shown that the evaluation of selective pressures acting on mitogenome proteins can provide key insights into the adaptive evolution of the mitogenome [28][29][30][31][32][33]. To date, various studies of mammalian [34,35], avian [36], frog [33], fish [37], and insect mitogenomes have indicated that adaptive evolution has occurred. Among insects, this has included studies of Hymenoptera [38], Orthoptera [39], Ephemeroptera [40], Diptera [41], and Lepidoptera [42]. Thus, the mitogenome can be used as a molecular tool to explore/assess adaptive evolution.
In the present study, we sequenced the mitogenomes of C. coreanus and C. songi, and compared these mitogenomes with those of other bristletails available in the GenBank database. The 13 protein-coding genes were used to construct phylogenetic relationships of Microcoryphia, in order to discuss the relationship of Coreamachilis. A positive selection analysis was also used to assess whether C. coreanus, with parthenogenic reproduction, and C. songi, with sexual reproduction, were under positive selection at the mitogenome level.

Specimen Collection and DNA Extraction
The specimens of C. songi were collected from Lianyungang city, Jiangsu Province, China, and the specimens of C. coreanus were collected from Gongchangling, Liaoyang city, Liaoning Province, China. All specimens were identified by JY Zhang and stored in 100% ethanol in a −40 • C freezer. Total DNA was extracted from individual female specimens using QIAGEN DNeasy blood and tissue kit (QIAGEN, Hilden, Germany).

Polymerase Chain Reaction (PCR) Amplification and Sequencing
After modifying the primers designed by Simon et al. [43] and Ma et al. [44] in combination with published sequences of Machilinae, we designed 15 pairs of universal primers for amplification of mitogenomes (Table 1). PCR was performed using a BioRAD MJ mini personal thermal cycler (California, USA). Based on the sequence information obtained from earlier PCR runs using universal primers, several pairs of specific primers were then designed to obtain the whole mitogenome. Both PCR amplifications and reaction volume were carried out using methods described previously [17]. All PCR products were sequenced in both directions by Sangon Biotech Company (Shanghai, China).

Phylogenetic Analyses of Microcoryphia
Since Microcoryphia is located at the base of the class Insecta, we chose Xibalbanus tulumensis (Nectiopoda: Speleonectidae) [52] and Daphnia magna (Anomopoda: Daphniidae) (MT199637) as the outgroups for phylogenetic analyses in this study. To date, only 10 mitogenomes of Microcoryphia were available for download from the NCBI [44,47,[53][54][55][56]. Therefore, we combined all mitogenomes of Microcoryphia as the ingroup including C. coreanus and C. songi along with the two outgroups in order to discuss the phylogenetic relationships in Microcoryphia. Accession numbers of all bristletail mitogenomes are listed in Table 2. Sequences of the thirteen PCGs were extracted from the mitogenomes using PhyloSuite v1.1.16 [49], and aligned in batches with MAFFT integrated into PhyloSuite v1.1.16 using codon-alignment mode. Gblocks in PhyloSuite v1.1.16 was used to remove ambiguous sites. Finally, the 13 protein-coding genes were linked into a single line using concatenate sequence in PhyloSuite v1.1.16. We used the program PartionFinder 1.1.1 [57] to estimate the best partitioning scheme and model according to the Bayesian information criterion (BIC). ML analysis was conducted in RAxML 8.2.0 [58] with the best model of GTRGAMMA and 1000 bootstrap replications. BI analysis was conducted using the Mr-Bayes 3.1 [59] with the model of GTR + I + G set for 10 million generations with sampling every 1000 generations. The first 25% of generations were discarded as burn-in and the average standard deviation of split frequencies in MrBayes 3.1 was below 0.01.

Selection Pressure Analyses of Coreamachilis
The ratio of non-synonymous-to-synonymous (d n /d s ) ω can indicate selection pressure at the protein level, with ω = 1, ω > 1 and ω < 1 meaning neutral selection, positive selection and negative selection, respectively [60]. The ω value was calculated by the codon-based maximum likelihood method using the CodeML algorithm [60] implemented in the EasyCodeML [61]. Four different models were used to test the selection pressure, including the site model, clade model, branch model and branch-site model. We chose C. coreanus and C. songi as the foreground branch, respectively, and others as the background branch using three different model selections (the clade model, branch model and branch-site model). The branch models were performed under the one-ratio model (M0) or the two-radio model, the former presuming that ω in specific branches was different from the rest of the tree. Also, the branch-site models were combined with heterogeneous ω across sites and branches (model A) and can be compared against a null model (model A null), which made it possible to find neutral evolution and negative selection. Likelihood ratio tests (LRTs) and Bayes empirical Bayes (BEB) were used to assess these models and evaluate the posterior probability of positive selection sites, respectively. The likelihood ratio test (LRT) was used to compare the statistical model in order to determine whether there were differences between them.

Specimen Collection and Mitogenome Structure
We gathered over three hundred C. coreanus female individuals in our collection area, Liaoyang city, Liaoning Province, China, between 2005 and 2015, finding no males among these specimens. Similarly, Mendes [3] collected C. coreanus species in North Korea, but did not find any male individuals. Hence, it can be inferred that parthenogenic reproduction has emerged in C. coreanus. Only two species of Meinertellidae [62], two species of Petrobiinae [63], and several species of Machilinae [1,64] were previously found to exhibit parthenogenesis. In addition, we collected male and female individuals of C. Songi, indicating that it was indeed using the sexual reproduction strategy as reported [4,5].
The complete mitogenomes of C. coreanus and C. songi were circular DNA molecules of 15,578 bp and 15,570 bp in length, respectively (Table 3). By comparison, the other 10 sequenced mitogenomes of Microcoryphia ranged from 15,473 bp [47] to 16,197 bp [56]. Both the C. coreanus and C. songi genomes encoded 13 PCGs, two rRNA genes, and 22 tRNA genes ( Figure 1), which is consistent with typical insect mitogenomes [65]. Among these, 23 genes were located on the heavy (H) strand and the other 14 genes were on the light (L) strand ( Table 3). The gene order of the C. coreanus and C. songi mitogenomes were the same as those of typical insects. Among the 12 mitogenomes of Microcoryphia, including the two mitogenomes in this study, ten of them had the same gene arrangement, except for Petrobius brevistylis [55] and Trigoniophthalmus alternatus [56]. In the mitogenomes of C. coreanus and C. songi, the tRNA Ala -tRNA Arg -tRNA Asn -tRNA Ser1 -tRNA Glu -tRNA Phe cluster was found between ND3 and ND5 genes (Figure 1), as also occurred in eight other mitogenomes of bristletails, except for P. brevistylis [55], which showed a cluster of tRNA Arg -tRNA Asn -tRNA Ser1 -tRNA Glu -tRNA Ala -tRNA Phe and T. alternatus with a tRNA Ala -tRNA Arg -tRNA Asn -tRNA Ser1 -tRNA Glu -tRNA Tyr -tRNA Phe cluster.
There were some overlaps (34 bp and 35 bp) and intergenic regions (37 bp and 60 bp) in the mitogenomes of C. coreanus and C. songi, respectively. Among other Microcoryphia mitogenomes, the overlaps ranged from 13 bp to 38 bp, and the intergenic regions ranged from 77 bp to 230 bp. The intergenic region of C. coreanus was the shortest among the sequenced Microcoryphia, largely due to the short intergenic region between ND1 and tRNA Ser2 , which was always longer than 20 bp in most of the Microcoryphia species, but was only 6 bp in C. coreanus.   The nucleotide composition of the C. songi mitogenome was A = 36.38%, T = 35.33%, C = 16.68%, and G = 11.61%, and was very similar to C. coreanus, which was as A = 36.67%, T = 35.37%, C = 16.66%, and G = 11.30%. There were strong A + T biases in both C. songi and C. coreanus, 71.7% and 72.1%, respectively, and these were within the range of other Microcoryphia mitogenomes (67.9-74.4%). The skew statistics showed that there was a positive AT skew and negative GC skew in both C. coreanus and C. songi (Table 4). Table 4. Base composition, the percent of A and T, AT skew and GC skew of C. coreanus and C. songi.

Region
C

Protein-Codon Genes and Codon Usages
There were 13 PCGs identified in the mitogenomes of C. coreanus and C. songi, and nine of them were located on the heavy strand (H-strand), whereas the others were on the light strand (L-strand) (Figure 1). Most of the PCGs used ATN (N represents A, G, C, or T) as the start codon, except for ND2 in C. coreanus, which used TTG as the start codon (Table 3). Although ATN has been accepted as the canonical mitochondrial start codon for insect mitogenomes [65], TTG is also regarded as a start codon [66]. In fact, before being found in C. coreanus, TTG had only been found for the COX1 gene of the T. alternatus mitogenome in the Microcoryphia order [56]. In the C. coreanus mitogenome, there were seven genes that used ATG as the start codon, which was the same as in C. songi. Four genes used ATT as the start codon in C. coreanus and five genes used ATT in C. songi. The COX3 gene of C. coreanus and C. songi used ATA as the start codon (Table 5). The stop codon usage was the same in both C. coreanus and C. songi. Eleven genes used ATT as the stop codon, whereas an incomplete stop codon (T) was used in COX1 and ND1 genes ( Table 5). It is common, in metazoan mitochondrial genomes, to see an incomplete stop codon, and these truncated stop codons are presumed to be completed by post-transcriptional polyadenylation [67].
The A + T content of the 13 PCGs was 69.1% and 69% in C. coreanus and C. songi, respectively. The AT skew and GC skew were negative in both the species (Table 4). The RSCU of C. coreanus and C. songi is shown in Figure 2. The most frequently encoded amino Insects 2021, 12, 795 9 of 17 acids were Leu (UUR), Phe, and Ile (>300), whereas the least frequently used amino acid was Cys (<45). The A + T content of the 13 PCGs was 69.1% and 69% in C. coreanus and C. songi, respectively. The AT skew and GC skew were negative in both the species (Table 4). The RSCU of C. coreanus and C. songi is shown in Figure 2. The most frequently encoded amino acids were Leu (UUR), Phe, and Ile (>300), whereas the least frequently used amino acid was Cys (<45).

Ribosomal RNAs, Transfer RNAs and Hairpin Structures
Of the 22 tRNA genes, 14 tRNAs were encoded on the heavy strand (H-strand), whereas the other eight tRNAs were encoded on the light strand (L-strand). Both rRNAs were on the light strand ( Table 3). As in other Microcoryphia mitogenomes, the 16S rRNA gene was located between tRNA Leu and tRNA Val , with a length of 1332 bp and 1322 bp in C. coreanus and C. songi, respectively. Located between tRNA Val and the CR, the 12S rRNA gene was 813 bp and 806 bp in C. coreanus and C. songi, respectively.
The total length of the tRNAs was 938 bp and 1001 bp, with an average A + T content of 74% and 73.4% in C. coreanus and C. songi, respectively. In both the mitogenomes, most of the tRNA genes displayed the common cloverleaf secondary structure, except for tRNA Ser (AGN), which had lost the dihydrouridine (DHU) arm ( Figure S1), as often found in other insect mitogenomes [44].
A novel hairpin structure was located in the ND1 and 16S rRNA genes, with a length of about 60 bp in the mitogenomes of both ( Figure 3C,D). In fact, we found a similar hairpin structure in 7 out of the 10 known bristletails mitogenomes ( Figure 3A,B,E-I). Some hairpin structures have been found in other species, but all of them were short (less than 10 bp) and mostly located between PCGs [20][21][22][23][24][25]. Such sequences were suggested to have a potential role as a secondary strand-replication origin [25]. In this case, the stop codon of one of the PCGs was incomplete (T or TA), and the 3 -end region of this gene has the potential to form a hairpin structure with the beginning region of the neighboring protein-coding gene, which is important for the precise cleavage of the mature protein-coding genes [67][68][69]. In the mitogenomes of C. coreanus and C. songi, the stop codon of ND1 was incomplete, but the stem of the hairpin structure in the 16S gene was not at the beginning. Also, there were similar hairpin structures in other bristletail mitogenomes, but all the stop codons of ND1 in those mitogenomes were complete (Table 4), and the stems of the hairpin structures were not at the 3 -end region or beginning region. As for other genes with incomplete stop codons, we did not find any other long hairpin structures like this in other insect orders. In the A + T-rich region, it was reported that stem-loop structures could be formed. Also, it has been mentioned that conserved regions of single stem-loop structure (the 5 flanking sequences with "TTATA", while 3 flanking sequences with a "G(A)nT" motif) have been observed in a variety of insect orders [24,26]. This stem-loop structure in the A + T-rich region was suggested as the site of the initiation of secondary strand synthesis [21,24,26]. However, the hairpin structures we found were not in the control region, and none of the conserved structure was found. This hairpin structure between the ND1 and 16S RNA genes, first reported in Microcoryphia, appears to be novel and about 60 bp in length, with certain homologous segments (Figure 4), unlike the short hairpin structures previously found. No relevant studies about such a structure were found during an online search of the literature, and if a similar hairpin structure can be found in other species, it would be interesting to explore its function. Insects 2021, 12, x FOR PEER REVIEW 11 of 16

Phylogenetic Analyses and Selection Pressure Analyses
We constructed BI and ML trees with the nucleotide sequences of the 13 mitochondrial protein-coding genes, including 12 Microcoryphia species and two outgroup species (X. tulumensis and D. magna), to perform phylogenetic analyses (Figure 4). The topologies of BI and ML analyses were the same. The monophyly of Machilidae failed to be supported because Meinertellidae was clustered into Machilidae and formed a sister clade to Petrobiellinae. The monophyly of Machilinae and Petrobiinae also failed to be supported, because T. alternatus, belonging to Machilinae, was clustered into the subfamily Petrobiinae ( Figure 4). Similar results were also presented in the study by Ma et al. [44]. The previous main diagnosis of the Machilinae subfamily was via morphological characters of plesiomorphies [1], but our phylogenetic tree-based result was not congruent with this morphological taxonomy. The classification of Trigoniophthalmus genus should be reconsidered in combination with phylogenetic trees and morphological characteristics. The clade consisting of (Songmachilis xinxiangensis + (Allopsontus baii + (A. helanensis + (C. coreanus + C. songi)))) strongly supported the non-monophyly of Allopsontus and the monophyly of Coreamachilis in BI and ML analyses. The genus Coreamachilis was proposed as a monophyletic group, which is consistent with the morphological taxonomy, as well as sequence analysis of COX1 from the mitogenome [70]. However, the paraphyly of Allopsontus suggested that the subgenera of Allopsontus should be further assessed, because A. baii and A. helanensis belong to the genus Allopsontus s. s. Silvestri, 1911 and Allopsonsus (Anisopsontus) Mendes, 1990 [71]. Our results suggest that the morphological classification at the genus level can be further improved by combining phylogenetic relationships. It was interesting that the hairpin structure between ND1 and 16S RNA genes is a synapomorphy of the clade including Machilinae and Petrobiinae. (Figure 3). In addition, in the clade with the long hairpin structure, the homologous segments of hairpin structures aggregated together are more alike (Figure 4). The five species of Machilinae were grouped into a branch, and their stem sequences in the hairpin structure showed a high similarity ( Figure 5A). At the same time, the remaining clades with a hairpin structure were also of

A + T Rich Region
The largest non-coding region of bristletail mitogenomes was the control region, which was located between the 12S rRNA and tRNA Ile genes. The length of the A + T-rich region of C. coreanus and C. songi mtDNA was 736 bp and 734 bp, respectively. Ten other Microcoryphia mitogenomes, published online, showed lengths of the A + T-rich region ranging from 538 bp (N. australica) to 1149 bp (T. alternatus).

Phylogenetic Analyses and Selection Pressure Analyses
We constructed BI and ML trees with the nucleotide sequences of the 13 mitochondrial protein-coding genes, including 12 Microcoryphia species and two outgroup species (X. tulumensis and D. magna), to perform phylogenetic analyses (Figure 4). The topologies of BI and ML analyses were the same. The monophyly of Machilidae failed to be supported because Meinertellidae was clustered into Machilidae and formed a sister clade to Petrobiellinae. The monophyly of Machilinae and Petrobiinae also failed to be supported, because T. alternatus, belonging to Machilinae, was clustered into the subfamily Petrobiinae ( Figure 4). Similar results were also presented in the study by Ma et al. [44]. The previous main diagnosis of the Machilinae subfamily was via morphological characters of plesiomorphies [1], but our phylogenetic tree-based result was not congruent with this morphological taxonomy. The classification of Trigoniophthalmus genus should be reconsidered in combination with phylogenetic trees and morphological characteristics. The clade consisting of (Songmachilis xinxiangensis + (Allopsontus baii + (A. helanensis + (C. coreanus + C. songi)))) strongly supported the non-monophyly of Allopsontus and the monophyly of Coreamachilis in BI and ML analyses. The genus Coreamachilis was proposed as a monophyletic group, which is consistent with the morphological taxonomy, as well as sequence analysis of COX1 from the mitogenome [70]. However, the paraphyly of Allopsontus suggested that the subgenera of Allopsontus should be further assessed, because A. baii and A. helanensis belong to the genus Allopsontus s. s. Silvestri, 1911 and Allopsonsus (Anisopsontus) Mendes, 1990 [71].
Our results suggest that the morphological classification at the genus level can be further improved by combining phylogenetic relationships. It was interesting that the hairpin structure between ND1 and 16S RNA genes is a synapomorphy of the clade including Machilinae and Petrobiinae. (Figure 3). In addition, in the clade with the long hairpin structure, the homologous segments of hairpin structures aggregated together are more alike (Figure 4). The five species of Machilinae were grouped into a branch, and their stem sequences in the hairpin structure showed a high similarity ( Figure 5A). At the same time, the remaining clades with a hairpin structure were also of relative similarity ( Figure 5B). It is well known that gene order and genome organization can provide useful genetic information to understand evolutionary relationships [72]. Hence, intergenic regions can be used as a synapomorphy for a genus [18,19]. The hairpin structure found in this study is located in the gene spacer between ND1 and 16S, and as a special structure, it is highly correlated with the results of the phylogenetic tree. Hence, we indicate that the long hairpin structure may be useful as a molecular characteristic to discuss phylogenetic relationships within Microcoryphia.
can provide useful genetic information to understand evolutionary relationships [72] Hence, intergenic regions can be used as a synapomorphy for a genus [18,19]. The hairpin structure found in this study is located in the gene spacer between ND1 and 16S, and as a special structure, it is highly correlated with the results of the phylogenetic tree. Hence we indicate that the long hairpin structure may be useful as a molecular characteristic to discuss phylogenetic relationships within Microcoryphia.
The phylogenetic tree developed from the current data was used for selection pres sure analyses. Four different models were used in the test, including the site model (M3 vs. M0, M1a vs. M2a, M7 vs. M8, and M8 vs. M8a), clade model (CmC vs. M2a_rel), branch model (two-ratio model vs. M0), and branch-site model (A vs. A null) (Table S1). We chos Coreamachilis as the foreground branch and other species as the background branch in the clade model, branch model, and branch-site model. However, the LRTs did not indicate any sign of positive selection (p > 0.05). Since C. coreanus and C. songi reproduced differ ently, we calculated the selective pressures on each of them, to see if there was an adaptive evolution at the mitochondrial gene level for the differences in reproductive patterns However, the results did not indicate any sign of positive selection either (Table S2 and S3). Neither C. coreanus (with parthenogenesis) nor C. songi (with sexual reproduction were subject to any positive selection. In a word, it was suggested that differences in Coreamachilis reproductive patterns may not result from adaptive evolution on the mito genome.

Conclusions
We successfully sequenced the complete mitogenomes of C. coreanus and C. songi which showed similar gene characteristics to the other 10 species of Microcoryphia with published mitogenomes. We found a special long hairpin structure (about 60 bp) located between ND1 and 16S RNA genes of C. coreanus and C. songi, which was also found in seven known bristletail mitogenomes belonging to Petrobiinae and Machilinae. However the function and formation of this long hairpin structure remain unclear, and still need to be explored. Since the hairpin structure is highly correlated with the phylogenetic tree the long hairpin structure can be used as the synapomorphy for Machilinae and Petrobi inae. Both BI and ML analyses supported the genus Coreamachilis as a monophyleti group, whereas the monophyly of Allopsontus was not recovered. Neither C. coreanus (with parthenogenesis) nor C. songi (with sexual reproduction) were subject to any positive se lection.  (Table S1). We chose Coreamachilis as the foreground branch and other species as the background branch in the clade model, branch model, and branch-site model. However, the LRTs did not indicate any sign of positive selection (p > 0.05). Since C. coreanus and C. songi reproduced differently, we calculated the selective pressures on each of them, to see if there was an adaptive evolution at the mitochondrial gene level for the differences in reproductive patterns. However, the results did not indicate any sign of positive selection either (Tables S2 and S3). Neither C. coreanus (with parthenogenesis) nor C. songi (with sexual reproduction) were subject to any positive selection. In a word, it was suggested that differences in Coreamachilis reproductive patterns may not result from adaptive evolution on the mitogenome.

Conclusions
We successfully sequenced the complete mitogenomes of C. coreanus and C. songi, which showed similar gene characteristics to the other 10 species of Microcoryphia with published mitogenomes. We found a special long hairpin structure (about 60 bp) located between ND1 and 16S RNA genes of C. coreanus and C. songi, which was also found in seven known bristletail mitogenomes belonging to Petrobiinae and Machilinae. However, the function and formation of this long hairpin structure remain unclear, and still need to be explored. Since the hairpin structure is highly correlated with the phylogenetic tree, the long hairpin structure can be used as the synapomorphy for Machilinae and Petrobiinae. Both BI and ML analyses supported the genus Coreamachilis as a monophyletic group, whereas the monophyly of Allopsontus was not recovered. Neither C. coreanus (with parthenogenesis) nor C. songi (with sexual reproduction) were subject to any positive selection.

Data Availability Statement:
The data supporting the findings of this study are openly available in National Center for Biotechnology Information (https://www.ncbi.nlm.nih.gov), accession numbers were MW752137 (Coreamachilis coreanus) and MW752138 (Coreamachilis songi).