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

Identification of Wolbachia-responsive microRNAs in the two-spotted spider mite, Tetranychus urticae



The two-spotted spider mite, Tetranychus urticae, is infected with Wolbachia, which have the ability to manipulate host reproduction and fitness. MicroRNAs (miRNAs) are small non-coding RNAs that are involved in many biological processes such as development, reproduction and host-pathogen interactions. Although miRNA was observed to involve in Wolbachia-host interactions in the other insect systems, its roles have not been fully deciphered in the two-spotted spider mite.


Small RNA libraries of infected and uninfected T. urticae for both sexes (in total four libraries) were constructed. By integrating the mRNA data originated from the same samples, the target genes of the differentially expressed miRNAs were predicted. Then, GO and pathway analyses were performed for the target genes. Comparison of libraries showed that Wolbachia infection significantly regulated 91 miRNAs in females and 20 miRNAs in males, with an overall suppression of miRNAs in Wolbachia-infected libraries. A comparison of the miRNA and mRNA data predicted that the differentially expressed miRNAs negatively regulated 90 mRNAs in females and 9 mRNAs in males. An analysis of target genes showed that Wolbachia-responsive miRNAs regulated genes with function in sphingolipid metabolism, lysosome function, apoptosis and lipid transporting in both sexes, as well as reproduction in females.


Comparisons of the miRNA and mRNA data can help to identify miRNAs and miRNA target genes involving in Wolbachia-host interactions. The molecular targets identified in this study should be useful in further functional studies.


Wolbachia is a maternal transmitted endosymbiont that infect about 40% of all arthropod species as well as many filarial nematode species and has profound effects on host biology [1, 2]. Wolbachia manipulates host reproductive systems through a variety of strategies, including cytoplasmic incompatibility (CI), male killing, feminization, and parthenogenesis [3]. In some cases, Wolbachia can also provide benefits [4] and develop mutualistic relationships with their hosts [5, 6]. However, little is known about the molecular mechanisms underlying the above Wolbachia-mediated phenotypes in hosts.

