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

Background 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. Results 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. Conclusion 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. Electronic supplementary material The online version of this article (doi:10.1186/1471-2164-15-1122) contains supplementary material, which is available to authorized users.


Background
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 Wolbachiamediated phenotypes in hosts.
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,[18][19][20]. 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 miR-NAs 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).
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 (http://www.mirbase.org/ftp.shtml) 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 miR-NAs 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.
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).

Expression analysis of Wolbachia-responsive miRNAs
In females, 32 miRNAs were observed to be upregulated and 59 down-regulated in the Wolbachia-infected line compared to the Wolbachia-uninfected line (Additional file 3: Table S3, Figure 4), while 11 miR-NAs 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. 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 miR-NAs 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). 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 femalespecific differential miRNAs were significantly downregulated 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.

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.
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 upregulation 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 tetur 03g00110), novel_50 (targets tetur26g00300), may also regulate apoptosis (Additional file 3: Table S3). The    [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 resistanceassociated 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 miR-NAs 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 genderspecific 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.

Conclusion
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.

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 (http://bioinformatics.psb. ugent.be/orcae/overview/Tetur; 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 (https://bioinformatics.psb.ugent. be/gdb/tetranychus/small_RNAs/). Tags originating from repeat sequences were filtered by using a repeat sequences database (http://www.repeatmasker.org/cgi-bin/ WEBRepeatMasker/), and tags originating from proteincoding 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 (foldchange)| > 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 (https://bioinformatics.psb. ugent.be/gdb/tetranychus/) 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 www.ebi. ac.uk/arrayexpress 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 downexpressed 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-dayold 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 Prime-Script 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.