Novel mitochondrial genes and gene reannotation: conservation of gene arrangements among available astigmatid mites

Mitochondrial genomes (mitogenomes) of metazoans typically contain 37 genes, comprising 13 protein-coding genes, two rRNA genes, and 22 tRNA genes. To date, complete mitogenome sequences of 15 species of Astigmatina are available, and they present variation in a number of features, such as gene arrangements, tRNA unconventional secondary structures, and the number and internal structures of control regions. Furthermore, 11 astigmatid mites from six superfamilies share the same gene arrangement. Two available species from the genus Histiostoma reportedly have different mitochondrial (mt) tRNA gene arrangements. We sequenced the mitogenomes of Lepidoglyphus destructor and Gohieria fusca, both from the superfamily Glycyphagoidea (Astigmatina). In total, 37 mt genes were identied in the two Glycyphagoidea species. Based on AT content and stem-loop structures, we divided the largest non-coding regions (LNRs) in L. destructor and G. fusca into two domains, respectively. The novel feature of two domains for the LNR was also found in Acalvolia sp. (Astigmatina, Hemisarcoptoidea). Using MITOS 2, tRNAScan, ARWEN, and manual approaches, we reannotated the mitogenomes of Histiostoma blomquisti, H. feroniarum, and Trouessartia rubecula. We reannotated six tRNA genes in H. blomquisti and four tRNA genes in H. feroniarum. We were able to identify all of the mt tRNA genes that were reported as lost in Tr. rubecula. The phylogenetic relationships found in our study were fairly consistent with previous studies of astigmatid mites phylogeny. Within Astigmatina, Glycyphagoidea was recovered as a monophyletic group. tRNA-D- T-arms. 21 sequences in other astigmatid mites based on their anticodon and secondary structure. fusca, identied 19 automated identied other three tRNA genes (trnN, trnV, and trnA)

rearrangements in the astigmatida investigated to date.

