Skip to main content
  • Research article
  • Open access
  • Published:

Complete mitochondrial genomes of the human follicle mites Demodex brevis and D. folliculorum: novel gene arrangement, truncated tRNA genes, and ancient divergence between species

Abstract

Background

Follicle mites of the genus Demodex are found on a wide diversity of mammals, including humans; surprisingly little is known, however, about the evolution of this association. Additional sequence information promises to facilitate studies of Demodex variation within and between host species. Here we report the complete mitochondrial genome sequences of two species of Demodex known to live on humans—Demodex brevis and D. folliculorum—which are the first such genomes available for any member of the genus. We analyzed these sequences to gain insight into the evolution of mitochondrial genomes within the Acariformes. We also used relaxed molecular clock analyses, based on alignments of mitochondrial proteins, to estimate the time of divergence between these two species.

Results

Both Demodex genomes shared a novel gene order that differs substantially from the ancestral chelicerate pattern, with transfer RNA (tRNA) genes apparently having moved much more often than other genes. Mitochondrial tRNA genes of both species were unusually short, with most of them unable to encode tRNAs that could fold into the canonical cloverleaf structure; indeed, several examples lacked both D- and T-arms. Finally, the high level of sequence divergence observed between these species suggests that these two lineages last shared a common ancestor no more recently than about 87 mya.

Conclusions

Among Acariformes, rearrangements involving tRNA genes tend to occur much more often than those involving other genes. The truncated tRNA genes observed in both Demodex species would seem to require the evolution of extensive tRNA editing capabilities and/or coevolved interacting factors. The molecular machinery necessary for these unusual tRNAs to function might provide an avenue for developing treatments of skin disorders caused by Demodex. The deep divergence time estimated between these two species sets a lower bound on the time that Demodex have been coevolving with their mammalian hosts, and supports the hypothesis that there was an early split within the genus Demodex into species that dwell in different skin microhabitats.

Background

Mites of the genus Demodex live in the hair follicles and sebaceous glands of mammalian skin [1]. They are extremely widespread among mammalian lineages, with species having been described from hosts in three of the seven marsupial orders, and in 11 of the 18 eutherian orders [2]. Their bodies exhibit specializations that make them well adapted to inhabiting the constricted spaces of the pilosebaceous complex—they are typically just 100–300 μM long and cylindrical in shape, with extremely reduced legs and setation (Figure 1). Together, these observations suggest an ancient, coevolutionary relationship between mammals and the Demodex species that inhabit their skin; yet surprisingly little is known about the evolutionary history or dynamics of this association.

Figure 1
figure 1

Organization of the mitochondrial genome of Demodex folliculorum . Protein-coding, ribosomal RNA, and transfer RNA genes are depicted as green, blue, and red arrows, respectively. The arrows represent the direction of transcription. The AT-rich, putative control region is depicted as a gray box (labeled “D-Loop”). Inside the circle is a scanning electron micrograph of an individual D. folliculorum (image credit: Power and Syred). The identical gene arrangement was observed in the D. brevis mitochondrial genome. This is a novel gene arrangement among the acariform mitochondrial genomes sequenced to date. Like other members of the Acariformes, Demodex proved to have a compact and AT-rich mitochondrial genome.

Demodex are ubiquitous in humans and thought to contribute to medically important skin disorders. Both recent molecular studies [3], and studies in which the skin of cadavers was sampled extensively (reviewed in [2]), suggest that the incidence of infestation approaches 100%. Despite the high incidence of Demodex itself, the frequency with which Demodex causes pathogenicity appears to be very low, which has led to their characterization as human commensals by many authors. Nevertheless, there are reports of high Demodex densities associated with two medically important disorders of the skin—marginal blepharitis [4] and acne rosacea [5]. Furthermore, mite densities are known to increase and cause skin diseases when the immune system is compromised [6]. Treatment with acaricides, such as ivermectin, can be curative for some of these diseases [7]. These results suggest that Demodex should be considered parasitic in some host individuals. Future studies examining the potential roles of Demodex in skin disorders, as well as interactions between Demodex and the host immune system, would be aided by an increase in the genetic markers available for distinguishing mite populations.

More generally, greater genetic information would facilitate the study of Demodex variation within and between host species. For example, mite genetic variation could be investigated within a host species in order to better understand the dynamics of movement among host individuals [8], or even to provide information about the migration patterns of the hosts themselves [9]. Similarly, the question of how often these mites migrate across host species boundaries could be addressed by making comparisons of host and mite molecular phylogenies [10].

Two species of Demodex inhabit the skin of humans, with histological studies suggesting that each occupies a different niche: Demodex folliculorum typically resides in the hair follicle nearer the skin surface, whereas D. brevis is generally found deeper in the sebaceous glands [1]. Although distinguished initially based on morphological differences, recent molecular work using both mitochondrial and nuclear gene markers has verified that these are distinct species [11, 12]. Indeed, the observed levels of sequence divergence suggest that they are not closely related species, although no estimate of divergence time is currently available.

Complete mitochondrial genomes have been determined for 685 arthropods; just 15 of these sequences are from the Acariformes, however, despite the more than 41,000 described species, and perhaps many times more undescribed species, within this division of the Acari [13]. The paucity of mitochondrial genomes available from among the Acariformes is especially unfortunate given their tendency to evolve unusual gene arrangements and truncated examples of tRNA genes [1419]. To improve our understanding of the evolution of these features, additional examples of mitochondrial genomes from this important group of organisms are needed.

Here we report the complete mitochondrial genomes of both D. brevis and D. folliculorum, the first available for any species within the genus Demodex. Both species shared a novel mitochondrial gene order that has diverged extensively from the ancestral chelicerate pattern, with tRNA genes apparently having moved much more often than other genes. Furthermore, the mitochondrial genomes of both Demodex species possessed extremely truncated tRNA genes, with several of these lacking the sequence necessary to code for both D- and T-arms. Finally, we estimated the divergence time between the two human-associated Demodex species by using the sequences of mitochondrial proteins to conduct relaxed molecular clock analyses among Acariformes.

Results and Discussion

Mitochondrial genome content and organization

Like other members of the Acariformes, both species of Demodex had compact and AT-rich mitochondrial genomes. The D. brevis genome [GenBank accession number: KM114225] length was 14,211 bp and AT-percentage was 69.0%; the D. folliculorum genome [Genbank accession number: KM114226] length was 14,150 bp and AT-percentage was 71.1%. For comparison, the 15 other species of Acariformes with complete mitochondrial genomes available in GenBank exhibited an average genome size of 14,240 bp (SE 281 bp) and an average AT-percentage of 74.9% (SE 1.7%). These values agree closely with those that we observed for both Demodex species, suggesting that the same combinations of mutation pressures and selective forces that determine these parameters for other Acariformes are also operating in the Demodex lineage.

Most metazoans have mitochondrial genomes that are at least several hundred base pairs larger than those of the Acariformes [20]. Having a small mitochondrial genome would seem to have the obvious benefit of being less costly and time-consuming to replicate. Nevertheless, it remains unclear why selection for reduced genome size would be especially effective in the Acariformes. One of the contributing factors, however, is the extreme truncation of mitochondrial tRNA genes that is widespread in these organisms (see below).

High AT content is a common feature of metazoan mitochondrial genomes [21]. It has been hypothesized that this bias results from being housed where there are high concentrations of reactive oxygen species that promote GC to AT mutations [22]. Alternatively, AT-richness could be an adaptation for metabolic efficiency, with the higher energy cost and more limited availability of G and C nucleotides driving this pattern [23].

Mitochondrial genomes of both Demodex species possessed the standard collection of 13 protein-coding, two ribosomal RNA (rRNA), and 22 tRNA genes (Figure 1) present in the mitochondrial genomes of most metazoans [24]. The arrangement of these genes was identical in the two Demodex species, but novel among Acariformes. The tendency to accumulate rearrangements of the mitochondrial genome appears to be a common feature among the Acariformes [1419]. Whether the accumulation of rearrangements in this taxonomic group is due to a higher structural mutation rate, or to relaxed natural selection on gene order, remains an open question.

