Mitochondrial genome evolution and tRNA truncation in Acariformes mites: new evidence from eriophyoid mites

The subclass Acari (mites and ticks) comprises two super-orders: Acariformes and Parasitiformes. Most species of the Parasitiformes known retained the ancestral pattern of mitochondrial (mt) gene arrangement of arthropods, and their mt tRNAs have the typical cloverleaf structure. All of the species of the Acariformes known, however, have rearranged mt genomes and truncated mt tRNAs. We sequenced the mt genomes of two species of Eriophyoidea: Phyllocoptes taishanensis and Epitrimerus sabinae. The mt genomes of P. taishanensis and E. sabinae are 13,475 bp and 13,531 bp, respectively, are circular and contain the 37 genes typical of animals; most mt tRNAs are highly truncated in both mites. On the other hand, these two eriophyoid mites have the least rearranged mt genomes seen in the Acariformes. Comparison between eriophyoid mites and other Aacariformes mites showed that: 1) the most recent common ancestor of Acariformes mites retained the ancestral pattern of mt gene arrangement of arthropods with slight modifications; 2) truncation of tRNAs for cysteine, phenylalanine and histidine occurred once in the most recent common ancestor of Acariformes mites whereas truncation of other tRNAs occurred multiple times; and 3) the placement of eriophyoid mites in the order Trombidiformes needs to be reviewed.


Results
General features of the mt genomes of the two eriophyoid mites, P. taishanensis and E. sabinae. The mt genomes of P. taishanensis and E. sabinae are 13, 475 bp and 13, 531 bp long respectively, are circular and have 37 genes: 13 protein-coding genes (PCG), two ribosomal RNA (rRNA) genes and 22 transfer RNA (tRNA) genes ( Fig. 1; Table 1). Genes are on both strands: one strand has 27 genes whereas the other strand has 10 genes. The start codons of the 13 PCGs were ATN in E. sabinae. In P. taishanensis, 11 of the 13 PCGs use ATN as start codons whereas cox1 and atp8 appear to start with CTG, which is a rare start codon for animal mitochondria (http://www.ncbi.nlm.nih.gov/Taxonomy/taxonomyhome.html/index.cgi ? chapter= cgencodes). We noticed that in the water mites, Unionicola parkeri and U. foili, another rare start codon, TTG, was used 34,35 . The stop codons were TAA or TAG in both eriophyoid species; incomplete stop codons, T, was found in protein-coding genes that precede a tRNA gene. In both species of eriophyoid mites, TAA was the most common stop codon and was used in nine of the 13 PCGs in P. taishanensis and 10 of the 13 PCGs in E. sabinae ( Table 1). The putative control region (CR) of P. taishanensis is short, only 47 bp in size, lying between rrnL and trnW. The putative CR of E. sabinae is 94 bp and is in a different location between trnK and atp6. No conserved regions were found between the CRs of the two eriophyoid mites. No other non-coding regions longer than 16 bp were found in the mt genomes of these two mites.
Eriophyoid mites have the least rearranged mt genomes among Acariformes mites. Like in other mites of the Acariformes, rearrangement of mt genes occurred in both of the eriophyoid mites. We calculated breakpoints with CREx as a measure of the extent of mt gene rearrangement from that of the hypothetical ancestor of arthropods 36 . Among the 12 different patterns of mt gene arrangement observed in the 27 species of Acariformes mites known, P. taishanensis has the least rearranged mt genome with 13 breakpoints and E. sabinae has the third least rearranged mt genome with 16 breakpoints (Table 2). In P. taishanensis, two rRNA genes (rrnS and rrnL) and seven tRNA genes (trnI, trnT, trnY, trnQ, trnV, trnW, trnM) changed their locations relative to their counterpart genes in the hypothetical ancestor of arthropods. In E. sabinae, two rRNA genes (rrnS and rrnL) and eight tRNA genes (trnK, trnI, trnT, trnY, trnQ, trnV, trnW, trnM) changed locations (Fig. 2). Four derived gene arrangements, trnS 1 -trnI-trnE, nad6-trnT-cob, trnY-trnQ-rrnS-trnV-rrnL, and trnW-nad2-trnM-trnC, were Protein-coding genes are color-coded (cox: blue; nad: green; atp: orange; cob: yellow); rRNA genes are in grey; control regions are in black; tRNA genes are in red or purple. Abbreviations of protein-coding genes are: atp6 and atp8 for ATP synthase subunits 6 and 8, cox1-3 for cytochrome oxidase subunits 1-3, cob for cytochrome b, nad1-6 and nad4L for NADH dehydrogenase subunits 1-6 and 4 L, rrnL and rrnS for large and small rRNA subunits. tRNA genes are indicated by the single letter IUPAC-IUB abbreviations for their corresponding amino acids. Arrows and arrowheads show the direction of gene transcription. Numbers at gene junctions indicate the length of non-coding sequences; negative numbers indicate overlap between genes.
Scientific RepoRts | 6:18920 | DOI: 10.1038/srep18920 found in both eriophyoid mites but not in any other mites; these novel gene arrangements were candidate synapomorphies (i.e. shared derived characters) for the eriophyoid mites.
Truncated tRNAs of the eriophyoid mites. Fifteen mt tRNA genes of P. taishanensis and 16 mt tRNA genes of E. sabinae were identified with tRNAscan-SE 37 or ARWEN 38 programs. The other seven tRNA genes of P. taishanensis (trnA, trnG, trnQ, trnR, trnS 1 , trnT, trnV) and the other six tRNA genes of E. sabinae (trnA, trnG, trnI, trnQ, trnR, trnS 1 ) were found manually based on conserved nucleotides and the anticodon sequences. The putative mt tRNAs were highly truncated in both P. taishanensis (47 to 61 bp) and E. sabinae (47 to 67 bp) (Figs 3 and 4). Sixteen of the 22 tRNAs have atypical secondary structures, missing either D-arm or T-arm in both mites. The majority of tRNAs also have mismatches on T-arm, D-arm, acceptor arm or anticodon arm.

Discussion
The ancestral pattern of mt gene arrangement of Acariformes mites. Mitochondrial gene arrangement is relatively stable within major animal lineages 8 . This is, however, not the case for the Acari (ticks and mites). For the superorder Parasitiformes, the ancestral pattern of mt gene arrangement to arthropods was found in 11 species of Argasidae, one species of Allothyridae and six Ixodes species of Ixodidae (Fig. 5). Four Mesostigmata species and 16 Ixodidae species have rearrangement of mt genes (Fig. 5, Table 2). In the other superorder Acariformes, however, rearrangement of mt genes was found in all of the 27 species whose mt genomes have been sequenced. Even the species in the same genus differ from one other in mt gene arrangement, such as the two water mites in the genus Unionicola 34,35 and three chigger mites in the genus Leptotrombidium 11,42 (Fig. 5). The two eriophyoid mites investigated in the current study were from the same tribe (Phyllocoptini) of the family Eriophyidae but differ from one another in mt gene arrangement. Scientific RepoRts | 6:18920 | DOI: 10.1038/srep18920 What was the ancestral pattern of mt gene arrangement for Acariformes mites ? Establishment of the ancestral pattern will lay the ground for understanding how mt genome organization evolved in Acariformes mites. Intriguingly, in comparison to other species of the Acariformes, the two eriophyoid mites have the least rearranged mt genomes and retained most of the ancestral pattern of mt gene arrangement of arthropods (Figs 2 and 5). The conserved mt gene-arrangement characters present in the two eriophyoid mites and other Acariformes mites allowed us to infer that the ancestral pattern of mt gene arrangement of the Acariformes retained likely the ancestral pattern of mt gene arrangement of arthropods with slight modifications (Fig. 2). The modifications include: 1) translocations of trnQ and trnY; these two tRNA genes are rearranged in all of the 27 species of Acariformes mites investigated to date and furthermore, their locations vary among different Acariformes lineages (Figs 2 and 5) inversion of rrnS-trnV-rrnL as a cluster, which is seen in the two eriophyoid mites and five Sarcoptiformes mites (Fig. 5); however, in seven of the 27 Acariformes species, rrnS and rrnL are translocated but not inverted, so there are alternative interpretations for the location and transcription orientation of rrnS-trnV-rrnL in the ancestor of Acariformes mites. tRNA truncation in Acariformes mites. Truncated mt tRNAs have been found in several orders of the class Arachnida including Araneae 25-27 , Acariformes 11 , Pseudoscorpiones 28 , Scorpiones 29,30 , and Thelyphonida 30 . In the superorder Acariformes, truncated tRNAs are common and have been observed in all of the 27 species whose mt genomes have been sequenced. Sixteen of the 22 tRNAs in the two eriophyoid mites investigated in the current study lack D-arm or T-arm (Figs 3 and 4). Nineteen tRNAs lack D-arm or T-arm in spider mites (Tetranychidae); in some species, the tRNAs for phenylalanine and glutamine lack both D-arm and T-arm (Table S1) 43,44 . Fifteen tRNAs lost D-arm or T-arm in Demodicidae, including five tRNAs losing both D-arm and T-arm 45 (Table S1).    Parasitiformes  --99  100  98  100  100  100  100  100  100  100  100   Acariformes  100  100  100  100  100  100  100  100  100  100  100  100  100   Sarcoptiformes  100  100  97  100  96  100  99  100  100  100  100  100  100   "Trombidiformes"  -58  100  100  100  100  100  100  100  100  100  100  100   Mesostigmata  41  -99  100  98  100  100  100  100  100  100  100  100   Ixodida  66  100  95  100  90  100  99  100  87  100  100  100  100   Acaridae  100  100  100  100  100  100  100  100  100  100  100  100  100   Argasidae  100  100  100  100  100  100  100  100  100  100  100  100  100   Ixodidae  100  100  100  100  100  100  100  100  100  100  100  100  100   Eriophyidae  100  100  100  100  100  100  100  100  100  100  100  100  100   Demodicidae  100  100  100  100  100  100  100  100  100  100  100  100  100   Phytoseiidae  100  100  100  100  100  100  100  100  100  100  100  100  100   Tetranychidae  100  100  100  100  100  100  100  100  100  100  100  100  100   Trombiculidae  100  100  100  100  100  100  100  100  100  100  100  100  100   Unionicolidae  100  100  100  100  100  100  100  100  100  100  100  100  100   Table 3 Large-scale tRNA truncation was first observed, and best studied in nematodes 24 . Twenty of the 22 tRNAs lack T-arm and the two tRNAs for serine lack D-arm in all nematodes except for several species in the class Enoplea. In Trichinella spiralis, eight tRNAs have the cloverleaf secondary structure with both D-and T-arms whereas several others lack both D-and T-arms 51 . tRNA truncation in the Acariformes mites has some similarity to that in nematodes but has their distinctive features. First, trnK has the typical cloverleaf structure in all Acariformes mites except for Steganacarus magnus 9 (but K was not identified in this mite, see Table S1), whereas the other 21 tRNAs lack either D-arm or T-arm or both arms in one or more species of Acariformes mites. Three tRNAs, for cysteine, phenylalanine and histidine respectively, lack T-arm in all known Acariformes mites (Table S1), indicating that T-arm loss in these three tRNAs is likely ancestral to Acariformes mites. The other 18 tRNAs, however, vary in their secondary structures among the Acariformes mites: either having a cloverleaf structure, or lacking D-arm, T-arm or both arms (Table S1). The pattern of D-arm or T-arm loss appears to be consistent within a family but differs between families (Table S1). Thus, the loss of D-arm or T-arm in these tRNAs may have occurred multiple times independently in different lineages of the Acariformes mites.
Limited evidence from nematodes indicates that truncation does not seem to disable mt tRNAs from function. Okimoto et al. hybridized mt tRNA gene-specific probes to the RNAs of two nematodes, Caenorhabditis elegans and Ascaris suum, and obtained evidence for the transcription of at least nine C. elegans and three A. suum mt tRNA genes 52 . Each tRNA transcript has the same size of its corresponding tRNA gene, to which three nucleotides (CCA) were added after transcription. Okimoto et al. concluded that the mt tRNAs without D-arm or T-arm were functional. Further, Suematsu et al. and Ohtsuki and Watanabe showed that the evolution of truncated mt tRNAs in nematodes was linked to EF-Tu protein: the paralogs of this protein acquired differential binding abilities to tRNAs with deleted domains 53,54 . Recently, Juhling et al. showed computational evidence that tRNAs of the nematodes of the class Enoplea that lack both D-arm and T-arm were functional 51 . The Acariformes mites provided another system for further experimental and computational investigation into the evolution and function of truncated mt tRNAs.
Should Eriophyoidea be placed in the order Trombidiformes?. The subclass Acari comprises two superorders, Acariformes and Parasitiformes 55 . The former includes two orders (Trombidiformes and Sarcoptiformes), while the later include four orders (Opilioacarida, Holothyrida, Ixodida, Mesostigmata) 55 . The monophyly of Acari, however, is not without controversy, when species of other major lineages of Arachnida were included into analyses 28,40,41,56 . Pseudoscorpiones was always grouped with Acariformes in our analyses, which is consistent with Ovchinnikov and Masta 28 . Several phylogenetic studies of Acariformes [39][40][41] and Parasitiformes 57 have been conducted previously using molecular data; however, eriophyoid mites were not included. The current taxonomic assignment of Eriophyoidea to the order Trombidiformes 55 is controversial, due to the distinctive morphology of eriophyoid mites such as having only two pairs of legs, lack of ontogenetic diversity and absence of respirator system 58 . Indeed, André placed Eriophyoidea outside Trombidiformes 59 Our analyses based on different partitions and inference methods showed consistently that Acariformes, Parasitiformes and several families (Acaridae, Argasidae, Demodicidae, Ixodidae, Phytoseiidae, Tetranychidae, Trombiculidae and Unionicolidae) were monophyletic, which was consistent with previous studies using mt genes 9,18,27,39,[44][45][46]50 or nuclear genes 18,39,57 . The monophyly of Trombidiformes, however, was rejected in our analyses ( Figure S1, Figure  S2). The monophyly of Trombidifromes was also rejected in a recent study by Gu et al. based on the mt genome sequences of 16 species of Acariformes mites due to spider mites (Tetranychidae) grouped with sarcoptiform mites 50 . In our analyses, when the two species of eriophyoid mites were excluded, the rest of the Trombidiformes species were always together in a monophyletic group. Our results thus raised the need for further phylogenetics studies on Acariformes mites with more taxa included especially from Tydeoidea and Eupodoidea, which were thought to be the sister groups of Eriophyoidea 60 . Further morphological studies are also necessary to elucidate the position of eriophyoid mites in Acariformes.
In conclusion, we sequenced the mt genomes of two species of eriophyoid mites, P. taishanensis and E. sabinae. These two mites have the least rearranged mt genomes seen in the Acariformes and have highly truncated mt tRNAs. Our comparison between the eriophyoid mites and other mites and ticks showed that the most recent common ancestor of Acariformes mites retained the ancestral pattern of mt gene arrangement of arthropods with slight modifications. The truncation of three tRNAs (for cysteine, phenylalanine and histidine, respectively) likely occurred once in the common ancestor of Acariformes mites. Truncation of other tRNAs, however, occurred multiple times independently in different lineages of Acariformes mites. Our phylogenetic analyses of mt genome sequences rejected the monophyly of the order Trombidiformes when eriophyoid mites were included. Further phylogenetics studies on Acariformes mites including more taxa from different lineages is needed to clarify the position of eriophyoid mites.

Methods
Collection of mites. E. sabinae and P. taishanensis were collected in May 2013 in Nanjing, China. E. sabinae was collected from Juniperus chinensis (Cupressaceae) (China savin), while P. taishanensis from Cedrus deodara (Pinaceae) (Deodar cedar). Mite samples were either used immediately for DNA extraction or were preserved in 100% ethanol at − 20 °C prior to DNA extraction. Samples of each eriophyoid species were also mounted to slides as voucher, using modified Berlese medium 61 for morphological check with Zeiss A2 (microphoto camera AxioCam MRc) microscope. All of the specimens and vouchers were deposited at the Arthropod Collection, Department of Entomology, Nanjing Agricultural University, China. DNA extraction, mt genome amplification and sequencing. Genomic DNA was extracted from both individual and pooled specimens for each species, using a DNeasy Blood and Tissue Kit (QIAGEN), following the modified protocol 62 . For E. sabinae, a 658-bp fragment of cox1 and a 413-bp fragment of rrnL were initially Scientific RepoRts | 6:18920 | DOI: 10.1038/srep18920 amplified by PCR with the primer pairs LCO1490-HCO2198 63 and LRJ12287-LRN13398 64 (see Additional file Table S2). PCR products were purified and sequenced directly using Sanger method at Majorbio (Shanghai, China). Specific primers for E. sabinae, ECOISR3 and E16SSR2, were designed from the sequences of the cox1 and rrnL fragments, respectively. PCR with these two primers produced a 1.7-kb amplicon, which was sequenced using Sanger method at Majorbio. Another pair of primers, PF2-PR2, were designed from the sequence of the 1.7-kb amplicon. An 11.2-kb amplicon was produced with PF2-PR2 primer pair and was sequenced with Illumina Hiseq 2000 platform at the Beijing Genome Institute, Hong Kong (BGI-HK).
For P. taishanensis, a 658-bp fragment of cox1 and a 407-bp fragment of rrnL were initially amplified by PCR with the primer pairs LCO1490-HCO2198 63 and LRJ12287-LRN13398 64 (see Additional file Table S2). The PCR products were purified and sequenced directly using Sanger method at Majorbio. Two pairs of specific primers, TYCB3R2-TY16sR1 and RTYCB3F5-RTY16sF4, were designed from the sequences of the cox1 and rrnL fragments. The PCR with TYCB3R2-TY16sR1 produced a 1.7-kb amplicon, which was sequenced using Sanger method at Majorbio. The PCR with RTYCB3F5-RTY16sF4 produced an 11.2-kb amplicon, which was sequenced with Illumina Hiseq 2000 platform at the BGI-HK.
The initial PCRs contained 12.5 μ L of PCR SuperMix (Transgen Biotech Co., Ltd., Beijing, China), 2 μ l of template DNA, and 1.25 μ M of each primer, in a total volume of 25 μ L. The PCR cycling conditions were: 3-min denaturation at 96 °C; 35 cycles of 10-sec denaturation at 95 °C, 30-sec annealing at 46 °C and 1.5-min extension at 72 °C; 5-min final extension at 72 °C; and then held at 4 °C. PCR products were checked on 1% agarose gel. PrimeSTAR GXL DNA polymerase (TAKARA) was used in the long PCRs with the cycling conditions: 98 °C for 10 sec, 68 °C for 2 to 10 min (depends on the length of regions between rrnL and cox1). The reaction mixture contained 0.5 μ l GXL DNA Polymerase, 5 μ l buffer, 2 μ l dNTP mixture, 0.75 μ l of each primer, 1 μ l of template DNA and Milli-Q water added to total volume of 25 μ l. Positive and negative controls were executed with each PCR. PCR products were checked on 1% agarose gel. PCR products were purified with QIAquick Spin PCR Purification Kit (QIAGEN).