Mitochondrial genomes of Lepidoglyphus destructor and Gohieria fusca
We determined the complete mitochondrial genomes of L. destructor (GenBank accession number: MT075728) and G. fusca (GenBank accession number: MN608156), with lengths of 14,663 bp and 14,741 bp, respectively. Both contained the typical set of 13 protein-coding genes (PCGs), 22 transfer RNAs (tRNAs), two ribosomal RNAs (rRNAs), and ( Table 1). The overall base composition of the entire L. destructor mitogenome consisted of 25.0% A, 39.2% T, 15.2% C, and 20.6% G. In G. fusca, the overall base composition of the entire mitogenome consisted of 25.8% A, 38.4% T, 15.8% C, and 20.1% G. In both species, the nucleotide composition had the same AT bias (64.2%). The combined protein-coding sequence lengths in the two species were 10,785 bp and 10,827 bp, respectively. In the G. fusca mitogenome, 11 protein-coding genes used standard start codons known as ATN. In addition, several unorthodox initiation codons were used. Atp8 initiated with GTG, and nad3 started with codon TTG (Table 1). In L. destructor, 12 protein-coding genes used standard start codons, and nad3 starts with codon TTG (Table 1). Table 1. Mitochondrial genome organization of Lepidoglyphus destructor and Gohieria fusca. Int = intergenic nucleotides; negative numbers indicate overlapping nucleotides between adjacent genes.
Mitogenomes of L. destructor and G. fusca have been predicted to have the full set of tRNA genes. The putative mt tRNA genes were highly truncated in both L. destructor (46 to 61 bp) and G. fusca (48 to 61 bp). Only trnK showed the typical cloverleaf in these two species (Fig. S1, S2). Other tRNAs showed the reduction of tRNA-D-and/or T-arms. In L. destructor, we identi ed 21 tRNAs using automated prediction. We determined the trnV manually by alignment with homologous sequences in other astigmatid mites based on their anticodon and secondary structure. In G. fusca, we identi ed 19 tRNA genes using automated prediction. We identi ed the other three tRNA genes (trnN, trnV, and trnA) manually.
The novel feature of the largest non-coding region Two large (>50 bp) non-coding regions (NCRs) could be identi ed in the L. destructor and G. fusca mitogenomes ( Table 2). In L. destructor, the largest noncoding region (LNR, 813 bp) had an overall AT content of 68.4% and could be divided into two domains, domain I (463 bp) and domain II (350 bp), based on AT content and stem-loop structures (Fig. 1A). Domain I was on the J-strand at the 3' end of trnF and had an AT content of 62.2%. A peculiar feature of this domain is the presence of a stem-loop structure with a length of 432 bp and short sequences (5'-GGGGGTAGGGG and CCCCTACCCCC-3') ( Fig. 1A I-box).
Domain II was on the J-strand at the 5' end of the trnS1 gene. The AT content was 78.9%, which was more AT-rich than domain I. Comparative sequence analysis of this domain identi ed conserved sequences, one microsatellite-like AT-repeat element, and several stem-loop structures (Fig. 1). In the mitochondrial genome of G. fusca, the LNR was 861 bp long and had an overall AT content of 70.7%. The LNR also could be divided into two domains. Domain I had a stem loop structure with a length of 435 bp (Fig. 1A). The short sequences (5'-GGGGGTAGGGG and CCCCTACCCCC-3') were not found in domain I. The AT contents of domain I and domain II were 57.3% and 86.6%, respectively. The conserved sequences, one microsatellite-like AT-repeat element, and several stem-loop structures, were found in domain II. In total, we found the feature of two domains in both species from Glycyphagoidea.  Table 6.
To investigate whether the feature of two domains existed in other astigmatid mites, we conducted a comparative analysis of non-coding regions from all available Astigmatid mitogenomes (Table 2). In astigmatid mites, most reported mitogenomes feature a compact structure that usually contains two conserved site-speci c NCRs (except for Ty. longior, H. blomquisti, and H. feroniarum) and several nonconserved NCRs. These conserved site-speci c NCRs are located between trnF-trnS1 and trnW-nad1 with high AT content (the average content is 84.8.0% and 71.9%, respectively). The LNR is usually located in trnF-trnS1, and the length of this region averages 420 bp. The three longest conserved NCRs are in L. destructor, G. fusca, and Acalvolia sp. (813, 861, and 753 bp, respectively). The shortest NCRs located in trnF-trnS1 is only 76 bp in Sarcoptes scabiei. Table 2. Distribution of NCRs in the astigmatid mites mitochondrial genomes.
Note: Non-coding regions marked with a star are assumed to be putative control regions. The data for Trouessartia rubecula, Histiostoma blomquisti, and H. feroniarum in the table are after our revision.
The other conservative non-coding region was located in trnW-nad1, and the length of this region averaged 72 bp ( Table 2). The longest reached 273 bp in Acalvolia sp.. The shortest was in Ardeacarus ardeae, with only 29 bp. A similar sequence [5'-(G) n TA(G) n -3'] was found in the NCR of most available astigmatid mites (including the Tr. rubecula and the H. feroniarum after reannotation), except for Ardeacarus ardeae. A similar sequence [5'-(G)nTA(G)n -3'] also was found in domain I of the LNR for L. destructor and Acalvolia sp., except domain I of the LNR for G. fusca. This sequence, [5'-(G)nTA(G)n -3'], seemed to be conserved in most available astigmatid mites.
Additionally, some astigmatid mites also exhibited other nonconservative NCRs (Table 2). Both H. blomquisti and H. feroniarum had one nonconservative NCR, and several stem-loop structures were found in these NCRs. The sequences and the stem-loop structures did not seem to be conserved. Compared with the other described astigmatid mite NCRs, only Acalvolia sp. (Astigmata, Hemisarcoptoidea) showed a similar feature of two domains in the LNR. In Acalvolia sp.,the LNR also could be divided into two domains: a 469 bp fragment of domain I and a 284 bp fragment of domain II. The AT contents of domain I and domain II were 71.6% and 84.9%, respectively. In domain I, we found short sequences (5'-GGGGGTAGGGG and CCCCTACCCCC-3') and a stemloop structure. The conserved sequences among astigmatid mites were found in domain II (Fig. 1A).
As described earlier, the feature of two domains for the LNR was found only in L. destructor, G. fusca, and Acalvolia sp. among all available astigmatid mites.