The ancestral chelicerate gene arrangement, one shared by members of at least five orders, is exemplified by Limulus polyphemus[25, 26]. If we consider just the protein-coding and rRNA genes, the gene arrangement shared by both Demodex species can be derived from the ancestral chelicerate gene arrangement via a single block interchange (Figure 2), requiring three breakpoints (one breakpoint denotes when the linear sequence of genes is broken and rearranged such that what were two adjacent genes in the original genome no longer appear consecutively in the rearranged genome). Mechanistically, this block interchange could have resulted from a tandem duplication of the entire interval, followed by differential loss of functional genes in the two tandem copies. Alternatively, it could have resulted from direct translocation, either of the 12S gene or of the interval that stretches from ND5 to the 16S gene.

Figure 2
figure 2

Evolution of the Demodex gene arrangement from the ancestral chelicerate. Circular mitochondrial genomes are depicted as linearized, starting at the 5′ end of the COX1 gene. The one letter amino acid code is used to designate the tRNA genes, with the exceptions that L1 = CUN; L2 = UUR; S1 = AGN; and S2 = UCN. The putative control region is designated as A + T to indicate that this has a higher AT% than the rest of the genome. Genes above the medial line are encoded on one strand, while those below the line are encoded on the other strand. As shown, the arrangement of protein-coding plus rRNA genes can be derived via one block interchange, with a minimum of three breakpoints. This could have occurred via a tandem duplication followed by differential loss of genes between the two tandem copies. Alternatively, it could have resulted from direct translocation, either of the 12S gene or of the interval that stretches from ND5 to the 16S gene. The tRNA genes seem to have moved independently of the other genes, since it would require a minimum of 15 breakpoints to explain the evolution of the Demodex arrangement if tRNA genes are included in the analysis.

If tRNA genes are also considered in the analysis, however, the extent of rearrangement necessary to derive the Demodex mitochondrial gene order from the chelicerate ancestral pattern increases dramatically, requiring a minimum of 15 breakpoints. If tRNA genes usually move together with the other genes, then the estimated number of breakpoints should have been similar in both analyses. Hence, this increase in the minimum number of breakpoints necessary suggests that mitochondrial tRNA genes have a tendency to move more readily than protein-coding and rRNA genes. Evidence that tRNA genes tend to be more mobile than other genes within the mitochondria of some metazoan groups has been reported previously [20, 27, 28]. One plausible hypothesis is that this rapid movement is simply a consequence of their smaller size—structural rearrangements involving smaller elements may just happen to be more likely. For example, all other things being equal, the probability of a mitochondrial gene rearrangement causing a selective disadvantage may be greater for a larger gene than it is for a smaller gene, simply because the chances are greater that the rearrangement will encompass the entire gene. Alternatively, tRNA genes may be more mobile because they tend to move via a different, and unknown, mechanism.To compare the degree of rearrangement among Acariformes, and to determine whether a similar trend towards greater movement of tRNA genes holds throughout this group, we determined the minimum number of breakpoints necessary to derive each available species’ gene sequence from that of every other species. The results can be depicted as a breakpoint distance tree, in which the number of breakpoints that separates two genome sequences is represented as a distance. We conducted this analysis separately for just the protein-coding and rRNA genes (Figure 3a), and for all of the mitochondrial genes together, including the tRNA genes (Figure 3b). For all lineages, the number of breakpoints necessary to derive the gene order from the ancestral arrangement went up substantially when tRNA genes were included. This result indicates that the more rapid movement of tRNA genes is a general feature of mitochondrial genomes in the Acariformes.

Figure 3
figure 3

Tree showing the breakpoint distances between gene arrangements. Distances on this tree represent the minimum numbers of breakpoints necessary to transform one gene arrangement into another for all pairwise species comparisons. Limulus polyphemus (the horseshoe crab) is considered to represent the ancestral gene arrangement, since it has a gene order that is shared broadly across the chelicerates. In panel (a), only the arrangements of protein-coding and rRNA genes were considered. In panel (b), all mitochondrial genes were considered, including tRNA genes. Both trees are drawn to the same scale. All of the Acariformes have gene arrangements that have diverged substantially from that of L. polyphemus, and the tRNA genes have experienced substantially more rearrangements than the protein-coding or rRNA genes in all lineages. Both Demodex species shared an arrangement that is closer to the chelicerate ancestral arrangement than any other acariform species sequenced to date.

Regardless of whether tRNA genes are included in the analysis, Demodex exhibited a lower degree of rearrangement than has been observed in any other member of the Acariformes to date (Figure 3). In other words, the Demodex arrangement was more similar to the ancestral chelicerate gene arrangement, whereas taxa such as Leptotrombidium or Tetranychus have apparently experienced many more rearrangements. It remains unclear whether different lineages within the Acariformes actually differ in the likelihood that they will experience rearrangements within a given period of time, or if this is simply the amount of variation among lineages that we should expect to see when considering the accumulation of low-frequency events. In support of the latter hypothesis, some taxa in different genera (i.e., Panonychus versus Tetranychus), or even different families (i.e., Aleuroglyphus versus Dermatophagoides), share identical gene arrangements, which suggests that the observed changes in gene order must happen rarely.

The putative control regions were identified in both species based on the following features: (1) they were by far the largest stretches of apparently noncoding DNA (length of 587 bp in D. brevis, 563 bp in D. folliculorum); (2) they had higher than average AT% (81.7% in D. brevis, and 75.7% in D. folliculorum); (3) they were in a conserved position in both species; and (4) the sequences of these regions have apparently been evolving more quickly than the rest of the mitochondrial genome, since they are so divergent that most of the segment could not even be aligned between the two species. All of these features are common observations of the control region in metazoan mitochondrial genomes [29]. Furthermore, the almost complete lack of similarity observed between control region sequences suggests that these lineages split a long time ago.

tRNA gene structures

We identified models for the standard suite of 22 tRNA genes within the mitochondrial genomes of each Demodex species (Figure 4). Importantly, each of these putative tRNA genes was located in the same relative position and orientation within the mitochondrial genome in both species. Furthermore, orthologous tRNA sequences could be aligned reasonably well between the two Demodex species. Finally, for all of the putative genes, the inferred tRNA secondary structures looked similar in both species; this can be seen by comparing the inferred structures from D. brevis (Figure 4a) and D. folliculorum (Figure 4b) for each orthologous pair of tRNA genes. Taken together, these observations provide substantial support for these gene models. If these were not tRNA genes, then it would seem extremely unlikely that the same relative stretch of DNA in both species would happen to have the potential to serve as the template for the same type of tRNA molecule, to exhibit substantial sequence similarity, and to have similar inferred secondary structures. We concluded that these gene models represent a plausible set of tRNA genes for both species, despite the fact that most of them possessed unusual structures.

Figure 4
figure 4

Inferred structures of the mitochondrial tRNAs in both Demodex species. Structures are arranged in alphabetical order. Each tRNA gene is named according to the one-letter amino acid abbreviation, except that L1 = CUN; L2 = UUR; S1 = AGN; and S2 = UCN. A solid line indicates a Watson-Crick bond, whereas a circle indicates a bond between G and uracil. All 22 of these tRNA genes were located in the same relative position and orientation within the mitochondrial genomes of both Demodex species. Furthermore, orthologous tRNA structures all looked similar between Demodex species; to see this, compare each tRNA in panel (a) versus panel (b). Most of the tRNA molecules lacked at least one of the side arms, and some lacked both, suggesting that Demodex must have the molecular machinery necessary for these extremely truncated tRNAs to function.

Interestingly, most of the Demodex tRNA genes were truncated, lacking the sequence necessary to encode arms of the canonical cloverleaf-shaped tRNA (Figure 4). Eight pairs of orthologous tRNAs had just D-arms; three had just T-arms; and five had neither D- nor T-arms. In contrast, just six pairs of orthologous tRNAs retained both D- and T-arms and could fold into the canonical cloverleaf-shaped tRNA. Because the same basic structure could be inferred for each gene in both Demodex species, we concluded that these tRNA structures were probably present in the common ancestor of these lineages. The average sizes of tRNA genes in D. brevis and D. folliculorum agreed closely with each other: 53.5 bp (SE 1.4 bp) and 53.3 bp (SE 1.2 bp), respectively. These collections of tRNA genes were about the same size as the average value of 54.8 bp (SE 1.0 bp) reported previously for the tRNAs within mitochondrial genomes of Acariformes [19]. It appears that truncated mitochondrial tRNA genes were already widespread in the common ancestor of the Acariformes.

