Sampling and DNA extraction

Samples of C. malaccensis were collected from Haikou, China, then cultured in Academy of National Food and Strategic Reserves Administration in China. The samples were identified based on their morphology characteristics (Shen 1997). Genomic DNA was extracted from 50 adult mites using TIANamp Micro DNA Kit following the instruction. The UV-Vis Spectrophotometer Q5000 (Q5000, Quawell Technology, Inc. USA) was used to quantify the extracted DNA (23.76ug) for library construction.

DNA Sequencing and mt genome assembly

Total DNA of C. malaccensis was sent to Berry Genomics Company for library preparation, and finally sequenced on an Illumina sequencer. The insert size was 250bp with 150bp pair-end sequencing.

After getting the raw data, quality assessment was conducted using FastQCv0.11.7 ( http://www.bioinformatics.babraham.ac.uk/projects/fastqc). Trimmomatic v0.36 (Lohse et al. 2012) was used to the remove adapter sequences. The mt genome assembly strategy followed Feng et al. 2018. Briely, to reconstruct the mt genome, cox1 gene fragments were chosen as “anchors”. We sequenced the cox1 gene fragments using universal primer pairs LCO1490–HCO2198 (Folmer et al. 1994). We applied “map to reference” strategy and mapped all cleaned reads to the “anchor” using Geneious (Kearse et al. 2012). Illumina sequence-reads were assembled using the sequence cox1 as initial references. The contigs were then extended using the assembly parameters: (1) minimum overlap 50 bp, and (2) minimum similarity 99%, until the full circular mt chromosome sequences were obtained.

The PCR amplification of cox1 gene was in 25 ul containing 12.5ul of 2×Taq Mix (Tiangen, Beijing, China), 1 ul of each primer (10 uM), 1ul genomic DNA and 9.5ul ddH₂O. PCR cycling conditions were: 95°for 3 min, followed by 35 cycles of 94°for 1 min, 53°for 1 min, 72°for 1 min, and finally 72°for 10 min. Bioinformatic analyses

The protein coding genes (PCGs), ribosome RNA (rRNA) genes and transfer RNA (tRNA) genes were identified using MITOS Web Server ( http://mitos.bioinf.uni-leipzig.de/index.py) (Bent et al., 2013). PCGs and RNA genes were confirmed by alignment against close species while tRNA genes were confirmed by ARWEN. During tRNA search, we alternatively searched tRNA genes based on the methods by Xue et al. (2018). The base composition was analyzed with muscle algothrim in MEGA 7.0 (Kumar et al., 2016). Nucleotide compositional skew was measured using the following formula: AT skew = A-T/A+T and GC skew = G-C/G+C.

Phylogenetic analysis

To better understand molecular phylogeny of C. malaccensis in Arachnida, a total of 30 Arachnida species were used in phylogenetic analysis, including 29 representative stored mite species (Table S1). All 13 PCGs were aligned in MEGA 7.0. Two rRNA genes were aligned using MAFFT v7.0 online sever (Katoh et al. 2019). Then ambiguous positions in the alignment of PCGs and rRNA genes were removed by using Gblocks v0.91b web server (Castresana 2000). The concatenated phylogenetic trees were reconstructed based on the two datasets: 1) dataset “PCGsrRNA”, in which there are concatenated 13 protein-coding genes (atp6, atp8, cox1, cox2, cox3, cob, nad1, nad2, nad3, nad4, nad4l, nad5, nad6) and two rRNA genes (rrnL and rrnS); 2) dataset “PCGs”, including only 13 protein-coding genes. The Bayesian methods (BI) and maximum likehood (ML) were used to reconstruct the phylogenetic tree. The ML method was performed with PhyML (Guindon et al. 2010) ( http:/www.atgc-montpellier.fr/phyml/) and the online Smart Model Selection for optimal model selection. General time reversible (GTR) model was finally chosen. We also applied Bayesian method using MrBayes 3.2.2 (Ronquist & Huelsenbeck 2003). The GTR+I+G model was used. The datasets were conducted with two simultaneous runs of 2 million generations, each with one cold and three heated chains. Samples were drawn every 1,000 Markov chain Monte Carlo (MCMC) steps, with the first 25% discarded as burn-in. The stationarity was considered reached and stopped run when the average standard deviation of split frequencies was below 0.01.

Mt genome organization and nucleotide composition