Reannotation of four tRNA genes for Histiostoma blomquisti
In terms of previous study for H. blomquisti, the smallest four tRNAs (trnA, trnS2, trnR, and trnV) could only be annotated manually, but the remaining tRNAs (trnC, trnF, and the other 16 tRNAs) were identi ed using more than a manual approach [27]. To nd more probable structures, we reannotated four tRNAs (trnC, trnA, trnF, and trnS2) based on the mitogenome from our analysis. We identi ed trnF and trnS2 manually, and veri ed trnC and trnA based on the minimum free energy (MFE).
Our trnF identi ed by manual sequence alignment was more conservative, showing fewer mismatches in the acceptor stem. This tRNA lies on the J-strand at the 3′ end of nad5. Without predictions from tRNA search programs, we manually retrieved a tentative trnS2 from the sequence of previous study [27]. Similarly, we manually inferred the putative trnS2 to be a D-loop with fewer mismatches on stems. As a common phenomenon in Astigmatina, this D-loop was extremely truncated (43 bp). Compared with H. feroniarum and other astigmatid species, trnS2 in our study was more conserved (Table 3). Table 3. The alignment of nucleotide sequences of four mitochondrial tRNA genes (trnS2, trnW, trnF, and trnI) in ve reported astigmatid species in four different superfamilies.
Whenever contradictory predictions occurred, we calculated the minimum free energy (MFE) as a proxy. According to the MFE, we determined trnC and trnA in our study. These were more probable as having fewer mismatches on stems and arms (Fig. 2). After being amended, trnC apparently shared 54bp nucleotides with the adjacent trnP, these two genes overlapped and were on opposite strands (Table S1). Because of the trnF occupancy, the boundaries of nad5 in H. blomquisti changed. Due to the tRNA reannotation, the NCRs also were changed. In a previous study, the NCRs were located in rrnS-trnV, trnF-trnA. After the revision, one conserved site-speci c NCR was identi ed between trnF-trnS1, with high AT content. The other NCR was still located in rrnS-trnV, with the length increased to 1598 bp ( Table 2).
In this study, we used manual annotation to reannotate trnF, trnW, and trnI of H. feroniarum. Our trnF for H. feroniarum was identi ed by manual alignment and was conserved with the reannotated trnF in H. blomquisti (Table 3). In our analysis, trnW was reannotated by sequence alignment. This sequence also had fewer mismatches (Table 3). In the previous study, the D-loop trnI with the less common anticodon sequence AAU was reported in H. feroniarum [18], but D-loop trnI had never been reported in other Astigmatina. Our trnI was reannotated manually based on anticodon sequences and secondary structures. In our study, we also observed the less common anticodon sequence AAU. We inferred this reannotated trnI to be a TV-loop structure. This sequence was aligned to H. blomquisti but was the only other remaining member of Histiostomatoidea. The reannotated trnI presented considerable similarity (Table 3).
To select the most probable tRNAs, minimum free energy (MFE) was calculated for three tRNA genes (trnC, trnV, trnQ). We found smaller MFE values than in previous annotations. In our study, trnC was retrieved between trnS2 and trnP, and trnV was retrieved between 12S and 16S, whereas the positions were coincidentally the opposite from those predicted by Xue et al. [18]. Our trnC had common anticodon sequences GCA among astigmatid mites, whereas the previously described trnC used ACA instead of GCA in anticodon sequences. When we calculated MFE, both trnC and trnV in our study were more thermodynamically stable (Fig. 2). After reannotation, the trnV had 52 bp overlapping with the contiguous rrnS on the same chain. Notably, the previously described trnQ was annotated on the J-strand at the 3' end of nad5 [18], but many mismatches on stems and arms were found in this position. In fact, most nucleotides of trnQ among Astigmatina were highly conserved, and when located on the N-strand downstream of trnS1, the reannotated trnQ had a smaller MFE.
After amendment of the tRNA genes, the boundaries of nad2 of H. feroniarum changed because of trnI occupancy. We also observed changes to the rrnL size of H. feroniarum on account of trnW occupancy (Table S1). With alignment by Clustal W 2.0, the rrnL boundaries of the two Histiostoma species were more conserved. In addition, we changed one position of the NCR because of the reannotation of trnW and trnI. After the revision, we could identify one conserved site-speci c NCR with high AT content between trnW and nad1.
Notably, the reannotations of H. blomquisti and H. feroniarum mitogenomes led to the mtDNA gene order for both Histiostoma species sharing the gene arrangement of a possible common ancestor of astigmatid mites, excluding the different locations of non-coding regions (Fig. 3). Fig. 3. Mitochondrial gene arrangements in the sarcoptiform mites. The arrow pointing to the right represents the (+)-strand, and the arrow to the left the (−)strand. Translocated or inverted genes are color-coded (green: inversion and translocation; pink: translocation; yellow: inversion). rRNA genes are in blue. The possible common ancestor among astigmatid mites is found in different genera and families. Species marked with a star: tRNAs reannotated in our study. Species marked with a triangle: the gene order in previous studies.
Retrieving the "lost" mt tRNA genes in Trouessartia rubecula In Trouessartia rubecula, Esteban et al. predicted tRNAs based solely on the MITOS WebServer, and ve tRNA genes (trnA, trnE, trnI, trnY, and trnV) were not identi ed [23]. In protein annotation, the common start codon is ATN, but the cox1 gene start codon of T. rubecula is the uncommon TCT. T. Rubecula was reported to contain three D-loop sequences. The lengths of these non-coding regions are abnormally short, and the location of these region is different from the common non-coding regions.
In the current study, we used tRNAscan-SE, ARWEN, and MITOS 2 to predict the missing trnE, trnI, and trnY, and we annotated trnA and trnV using manual comparison. The secondary structure with the smallest constrained MFE was considered to be the most likely one (Fig. 4). We veri ed the tRNA genes of the two sarcoptiform mites by comparison of secondary structures and conserved nucleotide sequences with several Astigmatina species ( Table 4). Nucleotides that pair at the arms (acceptor arm, D arm, anticodon arm, and T arm) are underlined in the gure. We analyzed the codon usage in the mt protein-coding genes of the 17 species of astigmatid mites to see whether and how the corresponding codons of the "lost" tRNA genes were used. Overall, the codon use was similar between Tr. rubecula and other astigmatid mites (Fig. 5).
Note: Bars indicate Watson Crick base pairings; dots and circles represent GU pairs and mismatched pairs, respectively. MFE indicated the minimum free energy value. Table 4 The alignment of nucleotide sequences of two mitochondrial tRNA genes (trnA and trnV) in six species: Conserved nucleotides are shaded in gray across these ve mites.  After revision, we reannotated the position of the cox1 gene. The reannotated cox1 gene had a length of 1536 bp and had ATA as the initiation codon rather than the uncommon start codon TCT identi ed in a previous study. The result from BLASTn supported our boundaries for cox1. The LNR (301 bp in length) was arranged in the same way as most LNRs in the available sequences for astigmatid mites.
New insight into astigmatid mitochondrial gene arrangement Compared with Limulus polyphemus as the hypothetical ancestor of arthropods, gene orders were rearranged in the mitogenomes of 17 astigmatid mites (Fig.  3). We found that the mitogenomes of 13 astigmatid mites from six superfamilies (Glycyphagoidea, Acaroidea, Analgoidea, Sarcoptoidea, Pterolichoidea, and Hemisarcoptoidea) shared a consistent gene order, indicating that this gene arrangement was from the possible common ancestor of astigmatid mites. The mitogenomes of Histiostoma blomquisti, Histiostoma feroniarum, and Acalvolia sp. represented three types of gene arrangements resulting from the different locations of the non-coding regions from the gene arrangement of the possible common ancestor of astigmatid mites (Fig 3). The mt gene order of Tyrophagus longior was rearranged compared with other astigmatid mites because of tRNA loss.

Phylogenetic analysis
We constructed a phylogenetic tree based on a nucleotide data set from 13 mt PCGs of 17 astigmatid mites. The ML and BI analyses showed fully consistent topologies, and we provided the percentage of the bootstrap support at each node. The phylogenetic tree showed that the monophyly of Glycyphagoidea recovered with strong support (BPP = 1 and BSP = 100%), as well as the monophyly of Acaroidea, Hemisarcoptoidea, and Histiostomatoidea (Fig. 6). The monophyly of two superfamilies, Analgoidea and Sarcoptoidea, was always rejected. The logical next steps would be the integration of more astigmatina mitogenomes, including species from each of the six superfamilies.

Discussion
For the rst time, we determined the complete L. destructor and G. fusca mitogenomes in the superfamily Glycyphagoidea. Like most astigmatid mites, we identi ed 37 mitochondrial genes and one CR in the two Glycyphagoidea species. For comparison, the 15 other species of Astigmatina with complete mitochondrial genomes available in GenBank exhibited an average genome size of 14,227 bp with high AT content. These values closely agreed with those we observed for both Glycyphagoidea species. The length variation of mitogenomes was usually due to length variation in the LNRs [28].
The LNRs in the mitochondrial genomes were less conserved than the coding genes, and they have not been well studied compared with other mitochondrial genes [29]. Astigmatid mites have two conservative non-coding regions and several nonconservative non-coding regions. The LNRs were usually located in trnF-trnS1, and several conserved sequences were found in the regions, including one microsatellite-like (AT)n element and several stem-loop structures. The conserved sequences were common and present in 12 astigmatid mites, including L. destructor and G. fusca, and these sequences were involved in the regulation of mitogenome replication [30].
In addition to abovementioned features, in L. destructor and G. fusca, based on AT content and hypothetical stem-loop structures, we divided the LNR into two domains and identi ed several conserved features in both species. Domain II was more AT-rich than domain I. Similar LNR features have not been reported for other Acari. A comparison of the LNR with available Astigmata mitochondrial sequences showed that similar unusual features of the LNR could be found only in Acalvolia sp. (Hemisarcoptoidea). Whether this LNR feature is general needs to be investigated by sequencing of additional Astigmata mitogenomes.
Truncated tRNAs have been found in all published mitogenomes of Acariformes [11,12]. In addition, mismatches in the stems occur regularly in astigmatid mites. Both such tRNAs pose challenges for tRNA prediction programs [31]. In previous studies with two species of the genus Histiostoma, tRNA gene arrangements were different from each other. In Tr. rubecula, ve tRNAs were reported to be lost. In the current study, we reannotated the mitogenomes of H. blomquisti, H. feroniarum, and Tr. rubecula. In Histiostoma, in addition to prediction programs, we used manual annotation by sequence comparison when alternative anticodons were considered [18]. We reannotated trnC, trnA, trnF, and trnS2 in H. blomquisti, and trnC, trnV, trnQ, trnF, trnW, and trnI in H. feroniarum. The reannotation of trnF (H. blomquisti) was based on a comparison among astigmatid mites. The trnS2 (H. blomquisti) and the trnF and trnW (H. feroniarum) were identi ed by alignment with the tRNA of the only other species in the same genus. The trnI was reannotated to have a TV-loop, which was a more conserved secondary structure than the D-loop of trnI identi ed previously. The D-loop of trnI was extremely rare in Astigmatina, only known in H. feroniarum [18]. We used the MFE to verify trnC and trnA (H. blomquisti) and trnC, trnV, and trnQ (H. feroniarum). Compared with the trnC from a previous study, our trnC (H. feroniarum) had the common anticodon among astigmatid mites. In conclusion, the tRNAs after our reannotation contained more conserved sequences among homologues in species of Histiostoma or Astigmatina, and had more stable secondary structures than previously identi ed. Among these manually retrieved tRNAs, the extreme situations were found in trnC of H. blomquisti and trnV of H. feroniarum. The genes overlapped between trnC-trnP on the opposite strand; thus, the transcripts did not overlap [32]. The overlap between trnV-rrnS in H. feroniarum was on the same chain. Sequence overlaps universally existed in astigmatid species and may be corrected later through post-transcriptional editing, alternate treatment, and relaxed restrictions on tRNA structure [33,34]. Whether the consistent gene arrangement for Histiostoma is general needs to be investigated by sequencing additional Histiostomatoidea mitogenomes.
After the revision, the Tr. rubecula had a full set of tRNAs. Codons for all 22 amino acids were present in the protein-coding genes in all astigmatid mites. No evidence has yet indicated that nuclear tRNA genes can be imported into mitochondria in mites or other animals [19]. If an mt tRNA gene is indeed lost and there is no nuclear replacement, then its corresponding codons in the mt protein-coding genes will not be translated. Thus, the "loss" of mt tRNAs may not exist in astigmatid mites. Rearranged gene orders in the mitogenomes not only can provide phylogenetic information but also are useful in resolving the phylogeny of insects [35,36]. Therefore, the correct annotation of genes is important for phylogenetic and evolutionary studies. Our study had four types of gene arrangements in Astigmatina. The mitogenomes of 13 astigmatid mites from six superfamilies shared one type of gene rearrangement. This gene rearrangement likely was the most conserved gene arrangement among the astigmatid mites. In addition, these four types of gene arrangements for 37 genes (13 PCGs, 22 tRNAs, and 2 rRNAs) were consistent. The gene arrangement of Ty. longior represented a unique type because of missing tRNAs [20]. Two species had complete mitogenomes reported in the genus Tyrophagus. The gene arrangement of Ty. putrescentiae was consistent with the possible common ancestor of astigmatid mites [22]. Murillo et al. used COI to analyze molecular phylogenetics and suspected that a misidenti cation existed in the Ty. longior [20] sequence, as it clustered with other sequences of Ty. putrescentiae [36]. Hence, our results did not support the loss of mt tRNA genes in these astigmatid mites. Our results did support the hypothesis that the conserved gene arrangement in most available species was likely the arrangement of the possible common ancestor of astigmatid mites [11].
Within the Astigmatina, we observed the monophyly of Glycyphagoidea, Hemisarcoptoidea, Acaroidea, and Histiostomatoidea. Consistent with a previous study, the monophyly of two superfamilies, Sarcoptoidea and Analgoidea, was rejected [11,23,31]. Our nding supported the suggestion that the phylogeny of Astigmata is far from resolved at the superfamily level [11,38,39].