Among metazoans, more than 90% of mitochondrial tRNAs are inferred to share the common cloverleaf-shaped secondary structure of nuclear-encoded tRNA sequences [30]. Nevertheless, tRNAs that have lost one arm—either a D-arm or a T-arm—are extremely common among Acariformes [1419]. Indeed, numerous examples of mitochondrial tRNAs that have lost one arm or the other can be found in several additional lineages of the chelicerates [28, 31]. Furthermore, it is common for metazoan mitochondrial genomes generally to possess one to a few tRNA genes that lack either D-arm or T-arm sequence [32]. These noncanonical tRNAs are thought to remain functional due to coevolution of interacting factors. For example, seryl-tRNA synthetase has evolved to be able to recognize the noncanonical tRNA-Ser in mammalian mitochondria [33]. In nematodes, one paralog of EF-Tu has evolved to bind tRNAs that lack the T-arm [34]. Presumably, similar forms of coevolution have occurred in order for the many truncated mitochondrial tRNAs of Acariformes, including Demodex, to remain functional.

Mitochondrial tRNA genes that have lost both D- and T-arms—resulting in minimal tRNAs with only acceptor and anticodon arms—are extremely unusual; nevertheless, these “armless” tRNAs have been inferred in some mitochondrial genomes previously [30]. Indeed, another member of the Tetranychidae (Panonychus ulmi) apparently has at least one such armless tRNA [19]. The trend towards the loss of these tRNA arms appears to have continued in Demodex, resulting in several examples of tRNAs that lack both D- and T-arms. Interestingly, a recent study used RT-PCR and 5′- and 3′-RACE to show that several of these miniaturized, armless tRNAs are indeed transcribed and correctly processed by CCA addition, at least in one species of nematode [35]. If these are actually not functional tRNA genes, then the functional copies must have been imported into the nuclear genome, and we are left with the puzzle of what could be maintaining these sequences.

The similar sequences and conservation of apparent structures for these putative tRNA genes, despite the fact that these are not closely related species, suggests that selection is maintaining these structures. If the truncated tRNA genes are actually no longer functional, then one hypothesis is that these sequences actually function as promoters, and just happen to fold into a remnant of a tRNA gene in both species. This seems implausible, however, given the high levels of divergence often seen in mitochondrial promoter sequences (e.g., [36]), which would lead us to expect the loss of anything resembling a folded tRNA.

Many of the inferred tRNA structures in both species had acceptor arms with fewer than seven paired bases (Figure 4). Indeed, the number of inferred structures with the full complement of seven paired bases in the acceptor arm was only six in D. brevis and 11 in D. folliculorum. To be functional, these molecules with truncated acceptor arms would presumably require either coevolved interacting molecules or RNA editing. Post-transcriptional editing of minimal tRNA molecules, including extensive modification of the acceptor arm, has been demonstrated in velvet worms [37].

The parallel evolution of truncated mitochondrial tRNA genes among chelicerates suggests that this group has evolved a novel mechanism that either permits these tRNAs to function in a truncated state and/or allows post-transcriptional editing to repair them. Determining the mechanisms responsible for keeping such minimal tRNA genes functional in Demodex will be an interesting subject of future research. These results could have implications for developing treatments that affect Demodex without harming the human host, since it has the potential to provide an avenue for the development of new acaricides that target the molecular machinery necessary for the function of such highly unusual tRNA genes.

Phylogenetic position and species divergence

Our phylogenetic analysis recovered traditional groupings within the Acari (Figure 5). For example, the maximum clade credibility tree recovered the Acariformes as monophyletic, with the expected deep split into Trombidiformes and Sarcoptiformes. Within the Trombidiformes, the Eleutherengona and Parasitengona were recovered as monophyletic groups. The posterior probabilities for most of the clades were high. When the Solifugae were included in the analysis, the Unionicola actually fell outside of a clade that included the Eleutherengona plus five members of the Trombiculidae (Additional file 1: Figure S1). Otherwise, tree topologies were completely stable across analyses, regardless of whether the Solifugae were included or not. Demodex always clustered with the Tetranychidae—the spider mites—as expected, since both are considered members of the Eleutherengona (Figure 5). This result agrees with both traditional taxonomy [38] and a recent phylogenetic analysis based on nuclear rRNA sequences [11].

Figure 5
figure 5

Phylogenetic analysis using a relaxed molecular clock based on mitochondrial protein sequences to estimate the age of divergence between Demodex species. Traditional groupings within the Acari were recovered in all phylogenies constructed. This phylogeny represents the results for one particular fossil calibration density, which is based on the minimum age of 410 mya for the time to the most recent common ancestor of the Acariformes (as inferred based on the four oldest Acariforme fossils identified in [39]); the particular fossil calibration density used in this analysis is depicted at the bottom of the tree, along with dashed arrow pointing to the common ancestor of the Acariformes. Amino acid sequences were used for all age-estimation analyses to minimize the effects of mutation saturation. For the particular analysis depicted in this phylogeny, the resulting estimated mean time to the most recent common ancestor of the two Demodex lineages was 134 mya, with a 95% highest probability density interval of 86 – 180 mya. Similar age estimates were obtained for two other fossil calibration densities, one that placed greater probability density near 410 mya, and one with a uniform distribution between 410–510 mya (Table 1). Furthermore, these estimates were not altered substantially when two species of Solifugae were included in the analysis (Table 1 and Additional file 1: Figure S1). Finally, these results are broadly overlapping with the results obtained based on an alignment of 18S sequences. Overall, these results suggest that Demodex have probably been coevolving with their mammalian hosts since before the placental radiation.

In the absence of either a fossil record or estimates of molecular divergence, it is unclear how long Demodex have inhabited the mammalian skin. For example, it is possible that Demodex were present on the mammalian common ancestor, and have been speciating with their hosts ever since. Since we do not have estimates of how often these mites can move between host lineages, however, it is also possible that Demodex evolved to parasitize mammals much more recently and are widespread across mammalian lineages simply because they have repeatedly colonized new host species.

Interestingly, the results of the relaxed molecular clock analysis based on the amino acid sequences of mitochondrial proteins suggest that the split between D. brevis and D. folliculorum happened a very long time ago (Figure 5 and Table 1). For example, the youngest average estimated time of divergence between these two Demodex lineages was 129 mya, with a 95% highest posterior density (HPD) interval of 87 – 173 mya. Across six different relaxed molecular clock analyses, the average lower bound of the 95% HPD interval was 87 mya for the time of the most recent common ancestor of the two Demodex lineages. Assuming that the transition to living on mammals happened once, these results suggest that the genus Demodex has been evolving and diversifying with their mammalian hosts for at least this long.

Table 1 Summary statistics of relaxed molecular clock estimates of the divergence time between Demodex species based on amino acid sequence alignments of mitochondrial proteins

The average age estimate across analyses for the D. brevis – D. folliculorum split was 136 mya (Table 1). For comparison, these same analyses resulted in average estimated times to the most recent common ancestor of the Acari of 483 – 520 mya. It is remarkable that the split between D. brevis and D. folliculorum appears to have happened about 26 – 28% as long ago as the original divergence between Acariformes and Parasitiformes within the Acari, especially given that these Demodex species are so similar morphologically [1].

To calibrate our relaxed molecular clock analyses, we used a date of ~410 mya as the minimum age of the Acariformes, based on the four oldest acariform fossils so far identified [39]. These fossils include three that are considered members of the Endeostigmata, and one that is considered a member of the Trombidiformes. Molecular phylogenies support the hypothesis that the Endeostigmata represent a basal lineage of the Sarcoptiformes [40]. If, however, the endeostigmatan fossils actually turn out to represent a basal lineage of the Acariformes (i.e., falling outside of the Sarcoptiformes-Trombidiformes clade), then our fossil calibration date of 410 mya would not accurately represent the common ancestor of the Sarcoptiformes and Trombidiformes, which would bias the estimated age of the D. brevisD. folliculorum split, making this split appear older than it is in reality.

To test the impact of different fossil calibration densities on our age estimates, we used three different distribution densities for the priors on age of the Acariformes (Table 1), all based on the ~410 mya time of appearance of members of this group in the fossil record [39]. Prior distributions used to represent the age of the Acariformes varied from highly restrictive (i.e., using an exponential distribution with most of the density near the 410 mya minimum age), to much less informative (i.e., using a uniform distribution with a range from 410 – 510 mya). Altering the prior distributions had no substantial effect on the estimated age of the Demodex divergence.

To test the impact of different taxon combinations on our age estimates, we varied the relaxed molecular clock analysis to either include or exclude two members of the Solifugae (Table 1). These taxa were chosen because two recent studies place them as the sister group to the Acariformes, which would make the Acari a paraphyletic grouping [40, 41]. We did not observe any sensitivity of the Demodex age estimates to whether Solifugae were included in the analysis. Furthermore, we did not see any evidence that the Solifugae are the sister taxon to Acariformes; in other words, when the Solifugae were included in the analyses, we always recovered a monophyletic Acari, with the posterior probability of the node uniting the Parasitiformes and Acariformes into the Acari always equal to one, providing strong support for the traditional grouping (to the exclusion of the Solifugae). Although this was by no means a comprehensive test of the higher order relationships among the Chelicerata, these results should nevertheless be added into the growing debate as to the correct topology of these relationships.

The common ancestor of the placental mammals has been estimated to have lived 72 – 108 mya, with the radiation of placental ordinal level crown groups occurring later [42]. Hence, our molecular clock results indicate that the two species found on humans may have last shared a common ancestor prior to the ordinal level radiation of the placental mammals, which is consistent with an ancient colonization of mammals by the Demodex lineage. Such an ancient colonization event would help to explain the wide distribution of Demodex across mammalian orders, as well as the clear morphological adaptations to living in the pilosebaceous complex. Furthermore, since Jurassic mammals had fur, these data are in agreement with the fossil record of mammals [43].

One potential shortcoming of our dating analyses is that mutation saturation could be influencing the results. In particular, if there is substantial saturation of the substitutions among the oldest splits in the phylogeny, then this will tend to cause the over-estimation of split ages near the tips of the tree, which could make our estimation of the D. brevis – D. folliculorum split look older than it actually is. For this reason, we used amino acid sequences rather than DNA sequences for our molecular dating analyses. Owing to their reduced state space (four possible bases), nucleotide sequences will saturate much more rapidly than protein sequences (20 possible amino acids). Furthermore, to assess whether there was substantial saturation among the amino acid substitutions, we also used the software AsaturA [44] to visualize the accumulation of amino acid substitutions with genetic distance, which showed that there was only very slight saturation of even the most frequent types of amino acid substitutions in these data. This suggests that our estimate of an ancient split between D. brevis and D. folliculorum was not being caused artificially by mutation saturation.

A second potential shortcoming of our dating analyses is that they are based on a single fossil calibration point. If there have been substantial changes over time in the rates of molecular evolution of the mitochondrial protein-coding genes used in the analysis, then this will tend to reduce the accuracy of our estimated D. brevis – D. folliculorum divergence time. On the other hand, the analysis method modeled evolution in each branch of the phylogeny with a relaxed molecular clock, meaning that rates of change were allowed to vary, which should have helped to ameliorate problems associated with slowdown or speedup of the rates. Furthermore, the age estimate is based on aligned sequences of four different proteins, each of which was allowed to slow down or speed up independently in the analysis. Finally, inspection of the rates estimated for each gene along different branches of the phylogeny did not indicate that there was a systematic change in evolutionary rates in the lineage leading to the Demodex species within the Acariformes. Nevertheless, because we based our divergence time estimates on a single fossil calibration date, the possibility remains that undetected changes in molecular evolutionary rates influenced the results.

One way to address these possible shortcomings is to ask whether the divergence estimate based on another gene is consistent with those based on the mitochondrial protein sequences. To this end, we conducted a relaxed molecular clock analysis based on a collection of 18S rDNA sequences for the Acariformes, with L. polyphemus as an outgroup. Based on the 18S alignment, the mean time to most recent common ancestor of D. brevis and D. folliculorum was 74 mya, with a 95% HPD interval of 12 – 150 mya. Unfortunately, this is a large time interval, which does not provide a strong test of our original age estimates. Nevertheless, this interval is broadly overlapping with the results of our analyses based on the mitochondrial protein sequences (i.e., compare this result to the intervals reported in Table 1), so our original estimates based on mitochondrial proteins are not contradicted by this estimate using 18S rDNA sequences. Furthermore, the average estimate of 74 mya based on the 18S genes would still make this an ancient split between D. brevis and D. folliculorum (i.e., still occurring prior to the radiation of placental mammals). If the area of overlap between the 95% HPD intervals based on both 18S genes and mitochondrial proteins is considered, then we would estimate that these taxa split somewhere between about 87 – 150 mya.

D. brevis and D. folliculorum live in distinct habitats within the skin: the former inhabits the sebaceous glands whereas the latter resides in the hair follicles nearer the skin surface, alongside the hair shaft. The ancient divergence observed between these lineages is consistent with the hypothesis that there was an early separation into mite lineages that specialize in either the hair follicle or the sebaceous gland. This hypothesis has been suggested previously [45], based on the observation that there are consistent morphological differences among mites from different skin microhabitats—mite species from different hosts but found in the same microhabitat tend to display more similarities than those living on the same host but in different microhabitats. These similarities are not restricted to overall shape, as might be expected to result from natural selection for squeezing into particular skin structures, but have also been observed in morphological characters that are regarded as diagnostic in taxonomy of Demodicidae. This hypothesis could be tested by comparing sequence divergence among additional Demodex species that reside either deep in glands within the skin versus near the surface of the follicle in other mammalian species. If correct, this hypothesis predicts that gland-dwelling versus follicle-dwelling species will tend to group into distinct clades. Although a small amount of sequence information is presently available for several other members of the genus, the classification of these species into primarily follicle- versus gland-dwelling taxa is unfortunately not clear, unlike the case for D. brevis versus D. folliculorum. A further test of this hypothesis will need to await further information about the primary habitats for taxa with available sequence information.

Conclusions

These represent the first determinations of the complete mitochondrial genome sequences from any member of the genus Demodex. The availability of this additional genetic information promises to open the way for studies of variation within and between Demodex species. This is important because there has apparently been an ancient radiation of these mites on their mammalian hosts. Furthermore, Demodex are known to be ubiquitous in humans, and to cause medically important skin disorders.

We found that Demodex mitochondrial genomes share the following features with other members of the Acariformes: (1) they are compact in size and AT-rich; (2) they have experienced multiple gene rearrangements, especially among the tRNA genes; and (3) they have many truncated tRNA genes.

It remains unclear what is driving the evolution of truncated tRNA genes. Perhaps these are favored simply because of their reduced size, resulting in less RNA needing to be transcribed in order for translation to occur, and also less mitochondrial genome to be replicated. This hypothesis does not explain why selection for smaller tRNA genes has been so much more effective in certain lineages, such as the Acariformes. To be functional, the many truncated tRNA genes found in the mitochondrial genome would seem to require some combination of coevolved interacting factors and/or extensive tRNA editing. The molecular machinery necessary for dealing with truncated tRNAs might provide avenues to treat the apparently harmful Demodex found in some patients.

Based on the four most slowly evolving protein-coding genes within the mitochondrial genomes of the Acari, we examined the phylogenetics of the Acariformes, and supported the traditional hypothesis that the Demodex are members of the Eleutherengona. Interestingly, we found that the two Demodex species found on humans apparently diverged more than 87 mya, which is prior to the estimated time of the radiation of the placental mammals. Assuming a single transition to living in mammalian skin, this estimate places a lower bound on how long Demodex have been living on their mammalian hosts. Furthermore, the deep divergence time estimated between these two species is consistent with the hypothesis that there was an early split into distinct forms that either live deep in the sebaceous glands or near the surface alongside the hair shaft. This niche-separation hypothesis can be tested by examining the phylogenetics of additional Demodex species that have been shown to live in different skin microhabitats.

Methods

Ethics statement

Prior to sampling, each participant to be sampled for Demodex was informed both verbally and in writing about the goals of the project and the sampling protocol. All participants signed a written Informed Consent form. All human Demodex sampling procedures and the Informed Consent form were approved by Bowdoin College’s Research Oversight Committee, Approval No. 2007–34.

Mite isolation and molecular techniques

Mites were isolated by drawing the curved end of a new bobby pin across the forehead of each participant. The resulting exudate was searched for mites in mineral oil under a stereomicroscope. Individual mites were verified as exhibiting the genital morphology of either D. folliculorum or a D. brevis (per [1]) using 600x magnification on a compound light microscope. To purify DNA, each mite was washed several times in fresh mineral oil, then the mineral oil was removed by washing 10x with 100% ethanol; the ethanol was evaporated by heating 2 min at 95°C, and the dried mite was then resuspended in 10 μL Lysis Buffer (10 μL PCR Buffer + 8 units Proteinase K in 10 μL H2O + 80 μL 1% Triton X), incubated 60 min at 65°C followed by 10 min at 95°C, then frozen for at least 1 hr. For initial PCR experiments the DNA prepared from multiple mites was pooled, whereas DNA from a single female mite from each species that had been diluted 100-fold in water was used as template for the PCR experiments necessary to generate the final mitochondrial sequences.

The first fragment successfully amplified from D. folliculorum was a fragment of the 12S rRNA gene, which was obtained using degenerate PCR based on an alignment of 12S rRNA from various arthropods. Primers were then designed within this small fragment so that long-range PCR could be used to amplify the rest of the mitochondrial genome. The resulting PCR product was ~14 kb in length, and the sequence of that large product was determined with a combination of primer walking and subcloning of fragments into plasmids using various restriction enzyme digestions and ligations, followed by plasmid sequencing, together with PCR experiments designed to link subcloned fragments. Once a draft of the mitochondrial genome had been assembled, a collection of PCR products that covered the entire sequence was generated from a single mite; these products were sequenced on both strands to determine the final mitochondrial genome sequence.

The first fragment successfully amplified from D. brevis was part of the ND5 gene, using PCR primers designed from the D. folliculorum mitochondrial genome. Otherwise, the same techniques were employed as above to generate a draft of the D. brevis genome, and then to determine the mitochondrial genome sequence from a single mite.

To verify the species identity of our final genome sequences, 20 individual mites from each species were isolated and identified based on the genital morphology described by Desch and Nutting [1]. DNA was prepared for each mite separately, and PCR amplifications using species-specific primer pairs (Additional file 2: Table S1), designed based on each mitochondrial genome sequence, were carried out to test species identity. Successful PCR amplifications with species-specific primers correlated perfectly with identifications based on morphology (P < 0.0001, Fisher’s Exact Test). Finally, we verified that our sequences agreed with fragments of the 16S [12] and COX1 [3] genes that have already been determined for both species. We concluded that the mitochondrial genome sequences correspond to the correct species.

Additional file 2: Table S1 lists the PCR primer pairs, annealing temperatures, and product sizes for all of the amplifications that were necessary to determine the mitochondrial genome sequences of both species. Additional primers were designed, as needed, for linking plasmid subclones during the compilation of the initial genome sequences, and for sequencing PCR products and plasmids. All primers were purchased from Invitrogen (Carlsbad CA, USA). Routine PCR was performed using GoTaq Hot Start DNA Polymerase (Promega, Madison WI, USA) or Phusion Hot-Start High-Fidelity DNA Polymerase (New England Biolabs, Ipswich MA, USA). Long-range PCR was accomplished using the Hot-Start version of TaKaRa LaTaq DNA Polymerase (Clontech Laboratories, Mountain View CA, USA). Both PCR products and plasmids were analyzed on 1% agarose gels. PCR was carried out using standard techniques in 50-μL volumes using an MJ Research PTC-200 thermocycler (Bio-Rad, Hercules CA, USA). PCR products were purified with the QIAquick PCR Purification Kit and plasmid DNA was purified using the QIAprep Spin Miniprep Kit (Qiagen Inc, Valencia CA, USA). Long-range PCR products were digested and subcloned into pBluescript II SK+ vector (Stratagene, La Jolla CA, USA) for plasmid sequencing using standard techniques. PCR products and plasmids were sent to either Geneway Research (Hayward, California, USA) or Mount Desert Island Biological Laboratory DNA Sequencing Core (Salisbury Cove, Maine, USA) for Sanger sequencing.

Sequence annotation and inferences of secondary structures

All sequencing results were compiled and aligned for the purposes of determining the mitochondrial genomes, and oligonucleotide primers for both PCR and sequencing were designed, using MacVector (version 10.6, North Carolina, USA). D. brevis and D. folliculorum protein-coding genes were mapped by aligning the new genome sequences against both translated and untranslated sequences from a collection of taxa within the Chelicerata (Additional file 3: Table S2). The beginnings and endings of these coding sequences were based on putative start and stop codons nearest the ends of the alignable sequences. Locations and orientations of rRNA genes were mapped by alignment with rRNA genes from the same collection of chelicerate taxa, although the beginnings and endings of these genes were simply based on the limits of neighboring genes and so must be considered extremely tentative. The tRNA genes were identified by first analyzing the mitochondial sequences with four separate software packages that have been designed to search for tRNA genes: ARWEN [46], DOGMA [47], tRNAscan-SE [48], and MITOS [49]. Where possible, these software packages were all run using the most relaxed possible settings for recognition of tRNA structures. The resulting potential tRNAs were also aligned with the orthologous tRNAs from the standard collection of chelicerate taxa in order to verify their identity. Finally, a small number of tRNA genes could only be identified manually, and these were also verified by alignment. tRNA genes that were in the same relative physical locations in both Demodex species, exhibited a high degree of sequence similarity between Demodex species, showed at least a moderate degree of sequence similarity with other chelicerate sequences, and could be folded into plausible secondary structures, were considered real. Inferred tRNA structures were drawn using VARNA [50] and imported into Illustrator CS5 (Adobe Systems, Carlsbad CA, USA) for final graphics production.

The putative control regions were inferred as likely to comprise the largest noncoding stretch in each genome, while verifying that each of these also had a higher-than-average AT% nucleotide composition, although these stretches were too divergent to be aligned reliably between species.

Gene order comparisons

Gene order comparisons were made using the SPRING web server [51], which takes a sample of chromosome gene orders and orientations as its input and then computes a minimum series of reversals and/or block-interchanges necessary for transforming each chromosome into every other chromosome. Codified gene orders and orientations for the 13 protein-coding genes plus the two rRNA genes were input into web server for the following: (1) gene arrangement of the horseshoe crab Limulus polyphemus, representing the ancestral chelicerate gene arrangement; (2) fourteen of the Acariformes species with available mitochondrial genome sequences (Leptotrombidium pallidum was excluded from the analysis because this species has a duplication within the mitochondrial genome); and (3) both D. folliculorum and D. brevis. The same analysis was completed including tRNA genes, except that Steganacarus magnus was excluded in this case because one tRNA gene is missing from that genome. The SPRING web server also estimates a distance tree based on the minimum number of breakpoints necessary to explain the gene arrangements; we downloaded the topology in Newick format and generated a drawing of this tree in FigTree, which was exported to Illustrator for final graphics generation.

Phylogenetics and divergence time estimation

We chose to focus our phylogenetic and molecular clock analysis on four protein-coding genes that exhibit the slowest rates of nonsynonymous substitution, at least among the Acari [19]—COX1, COX2, COX3, and CYTB. Based on the protein sequences of these four genes, we inferred the phylogenetic position and divergence time of D. brevis and D. folliculorum based on a concatenated multiple sequence alignment that was generated using the GUIDANCE webserver [52], and the resulting alignment was used to estimate phylogenetic relationships and divergence times using the Markov chain Monte Carlo search of parameter space implemented in BEAST (v. 1.7.4) [53].

We extracted from the NCBI database the sequences for the COX1, COX2, COX3, and CYTB genes from the complete mitochondrial genomes of 21 Acari, plus the two available species of Solifugae, plus the horseshoe crab (Limulus polyphemus) as an outgroup (Additional file 3: Table S2). All of the available Acariformes in the database were included in the analysis (15 altogether), whereas only a small subset of the available Parasitiformes taxa were included, although these were chosen in order to sample a wide diversity within this clade. The Solifugae were included because both Dabert et al. [40] and Pepato et al. [41] provide support for the hypothesis that the Acari are actually paraphyletic, suggesting that the Solifugae actually group with the Acariformes to the exclusion of the Parasitiformes. We conducted the same analyses separately with the Solifugae included or excluded in order to compare the results of different taxon combinations on the phylogenetic positions and estimated divergence time of the two Demodex species. All DNA sequences were translated into amino acid sequences to minimize the influence of mutation saturation on the results.

To generate a multiple sequence alignment that consisted of only reliably aligned amino acids, thereby minimizing possible artifacts introduced due to alignment error, we used the online web server GUIDANCE [52]. Each individual gene was aligned separately, and the separate gene alignments were combined into a concatenated alignment (totaling 1279 or 1293 amino acids for the alignment including or excluding the Solifugae, respectively). The GUIDANCE alignment analyses were conducted using the MAFFT algorithm and 100 bootstrap replicates, with variable amino acid position columns identified using the GUIDANCE method and then removed based on a column-score cutoff of <0.93. The overall GUIDANCE scores for the each gene was always above 0.94, so most of each gene was included in the final alignment. The alignment file was converted to NEXUS format to be used as input for the software BEAUTi (version 1.7.4), which was used to set the Bayesian analysis parameters and to generate the input file for BEAST [53].

All BEAST runs were conducted using a chain length of 100,000,000, lognormal relaxed clock (uncorrelated), Yule process speciation tree prior, and an uninformative ucld.mean prior (exponential, mean = 1). The substitution and clock models for the four genes were unlinked, whereas the tree models were linked. Two separate BEAST analyses were conducted for each set of conditions, and the results inspected to verify that they converged on very similar parameter values and identical tree topologies, although the results of just one run for each set of conditions is reported. To calibrate the molecular clock, we used three different fossil calibration densities as priors for the age of the Acariformes, all based on the ~410mya time of appearance of members of this group in the fossil record [39]: (1) Exponential [10] with offset 410—this is the most restrictive prior, with most of the density very close to the 410 mya minimum age based on the fossil record, and 95% of the distribution density falling between 410 – 440 mya; (2) Gamma [2, 15] with offset 410—this distribution is somewhat less restrictive, and has greater density concentrated near the estimated Acariformes age of ~435mya reported by Dabert et al. [40], with 95% of the distribution density falling between 414 – 494 mya; and (3) Uniform [410, 510]—this is a relatively uninformative prior stretching over the 100-my interval starting from the oldest known Acariformes fossils. This allowed us to determine the effects of assuming different fossil calibration densities on the estimated age of the Demodex species divergence. All other settings were left in the default state or set to the uniform distribution with an arbitrarily large range. Tracer (version 1.5, [54]) was used to verify the convergence of runs by discarding the first 10% of the samples as burn-in, and examining the effective sample size of all parameters—the results were considered reliable only when the effective sample size of all parameters was above 100—as well as by visual inspection of the trace graphs for good mixing. We summarized the trees using TreeAnnotator ver. 1.7.4 [53], with the first 10% of trees discarded as burn-in; samples from the posterior were summarized on the maximum credibility tree for both the probability of each node as well as the 95% HPD of node heights. Finally, FigTree ver. 1.4.0 was used for tree visualization. The scalable vector graphic tree output from FigTree was imported into Adobe Illustrator (CS5) for final graphics generation.

Finally, we conducted a similar analysis to estimate time of divergence based on 18S rDNA gene sequences of various Acariformes that are available in GenBank. We conducted an analysis of estimated divergence based on a collection of 18S sequences, including D. brevis, D. folliculorum, 11 other Acariformes from both Sarcoptiformes and Trombidiformes, and Limulus polyphemus as the outgroup (Additional file 4: Table S3). The 18S rDNA sequences were aligned using the GUIDANCE web server, and the resulting alignment was inspected visually to insure that the inserted gaps resulted in a reasonable alignment across these taxa. The alignment was used to conduct a relaxed molecular clock analysis in BEAST using the same fossil calibration point as was used previously (i.e., the tmrca of Acariformes was assumed to be at least 410 mya) and the intermediate level of restriction on the prior for the time to most recent common ancestor of the Acariformes (i.e., assuming the Gamma [2, 15] distribution with an offset of 410 mya); we also assumed a general-time-reversible model of nucleotide substitution, with site heterogeneity of four gamma categories plus invariant sites, a lognormal relaxed clock (uncorrelated), and speciation was assumed to follow a Yule Process. Other settings were left as defaults.

Availability of supporting data

The new mitochondrial genome sequences reported in this article are available in the GenBank repository, accession numbers KM114225 and KM114226. The accession numbers for the other complete mitochondrial genomes used in the analyses are listed in Additional file 3: Table S2. The accession numbers for the 18S rDNA gnes used for comparison are listed in Additional file 4: Table S3.

Abbreviations

A:

Adenine

T:

Thymine

T-arm:

Pseudouridine-arm of a tRNA secondary structure

D-arm:

Dihydrouridine-arm of a tRNA secondary structure

rRNA:

Ribosomal RNA

tRNA:

Transfer RNA

12S:

Gene for small subunit ribosomal RNA

16S:

Gene for large subunit ribosomal RNA

Atp6 and −8:

ATPase subunit 6 and 8

CO1-3 or COX1-3:

Cytochrome oxidase subunits I-III

Cytb:

Cytochrome b

ND1-6 or NAD1-6 and ND4L:

NADH dehydrogenase subunits 1–6 and subunit 4 L

A+T:

Large AT-rich non-coding region/control region

Bp:

Base pairs.

References

  1. Desch CE, Nutting WB: Demodex folliculorum (Simon) and D. brevis Akbulatova of man: redescription and reevalution. J Parasitol. 1972, 58 (1): 169-177. 10.2307/3278267.

    Article  CAS  PubMed  Google Scholar 

  2. Desch CE: Human hair follicle mites and forensic acarology. Exp Appl Acarol. 2009, 49 (1–2): 143-146.

    Article  PubMed  Google Scholar 

  3. Thoemmes ME, Fergus DJ, Urban J, Trautwein M, Dunn RR: Ubiquity and diversity of human-associated Demodex mites. PLoS One. in press

  4. Zhao YE, Wu LP, Hu L, Xu JR: Association of blepharitis with Demodex: a meta-analysis. Ophthalmic Epidemiol. 2012, 19 (2): 95-102. 10.3109/09286586.2011.642052.

    Article  PubMed  Google Scholar 

  5. Jarmuda S, O’Reilly N, Zaba R, Jakubowicz O, Szkaradkiewicz A, Kavanagh K: Potential role of Demodex mites and bacteria in the induction of rosacea. J Med Microbiol. 2012, 61 (Pt 11): 1504-1510.

    Article  PubMed  Google Scholar 

  6. Nara T, Katoh N, Inoue K, Yamada M, Arizono N, Kishimoto S: Eosinophilic folliculitis with a Demodex folliculorum infestation successfully treated with ivermectin in a man infected with human immunodeficiency virus. Clin Exp Dermatol. 2009, 34 (8): e981-e983. 10.1111/j.1365-2230.2009.03621.x.

    Article  CAS  PubMed  Google Scholar 

  7. Chen W, Plewig G: Human demodicosis: revisit and a proposed classification. Br J Dermatol. 2014, 170 (6): 1219-1225. 10.1111/bjd.12850.

    Article  CAS  PubMed  Google Scholar 

  8. Escobar-Páramo P, Grenet K, Le Menac’h A, Rode L, Salgado E, Amorin C, Gouriou S, Picard B, Rahimy MC, Andremont A, Denamur E, Ruimy R: Large-scale population structure of human commensal Escherichia coli isolates. Appl Environ Microbiol. 2004, 70 (9): 5698-5700. 10.1128/AEM.70.9.5698-5700.2004.

    Article  PubMed Central  PubMed  Google Scholar 

  9. Jones EP, Eager HM, Gabriel SI, Jóhannesdóttir F, Searle JB: Genetic tracking of mice and other bioproxies to infer human history. Trends Genet. 2013, 29 (5): 298-308. 10.1016/j.tig.2012.11.011.

    Article  CAS  PubMed  Google Scholar 

  10. Reed DL, Light JE, Allen JM, Kirchman JJ: Pair of lice lost or parasites regained: the evolutionary history of anthropoid primate lice. BMC Biol. 2007, 5: 7-10.1186/1741-7007-5-7.

    Article  PubMed Central  PubMed  Google Scholar 

  11. Zhao YE, Wu LP, Hu L, Xu Y, Wang ZH, Liu WY: Sequencing for complete rDNA sequences (18S, ITS1, 5.8S, ITS2, and 28S rDNA) of Demodex and phylogenetic analysis of Acari based on 18S and 28S rDNA. Parasitol Res. 2012, 111 (5): 2109-2114. 10.1007/s00436-012-3058-8.

    Article  PubMed  Google Scholar 

  12. Zhao YE, Hu L, Ma JX: Molecular identification of four phenotypes of human Demodex mites (Acari: Demodicidae) based on mitochondrial 16S rDNA. Parasitol Res. 2013, 112 (11): 3703-3711. 10.1007/s00436-013-3558-1.

    Article  PubMed  Google Scholar 

  13. Walter DE, Proctor H: Mites: Ecology, Evolution & Behaviour. 2013, Netherlands: Springer, 2

    Book  Google Scholar 

  14. Shao R, Barker SC, Mitani H, Takahashi M, Fukunaga M: Molecular mechanisms for the variation of mitochondrial gene content and gene arrangement among chigger mites of the genus Leptotrombidium (Acari: Acariformes). J Mol Evol. 2006, 63 (2): 251-261. 10.1007/s00239-005-0196-y.

    Article  CAS  PubMed  Google Scholar 

  15. Domes K, Maraun M, Scheu S, Cameron SL: The complete mitochondrial genome of the sexual oribatid mite Steganacarus magnus: genome rearrangements and loss of tRNAs. BMC Genomics. 2008, 9: 532-10.1186/1471-2164-9-532.

    Article  PubMed Central  PubMed  Google Scholar 

  16. Dermauw W, Van Leeuwen T, Vanholme B, Tirry L: The complete mitochondrial genome of the house dust mite Dermatophagoides pteronyssinus (Trouessart): a novel gene arrangement among arthropods. BMC Genomics. 2009, 10: 107-10.1186/1471-2164-10-107.

    Article  PubMed Central  PubMed  Google Scholar 

  17. Ernsting BR, Edwards DD, Aldred KJ, Fites JS, Neff CR: Mitochondrial genome sequence of Unionicola foili (Acari: Unionicolidae): a unique gene order with implications for phylogenetic inference. Exp Appl Acarol. 2009, 49 (4): 305-316. 10.1007/s10493-009-9263-1.

    Article  CAS  PubMed  Google Scholar 

  18. Klimov PB, OConnor BM: Improved tRNA prediction in the American house dust mite reveals widespread occurrence of extremely short minimal tRNAs in acariform mites. BMC Genomics. 2009, 10: 598-10.1186/1471-2164-10-598.

    Article  PubMed Central  PubMed  Google Scholar 

  19. Yuan M-L, Wei D-D, Wang B-J, Dou W, Wang J-J: The complete mitochondrial genome of the citrus red mite Panonychus citri (Acari: Tetranychidae): high genome rearrangement and extremely truncated tRNAs. BMC Genomics. 2010, 11: 597-10.1186/1471-2164-11-597.

    Article  PubMed Central  PubMed  Google Scholar 

  20. Gissi C, Iannelli F, Pesole G: Evolution of the mitochondrial genome of Metazoa as exemplified by comparison of congeneric species. Heredity (Edinb). 2008, 101 (4): 301-320. 10.1038/hdy.2008.62.

    Article  CAS  Google Scholar 

  21. Min XJ, Hickey DA: DNA barcodes provide a quick preview of mitochondrial genome composition. PLoS One. 2007, 2 (3): e325-10.1371/journal.pone.0000325.

    Article  PubMed Central  PubMed  Google Scholar 

  22. Martin AP: Metabolic rate and directional nucleotide substitution in animal mitochondrial DNA. Mol Biol Evol. 1995, 12 (6): 1124-1131.

    CAS  PubMed  Google Scholar 

  23. Rocha EP, Danchin A: Base composition bias might result from competition for metabolic resources. Trends Genet. 2002, 18 (6): 291-294. 10.1016/S0168-9525(02)02690-2.

    Article  CAS  PubMed  Google Scholar 

  24. Boore JL: Animal mitochondrial genomes. Nucleic Acids Res. 1999, 27 (8): 1767-1780. 10.1093/nar/27.8.1767.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  25. Fahrein K, Talarico G, Braband A, Podsiadlowski L: The complete mitochondrial genome of Pseudocellus pearsei (Chelicerata: Ricinulei) and a comparison of mitochondrial gene rearrangements in Arachnida. BMC Genomics. 2007, 8: 386-10.1186/1471-2164-8-386.

    Article  PubMed Central  PubMed  Google Scholar 

  26. Masta SE: Mitochondrial rRNA secondary structures and genome arrangements distinguish chelicerates: comparisons with a harvestman (Arachnida: Opiliones Phalangium opilio). Gene. 2010, 449 (1–2): 9-21.

    Article  CAS  PubMed  Google Scholar 

  27. Dowton M, Cameron SL, Dowavic JI, Austin AD, Whiting MF: Characterization of 67 mitochondrial tRNA gene rearrangements in the Hymenoptera suggests that mitochondrial tRNA gene position is selectively neutral. Mol Biol Evol. 2009, 26 (7): 1607-1617. 10.1093/molbev/msp072.

    Article  CAS  PubMed  Google Scholar 

  28. Ovchinnikov S, Masta SE: Pseudoscorpion mitochondria show rearranged genes and genome-wide reductions of RNA gene sizes and inferred structures, yet typical nucleotide composition bias. BMC Evol Biol. 2012, 12: 31-10.1186/1471-2148-12-31.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  29. Lee WJ, Conroy J, Howell WH, Kocher TD: Structure and evolution of teleost mitochondrial control regions. J Mol Evol. 1995, 41 (1): 54-66.

    Article  CAS  PubMed  Google Scholar 

  30. Jühling F, Pütz J, Bernt M, Donath A, Middendorf M, Florentz C, Stadler PF: Improved systematic tRNA gene annotation allows new insights into the evolution of mitochondrial tRNA structures and into the mechanisms of mitochondrial genome rearrangements. Nucleic Acids Res. 2012, 40 (7): 2833-2845. 10.1093/nar/gkr1131.

    Article  PubMed Central  PubMed  Google Scholar 

  31. Masta SE, Boore JL: Parallel evolution of truncated transfer RNA genes in arachnid mitochondrial genomes. Mol Biol Evol. 2008, 25 (5): 949-959. 10.1093/molbev/msn051.

    Article  CAS  PubMed  Google Scholar 

  32. Watanabe YI, Suematsu T, Ohtsuki T: Losing the stem-loop structure from metazoan mitochondrial tRNAs and co-evolution of interacting factors. Front Genet. 2014, 5: 109-

    Article  PubMed Central  PubMed  Google Scholar 

  33. Chimnaronk S, Jeppesen MG, Suzuki T, Nyborg J, Watanabe K: Dual-mode recognition of noncanonical tRNAsSer by seryl-tRNA synthetase in mammalian mitochondria. EMBO J. 2005, 24: 3369-3379. 10.1038/sj.emboj.7600811.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  34. Ohtsuki T, Watanabe Y, Takemoto C, Kawai G, Ueda T, Kita K, Kojima S, Kaziro Y, Nyborg J, Watanabe K: An “elongated” translation elongation factor Tu for truncated tRNAs in nematode mitochondria. J Biol Chem. 2001, 276 (24): 21571-21577. 10.1074/jbc.M011118200.

    Article  CAS  PubMed  Google Scholar 

  35. Wende S, Platzer EG, Jühling F, Pütz J, Florentz C, Stadler PF, Mörl M: Biological evidence for the world’s smallest tRNAs. Biochimie. 2014, 100: 151-158.

    Article  CAS  PubMed  Google Scholar 

  36. Fisher RP, Parisi MA, Clayton DA: Flexible recognition of rapidly evolving promoter sequences by mitochondrial transcription factor 1. Genes Dev. 1989, 3 (12B): 2202-2217. 10.1101/gad.3.12b.2202.

    Article  CAS  PubMed  Google Scholar 

  37. Segovia R, Pett W, Trewick S, Lavrov DV: Extensive and evolutionarily persistent mitochondrial tRNA editing in Velvet Worms (phylum Onychophora). Mol Biol Evol. 2011, 28 (10): 2873-2881. 10.1093/molbev/msr113.

    Article  CAS  PubMed  Google Scholar 

  38. Mironov SV, Bochkov AV: Modern conceptions concerning the macrophylogeny of acariform mites (Chelicerata, Acariformes). Entomol Rev. 2009, 89 (8): 975-992. 10.1134/S0013873809080120.

    Article  Google Scholar 

  39. Dunlop JA, Selden PA: Calibrating the chelicerate clock: a paleontological reply to Jeyaprakash and Hoy. Exp Appl Acarol. 2009, 48 (3): 183-197. 10.1007/s10493-009-9247-1.

    Article  PubMed  Google Scholar 

  40. Dabert M, Witalinski W, Kazmierski A, Olszanowski Z, Dabert J: Molecular phylogeny of acariform mites (Acari, Arachnida): strong conflict between phylogenetic signal and long-branch attraction artifacts. Mol Phylogenet Evol. 2010, 56 (1): 222-241. 10.1016/j.ympev.2009.12.020.

    Article  PubMed  Google Scholar 

  41. Pepato AR, da Rocha CEF, Dunlop JA: Phylogenetic position of the acariform mites: sensitivity to homology assessment under total evidence. BMC Evol Biol. 2010, 10: 235-10.1186/1471-2148-10-235.

    Article  PubMed Central  PubMed  Google Scholar 

  42. dos Reis M, Donoghue PC, Yang Z: Neither phylogenomic nor palaeontological data support a Palaeogene origin of placental mammals. Biol Lett. 2014, 10 (1): 20131003-10.1098/rsbl.2013.1003.

    Article  PubMed Central  PubMed  Google Scholar 

  43. Zhou CF, Wu S, Martin T, Luo ZX: A Jurassic mammaliaform and the earliest mammalian evolutionary adaptations. Nature. 2013, 500 (7461): 163-167. 10.1038/nature12429.

    Article  CAS  PubMed  Google Scholar 

  44. Van de Peer Y, Frickey T, Taylor JS, Meyer A: Dealing with saturation at the amino acid level: a case study involving anciently duplicated zebrafish genes. Gene. 2002, 295: 205-211. 10.1016/S0378-1119(02)00689-3.

    Article  CAS  PubMed  Google Scholar 

  45. Izdebska JN: Selected aspects of adaptations to the parasitism of hair follicle mites (Acari: Demodecidae) from hoofed mammals. Eur Bison Cons News. 2009, 2: 80-88.

    Google Scholar 

  46. Laslett D, Canbäck B: ARWEN: a program to detect tRNA genes in metazoan mitochondrial nucleotide sequences. Bioinformatics. 2008, 24 (2): 172-175. 10.1093/bioinformatics/btm573.

    Article  CAS  PubMed  Google Scholar 

  47. Wyman SK, Jansen RK, Boore JL: Automatic annotation of organellar genomes with DOGMA. Bioinformatics. 2004, 20 (17): 3252-3255. 10.1093/bioinformatics/bth352.

    Article  CAS  PubMed  Google Scholar 

  48. Lowe TM, Eddy SR: tRNAscan-SE: a program for improved detection of transfer RNA genes in genomic sequence. Nucleic Acids Res. 1997, 25 (5): 955-964. 10.1093/nar/25.5.0955.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  49. Bernt M, Donath A, Jühling F, Externbrink F, Florentz C, Fritzsch G, Pütz J, Middendorf M, Stadler PF: MITOS: improved de novo metazoan mitochondrial genome annotation. Mol Phylogenet Evol. 2013, 69 (2): 313-319. 10.1016/j.ympev.2012.08.023.

    Article  PubMed  Google Scholar 

  50. Darty K, Denise A, Ponty Y: VARNA: Interactive drawing and editing of the RNA secondary structure. Bioinformatics. 2009, 25 (15): 1974-1975. 10.1093/bioinformatics/btp250.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  51. Lin YC, Lu CL, Liu YC, Tang CY: SPRING: a tool for the analysis of genome rearrangement using reversals and block-interchanges. Nucleic Acids Res. 2006, 34 (Web Server issue): W696-W699.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  52. Penn O, Privman E, Ashkenazy H, Landan G, Graur D, Pupko T: GUIDANCE: a web server for assessing alignment confidence scores. Nucleic Acids Res. 2010, 38 (Web Server issue): W23-W28.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  53. Drummond AJ, Suchard MA, Xie D, Rambaut A: Bayesian phylogenetics with BEAUTi and the BEAST 1.7. Mol Biol Evol. 2012, 29 (8): 1969-1973. 10.1093/molbev/mss075.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  54. Rambaut A, Drummond AJ: Tracer v1.5 2007. Available from http://tree.bio.ed.ac.uk/software/tracer