The length of the mt genome of C. malaccensis is 14,732bp (Figure 1), including 13 PCGs, two rRNAs and 22 tRNAs (Table 1). The symmetric nucleotide compositions are A:47.3%, C:14.1%, G:6.8% and T:31.8%. The content of G is small and concentrated in cox1 (10.5%) and cox2 (8.7%) genes, which shows that these two genes are more conservative and stable. The AT skew is 0.196, which is similar with most AT skews of the arthropod animals (Dermauw et al. 2009; Sun et al. 2014a). Four PCG genes (nad5, nad4, nad4l, nad1), two rRNA genes and seven tRNA genes (trnL1, F, P, Q, L2, Y, C) of C. malaccensis were encoded by the minor strand and the other were encoded by the major strand. Two tRNA genes were contained by protein coding gene, trnH gene located in nad5 and trnV was in 16S rRNA. In addition, trnA was in control region. Compared with the Arthropod ancestral pattern (Palopoli et al. 2014), the mt genome of C. malaccensis has different arrangements. Rearrangement was found in 12 genes, including 11 tRNA genes (trnE, S1, N, M, Y, I, L2, Q, R, V, L1) and 12S rRNA gene. Except trnR and trnV were in different location, C. malaccensis has the similar gene rearrangement with Demodex brevis and D. folliculorum (Palopoli et al. 2014), belonging to Cheyletoidea. To fully understand the mitochondrial arrangement of Cheyletoidea, and analyze its characteristics, more species need to be tested.

TABLE 1

Annotation of the mitochondrial genome of Cheyletus malaccensis.

FIGURE 1.

The mitochondrial genome arrangements of Cheyletus malaccensis.

Note: the blue is PCGs, the red is rRNA genes, the purple is tRNA genes.

The total length of intergenic sequences in C. malaccensis was 548 bp, including two intergenic sequences longer than 30 bp, 366 bp (between nad4 gene and nad4l gene) and 47 bp (between nad4 gene and nad5 gene). The total length of C. malaccensis overlapping regions was 195 bp with three overlapping regions longer than 30 bp: 53 bp (between trnF gene and nad5 gene), 43 bp (between trnQ gene and trnV gene) and 37 bp (between trnS2 gene and cob gene), respectively.

Protein coding genes

The mt genome of C. malaccensis has three kinds of start codon (ATT, ATA, ATG) and two kinds of stop codon (TAA, TAG). ATT is the start codon of cox1 and nad5, ATA is the start codon of atp8, atp6, nad3, nad4l, nad1 and nad2. ATG is the start codon for other genes. TAG is the stop codon of atp8, nad5 and nad4, while the other genes end with TAA.

The tRNA and rRNA genes

Twenty-two tRNA genes were found in the mt genome of C. malaccensis. All present tRNA genes were extremely truncated, ranging in size from 47 bp to 67 bp. Four tRNA (trnM, W, N and K) have complete cloverleaf secondary structures, most of predicted tRNA genes are folded into the atypical cloverleaf secondary structures with the absence of either D-arm or TψC-arm (Figure 2). As reported in some species in Arachnida, which also found the tRNA lost one arm either the D-arm or TψC-arm such as Panonychus citri (Yuan et al. 2010), D. brevis and D. folliculorum (Palopoli et al. 2014), Caloglyphus berlesei (Sun et al. 2014b) and Tyrophagus longior (Yang & Li 2015). It seems that lose arms is a normal situation in Arachnida and it does not affect mitochondrial protein synthesis.

FIGURE 2.

The 22 tRNAs of Cheyletus malaccensis

For C. malaccensis, rearrangements of tRNA genes occur much more than other genes. The truncated tRNA genes observed in C. malaccensis would seem to require the evolution of extensive tRNA editing capabilities. The molecular machinery necessary for these unusual tRNAs to function might provide an explanation for C. malaccensis adaptability and becoming the dominant predator mite in grain depot.

Phylogeny

Leveraging two datasets (PCG123rRNA and PCG) as well as two methods (BI and ML), we reconstructed the phylogeny of Acariformes which had been reported on NCBI (Figure 3). We uncovered for the first time the mt genome belonging to Cheyletidae, which provided a good opportunity to resolve its phylogenetic position. Cheyletidae had the closest relationship to Tetranychidae (PCGsrRNA Bayesian posterior probability / PCGs Bayesian posterior probability / PCGsrRNA ML bootstrap value/PCGs ML bootstrap value = 1/1/89/87) and they formed a clade with sister relationship with Demodicidae (support values = 1/1/96/94). The phylogenetic relationship among the three families would be clear with the implement of our mt genome data. Phylogenetic relationships among these families could be more clear when supplementing enough molecular information of Acariformes.

FIGURE 3.

Phylogenetic tree of 30 species from Acariformes.

Note: The values in each node represented posterior probabilities of Bayesian methods and bootstrap value from maximum likehood method: PCGsrRNA Bayesian/PCGs Bayesian/PCGsrRNA ML/PCGs ML.