Conclusion
In this study, we determined the complete L. destructor and G. fusca mitogenomes. A novel feature of the LNR was found in L. destructor, G. fusca, and Acalvolia sp., the feature that was not found in other available Astigmatina mitochondrial sequences. We reannotated the mitogenomes of Tr. rubecula, H. blomquisti, and H. feroniarum. The gene order of 13 astigmatid mites from six superfamilies shared the gene order that was identical to the possible common ancestor of astigmatid mites. This information from our study has important rami cations for our understanding of mitogenome evolution in astigmatid mites.

Mite collection and DNA extraction
We collected samples of Gohieria fusca from a our factory in Wuhu, southeast China, in July 2016. We collected samples of L. destructor from a our shop in Wuhu in spring 2018. L. destructor and G. fusca were reared at 25°C and 85% relative humidity in the dark using the following medium: wheat germ, corn our, and baker's yeast at a ratio of 10:10:1 by weight [40]. Mites were stored in 100% ethanol at -20°C until use. We identi ed specimens based on their morphological characteristics [41]. The cox1 sequences of L. destructor and Gohieria fusca were obtained with the primers LWCOI_U and LWCOI_L, as suggested by Webster et al. [42]. After performing BLASTn searches of the nucleotide collection (nr/nt) database of the NCBI, the results showed that the cox1 sequence from two specimen were very similar to several sequences from L. destructor and Gohieria fusca. In L. destructor, the identities were 98.91% (accession number: EU078972) and 98.64% (accession number: AY525569). In G. fusca, the identities were 98.64% (accession number: MG279724) and 98.38% (accession number: MG279719). Therefore, the two specimens collected in the current study were identi ed as L. destructor and G. fusca. We extracted whole-genomic DNA by standard phenol-chloroform extraction [43].

Sequence assembly, annotation, and analysis
We obtained the mitogenomes of L. destructor and G. fusca by the Shanghai Majorbio Bio-pharm Biotechnology Company (Shanghai, China) and Shanghai BIOZERON Company (Shanghai, China), respectively. Both were sequenced with an Illumina HiSeq sequencer [44]. We annotated the assembled genomes using the MITOS WebServer [15] using the mitochondrial invertebrate genetic code for invertebrates. We calculated relative synonymous codon usage and amino acid frequencies using MEGA X [45]. The PCGs boundaries were con rmed manually in MEGA X software and BLASTp searches conducted at the NCBI website [46]. We identi ed secondary structures of the LNR using the Mfold web server [47]. The two rRNA genes, rrnL and rrnS, were curated by BLAST searches based on sequence similarity. The boundaries of rRNAs were assumed to be immediately after their upstream genes and before their downstream genes [18]. To predict tRNAs, we used tRNAscan-SE [13], ARWEN [14], MITOS, and MITOS2 [15]. The tRNAs that could not be predicted by programs were determined manually by alignment with their homologues in other astigmatid mites [11]. We also used the manual comparative approach to determine the results of tRNA search programs.

The largest non-coding region of mitogenomes ampli cation and sequencing
To verify the number of variable length AT-repeats in the LNR, we extracted DNA from 10 individual mite offspring from a single female using the DNeasy Blood and Tissue Kit (QIAGEN, Hilden, Germany). According to the obtained sequences, we designed speci c polymerase chain reaction (PCR) primers for each species (Table 5) to amplify the remaining genome by long PCR as a single fragment using the manufacturer's rapid PCR protocol.
We ampli ed the LNRs for G. fusca and L. destructor using four anking primers. The primer speci cations and the length of each ampli ed fragment are given in Table 5. We used prime STAR GXL DNA polymerase (Takara) for PCR under the following cycling conditions: 30 cycles of 98°C for 10 s, 52°C for 16 s, and 68°C for 1 min 40 s. We conducted PCR ampli cation reactions in 12.5-μl volumes containing 2.5 μl 5 × buffer, 1 μl dNTP, 0.25 μl GXL DNA Polymerase, 0.25 μl of each primer, 0.5 μl of template DNA, and double-distilled water to a nal volume of 12.5 μl. Finally, PCR products were separated by electrophoresis on a 1% agarose gel. PCR products were puri ed with a QIAquick Spin PCR Puri cation Kit (Qiagen). PCR products were cloned into a Blunting Kination Ligation Kit (BKL) (Takara). After heat-shock transformation of Escherichia coli (DH5a cells), we aligned their similarity against the sequence by using Clustal W 2.0 [48]. Table 5 The primers used for the ampli cation of Gohieria fusca and Lepidoglyphus destructor mitochondrial genomes, with the size of each fragment that was ampli ed in bp.