Download references

Acknowledgements

We thank Haley Bridger and Van Tra for help in the research laboratory. This project was supported by grants from the National Center for Research Resources (5P20RR016463-12) and the National Institute of General Medical Sciences (8 P20 GM103423-12) from the National Institutes of Health, the HHMI Undergraduate Science Program, and Bowdoin College.

Author information

Authors and Affiliations

Authors

Corresponding author

Correspondence to Michael F Palopoli.

Additional information

Competing interests

The authors declare that they have no competing interests.

Authors’ contributions

All authors participated in obtaining the sequence data and writing the manuscript. MFP participated in annotating the genomic data, comparing gene order rearrangements, inferring the secondary structures of tRNAs, estimating the age of divergence between species, and wrote the first draft of the manuscript. JE, SM, and DP participated in annotating the protein-coding genes. AS participated in annotating the tRNA genes and inferring the secondary structures of tRNAs. All authors read and approved the final manuscript.

Electronic supplementary material

12864_2014_6923_MOESM1_ESM.png

Additional file 1: Figure S1: Phylogenetic analysis including two species of Solifugae. This phylogeny represents the same type of analysis as represented in Figure 5, except that in this case two species of Solifugae are included. The fossil calibration density is based on the minimum age of 410 mya for the time to the most recent common ancestor of the Acariformes, and amino acid sequences were utilized to minimize the effects of mutation saturation. For this analysis, the resulting estimated mean time to the most recent common ancestor of the two Demodex lineages was 140 mya, with a 95% highest probability density interval of 88 – 191 mya. Similar age estimates were obtained for two other fossil calibration densities, and these estimates were not altered substantially whether or not the Solifugae were included in the analysis (see Table 1 and Figure 5). The Solifugae were always recovered as a clade outside of the Acari. (PNG 902 KB)

12864_2014_6923_MOESM2_ESM.docx

Additional file 2: Table S1: Primer pairs used in PCR experiments for determination of Demodex mitochondrial genome sequences. (DOCX 163 KB)

12864_2014_6923_MOESM3_ESM.docx

Additional file 3: Table S2: Alphabetical list of taxa used for alignments to determine mitochondrial annotations and to estimate divergence time. (DOCX 113 KB)

12864_2014_6923_MOESM4_ESM.docx

Additional file 4: Table S3: Alphabetical list of taxa used for estimation of divergence time based on 18S rRNA genes. (DOCX 83 KB)

Authors’ original submitted files for images

Rights and permissions

This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly credited. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Palopoli, M.F., Minot, S., Pei, D. et al. Complete mitochondrial genomes of the human follicle mites Demodex brevis and D. folliculorum: novel gene arrangement, truncated tRNA genes, and ancient divergence between species. BMC Genomics 15, 1124 (2014). https://0-doi-org.brum.beds.ac.uk/10.1186/1471-2164-15-1124

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://0-doi-org.brum.beds.ac.uk/10.1186/1471-2164-15-1124

Keywords