Assembly of Illumina sequence-reads, gene identification and gene rearrangement analysis.
Illumina sequence-reads obtained from the mt genome amplicons of P. taishanensis and E. sabinae were assembled into contigs with Geneious 6.1.6 (Biomatters Ltd.). The transfer RNA (tRNA) genes were identified using tRNAscan-SE 37 and ARWEN 38 or identified manually based on anticodons and secondary structures. tRNA genes of the two eriophyoid mites were verified by comparison of secondary structures and conserved nucleotide sequences with those of the Acari species reported in published literature. PCGs were identified by open reading frame search in Geneious and BLAST searches of GenBank 65 . The two rRNA genes, rrnL and rrnS, were also identified by BLAST searches of GenBank based on sequence similarity and conserved sequence motifs. The start and stop nucleotides of rrnL and rrnS cannot be determined exactly and were assumed to be immediately after their upstream genes and before their downstream genes. Breakpoints were calculated with CREx 36 web server (http://pacosy.informatik.uni-leipzig.de/crex) as a measure of the extent of mt gene rearrangement; the two eriophyoid mites were compared with the hypothetical ancestor of arthropods (see Table 3 for details). The nucleotide sequences of mt genomes of P. taishanensis and E. sabinae have been deposited in GenBank under accession numbers KR604967 and KR604966.
We analyzed the datasets of mt genome sequences as two types of matrix: amino acid sequences of PCGs and nucleotide sequences of all genes (13 PCGs, 2 rRNA genes and 22 tRNA genes). The datasets were partitioned by genes, by codon positions and optimal partitioning as determined by PartitionFinder (Table S4). Amino acid sequences were partitioned by 13 genes. Nucleotide sequences partitioned by genes resulted in 16 datasets: 13 PCGs, 2 rRNA genes and concatenated tRNA genes. Partitioning by codon positions resulted in six datasets: one for each base of codons, two for each rRNA gene and one for the combined tRNAs. Both types of partitioned nucleotide sequences were run twice independently, once with the third codon positions of PCGs included and once without the third codon positions. To avoid redundant sampling of each genus, which may potentially affect a complex lower-level phylogeny, two representative species from each genus (Argas, Ornithodoros, Amblyomma, Haemaphysalis, Ixodes, Rhipicephalus, Tetranychus, Leptotrombidium) were included in the analyses. The reduced datasets included 43 taxa and only Baysian inference was run based on the best models found by PartitionFinder. To check if RNA genes could affect the topology, we also ran Baysian inference with the nucleotide sequences of the 13 PCGs based on the best models found by PartitionFinder. The datasets with 18 additional Arachnid mt genomes were ran in ML and Baysian analyses with nucleotide sequences of the 13 PCGs by 13 partitions.
The best models of the datasets were predicted by jModelTest 2.1.1 75 and PartitionFinder v1.1.1 76 , using the Bayesian Information Criterion (BIC). PartitionFinder was set using unlinked branch lengths, greedy search for nucleotide sequences and amino acid sequences. MtRev + I + G + F was chosen as the best amino acid substitution Scientific RepoRts | 6:18920 | DOI: 10.1038/srep18920 model for most PCGs, except for nad4, for which JTT + I + G + F was chosen as the best model. Nucleotide substitution model GTR + I + G was chosen as the best of 16 partitions. The ML analyses were performed using the GTRGAMMA model for nucleotide partitions and MtRev + I + G + F (JTT + I + G + F for nad4) for amino acid partitions in the program raxmlGUI1.3 77,78 . For nodal support evaluation, a nonparametric bootstrap with 1,000 replicates was used. Mixed-model Bayesian analyses were performed with MrBayes 3.2.2 79 using separate data partitions. For BI, one cold chain and three heated chains were run with the combined dataset for 2 million generations. The average standard deviation of split frequencies fell down quickly. After 0.2 million generations, the average standard deviation number was below 0.01 in most of the BI trees. The convergence of the parameter estimates was performed with Tracer v1.6. The consensus three was edited with FigTree1.4.0.