MicroRNAs (miRNAs) are endogenous non-coding small RNAs (~22nt) that play significant roles in regulating a range of cellular processes, including development, differentiation, apoptosis, and immunity [710]. Originally, miRNAs were considered to control gene expression post-transcriptionally by repressing mRNA translation or degrading mRNA in the cytoplasm [11, 12]. Recently, miRNAs was observed to induce gene expression [12, 13] and even regulate pre-mRNA processing in the nucleus, which has a significant role in transcriptional gene silencing [1416]. To date, more than 30,000 miRNAs from over 100 organisms, including two ticks and one mite, have been discovered and deposited in miRNA database (miRBase v. 21.0; [17].

Recent studies have shown that miRNAs are important players in response to bacterial and viral infections in animals and plants [8, 15]. In the mosquito Aedes aegypti, the wMelPop-CLA strain of Wolbachia significantly alters the miRNA profile of A. aegypti, thereby manipulating a broad array of host genes to facilitate self-maintenance [13, 16, 1820]. These observations indicate that host miRNAs may play essential role in Wolbachia-host interactions.

The two-spotted spider mite Tetranychus urticae is a cosmopolitan agricultural pest that feeds on hundreds of plant species [21]. We previously demonstrated that Wolbachia induces strong CI and increases host fecundity in T. urticae[22]. In order to explore whether and how host miRNAs are involved in T. urticae-Wolbachia interactions, we constructed four small RNA (sRNA) libraries representing female mites infected (FI) and uninfected (FU) with Wolbachia, and male mites infected (MI) and uninfected (MU) with Wolbachia. These libraries also offer the opportunity to identify protein-encoding genes that were differentially expressed in response to Wolbachia infection. By integrating the miRNA and mRNA data, the target genes of the differentially expressed miRNAs were identified. Furthermore, several altered miRNAs associated with reproduction were analyzed. Our results provide new insights on how Wolbachia affects its hosts through mediating host miRNAs.

Results and discussion

Identification of known and novel miRNAs in T. urticae

Four sRNA libraries were constructed with a total of about 49 million raw reads. After sequencing process, 48 million clean reads were generated (Table 1). About 80% of the sRNAs were in the range 18-25nt and the peak sizes were 21nt for the FI and FU libraries and 22nt for the MI and MU libraries (Figure 1). The sRNAs were divided into ten classes, such as known miRNAs, novel miRNAs and tRNAs (Figure 2). The percentages of known and novel miRNAs were highest in the MU library (23.73% and 10.17%, respectively), but lowest in the FI library (1.55% and 0.66%, respectively).

Table 1 Statistics of small RNA sequences of the four libraries
Figure 1
figure 1

The length distribution of small RNA (sRNA) in the four libraries. FI, female infected; FU, female uninfected; MI, male infected; MU, male uninfected; nt, nucleotides.

Figure 2
figure 2

Classification of sRNAs in the four libraries. ‘rRNAetc’ includes rRNA, snRNA (small nuclear RNA) and snoRNA (small nucleolar RNA); others, unannotated sRNA.

A total of 83 known miRNAs and 112 novel miRNAs were identified in the four libraries (Table 1). Most of the predicted novel miRNAs had the guide strand retained and lost the passenger strand (miRNA*) (data not shown), so only the guide strand of a novel miRNA was analyzed further in this study. These novel miRNAs showed all characteristic signatures of a miRNA as described previously [23]. The predicted secondary structures of nine novel miRNAs were shown in Additional file 1: Figure S1. About 45% of novel miRNAs were 25nt and 77% of the novel miRNAs were below 100 readcount (Additional file 2: Table S2). According to miFam.dat ( and Rfam [24], forty conserved miRNAs and five novel miRNAs (novel_10, novel_62, novel_65, novel_67 and novel_112) were classified into 31 miRNA families (Additional file 2: Table S2).

There were significant overlaps on identified miRNA between different groups (Figure 3A). Specifically, 193 miRNAs were identified in female libraries. Among them, 5 miRNAs were detected only in FI and 29 miRNAs were found only in FU, while 159 miRNAs were found in both libraries (Figure 3A). 189 miRNAs were identified in the male libraries, with 158 miRNAs shared by MI and MU libraries. Similar to the observation that less miRNAs present in FI than FU, the MI and MU libraries had 4 and 25 library-specific miRNAs, respectively (Figure 3A). These results indicate that the presence of Wolbachia may inhibit the synthesis of miRNAs in both female and male mites.

Figure 3
figure 3

Venn chart of miRNA expression in the four libraries. A. Unique and commonly expressed miRNAs among four libraries (FI, FU, MI and MU). Blue oval represents FU; magenta oval represents MU; green oval represents MI; red oval represents FI. B. Unique and commonly differentially expressed miRNAs among four comparisons (FI vs FU, MI vs MU, FI vs MI, FU vs MU). Blue oval represents FI vs FU; magenta oval represents MU vs FU; green oval represents MI vs FI; red oval represents FI vs MU.

Abundance of miRNAs

miRNAs with more than 1,000 transcripts per million (TPM) were classified as abundant while those with less than 10 TPM were classified as rare. The 20 most abundant miRNAs in each of the libraries (accounting for ~90% of total miRNA reads) are listed in Table 2. Three of them (novel_85, novel_1 and novel_2, shown in bold) were novel miRNAs. The number of rare miRNAs in FI (52) was twice the number of FU (26), while showed no significant difference between MI (91) and MU (90), indicating that many miRNAs may be down-regulated in infected females (See Additional file 2: Table S2).

Table 2 Top 20 most abundant miRNAs expressed in the four libraries (readcount were shown)

Expression analysis of Wolbachia-responsive miRNAs

In females, 32 miRNAs were observed to be up-regulated and 59 down-regulated in the Wolbachia-infected line compared to the Wolbachia-uninfected line (Additional file 3: Table S3, Figure 4), while 11 miRNAs were up-regulated and 9 down-regulated in the Wolbachia-infected line in males (Additional file 4: Table S4, Figure 4). Among the up-regulated miRNAs, tur-miR-1-5p and tur-miR-2-3p were strongly expressed in MI and FI, respectively, while they were not among the top 20 miRNAs in MU and FU (Table 2). tur-miR-2-3p belongs to the mir-2 miRNA family, which is widespread in invertebrates, and is the largest family of miRNAs in the model species Drosophila melanogaster[17, 25]. Based on their predicted targets, mir-2 miRNAs have been suggested to have a role in neural development and maintenance [26]. Leaman et al.[27] showed that the miR-2 family regulates cell survival by translational repression of proapoptotic factors. Thus, the up-regulation of tur-miR-2-3p in Wolbachia-infected females may be related to changes in neural or apoptotic processes.

Figure 4
figure 4

Expression of miRNAs in FI vs FU and MI vs MU. The x-axis and y-axis show the expression levels of miRNAs in each library. Each point represents a miRNA. Red points represent the up-regulated miRNAs. Blue points indicate equally-expressed miRNAs. Green points represent down-regulated miRNAs. FU or MU was set as control.

Thirteen differentially expressed miRNAs were common between the female and male (overlap of blue and red ovals in Figure 3B). Among them, novel_23 was up-regulated in Wolbachia-infected males but down-regulated in Wolbachia-infected females (Table 3). Four other miRNAs were up-regulated in Wolbachia-infected lines, and eight miRNAs were down-regulated in Wolbachia-infected lines (Table 3). Interestingly, novel_126, novel_23 and novel_32 were detected only in Wolbachia-uninfected females (See Additional file 3: Table S3). Similarly, novel_105, novel_109, novel_126, novel_16 and novel_32 were found only in Wolbachia-uninfected males (Additional file 4: Table S4).

Table 3 Predicted target genes of the 13 common Wolbachia -responsive miRNAs in the female and male comparisons

We also identified 74 female-specific differentially expressed miRNAs, including 18 substantially (decreased by more than 16-fold) down-regulated miRNAs and 5 substantially up-regulated miRNAs (increased by more than 16-fold) (Additional file 3: Table S3). In Wolbachia-infected females, the most down-regulated miRNA was novel_126, which had a log2 (fold-change) value of 9.603. In contrast, novel_59 was most up-regulated with a log2 (fold-change) value of 7.332 (Additional file 3: Table S3). Furthermore, several miRNAs like novel_59, novel_117, novel_45, novel_120 and tur-miR-5729a-3p were found only in infected females (Additional file 3: Table S3).

miRNA expression patterns

To further validate the expression patterns of miRNAs identified in this work, we measured the expression levels of nine differential miRNAs (novel_126, novel_32, novel_105, novel_36, novel_16, novel_89, novel_109, novel_43 and one novel_23) shared by the female and male and five female-specific differential miRNAs (novel_101, novel_18, novel_80, novel_21 and novel_87) using quantitative real time polymerase chain reaction (qRT-PCR). Except for novel_23, all other 8 miRNAs were significantly (P < 0.05) down-regulated in Wolbachia-infected female and male lines, which was in agreement with our sequencing data (Figure 5A,B). novel_126, novel_32, novel_36 and novel_89 were decreased by more than five times (P < 0.001). Three of the five female-specific differential miRNAs were significantly down-regulated in infected females, which was consistent with our sequencing data (Figure 5C). As a whole, these results supported that the Illumina sequencing data accurately reflected the abundance of the miRNAs.

Figure 5
figure 5

Expression patterns of differentially expressed miRNAs and one target gene (tetur0906680) in the female and male comparisons. A. and B. Nine common differentially expressed miRNAs in the comparisons FI vs FU (control) and MI vs MU (control) identified by using Illumina small RNA deep sequencing (log2 fold-change) and validated by qRT-PCR. C. Validation of five female-specific differentially expressed miRNAs by qRT-PCR. B and C. Error bars represent standard errors of averages from three biological replicates (NS: not significant; *P < 0.05; **P < 0.01; ***P < 0.001; Fisher's LSD test). D. qRT-PCR and Illumina small RNA deep sequencing both showed that candidate target gene tetur09g06680 of novel_16 was up-regulated in Wolbachia-infected mites. All the experiments had three biological replicates.

Target gene prediction of differentially expressed miRNAs

To further determine the biological functions of the differentially expressed miRNAs, we integrated the miRNA and the mRNA transcriptome data to predict the candidate target genes. At first, miRanda was used to predict the target genes of the 195 miRNAs, which resulted in the identification of 41,914 miRNA-target pairs. Then, the predicted target genes of differentially expressed miRNAs were aligned with differentially expressed mRNAs. This approach revealed possible interactions between 115 mRNAs and 79 miRNAs in the female comparison (Additional file 5: Table S5), and between 92 mRNAs and 15 miRNAs in the male comparison (Additional file 6: Table S6). In addition, 65 differential miRNAs found in the female comparison were predicted to negatively regulate 90 candidate target genes (Additional file 3: Table S3). While in the male comparison, 9 candidate target genes were suggested to be negatively regulated by the 11 differentially expressed miRNAs (Additional file 4: Table S4).

The biological functions of the predicted target genes were then analyzed using the gene ontology (GO) annotations of the T. urticae genes. The GO annotation enrichment results showed that genes related to binding, catalytic activity, and metabolic and cellular processes were most enriched in the two comparisons (Figure 6A and B). In both female and male mites, GOs corresponding to the up-regulated genes included establishment of localization, membrane and transporter activity (data not shown). However, several GOs, such as death, enzyme regulator activity and structural molecular activity, were exclusively enriched in the female mites (Figure 6A and B).Pathway analysis resulted in observation of 28 and 5 different pathways corresponded to the differentially expressed genes in the female and male comparison, respectively. Notably, the five KEGG pathways (lysosome, sulfur metabolism, glycan degradation, sphingolipid metabolism and metabolic pathways) and their related three candidate genes (tetur08g05010, tetur09g06680 and tetur03g03470) of the male comparison were also highlighted in the female comparison. The 20 most enriched KEGG pathways of females included the degradation of valine, leucine and isoleucine, phenylalanine metabolism, lysosome function and so on (Figure 7). The target gene tetur09g06680 was the only enriched gene in both female’s and male’s KEGG pathways which may be negatively regulated by a differential miRNA (novel_16). Sequencing and qRT-PCR results showed that tetur09g06680 was up-regulated in both infected female and male mites (Figure 5D). The up-regulation of tetur09g06680 indicated that novel_16 potentially represses the expression of this gene.

Figure 6
figure 6

GO analysis results of the target genes of the two comparisons. A. FI vs FU (control); B. MI vs MU (control). The x-axis is the GO category and the y-axis is the percent and number of genes.

Figure 7
figure 7

The 20 most enriched KEGG pathways based on target genes of differentially expressed miRNAs in the female comparison. The x-axis shows the rich factor. The y-axis shows the pathway names. The size of each point represents the number of genes enriched in a particular pathway. The bigger the value of rich factor and the smaller the value of Q-value indicate the degree of enrichment is more significant.

Wolbachia-responsive miRNAs are related to regulation of apoptosis in mites of both sexes

The target genes that may be negatively regulated by the 13 common differentially expressed miRNAs in the female and male comparisons were presented in Table 3. In general, the target genes of some miRNAs revealed multiply. Despite the similarities in miRNAs regulation, we found that only novel_16 and novel_89 targeted the same genes in both comparisons (Table 3). The miRNA novel_16 targets two genes potentially associated with apoptosis: tetur09g06680 encoding glucosylceramidase and tetur26g00300 encoding protein lifeguard 1, which were up-regulated in Wolbachia infected mites with more than 4-fold (Table 3). The gene tetur09g06680 is involved in several KEGG pathways, including lysosome function, sphingolipid metabolism and glycan degradation (data not shown). Deficiency of GCase causes several diseases in mice, including severe neurodegeneration and apoptosis in the brain [28, 29]. Therefore, the observed up-regulation of tetur0906680 in both infected lines may indicate its role in inhibiting apoptosis. Moreover, some other miRNAs, such as novel_100 (targets tetur03g00110), novel_84 (targets tetur26g00300), novel_92 (targets tetur03g00110), novel_50 (targets tetur26g00300), may also regulate apoptosis (Additional file 3: Table S3). The expression level of tetur03g00110, encoding Protein lethal (2) essential for life, was elevated in Wolbachia-infected females. Previous studies have proved that Wolbachia can inhibit apoptosis in several ways. For example, Wolbachia surface protein (Wsp) inhibits apoptosis of Polymorphonuclear (PMN) cells in vitro[30]. In the parasitic wasp Asobara tabida, Wolbachia helps its host to avoid becoming reproductively sterile by repressing massive apoptosis of nurse cells in the ovary [31]. Similarly, Wolbachia infection leads to a decrease in programmed cell death in the germarium of D. mauritiana[32]. However, a virulent Wolbachia strain wMelPop increase the frequency of apoptosis in the female Drosophila[33]. Regardless of the strategy that Wolbachia uses to interfere with host, the apoptosis inhibition will in turn benefit to maintenance of Wolbachia itself. The down-regulation of the anti-apoptotic miRNAs (novel_16, novel_100, novel_84, novel_92 and novel_50) reinforces the importance of investigating the apoptotic signaling cascade upon Wolbachia infection (Additional file 3: Table S3 and Additional file 4: Table S4).

Wolbachia-responsive miRNAs may regulate reproduction in female mites

Since Wolbachia causes strong CI in this spider mite population, an important aim of our study was to clarify the underlying mechanism of Wolbachia-induced CI. Some remarkable features observed in CI embryos include delayed paternal nuclear envelope breakdown and activation of Cdk1 [34], a failure of the maternal histones H3.3 to deposit in the paternal genome, and slowed replication of sperm DNA [35]. In this study, we found several target genes that may be related to females’ reproduction. For instance, target gene tetur13g00510 (targeted by novel_126) codes for Histone H2B was induced by 6.9-fold in Wolbachia-infected females (Table 3). Histone play important roles in transcription regulation, DNA repair, DNA replication and chromosomal stability [36]. Target gene tetur11g00940 of novel_36 was exclusively expressed in Wolbachia-uninfected females (Additional file 5: Table S5), which codes for meiosis arrest female protein 1 and involved in oogenesis (Table 3). Target gene tetur11g00600 and tetur01g12840 (targeted by novel_105) code for cuticle protein, which may be related to egg formation (Table 3). In addition, target gene tetur04g05770 (potentially regulated by novel_80) codes for sorbitol dehydrogenase, an enzyme which possibly play a role in sperm motility by providing a source of energy for sperm [37] (Additional file 3: Table S3). Taken together, it supports the hypothesis that Wolbachia can manipulate female mites’ reproduction via host miRNAs and further studies are needed to test this hypothesis.

Other function genes mediated by Wolbachia-responsive miRNAs

Among the seven male-specific differential miRNAs, only two had target genes with description (Additional file 4: Table S4). The target gene exclusively found in the male comparison was tetur09g04610 (targeted by novel_109), which codes for multidrug resistance-associated protein 1 (Table 3). The target gene of tur-miR-9-3p (tetur09g04720) codes for apolipoprotein D, which is related to lipid transporting. We also found target genes related to drug resistance (gene tetur10g01260 potentially regulated by novel_80) and apolipophorins biogenesis (gene tetur20g01230 potentially regulated by novel_100) were up-regulated in the infected females (Additional file 3: Table S3), which indicated that these two biological processes may be regulated by different miRNA-target pairs in the female and male mites. Apolipophorins may be up-regulated in female mites because they are essential for reproduction, but they may also be up-regulated by Wolbachia because Wolbachia needs them to produce lipopolysaccharide, the major component of its outer membrane [38]. Another observation worth to note here was that Wolbachia also enhances the expression of the genes that codes for riboflavin transporter in female mites (Additional file 3: Table S3). Riboflavin that Wolbachia provided to nematodes was found to be crucial for worm health and fertility (6), and in cultured mosquito cells, depletion of host cell riboflavin decreased Wolbachia abundance [39]. Taken together, these results raise the possibility that Wolbachia and T. urticae have a symbiotic relationship.

Sex-specific interaction mechanisms

The divergences of the differentially expressed miRNAs and target genes found in the female and male comparisons suggest that Wolbachia used different strategies to regulate miRNA and mRNA in females and males. We previously showed that Wolbachia mainly localized in the ovary of female mites and its density increased during the reproduction process [22]. In contrast, the density of Wolbachia in male mites decreased with age [22]. These phenomena indicated that Wolbachia replicate massively in female reproductive tissues to ensure its successful vertical transmission. In the present study, Wolbachia infection led to the down-regulation of 65% of the differential miRNAs and up-regulation of 69% of the predicted target genes in the female mites (Additional file 3: Table S3). GO analysis showed that the up-regulate genes were related to establishment of localization, membrane, structural molecular activity and transporter activity (data not shown), which also suggested that Wolbachia replication is increasing in young females. Besides Wolbachia replication, the reproductive tissues of female mites probably experience high metabolic activity during egg production. That could be two reasons why there are more specific- and up- regulated mRNAs in females. Some studies suggested that Wolbachia interacts with female hosts in a broader manner than it does with male hosts. In some cases, Wolbachia acts as an obligate mutualism, i.e., it is required for normal reproduction of the host [5, 40]. In other cases, Wolbachia was found to benefit to its female host by increasing fecundity [4, 22]. Sex-specific regulations of miRNAs were also found in flour beetles, in which 245 (54%) of the miRNAs exhibited gender-specific expression patterns upon exposure to environmental stress [41]. This discrepancy in gender- specific miRNA expression is consistent with Bateman’s principle, which states that males gain fitness by increasing their mating success whereas females gain fitness through increasing longevity because their reproductive effort is higher [42]. If the female and male mites gain fitness through different ways following Wolbachia infection, the expression of miRNAs and mRNAs may be diversely influenced.


Four libraries of young female and male mites with and without Wolbachia infection, were constructed, amplified and sequenced. As a result, we detected 91 and 20 miRNAs exhibiting differential expression in response to Wolbachia infection in female and male mites, respectively. There was an overall decrease of miRNAs in Wolbachia-infected lines. Integrating the miRNA and mRNA data uncovered many target genes, and enabled us to develop new hypotheses for Wolbachia-regulated reproduction and apoptosis mechanisms. The discovery of these Wolbachia-responsive miRNAs and their target genes provides a basis for understanding Wolbachia-host interaction in a CI phenotype host. In addition, we described a complex interaction network of miRNAs and target mRNAs that will encourage future studies to examine the contribution of the specific newly identified miRNAs on the regulation of biological processes in response to Wolbachia infection in more detail.

Eventually, this study reports the discovery of sexually differentially expressed miRNAs and miRNA- target gene pairs. Wolbachia appears to affect a large number of cellular processes in female mites. Those findings lay the foundation for future studies to identify miRNAs or genes responsible for not only the Wolbachia maintenance, but also the sexually host-endosymbiont interaction mechanism.


Two-spotted spider mite

Mites used in this study were originally collected from Hohhot, Inner Mongolia, northeast China in July 2010 and reared on leaves of the common bean (Phaseolus vulgaris L.) placed on a water-saturated sponge mat in Petri dishes (dia. 9) at 25 ± 1°C, 60% r. h. and under L16-D8 conditions. 100% infected and 100% uninfected Wolbachia lines with identical genetic backgrounds were prepared according to a method described by Xie et al.[43]. Through PCR assays, neither line was found to be infected with Cardinium (CLO-f1: 5’-GGAACCTTACCTGGGCTAGAATGTATT-3’, CLO-r1: 5’-GCCACTGTCTTCAAGCTCTACCAAC-3’) or Rickettsia (R1: 5’-GCTCTTGCAACTTCTATGTT-3’, R2: 5’-CATTGTTCGTCAGGTTGGCG-3’) [44], which can manipulate host reproduction. Within the two lines, 1–3 day old adult females (Turt_FI and Turt_FU, ‘I’ indicates infection and ‘U’ indicates uninfection) and 1 day-old adult virgin males (Turt_MI and Turt_MU) were respectively collected. The samples were stored in liquid nitrogen until required for RNA isolation.

Small RNA library construction for Illumina sequencing

Total RNA was extracted by using Trizol reagent (Invitrogen catalog no. 15596–026) according to the manufacturer’s protocol with small modifications. Total RNA (>3 μg) with good quality was used to construct a sRNA library for each sample by using a TruSeq small RNA Sample Pre Kit (Illumina). Briefly, total RNA was ligated with 5’ and 3’ adaptors followed by reverse transcription using RT primers. After PCR amplification of the cDNAs, amplified PCR products within 130–160 bp were purified on an 8% polyacrylamide gel (100 V, 80 min). The library quality was assessed on the Agilent Bioanalyzer 2100 system using DNA High Sensitivity Chips and then sequenced on a HiSeq2000 sequencer (Illumina).

Bioinformatics analysis

After Illumina sequencing, raw data were processed through Novogene Company’s Perl and Python scripts. In this step, clean data were obtained by removing reads containing more than three N (undetermined bases), with 5’ adapter contaminants, without 3’ adapter or the insert tag, containing poly A or T or G or C and low quality reads from the raw data. Then, sRNAs with lengths of 18–30 nt were selected for further analyses. In the alignment and annotation procedure, we used the following priority rule: known miRNA > rRNA > tRNA > snRNA > snoRNA > repeat > gene > novel miRNA so that every unique sRNA mapped to only one annotation. Briefly, the first step was mapping the sRNA tags to the T. urticae genome sequence (; released in Sep, 2009) by Bowtie [45] without mismatch to analyze their expression and distribution on the reference sequence. Next, the mappable sRNA tags were aligned to the miRNA precursor of T. urticae in the miRNA database (miRBase v. 20.0; released in June, 2013) to obtain the known miRNA count. Then, tags originating from rRNAs, tRNAs, snRNAs, and snoRNAs were removed by mapping the remained sRNA tags to the ncRNA annotation database of T. urticae ( Tags originating from repeat sequences were filtered by using a repeat sequences database (, and tags originating from protein-coding genes were discarded by mapping to the exon and intron of mRNAs of T. urticae. Finally, novel miRNAs were predicted by exploring the secondary structure, the Dicer cleavage site and the minimum free energy of the former unannotated sRNA tags which could be mapped to the reference sequence by integrating two available software miREvo [46] and mirdeep2 [47].

Comparison of miRNA libraries

The readcounts of miRNAs were transformed into TPM (transcript per million) through the following criteria: Normalization formula: Normalized expression = (mapped readcount/Total reads)*1000000 [48] and then analyzed the abundance of miRNAs. Since we didn’t set up biological replicates for each sample, when analyzing differentially expressed miRNAs between libraries, we first normalized the readcount data by using TMM (trimmed mean of M values) [49], then used the DEGseq R package to analyze the differences [50]. The P-value was adjusted using the q value [51]. q Value <0.01 and |log2 (fold-change)| > 1 was set as the threshold for significant differential expression by default. We compared the expression level of miRNAs between FI and FU (control), MI and MU (control), MI and FI (control), MU and FU (control).

Target prediction

The 3’ UTR annotation information originated from the genome database of T. urticae ( was used by miRanda [52] to predict target genes of miRNAs. As we used the same sample to get the RNA transcriptome data (Transcriptome raw data are available in the ArrayExpress database at under accession number E-MTAB-2491), we integrated the target genes of differentially expressed miRNAs with the differentially expressed genes found by comparing the four transcriptomes to obtain more precise miRNA-mRNA correlation information. Since most miRNAs were found to repress mRNA translation or induce mRNA degradation, the datasets in our study were analyzed through intersecting elements of target genes of 1) differentially expressed miRNAs and differentially expressed mRNAs, 2) differentially down-expressed miRNAs and differentially up-expressed mRNAs and 3) differentially up-expressed miRNAs and differentially down-expressed mRNAs. The datasets obtained above were used for GO enrichment analysis [53] and KEGG enrichment analysis [54]. KOBAS [55] software was used to test the statistical enrichment of the target gene candidates in the KEGG pathways.

miRNA expression patterns

miRNA expression patterns were established in 1-day-old virgin Wolbachia-infected and Wolbachia-free females and males. The experiment was run on an ABI 7500 thermal cycler (Applied Biosystems) with three technical replicates and three biological replicates. About 500 females and 1000 males were prepared for each biological replicate. For qRT-PCR, RNA was extracted with a miRNeasy extraction kit (Qiagen catalog no. 217004) together with RNase-Free DNase Set (Qiagen catalog no. 79254). RNA was extracted according to the manufacturer’s protocol except that after grinding with a homogenizer, the homogenate was centrifuged at 12,000 × g to remove the precipitate. RNA quality was estimated with a Nanodrop spectrophotometer and by agarose gel electrophoresis. The expression of miRNA was checked by using SYBR PrimeScript miRNA RT-PCR Kit (Takara catalog no. RR716). Total RNA (1.5 μg) was reverse transcribed and 30 ng cDNA product was put into a 20 μL quantification system. The primers used in this study are shown in Additional file 7: Table S1. The reactions were incubated at 95°C for 30 s, followed by 40 cycles of 95°C for 5 s, 60°C for 34 s. A dissociation curve was obtained to ensure that only one product was amplified after the amplification phase. The 2-ΔΔCt method for relative quantification of gene expression was used to determine the level of miRNA expression and U6 snoRNA was used for normalizing the data.

Candidate target gene expression pattern

The same RNA used for miRNA expression analysis was used for target gene expression analysis. Total RNA (1.5 μg) was reverse transcribed with SuperScript III Reverse Transcriptase (Invitrogen catalog no.18080-044) according to the instructions. SYBR Premix Ex Taq™ II (Takara catalog no. RR820A) was used to quantify the target genes by putting 30 ng cDNA products into 20 μL quantification system. The qPCR cycling parameters included 95°C for 30 s, and 40 cycles of 95°C for 5 s, and 60°C for 34 s. At the end of the PCR reaction, a melting curve was generated. β-Tubulin was used as an endogenous reference gene. The experiment was conducted three times.

Availability of supporting data

Small RNA sequencing raw data are available in the ArrayExpress database ( under accession number E-MTAB-2606.


  1. Zug R, Hammerstein P: Still a host of hosts for Wolbachia: analysis of recent data suggests that 40% of terrestrial arthropod species are infected. PLoS One. 2012, 7: e38544-10.1371/journal.pone.0038544.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  2. Ferri E, Bain O, Barbuto M, Martin C, Lo N, Uni S, Landmann F, Baccei SG, Guerrero R, de Souza LS, Bandi C, Wanji S, Diagne M, Casiraghi M: New insights into the evolution of Wolbachia infections in filarial nematodes inferred from a large range of screened species. PLoS One. 2011, 6: e20843-10.1371/journal.pone.0020843.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  3. Werren JH, Baldo L, Clark ME: Wolbachia: master manipulators of invertebrate biology. Nat Rev Microbiol. 2008, 6: 741-751. 10.1038/nrmicro1969.

    Article  CAS  PubMed  Google Scholar 

  4. Weeks AR, Turelli M, Harcombe WR, Reynolds KT, Hoffmann AA: From parasite to mutualist: rapid evolution of Wolbachia in natural populations of Drosophila. PLoS Biol. 2007, 5: e114-10.1371/journal.pbio.0050114.

    Article  PubMed Central  PubMed  Google Scholar 

  5. Hosokawa T, Koga R, Kikuchi Y, Meng XY, Fukatsu T: Wolbachia as a bacteriocyte-associated nutritional mutualist. Proc Natl Acad Sci U S A. 2010, 107: 769-774. 10.1073/pnas.0911476107.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  6. Li Z, Carlow CK: Characterization of transcription factors that regulate the type IV secretion system and riboflavin biosynthesis in Wolbachia of Brugia malayi. PLoS One. 2012, 7: e51597-10.1371/journal.pone.0051597.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  7. Smibert P, Lai EC: A view from Drosophila: Multiple biological functions for individual microRNAs. Semin Dev Biol. 2010, 21: 745-753. 10.1016/j.semcdb.2010.03.001.

    Article  CAS  Google Scholar 

  8. Asgari S: Role of microRNAs in insect host-microorganism interactions. Front Physiol. 2011, 2: 48-

    Article  PubMed Central  PubMed  Google Scholar 

  9. Wang XH, Aliyari R, Li WX, Li HW, Kim K, Carthew R, Atkinson P, Ding S-W: RNA interference directs innate immunity against viruses in adult Drosophila. Sci Signal. 2006, 312: 452-

    CAS  Google Scholar 

  10. Lourenço AP, Guidugli-Lazzarini KR, Freitas FC, Bitondi MM, Simões ZL: Bacterial infection activates the immune system response and dysregulates microRNA expression in honey bees. Insect Biochem Mol Biol. 2013, 43: 474-482. 10.1016/j.ibmb.2013.03.001.

    Article  PubMed  Google Scholar 

  11. Bartel DP: MicroRNAs: target recognition and regulatory functions. Cell. 2009, 136: 215-233. 10.1016/j.cell.2009.01.002.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  12. Asgari S: MicroRNA functions in insects. Insect Biochem Mol Biol. 2013, 43: 388-397. 10.1016/j.ibmb.2012.10.005.

    Article  CAS  PubMed  Google Scholar 

  13. Hussaina M, Frentiua FD, Moreiraa LA, O'Neill SL, Asgaria S: Wolbachia uses host microRNAs to manipulate host gene expression and facilitate colonization of the dengue vector Aedes aegypti. Proc Natl Acad Sci U S A. 2011, 108: 9250-9255. 10.1073/pnas.1105469108.

    Article  Google Scholar 

  14. Fagegaltier D, Bougé AL, Berry B, Poisot É, Sismeiro O, Coppée JY, Théodore L, Voinnet O, Antoniewski C: The endogenous siRNA pathway is involved in heterochromatin formation in Drosophila. Proc Natl Acad Sci U S A. 2009, 106: 21258-21263. 10.1073/pnas.0809208105.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  15. Fabian MR, Sonenberg N, Filipowicz W: Regulation of mRNA translation and stability by microRNAs. Annu Rev Biochem. 2010, 79: 351-379. 10.1146/annurev-biochem-060308-103103.

    Article  CAS  PubMed  Google Scholar 

  16. Hussain M, O'Neill SL, Asgari S: Wolbachia interferes with the intracellular distribution of Argonaute 1 in the dengue vector Aedes aegypti by manipulating the host microRNAs. RNA Biol. 2013, 10: 1868-1875. 10.4161/rna.27392.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  17. Kozomara A, Griffiths-Jones S: miRBase: annotating high confidence microRNAs using deep sequencing data. Nucleic Acids Res. 2014, 42: D68-D73. 10.1093/nar/gkt1181.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  18. Osei-Amo S, Hussain M, O'Neill SL, Asgari S: Wolbachia-induced aae-miR-12 miRNA negatively regulates the expression of MCT1 and MCM6 genes in Wolbachia-infected mosquito cell line. PLoS One. 2012, 7: e50049-10.1371/journal.pone.0050049.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  19. Zhang G, Hussain M, O’Neill SL, Asgari S: Wolbachia uses a host microRNA to regulate transcripts of a methyltransferase, contributing to dengue virus inhibition in Aedes aegypti. Proc Natl Acad Sci U S A. 2013, 110: 10276-10281. 10.1073/pnas.1303603110.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  20. Mayoral JG, Etebari K, Hussain M, Khromykh AA, Asgari S: Wolbachia infection modifies the profile, shuttling and structure of microRNAs in a mosquito cell line. PLoS One. 2014, 9: e96107-10.1371/journal.pone.0096107.

    Article  PubMed Central  PubMed  Google Scholar 

  21. Jeppson LR, Keifer HH, Baker EW: Mites Injurious to Economic Plants. 1975, Berkeley: University of California Press

    Google Scholar 

  22. Zhao DX, Zhang XF, Chen DS, Zhang YK, Hong XY: Wolbachia-host interactions: host mating patterns affect Wolbachia density dynamics. PLoS One. 2013, 8: e66373-10.1371/journal.pone.0066373.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  23. Ambros V, Bartel B, Bartel DP, Burge CB, Carrington JC, Chen X, Dreyfuss G, Eddy SR, Griffiths-Jones S, Marshall M, Matzke M, Ruvkun G, Tuschl T: A uniform system for microRNA annotation. RNA. 2003, 9: 277-279. 10.1261/rna.2183803.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  24. Griffiths-Jones S, Moxon S, Marshall M, Khanna A, Eddy SR, Bateman A: Rfam: annotating non-coding RNAs in complete genomes. Nucleic Acids Res. 2005, 33: D121-124.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  25. Lagos-Quintana M, Rauhut R, Lendeckel W, Tuschl T: Identification of novel genes coding for small expressed RNAs. Science. 2001, 294: 853-858. 10.1126/science.1064921.

    Article  CAS  PubMed  Google Scholar 

  26. Marco A, Hooks K, Griffiths-Jones S: Evolution and function of the extended miR-2 microRNA family. RNA Biol. 2012, 9: 242-248. 10.4161/rna.19160.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  27. Leaman D, Chen PY, Fak J, Yalcin A, Pearce M, Unnerstall U, Marks DS, Sander C, Tuschl T, Gaul U: Antisense-mediated depletion reveals essential and specific functions of microRNAs in Drosophila development. Cell. 2005, 121: 1097-1108. 10.1016/j.cell.2005.04.016.

    Article  CAS  PubMed  Google Scholar 

  28. Enquist IB, Bianco CL, Ooka A, Nilsson E, Månsson JE, Ehinger M, Richter J, Brady RO, Kirik D, Karlsson S: Murine models of acute neuronopathic Gaucher disease. Proc Natl Acad Sci U S A. 2007, 104: 17483-17488. 10.1073/pnas.0708086104.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  29. Jmoudiak M, Futerman AH: Gaucher disease: pathological mechanisms and modern management. Br J Haematol. 2005, 129: 178-188. 10.1111/j.1365-2141.2004.05351.x.

    Article  CAS  PubMed  Google Scholar 

  30. Bazzocchi C, Comazzi S, Santoni R, Bandi C, Genchi C, Mortarino M: Wolbachia surface protein (WSP) inhibits apoptosis in human neutrophils. Parasite Immunol. 2007, 29: 73-79.

    Article  CAS  PubMed  Google Scholar 

  31. Kremer N, Voronin D, Charif D, Mavingui P, Mollereau B, Vavre F: Wolbachia interferes with ferritin expression and iron metabolism in insects. PLoS Pathog. 2009, 5: e1000630-10.1371/journal.ppat.1000630.

    Article  PubMed Central  PubMed  Google Scholar 

  32. Fast EM, Toomey ME, Panaram K, Desjardins D, Kolaczyk ED, Frydman HM: Wolbachia enhance Drosophila stem cell proliferation and target the germline stem cell niche. Science. 2011, 334: 990-992. 10.1126/science.1209609.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  33. Zhukova MV, Kiseleva E: The virulent Wolbachia strain wMelPop increases the frequency of apoptosis in the female germline cells of Drosophila melanogaster. BMC Microbiol. 2012, 12 Suppl 1: S15-

    Article  PubMed  Google Scholar 

  34. Tram U, Sullivan W: Role of delayed nuclear envelope breakdown and mitosis in Wolbachia-induced cytoplasmic incompatibility. Science. 2002, 296: 1124-1126. 10.1126/science.1070536.

    Article  CAS  PubMed  Google Scholar 

  35. Landmann F, Orsi GA, Loppin B, Sullivan W: Wolbachia-mediated cytoplasmic incompatibility is associated with impaired histone deposition in the male pronucleus. PLoS Pathog. 2009, 5: e1000343-10.1371/journal.ppat.1000343.

    Article  PubMed Central  PubMed  Google Scholar 

  36. Peterson CL, Laniel MA: Histones and histone modifications. Curr Biol. 2004, 14: R546-R551. 10.1016/j.cub.2004.07.007.

    Article  CAS  PubMed  Google Scholar 

  37. Cao W, Aghajanian HK, Haig-Ladewig LA, Gerton GL: Sorbitol can fuel mouse sperm motility and protein tyrosine phosphorylation via sorbitol dehydrogenase. Biol Reprod. 2009, 80: 124-133. 10.1095/biolreprod.108.068882.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  38. Wu M, Sun LV, Vamathevan J, Riegler M, Deboy R, Brownlie JC, McGraw EA, Martin W, Esser C, Ahmadinejad N: Phylogenomics of the reproductive parasite Wolbachia pipientis wMel: a streamlined genome overrun by mobile genetic elements. PLoS Biol. 2004, 2: e69-10.1371/journal.pbio.0020069.

    Article  PubMed Central  PubMed  Google Scholar 

  39. Fallon AM, Baldridge GD, Carroll EM, Kurtz CM: Depletion of host cell riboflavin reduces Wolbachialevels in cultured mosquito cells. In Vitro Cell Dev Biol Anim. 2014, 50 (8): 707-713. 10.1007/s11626-014-9758-x.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  40. Darby AC, Armstrong SD, Bah GS, Kaur G, Hughes MA, Kay SM, Koldkjær P, Rainbow L, Radford AD, Blaxter ML: Analysis of gene expression from the Wolbachia genome of a filarial nematode supports both metabolic and defensive roles within the symbiosis. Genome Res. 2012, 22: 2467-2477. 10.1101/gr.138420.112.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  41. Freitak D, Knorr E, Vogel H, Vilcinskas A: Gender- and stressor-specific microRNA expression in Tribolium castaneum. Biol Lett. 2012, 8: 860-863. 10.1098/rsbl.2012.0273.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  42. Rolff J: Bateman's principle and immunity. Proc Biol Sci. 2002, 269: 867-872. 10.1098/rspb.2002.1959.

    Article  PubMed Central  PubMed  Google Scholar 

  43. Xie RR, Chen XL, Hong XY: Variable fitness and reproductive effects of Wolbachia infection in populations of the two-spotted spider mite Tetranychus urticae Koch in China. Appl Entomol Zool. 2011, 46: 95-102. 10.1007/s13355-010-0014-x.

    Article  Google Scholar 

  44. Duron O, Bouchon D, Boutin S, Bellamy L, Zhou L, Engelstädter J, Hurst GD: The diversity of reproductive parasites among arthropods: Wolbachia do not walk alone. BMC Biol. 2008, 6 (1): 27-10.1186/1741-7007-6-27.

    Article  PubMed Central  PubMed  Google Scholar 

  45. Langmead B, Trapnell C, Pop M, Salzberg SL: Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. 2009, 10: R25-10.1186/gb-2009-10-3-r25.

    Article  PubMed Central  PubMed  Google Scholar 

  46. Wen M, Shen Y, Shi S, Tang T: miREvo: an integrative microRNA evolutionary analysis platform for next-generation sequencing experiments. BMC Bioinformatics. 2012, 13: 140-10.1186/1471-2105-13-140.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  47. Friedländer MR, Mackowiak SD, Li N, Chen W, Rajewsky N: miRDeep2 accurately identifies known and hundreds of novel microRNA genes in seven animal clades. Nucleic Acids Res. 2012, 40: 37-52. 10.1093/nar/gkr688.

    Article  PubMed Central  PubMed  Google Scholar 

  48. Zhou L, Chen J, Li ZZ, Li XX, Hu XD, Huang Y, Zhao X, Liang C, Wang Y, Sun L: Integrated profiling of microRNAs and mRNAs: microRNAs located on Xq27. 3 associate with clear cell renal cell carcinoma. PLoS One. 2010, 5: e15224-10.1371/journal.pone.0015224.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  49. Robinson MD, Oshlack A: A scaling normalization method for differential expression analysis of RNA-seq data. Genome Biol. 2010, 11: R25-10.1186/gb-2010-11-3-r25.

    Article  PubMed Central  PubMed  Google Scholar 

  50. Wang L, Feng Z, Wang X, Wang X, Zhang X: DEGseq: an R package for identifying differentially expressed genes from RNA-seq data. Bioinformatics. 2010, 26: 136-138. 10.1093/bioinformatics/btp612.

    Article  PubMed  Google Scholar 

  51. Storey JD: The positive false discovery rate: A Bayesian interpretation and the q-value. Ann Stat. 2003, 31: 2013-2035. 10.1214/aos/1074290335.

    Article  Google Scholar 

  52. Enright AJ, John B, Gaul U, Tuschl T, Sander C, Marks DS: MicroRNA targets in Drosophila. Genome Biol. 2004, 5: R1-R1.

    Article  PubMed Central  Google Scholar 

  53. Young MD, Wakefield MJ, Smyth GK, Oshlack A: Gene ontology analysis for RNA-seq: accounting for selection bias. Genome Biol. 2010, 11: R14-10.1186/gb-2010-11-2-r14.

    Article  PubMed Central  PubMed  Google Scholar 

  54. Kanehisa M, Araki M, Goto S, Hattori M, Hirakawa M, Itoh M, Katayama T, Kawashima S, Okuda S, Tokimatsu T: KEGG for linking genomes to life and the environment. Nucleic Acids Res. 2008, 36: D480-D484.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  55. Mao X, Cai T, Olyarchuk JG, Wei L: Automated genome annotation and pathway identification using the KEGG Orthology (KO) as a controlled vocabulary. Bioinformatics. 2005, 21: 3787-3793. 10.1093/bioinformatics/bti430.

    Article  CAS  PubMed  Google Scholar 

Download references


We thank Dr. Dong-Xiao Zhao for providing the spider mites. We thank Novogene Company (Beijing) for helping with small RNA Sequencing. We also thank the reviewers for insightful comments and suggestions on the early version of manuscripts. This research was supported in part by a Grant-in-aid from the Science and Technology Program of the National Public Welfare Professional Fund (201103020) from the Ministry of Agriculture of China, and a Grant-in-aid for Scientific Research (31172131, 31371944) from the National Natural Science Foundation of China.

Author information

Authors and Affiliations


Corresponding author

Correspondence to Xiao-Yue Hong.

Additional information

Competing interests

The authors declare that they have no competing interests.

Authors’ contributions

XR conceived and designed the experiments. XR, YZ carried out the data analysis and drafted the manuscript. XH directed and coordinated this study and revised the manuscript. KZ helped with the experiments. All authors read and approved the final manuscript.

Xia Rong, Yan-Kai Zhang contributed equally to this work.

Electronic supplementary material

Additional file 1: Figure S1: Predicted secondary structures of nine novel miRNAs. (TIFF 5 MB)

Additional file 2: Table S2: List of all the miRNAs identified in the four libraries. (XLS 72 KB)


Additional file 3: Table S3: Predicted target genes which may be negatively-regulated by differential miRNAs in the female comparison. (XLS 60 KB)


Additional file 4: Table S4: Predicted target genes which may be negatively-regulated by differential miRNAs in the male comparison. (XLS 29 KB)

Additional file 5: Table S5: Analysis of mRNA-miRNA pairs in the comparison FI vs FU. (XLS 90 KB)

Additional file 6: Table S6: Analysis of mRNA-miRNA pairs in the comparison MI vs MU. (XLS 68 KB)

Additional file 7: Table S1: Primers used in this study. (XLS 26 KB)

Authors’ original submitted files for images

Rights and permissions

Open Access  This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made.

The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

To view a copy of this licence, visit

The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Rong, X., Zhang, YK., Zhang, KJ. et al. Identification of Wolbachia-responsive microRNAs in the two-spotted spider mite, Tetranychus urticae. BMC Genomics 15, 1122 (2014).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: