Transcriptome of Small Regulatory RNAs in the Development of the Zoonotic Parasite Trichinella spiralis

Background Trichinella spiralis is a parasite with unique features. It is a multicellular organism but with an intracellular parasitization and development stage. T. spiralis is the helminthic pathogen that causes zoonotic trichinellosis and afflicts more than 10 million people worldwide, whereas the parasite's biology, especially the developmental regulation is largely unknown. In other organisms, small non-coding RNAs, such as microRNAs (miRNA) and small interfering RNAs (siRNA) execute post-transcriptional regulation by translational repression or mRNA degradation, and a large number of miRNAs have been identified in diverse species. In T. spiralis, the profile of small non-coding RNAs and their function remains poorly understood. Methodology and Principal Findings Here, the transcriptional profiles of miRNA and siRNA in three developmental stages of T. spiralis in the rat host were investigated, and compared by high-throughput cDNA sequencing technique (“RNA-seq”). 5,443,641 unique sequence tags were obtained. Of these, 21 represented conserved miRNAs related to 13 previously identified metazoan miRNA families and 213 were novel miRNAs so far unique to T. spiralis. Some of these miRNAs exhibited stage-specific expression. Expression of miRNAs was confirmed in three stages of the life cycle by qRT-PCR and northern blot analysis. In addition, endogenous siRNAs (endo-siRNAs) were found mainly derived from natural antisense transcripts (NAT) and transposable elements (TE) in the parasite. Conclusions and Significance We provide evidence for the presence of miRNAs and endo-siRNAs in T. spiralis. The miRNAs accounted for the major proportion of the small regulatory RNA population of T. spiralis, while fewer endogenous siRNAs were found. The finding of stage-specific expression patterns of the miRNAs in different developmental stages of T. spiralis suggests that miRNAs may play important roles in parasite development. Our data provide a basis for further understanding of the molecular regulation and functional evolution of miRNAs in parasitic nematodes.


Introduction
Parasites of the genus Trichinella are a group of pathogens with diverse biological and pathological features. Phylogenetic analysis of the mitochondrial DNA of the species identified so far indicated that there are two genetic clades that form unique monophyletic lineages [1]. One clade is represented by T. pseudospiralis, which does not capsulate in the muscle cells, while the other clade is represented by T. spiralis while does encapsulate in muscle cells, named nurse cells, with a parasite inside the cell, surrounded by a sick layer of collagen. Parasites of Trichinella genus are unique intracellular pathogens. Interestingly, their entire life cycle can be completed within an animal, and trichinellosis in human and other mammals was caused through the ingestion of parasite-contaminated meat. This is a typical zoonotic disease that affects more than 10 million people worldwide [2].
After being ingested with the infected muscle tissue, L1 larvae are released and activated in the small intestine, enter the epithelial layer and undergo four times of moultings before maturation into adult worms. Mating is initiated on day 2 after infection and newborn larvae are released by the females into the mucosa as early as 4 days post infection (dpi) [3]. The larvae migrate through the lymphatic and blood vessels, invade striated muscle cells and develop into the infective Ll stage over a period of 2-3 weeks to the next host which complete the life cycle [4]. Thus, unlike other nematodes, T. spiralis does not have an embryonic developmental stage in the egg, which differs markedly in biological and molecular characteristics from other nematodes, especially the well-characterized free-living worm Caenorhabditis elegans. Intriguingly, compared to other nematodes, T. spiralis has a much smaller genome with 64 Mb in nuclear DNA, which contains <15,808 genes [5]. The availability of genome sequence information has made it possible to dissect parasite biology.
Small non-coding RNAs (sncRNAs) are a large group of small endogenous RNAs that have been widely identified in animals [6], plants [7], fungi [8,9] and some viruses [10]. They are generally 21-23 nucleotides in length, which guide various processes involving sequence-specific silencing through chromatin modification, mRNA degradation, and translational repression [11][12][13]. Based on their origins, structures, associated proteins and biological roles, sncRNAs are divided into three general categories: microRNAs (miRNAs), endogenous small interfering RNAs (endo-siRNAs), and piwi-interacting RNAs (piRNAs) [14]. MiRNA and endo-siRNAs have been discovered in diverse animals and plants and fungi [8] [12], while piRNAs are found only in animals [15]. miRNAs are generated from precursor transcripts by two RNase III-type enzymes, Drosha and Dicer. In animals, single-stranded miRNA is incorporated into the argonaute (Ago) protein complexes (Ago) known as RNA induced silencing complexes (RISC) and binds by partially or completely complementary to the 39 untranslated region (39UTR) of a target mRNA [6], this results negative control of gene expression by cleavage or inhibition of translation or other regulatory functions [16]. The biological function of miRNAs was first demonstrated in C. elegans, where two miRNAs (Let-7 and Lin-4) were shown to be regulators for stage-specific differentiation of the worm [17,18]. Endogenous siRNAs are mainly derived from three sources: transposable elements (TEs), complementary annealed transcripts, also called natural antisense transcripts (NAT) and long 'fold-back' transcripts called hairpin RNAs (hpRNAs) [19]. Unlike miRNAs, endo-siRNA's function requires perfect match with the target mRNA [6,20]. MiRNAs and endogenous siRNAs play important roles in the regulation of fundamental cellular processes, including cell differentiation, stress response, apoptosis, proliferation [12] and metabolism [21,22]. Recent studies further differentiated the functions of endo-siRNAs based on the sources they were generated from [14]. The TE-derived siRNAs are likely functional in the germline cells by repression of transposon activity and consequently keep genome stability, while NAT-derived siRNAs are more functional in somatic cells [14].
In this study, the transcriptional profiles of both miRNA and endo-siRNA in three developmental stages, namely adult (Ad), new born larvae (NBL) and muscle larvae (ML) of T. spiralis were systematically investigated and compared by high throughput cDNA sequencing technology (''RNA-seq''). We found that miRNAs were mainly expressed in the adult worm stage and endo-siRNAs were predominantly derived from transposable elements in the genome.

Summary of small RNA sequencing
To identify miRNA and endo-siRNA involved in the development of T. spiralis, small RNA libraries were generated from three life cycle stages, the new-born larvae, musclar larvae and adult worms. The libraries were directly sequenced using Solexa sequencing technology. A total of about 40 million unfiltered sequence reads (12,368,833 reads from adult worms, 12,867,246 reads from newborn larvae and 14,495,649 reads from muscular larvae) with sizes between 18 and 30 nucleotides were obtained, respectively. After removal of low quality reads and adaptor sequences, the clean reads obtained in Ad, NBL and ML were 11,878,917 (96.0% of total reads), 12,357,960 (96.0% of total reads) and 14,078,375 (97.1% of total reads) respectively. Sequence characterization suggested that the small RNA pools contained miRNAs (20.63%), other non-coding RNAs (rRNA, tRNA, snoRNA) (0.55%), mRNA-related small RNAs (32.13%), TE-related small RNAs (0.02%) and unknown small RNA transcripts (46.66%) (Fig. 1A, Table S1). The GC content of the small RNAs was around 38.7% and the distribution of the small RNA populations in the libraries generally followed the patterns identified in other organisms, such as schistosomal parasites [23][24][25][26][27][28][29]. The libraries generated in this study have significantly high coverage than those in a recent report, which only showed a proportion of miRNAs in one single stage of the parasite [30].
As shown in Fig. 1B, the unique clean reads obtained in Ad, NBL and ML were 1,816,014, 1,691,015 and 1,936,612, respectively (Table S2). The number of small RNAs existed in all three life cycle stages was 393,026 reads, which accounted for 21.6%, 20.3% and 23.2% of the unique clean reads of Ad, NBL and ML stage, respectively (Fig. 1B). The small RNAs commonly expressed in Ad stage and ML stage were 3 times more than that between Ad stage and NBL stage (Fig. 1B), indicating that many genes of Ad stage were activated already in ML and NBL had more stage-specific small RNAs.

Length variation of small RNAs in T. spiralis
The length of predicted small RNAs varied from 18 to 30 nt in the three life stages, with most of them between from 18 to 27 nt. Here the size and number of sequences in Ad and ML almost followed the same pattern, whereas the number of sequences between 18 and 22 nt was enriched in NBL ( Fig. 2 and Table S3). This implied that, apart from stage-specific activation in NBL, small RNAs may undergo size differentiation during early development of the parasite. A previous study indicated that the length of the non-conserved connector helix in Dicer was the main determinant of product size of the mature small RNAs [31], providing a possible explanation for the stage-specific variation we observed. The mechanism behind the stage-specific small RNA processing remains unknown. Since both piRNA and endo-siRNA were more common in the germline, it is likely that the small RNAs found in NBL were predominantly miRNAs, which may play fundamental roles in stage-specific differentiation. Longer RNAs may be more important for parasite reproduction [32][33][34][35][36][37].

Identification of miRNAs in T. spriralis
To identify candidate miRNAs of T. spiralis, the clean unique reads obtained were mapped to the draft T. spiralis genome sequences (http://genome.wustl.edu/pub/organism/Invertebrates/Trichinella_ spiralis/assembly/Trichinella_spiralis-1.0/) [5] sing SOAP. A total of 72,110 (out of 2,099,966) unique small RNA reads that perfectly matched are referred to as miRNA candidates. After screening for secondary structure of the inverted repeats (found with Einverted of Emboss) with RNAfold and evaluation by MirCheck, a total of 240 T. spiralis predicted miRNAs were identified.
To further characterize miRNAs in T. spiralis, all predicted miRNAs obtained above were compared against a miRNA database, miRBase (Release 15.0). In total, 21 miRNAs were found that had been identified in other species, and belong to 13 different miRNA families ( Table 1, Table S4 and S5). Most miRNAs contain a 7 nt region (typically positioned at the 59 side 2-8 nt) known as the ''miRNA seed'' sequence ( Fig. 3). It has been suggested that the seed regions serv to anchor miRNAs to their mRNA targets [14].
The sequencing data showed that, of the 21 conserved miRNAs, tsp-let-7 and tsp-miR-87 were found to locate only in the 39 arm of their pre-miRNAs, and tsp-miR-31 was located only in the 59 arm of the hairpin structures. The remaining 18 miRNAs were derived from both arms of the pre-miRNAs (Table S4 and S5). According to the current model of miRNA maturation, Dicer recognizes the double-strand of pre-miRNAs, cleaves away the loop structure and de-associates the duplex to generate antisense and sense strands, i.e. the mature miRNA and miRNA* [38]. Although the antisense and sense strands may have different thermostability, which also determined there relative abundance in different tissues, both strands could be functional in miRNA-mediated regulatory pathways [39]. For instance, tsp-miR-1-3p was detected at 1,515 TPM (transcripts per million) in the Ad stage, 9,759 in the NBL stage and 3,351 in the ML stage, respectively, whereas its counterpart tsp-miR-1-5p was much less abundant ( Fig. 4, Table  S5). The TPM of tsp-miR-100-5p and tsp-miR-100-3p, tsp-miR-125-5p and tsp-miR-125-3p, tsp-miR-9-1-5p and tsp-miR-9-1-3p, tsp-miR-9-2-3p and tsp-miR-9-2-5p behaved in similar fashion (Fig. 4, Fig. 5A and Table S5). Thus the anti-sense strands complementary to mRNA targets may play a main regulatory role. However, it cannot be ruled out that both strands generated from the same pre-miRNA execute similar functions. This is because, in some miRNAs, both strands were either similarly expressed or showed stage-specific expression patterns ( Fig. 5 A and B; Table S5 and S6).
We found 213 novel miRNAs ( Table 2, Table S7). Comparison of novel miRNAs expressed in the three developmental stages of T. spiralis revealed stage-specific expression patterns (Fig. 5). A majority of the novel miRNAs in NBL stage has a lower expression level than those in Ad stage and ML stage (Fig. 5B, Table S6 and S7). For instance, tsp-novel-21, tsp-novel-50a-3p, tsp-novel-50b-3p were relatively enriched in adult worm, tsp-novel-46 were abundant in muscular larvae, whereas tsp-novel-108 and tspnovel-83 exhibited a high abundance in all the three stages (Table  S7). Of the novel miRNAs, the expression patterns of tsp-novel-93-5p and tsp-novel-93-3p, tsp-novel-50a-3p and tsp-novel-50b-3p, tspnovel-101-3p and tsp-novel-101a-5p and tsp-novel-101b-5p were interesting. Tsp-novel-93-5p and tsp-novel-93-3p were derived from the same pre-miRNA, but tsp-novel-93-5p was mainly expressed in the NBL stage, less in ML and least in Ad. While tsp-novel-93-3p was much less expressed than tsp-novel-93-5p in all three stages. (Fig. 5B, Table S7). Tsp-novel-50a-3p and tsp-novel-50b-3p were derived from the 39 arm of two gene copies with a similar sequence. Tsp-novel-50a-3p was dominantly expressed in Ad, while tsp-novel-50b-3p was mainly expressed in NBL and Ad but less in ML. The genes coding for tsp-novel-101-3p and tspnovel-101a-5p and tsp-novel-101b-5p and their expression patterns were even more complicated. They were encoded by three genes in the genome of T. spiralis. The sequences of tsp-novel-101-3p derived from the 3 genes were the same, while tsp-novel-101b-5p was derived from 1 gene which had a single nucleotide change from tsp-novel-101a-5p ( Fig. 5B and Table S6 and S7). The later was derived from two gene copies. All these data suggested that there may be multiple layers of control for the stage-specific expression of miRNA genes, as well as their potential regulatory function in the development of the parasite [40]. Although the mRNA targets of these regulatory miRNAs have not been identified, with the availability of the genome sequence of the parasite, genome-wide association study can be pursued. For instance, miR-1 has been found in many species, from Drosophila to human [41][42][43], suggesting that they are evolutionary conserved, which has been characterized to play essential functions in regulating proliferation and differentiation of muscle cell [41,44]. One homologous to mir-1, tsp-miR-1, was observed to express in  three developmental stages of T. spiralis and may have similar functions. Further, MiR-100 and let-7, the two conserved miRNA in metazoa, play a role in regulation of developmental timing [18,42,45]. Their homologs, tsp-miRNA-100 and tsp-let-7, were found throughout the life cycle of T. spiralis. Tsp-let-7 showed very low expression in NBL stage, whereas tsp-miRNA-100 was detected in rather high abundance at the same development stage. The abundance of tsp-miRNA-100 was almost identical with that of tsp-let-7 in both Ad and ML stage, indicating that miRNA-100 may be more functional in NBL stage.
Characterization of stage-associated miRNAs in T. spiralis Quantitative real-time PCR and northern blots were performed to validate miRNAs identified in T. spiralis and their relative expression levels at different developmental stages. Five conserved miRNAs (tsp-miR-228, tsp-miR-100, tsp-let-7, tsp-miR-1 and tsp-miR-31) and five novel miRNAs (tsp-novel-108, tsp-novel-83, tspnovel-46, tsp-novel-86 and tsp-novel-21) with relatively higher TPM values identified by sequencing were validated by qRT-PCR and Northern blot. First, the above miRNAs identified in sequencing were all amplified by qRT-PCR (Fig. 6A) and except the tsp-let-7 which was found less expressed in NBL, the qRT-PCR results were all consistent with the TPM values of sequencing results. The reason for tsp-let-7 being found less common in NBL is not known. This may be due to preferential amplification in the tag-generation step before RNA-seq or other unidentified factors. Further, all qPCR amplicons were cloned into the pMD-18T vector (Takara, Dalian, China) and sequenced. The sequencing results demonstrated that the amplicons were identical to the miRNAs sequences (data not shown). Four miRNAs were selected for characterization by northern blots. The results confirmed the qPCR data. However, the probe for tsp-mir-100 did not show any signal with total RNA purified from the Ad stage; this was likely due to probe modification or other unidentified factors that prevented annealing of the probes with the targets (Fig. 6B). Since this miRNA was found by both  sequencing and qPCR in the Ad stage, we have no doubt that tspmir-100 was indeed expressed at this stage.
TEs are major components of the intergenic regions in the genomes of eukaryotic organisms [47,48]. Based on their structures and modes of integration, TEs are comprised of two main classes [51]. One class includes retrotransposons and retrotransposon-like elements such as Long Interspersed Nuclear Elements (LINE), Long Terminal Repeat Elements (LTR), and Short Interspersed Elements (SINE). The other group includes DNA transposons. We predicted TE structures in the T. spiralis genome with RepeatMasker (http:// www.repeatmasker.org/), and endo-siRNAs matched to TEs were further analyzed. As reported recently [5], the repetitive sequences in the T. spiralis genome estimated to be around 18%, which is much lower than in schistosomes and Drosophila melanogaster [52,53]. Thus, TE-derived siRNAs accounted for a minor portion of the sncRNAs identified here (Fig. 1A).
The numbers for TE-derived siRNAs identified in the Ad, NBL and ML stages were 2,055, 2,234 and 2908 respectively. No obvious stage-related variations in endo-siRNA expression were found, suggesting that TE components in the genome might not be very active during embryogenesis or during other development stages. Further, most of these siRNAs were derived from transposons of Long Interspersed Nucleotide Elements (LINE) family and the transposon Charlie 24 accounted for a major portion (Table S8). Sequencing analysis showed that more siRNA are derived from the antisense strand than from the sense strands (data not shown). This suggested that antisense-derived siRNA may execute regulatory functions through hybridization with the mRNA template generated from the sense strand.
Since NAT-derived siRNAs were generated from doublestranded RNAs formed by complementary annealed transcripts (NAT) and long 'fold-back' transcripts (called hairpin RNAs, hpRNAs) [19], T. spiralis NAT-siRNAs were identified from predicted overlapping genes. The number of NAT-derived siRNAs identified in Ad, NBL and ML was 21,157, 16,243 and 21,251, respectively (Table S9). Thus the number of NAT-derived endo-siRNAs was much higher than that from TE, which suggests that these siRNAs played more regulatory roles in the development and parasitization of T. spiralis. Further, all NAT-derived siRNAs were trans-NAT siRNAs and no cis-NAT siRNAs were found ( Fig. 7 and data not shown), which differs from observations with D. melanogaster and mice [33,[47][48][49]. This also suggests that there are fewer coding sequences with internal inverted repeat sequences in the genome of T. spiralis.
In summary, we identified and analyzed the expression of miRNAs and endo-siRNAs in three development stages of T. spiralis through high through-put RNA sequencing techniques. We found vastly more transcripts of miRNAs than that of endo-siRNAs. A total of 21 conserved miRNAs in 13 metazoan miRNA familis and 213 novel miRNAs were identified in the parasite. Some showed clear stage-specific expression patterns, suggesting a potential regulatory function in the corresponding developmental stages. Endo-siRNAs were mainly derived from natural antisense transcripts with TEderived siRNAs accounting for only minor proportion of the small RNA population. Thus, the function of endo-siRNAs in T. spiralis is likely to regulate gene expression instead of maintaining genome stability. The data of this study provide insight information for further dissection of the parasite's biology.

Parasites
Muscle larvae (ML) of T. spiralis (strain ISS534) were obtained from rats 35 days post infection (dpi) by digestion of minced skeletal muscle in 1% pepsin, 1% HCl for 3 h at 37uC with agitation as previously described [54]. To purify adult worms and newborn larvae, Wistar rats at 6 weeks of age were orally inoculated with T. spiralis (strain ISS534) with a dose of 8000 larvae per rat. At 30 h or 6 dpi, all rats were killed and the entire intestines were removed, opened longitudinally and cut into small pieces (about 0.5-1 cm). The fragmented intestine was put on a layer of gauze which was immersed into 0.9% sodium chloride solution at 37uC and incubated for 3 h. Adult T. spiralis worms migrate into the liquid phase, which was harvested by centrifugation. To obtain new-born larvae, adult worms collected at 6 dpi were incubated in Iscove's Modified Dulbecco's Medium (IMDM) in 75-cm 2 cell culture plates at 37uC, and the newborn larvae were harvested every 12 h. Our study was reviewed and approved by the Ethics Committee of Jilin University (Ethical clearance application number IZ-2009-III). All animal work was conducted according to Chinese and international guidelines.

RNA isolation
Total RNA of T. spiralis (Ad, NBL and ML) was extracted using Trizol reagent (Invitrogen, CA, USA) according to the manufacturer's instructions. RNAs were dissolved in diethylpyrocarbonate (DEPC)-treated water, aliquoted and stored at 280uC. RNA was quantified by measuring the absorbance at 260 nm with a Nanodrop 1000 machine (Thermo Scientific CA, USA).

Construction of small RNA libraries and sequencing
For small RNA library construction and deep sequencing, the 15-30 nt size range base-pair fraction of each RNA sample from the three life cycle stages was first enriched by 15% TBE urea polyacrylamide gel electrophoresis and the Illumina's proprietary adaptors (UCAGAGUUCUACAGUCCGACGAUC and UC-GUAUGCCGUCUUCUGCUUGUidT) were ligated to the 59 and 39 termini of the purified small RNAs which were converted into single-stranded cDNA with Superscript II reverse transcriptase (Invitrogen, CA, USA) and Illumina's small RNA RT-Primer. The cDNAs were PCR-amplified with high fidelity Phusion DNApolymerase (Finnzymes Oy, Finland) in 18 PCR cycles using Illumina's small RNA primer sets. After purification using 6% TBE urea PAGE gels, the PCR products were sequenced by Solexa's sequencing-by-synthesis method (Fig. S1).

Mapping the sequence reads onto the reference genome
The bioinformatic analysis and work flow is shown in Figure S1. After removing the low quality sequence reads and the adapter   tags of more than 5 reads that matched perfectly or near-perfectly (no more than 1 mismatch and the mismatch not positioned in the seed region) were regarded to be conserved miRNA candidates. For reads that did not match to the miRNA database, we used the software Einverted of Emboss [57] to find the inverted repeats (stem loops or hairpin structure). Each inverted repeat was extended 10 nt on each side, and the secondary structure of the inverted repeat was predicted by RNAfold [58]. Unique reads with a folding free energy of at least 25 kcal/mole( gGu folding #225 kcal/mol) were evaluated by MirCheck [59] with modified parameters. Finally, precursors (hairpins) of miRNA that passed MirCheck were inspected manually in order to remove false predictions. The reads passing the inspections were regarded as novel miRNAs.
Similar to credibility interval approaches reported for the analysis of SAGE data [60], we employed IDEG6 [61] to identify miRNAs showing statistically significant differences in relative abundance (as reflected by the total count of individual sequence reads) between the three libraries (corresponding to the three developmental stages of the parasite). We used the general Chisquare method for comparison analysis, which has been commonly applied by others [27,36]. Finally, miRNA with a P value#0.01 were deemed to be significantly different between the samples of the three developmental stages of the parasite.
The repeated sequences, i.e. potential transposable elements, in the T. spiralis genome were predicted by using RepeatMasker (http://www.repeatmasker.org/) and the sequences of typical transposons were annotated (data not shown). The sequencing reads that perfectly matched the T. spiralis genome were aligned to repeats (TE) using SOAP. The reads that perfectly matched TEs were considered TE-derived siRNAs.
Natural antisense transcripts (NATs) were detected by aligning the T. spiralis predicted genes to each other. If a pair of overlapping genes were matched on opposite strands with an E value of #1e-9 [62], this pair of overlapping genes was defined as a NAT pair. The reads that perfectly matched the T. spiralis genome were aligned to overlapped sequences of NAT pairs with SOAP. The reads that perfectly matched the overlapped regions were considered NAT-derived siRNAs.

miRNAs quantification by real-time PCR
Total RNA purified from the three stages were polyadenylated with E. coli poly(A) polymerase (E-PAP) following the manufacturer's protocol of the poly(A)-tailing kit (Ambion, CA, USA). The RNA samples were purified separately from the parasites of the same developmental stages for sequencing purpose. After phenol-chloroform extraction and ethanol precipitation, the polyadenylated products were dissolved in DEPC-treated water and reversetranscribed with 200 U of SuperScript TM III Reverse Transcriptase (Invitrogen, CA) and 39 RACE Adapter (59-GCGAGCACAGAAT-TAATACGACTCACTATAGGT12VN-39) in the FirstChoice RLM-RACE Kit (Ambion) according to the manufacturer's protocol. The 20 ml RT reaction contained 1 mg total RNA, 2 ml 39 RACE Adapter, 0.5 mM dNTP mix (Takara, Dalian, China), 10 U RNase inhibitor and 200 U SuperScript TM III Reverse Transcriptase.
In quantitative RT-PCR reactions, GAPDH (glyceraldehyde-3phosphate dehydrogenase) was chosen as an endogenous reference. The forward and reverse primers for GAPDH were 59-GTGCTGATTACGCTGTTG-39 and 59-CTAAGCCATTGG-TAGTGC-39. PCR was done on an Applied Biosystems 7500 system. The following forward primers were designed to confirm the sequencing results of miRNAs that showed differential expression patterns: tsp-miR-  [63]. ANOVA and Tukey's HSD posthoc test were used to analyze significant differences among three stage; P,0.05 was considered significant.

(TIF)
Table S1 General information of the small RNA libraries. (DOC)