The reannotations of Histiostoma blomquisti, H. feroniarum, and Trouessartia Rubecula mitogenomes
The assembled genome was annotated using the MITOS WebServer [15]. The boundaries of PCGs were con rmed manually by MEGA X version 10.1.5 [45] software and BLASTp searches conducted in the NCBI database [46]. Based on highly conserved sequence motifs, rrnL and rrnS were identi ed by BLASTn searches of the NCBI database. To predict tRNAs, we used tRNAscan-SE [13], ARWEN [14], and MITOS 2. The tRNA genes that could not be identi ed by the tRNAscan-SE, ARWEN, and MITOS 2 were determined manually by alignment with their homologues in related species in the Astigmatid mites based on their anticodon and secondary structures [11]. If a tRNA can not be determined by automated prediction or manual comparative approach, such as predictions of the different tools were in some cases contradictory, the Mfold Server [47] and RNAeval (ViennaRNA package v.2.3.3) [49] were used for calculating the minimum free energy (MFE) to select the most probable among alternative structures (constrained analysis). The secondary structure with the smallest constrained MFE was considered as the most stable one [6]. We also used the manual comparative approach to determine the results of tRNA prediction programs.

Phylogenetic analysis
To infer the phylogenetic positions of L. destructor and G. fusca within the Astigmata, we generated a dataset of 19 mite taxa (17 astigmatid mites and two oribatid mites) ( Table 6) that included only those species with a complete set of 13 PCGs. First, the nucleotide sequences of each PCG was translated under the invertebrate mt genetic code. Then the amino acid sequences of 13 PCGs were aligned individually by MAFFT v7 [50]. Additionally, large gaps and ambiguous sites were deleted by Gblocks v.0.91b [51]. Finally, we conducted phylogenetic analyses using maximum likelihood (ML) and Bayesian inference (BI) methods. Table 6. Mitochondrial genomes employed in this study.
We analyzed the dataset of amino acid sequences of 13 PCGs as two types of matrices: combined algorithm and partitioned algorithm. Dataset partitioning was performed by PartitionFinder v.2.1.1 [52] based on an initial total of 13 data blocks (amino acid sequences, 13 protein-coding genes). We predicted the models by PartitionFinder using the Akaike information criterion (AICc). PartitionFinder used unlinked branch lengths, and the greedy search algorithm was used for amino acid sequences and the MrBayes model. For the dataset of partitioned amino acid sequences based on genes, two substitution models [GTR+I+G (COI, COII, ATP6, NAD3, NAD1, NAD2, NAD4, NAD4L, NAD5, NAD6, CYTB), GTR+G (ATP8, COIII)] were chosen by PartitionFinder. For the dataset of combined amino acid sequences, the substitution model GTR+I+G was chosen by PartitionFinder.
For the dataset of amino acid sequences, we performed BI analyses with MrBayes v.3.2.2 [53]. For MrBayes v.3.2.2, we used separate data partitions plus mixed models and conducted two independent runs each with four Markov chain Monte Carlo runs (one cold chain and three heated chains). The two datasets were run for 2 million generations, with trees sampled every 1000 generations. We then applied a conservative burn-in of 25%. We evaluated the convergence of the parameter estimates with TRACER v.1.6. All estimated parameters showed ESS values above 200. We edited the consensus tree with

Declarations
Ethics approval and consent to participate No speci c permits were required for the mites collected for this study in China. The mite specimens were collected from mushrooms, and the eld studies did not involve endangered or protected species. The species in our study are common mites and are not included in the "List of Protected Animals in China".

Consent for publication
Not applicable.

Competing interests
The authors declare that they have no competing interests.

Availability of data and materials
The datasets generated and/or analysed during the current study are not publicly available due [the mt genomes of Lepidoglyphus destructor and Gohieria fusca submitted on GenBank under the accession numbers MT075728 and MN608156, and not released yet] but are available from the corresponding author on reasonable request.   Note: Non-coding regions marked with a star are assumed to be putative control regions. The data for Trouessartia rubecula, Histiostoma blomquisti, and H. feroniarum in the table are after our revision. Table 3. The alignment of nucleotide sequences of four mitochondrial tRNA genes (trnS2, trnW, trnF, and trnI) in ve reported astigmatid species in four different superfamilies.