The specimen of Apachyus feae de Bormans, 1894 was collected from Laibin, Guangxi Province of China (24.1402°N, 110.1844°E) in October of 2019; the specimen of Diplatys flavicollis Shiraki, 1907 was collected from Jurong, Jiangsu Province of China (32.1325°N, 119.0743°E) in February of 2020. The specimens were identified by the author, preserved in 100% ethanol, and stored in the Insect Collection of Jiangsu University of Science and Technology (ICJUST). The total genomic DNA of the two earwigs was isolated using the E.Z.N.A. Tissue DNA Kit (Omega Bio-Tek, Inc.) and preserved at −20 °C before the sequencing process.

The Illumina TruSeq short-insert libraries (size = 450 bp) were constructed using 1 μg of purified DNA fragments and were sequenced by Illumina Hiseq 4000 (Shanghai Biozeron Biotechnology Co., Ltd). Raw reads were filtered prior to assembly; high-quality reads were retained and assembled into contigs by SOAPdenovo2.04 (Luo et al. 2012). The assembled contigs were then aligned to the reference mitogenomes of Dermaptera using BLAST. Subsequently, the aligned contigs (≥80% similarity and query coverage) were arranged according to the reference mitogenomes. Finally, the clean reads were mapped to the assembled draft mitogenomes to fix the wrong bases; gaps were filled using GapFiller v. 2.1.1 (https://sourceforge.net/projects/gapfiller/). The mitogenome sequences of A. feae and D. flavicollis were deposited in GenBank under the accession numbers MW291948 and MW291949, respectively. Mitochondrial gene analyses of A. feae and D. flavicollis were compared to four additional species of Dermaptera with available mitogenomes (Table 1). The gene order was compared with Drosophila yakuba Burla, 1954, which was considered to possess the ancestral arthropod mitochondrial gene arrangement (Clary and Wolstenholme 1985).

All protein-coding genes (PCGs) and ribosomal RNA (rRNA) genes were identified by homology alignments. Gene boundaries of each PCG were further confirmed by ORF finder (https://www.ncbi.nlm.nih.gov/orffinder/). All transfer RNA (tRNA) genes were predicted and illustrated by the MITOS online server (Bernt et al. 2013). The visual structure of the mitogenomes were depicted using CGView Server (http://stothard.afns.ualberta.ca/cgview_server/) (Grant and Stothard 2008). Nucleotide composition of each gene and codon usage of PCGs were calculated using MEGA v. 6.0 (Tamura et al. 2013). The composition skew values were calculated by AT-skew = [A – T] / [A + T] and GC-skew = [G – C] / [G + C] formulas (Perna and Kocher 1995). The synonymous substitution rate (Ks) and nonsynonymous substitution rate (Ka) were computed by DnaSP v. 5.10 (Librado and Rozas 2009). Presumed secondary structures in the control regions were predicted by the online tool Tandem Repeats Finder (http://tandem.bu.edu/trf/trf.advanced.submit.html), DNAMAN v. 6.0.3 and ARWEN (http://mbio-serv2.mbioekol.lu.se/ARWEN/) (Laslett and Canbäck 2008).

Nucleotide sequences of PCGs derived from six species of Dermaptera, including A. feae and D. flavicollis sequenced in this study, were used in the phylogenetic analysis (Table 1). The stonefly Kamimuria chungnanshana Wu, 1938 (Plecoptera, Perlidae; GenBank accession no. NC_028076) was used as the outgroup. The 13 PCGs were respectively aligned by MAFFT and concatenated as a combined dataset using SequenceMatrix v. 1.7.8 (Katoh and Standley 2013). The optimal nucleotide substitution models and partitioning schemes for the dataset was determined by PartitionFinder v. 2.1.1 using the Bayesian Information Criterion (BIC) and a greedy search algorithm (Lanfear et al. 2016). Bayesian inferences (BI) and Maximum likelihood (ML) analyses were conducted with the optimal partition schemes. The BI analysis was conducted by MrBayes v. 3.2.7, with 20 million generations sampling every 1000 generations, running one cold chain and three hot chains with a burn-in of 25% trees (Ronquist and Huelsenbeck 2003). TRACER v. 1.5 was used to examine the stability of the results of the BI analysis. The ML analysis was performed by RAxML v. 8.2.12 with 1000 bootstrap replicates (Stamatakis 2014). FigTree v. 1.4.2 was used to adjust and visualize the tree files generated by both BI and ML inferences.

The complete mitogenome of A. feae is a typical double-strand circular molecule with a length of 19,029 bp (Fig. 1). The obtained partial mitogenome of D. flavicollis is 12,950 bp in length (Fig. 1). The completely sequenced three mitogenomes of Dermaptera range in size from 16,087 bp in E. arcanum to 20,456 bp in C. fletcheri. In the mitogenome of A. feae, an extra trnV gene and two control regions are found in addition to the standard set of 37 genes (13 PCGs, 22 tRNA genes and two rRNA genes) (Table 2). In the partial mitogenome of D. flavicollis, 10 and a partial fragment of PCGs, 18 tRNA genes, two rRNA genes, and a control region are annotated (Table 3). In A. feae, there are 56 overlapping nucleotides located in three pairs of neighboring genes, and the longest overlap is 41-bp long and located between trnT and ND4L (Table 2). A total of 296 intergenic nucleotides (IGNs) are dispersed in 19 locations for A. feae. In D. flavicollis, 17 overlapping nucleotides and 504 IGNs are found, including a 227-bp long IGN between trnS2 (UCN) and ND1 (Table 3).

Figure 1. 

Mitochondrial maps of Apachyus feae and Diplatys flavicollis. Genes outside the map are transcribed clockwise, whereas those inside the map are transcribed counterclockwise. Names and other details of the genes are listed in Tables 2 and 3. The inside circles show the GC content and the GC skew. GC content and GC skew are plotted as the deviation from the average value of the entire sequence.

The mitogenomes of A. feae and D. flavicollis are biased toward A and T nucleotides (61.2% and 73.5%, respectively), which is consistent with other earwigs (Table 1). The A+T contents were also rich in each mitochondrial gene, showing the highest in trnD of A. feae and trnF of D. flavicollis.

In the sequenced earwigs, no PCG rearrangement are found (Fig. 2). In A. feae, most tRNA genes in the gene cluster trnA-R-N-S1-E-F are rearranged, and an extra trnV is present in the gene cluster. In D. flavicollis, the gene cluster trnA-R-N-S1-E-F is also rearranged and incorporates trnY, trnC and trnQ from other locations. In C. fletcheri, trnI, trnC, trnY, trnQ, and trnE are rearranged (Wan et al. 2012). In E. arcanum, trnQ, trnC, trnY, trnR, and trnS1 are rearranged, and trnY is lost. In E. metallica and P. flavocapitatus, both trnR and trnS1 are absent. These tRNA rearrangements mainly occur in the trnA-R-N-S1-E-F gene cluster. The two rRNA genes are located in the same location for all sequenced earwigs; however, they are variable in size interspecifically. In addition to the tRNA rearrangements, the control region of D. flavicollis transfers to the new location between ND3 and ND5; an extra control region is also found in A. feae and C. fletcheri (Wan et al. 2012).

Figure 2. 

Mitochondrial gene arrangement of six earwigs in comparison with Drosophila yakuba.

All PCGs of A. feae are annotated, whereas ND2, COX2, and partial COX1 of D. flavicollis are not sequenced. The PCGs of A. feae are similar in size to those of D. flavicollis and other earwigs. Most PCGs of A. feae and all PCGs of D. flavicollis utilize the standard ATN start codon (ATT, ATC, and ATG), whereas ATP8 and COX3 of A. feae start with special start codons (GTG and TTG, respectively) (Tables 2, 3). Most PCGs of A. feae and all PCGs of D. flavicollis have the complete termination codon TAN (TAA or TAG), whereas ND4 and CYTB of A. feae end with an incomplete stop codon T (Tables 2, 3). The relative synonymous codon usage (RSCU) values were calculated for the six earwig mitogenomes (Fig. 3). The most frequently used codon is TCT (Ser) for A. feae, TTG (Leu) for E. metallica, TTA (Leu) for D. flavicollis, C. fletcheri, E. arcanum, and P. flavocapitatus.

Figure 3. 

Relative synonymous codon usage (RSCU) of PCGs in six species of earwigs.

The ratio of Ka/Ks was calculated for each PCG of the six earwig mitogenomes to evaluate the evolutionary rates of the PCGs (Fig. 4). The results showed that COX1 of E. metallica has the highest evolutionary rate, followed by ND5 of A. feae and ND2 of P. flavocapitatus, whereas COX1 of A. feae and E. arcanum appear to be the lowest. The genes with ratios of Ka/Ks above 1 are evolving under positive selection. Other genes with ratios of Ka/Ks below 1 are expected to evolve under purifying selection.

Figure 4. 

Evolutionary rates of PCGs in six species of earwigs. The bar indicates each gene’s Ka/Ks value.

The typical set of 22 tRNA genes and an extra trnV gene are detected in the mitogenome of A. feae (Fig. 5). In D. flavicollis, 18 tRNA genes are recognized and the four tRNA genes trnI, trnM, trnW, trnL are absent due to the incomplete sequencing of 5´ end (Fig. 6). In other sequenced earwigs, C. fletcheri has all 22 tRNA genes (Wan et al. 2012), E. arcanum lacks trnY, and E. metallica and P. flavocapitatus lack trnR and trnS1. Individual tRNA gene of the two newly sequenced mitogenomes range in size from 60 to 75 bp; the longest tRNA gene is trnY in D. flavicollis, and the shortest tRNA gene is the extra trnV in A. feae. In the mitogenomes of A. feae and D. flavicollis, most of the tRNA genes exhibit cloverleaf secondary structures, but the dihydrouridine (DHU) arm is lost for the extra trnV of A. feae and is reduced for trnS1 of both species. The anticodons of the tRNA genes were identical among the earwigs. In the tRNA genes of A. feae and D. flavicollis, a total of 48 and 25 mismatched base pairs are respectively recognized and all of them are G-U pairs.

Figure 5. 

Secondary structures of tRNA genes in the mitogenome of Apachyus feae. Mismatched base pairs are indicated by red circles; reduced arms are indicated by red arrowheads.

Figure 6. 

Secondary structures of tRNA genes in the mitogenome of Diplatys flavicollis. Mismatched base pairs are indicated by red circles; reduced arms are indicated by red arrowheads.

Two rRNA genes are consistently found in all sequenced mitogenomes. Locations of the two rRNA genes are conserved among earwig species and similar to D. yakuba, but the lengths are variable. In A. feae, the large ribosomal RNA (rrnL) gene is 1280 bp in length with an A+T content of 67.9%; the small ribosomal RNA (rrnS) gene is 739 bp with an A+T content of 66.0%. In D. flavicollis, the rrnL gene is 1364 bp with an A+T content of 76.1%; the rrnS gene is 960 bp with an A+T content of 76.9%.

Two putative control regions (CRs) are found in the mitogenomes of A. feae, E. metallica and P. flavocapitatus. The CR1 of A. feae is 756 bp and located after rrnS, containing a stem-loop (SL) structure and a poly-[TA]n like stretch (Fig. 7). The CR2 of A. feae is 3285-bp long and located between trnS2 (UCN) and ND1, being composed of five SL structures and three copies of tandem repeats. The CR of D. flavicollis is 592 bp and located between trnQ and trnS1, comprising two and partial copies of tandem repeats, two tRNA-like structures, and a poly-[T]n stretch (Fig. 8). In C. fletcheri, the 1816-bp long CR1 contains a SL structure and two regions of tandem repeats; the entire 2856-bp long CR2 comprises 21.1 copies of tandem repeats (Fig. 9). The CR of E. arcanum is 686 bp in size, containing a SL structure, a poly-[TA]n stretch and a tandem repeats region (Fig. 9). The 891-bp long CR of E. metallica comprises four SL structures (Fig. 9). The CR of P. flavocapitatus is short, 227-bp in size, and contains one SL structure (Fig. 9).

Figure 7. 

Predicted structural elements in the control regions of Apachyus feae.

Figure 8. 

Predicted structural elements in the control region of Diplatys flavicollis.

Figure 9. 

Predicted structural elements in the control regions of Challia fletcheri, Euborellia arcanum, Eudohrnia metallica, and Paratimomenus flavocapitatus.

The phylogenetic analyses use the nucleotide sequences of six available earwig mitogenomes to investigate the mitochondrial phylogenetic relationships within Dermaptera. The two phylogenetic trees using BI and ML analyses generated identical topological structures for Dermaptera (Fig. 10). The monophyly of Forficulidae is supported with high values. Diplatyidae is recovered as the sister group of Anisolabididae and their combined clade is grouped with Pygidicranidae. Apachyidae is supported as the sister group to other sequenced families. Monophyly of the two infraorders Protodermaptera and Epidermaptera cannot be supported by either analysis. The three parvorders Paradermaptera, Metadermaptera, and Eteodermaptera are each represented by single family and their relationship was recovered as Paradermaptera + (Eteodermaptera + Metadermaptera).

Figure 10. 

Phylogenetic relationships within Dermaptera inferred by Bayesian inference and maximum likelihood analysis. Numbers at the nodes are posterior probabilities (left) and bootstrap values (right). The family names are listed after the species. Infraorders and parvorders are indicated below each family name.