A Deep Analysis of the Small Non-Coding RNA Population in Schistosoma japonicum Eggs

Background Schistosoma japonicum is a parasitic flatworm that causes zoonotic schistosomiasis. The typical outcome of schistosomiasis is hepatic granuloma and fibrosis, which is primarily induced by soluble egg-derived antigens. Although schistosomal eggs represent an important pathogenic stage to the host, the biology of this critical stage is largely unknown. We previously investigated the expression profiles of sncRNAs during different developmental stages of this parasite. However, using small RNA extracted from egg-deposited liver tissues generated limited information about sncRNAs in eggs. Here, we characterized the complete small RNAome in this stage of the parasite after optimization of RNA purification. Methodology and Principal Findings A library, SjE, was constructed with the small RNA extracted from S. japonicum eggs and subjected to high-throughput sequencing. The data were depicted by comprehensive bioinformatic analysis to explore the expression features of sncRNAs in the egg stage. MicroRNAs accounted for about one quarter of the total small RNA population in this stage, with a strongly biased expression pattern of certain miRNA family members. Sja-miR-71, sja-miR-71-5p, and sja-miR-36-3p were suggested to play important roles in embryo development. A panel of transfer RNA fragments (tRFs) precisely processed from the 5′ end of mature tRNAs was identified for the first time, which represented a strong egg stage-biased expression. The tRNA-Ala derived small RNAs were the most highly expressed Sj-tRFs in eggs. Further, the expression of siRNAs from 29 types of well-defined transposable elements (TEs) was observed to be relatively stable among different developmental stages. Conclusions and Significance In this study, we characterized the sncRNA profile in the egg stage of S. japonicum. Featured expression of sncRNAs, especially the tRNA-derived small RNAs, was identified, which was further compared with that of other developmental stages. These novel findings would facilitate a deeper understanding of the biology of schistosomal parasites.


Introduction
Schistosomiasis, a debilitating disease, caused by blood flukes of the genus Schistosoma afflicts more than 230 million people worldwide (http://www.who.int/mediacentre/factsheets/fs115/ en/index.html). The three major species infecting humans are Schistosoma haematobium, S. mansoni, and S. japonicum. The pathology of chronic infection with S. japonicum or S. mansoni is well known as hepatosplenic schistosomiasis, with clinical symptoms of granulomatous inflammation, periportal fibrosis, portal hypertension, hepatosplenomegaly, ascites, and the formation of vascular shunts [1,2]. The granulomatous responses induced by schistosome soluble egg antigens (SEA) released from the eggshell-enclosed miracidium are regarded as an evolutionary compromise, that is critical for the survival of the infected host, and but also beneficial for the transmission of eggs [3]. As a classical immune regulatory model, the host immune responses induced by SEA were intensively investigated [4][5][6][7]; however, the gene expression regulatory mechanism during schistosomal embryonic development is still poorly understood.
Small non-coding RNAs (sncRNAs) with a size of 18,30 nt have been found in most eukaryotes, and are increasingly recognized as powerful regulators of gene expression and genome stability [8,9]. Among them, microRNAs (miRNAs), small interfering RNAs (siRNAs), and piwi-interacting RNAs (piRNAs) are the three major categories. So far, numerous miRNAs have been extensively identified in animals [10], plants [11], fungi [12], and viruses [13]. In mammals, mutation or deletion of enzymes involved in miRNA biogenesis has been observed to lead the defects in germ-line division and differentiation, and embryonic morphogenesis [14,15]. In the nematode of Meloidogyne incognita, knockdown of drosha and pasha in undifferentiated eggs led to irregular growth and embryonic lethality [16]. Recent advances have also proved that non-coding RNAs (ncRNAs) are key regulators of embryogenesis, including miRNA-induced degradation of mRNAs and long ncRNA-mediated modification of chromatin [17]. Two other classes of sncRNAs, siRNAs and piRNAs, have known involvement in the defense against parasitic DNA elements to maintain genome stability. PiRNAs have been proposed to ensure germline stability in germ-line cells, whereas siRNAs were observed to play roles in both somatic and germ-line cells. These sncRNAs are loaded into the RNA-induced silencing complex (RISC) [18] or RNA-induced transcriptional silencing complex (RITS) [19] to function in chromatin architecture modelling, post-transcriptional repression and mRNA destabilization, mobile genetic elements suppression, and virus defence [8,[20][21][22]. In schistosome, several Argonaute orthologues were identified in both S. japonicum [23,24] and S. mansoni [25], and one of the three Argonaute proteins in S. japonicum, SjAgo2 was suggested to maintain genome stability via suppression of retroelements [26].
Since the tissue-trapped eggs are the major agents causing the severe pathology of schistosomiasis and those released from the host are relevant for the prevalence of the disease, it is indispensable to explore a complete repertoire of sncRNAs in schistosomal eggs, which will assist the discovery of novel intervention targets. In this study, small RNA extracted from purified S. japonicum eggs was subjected to high-throughput sequencing and deep analysis. The data provide a unique expression feature of egg sncRNAs, at a comparable level to those from other developmental stages of S. japonicum, which will shed light on the gene regulatory mechanisms during embryonic morphogenesis of the schistosomal parasite.

Results and Discussion
Isolation and purification of S. japonicum eggs Rapid isolation of viable S. japonicum eggs from host hepatic tissue was a critical step for extracting intact total RNA. In our previous study [38], we analyzed the sncRNAs by directly sequencing the total small RNA from infected liver tissues. However, the host small RNA population, which overwhelmed that originating from the eggs, significantly reduced the resolution of the egg-derived small RNAome. In this study, an improved sieving and enzymatic methodology was applied. The purified egg samples were examined under a light microscope, and we found that most of the eggs contained a developing embryonic larva (miracidium) ( Figure S1), whereas small-sized ones with immature embryos had either passed through the nylon mesh screens or were removed from the suspension in the Percoll column after centrifugation. Therefore, the data obtained here reflect mainly the sncRNA of mature S. japonicum eggs, which were able to release SEA to elicit the host hepatic granulomatous responses.

General features of the two small RNA libraries of S. japonicum eggs
To investigate the small RNA profiles in the S. japonicum eggs, three libraries were constructed with small RNAs extracted from purified eggs and sequenced separately. Preliminary analysis indicated that one tRNA-derived small RNA fragment was preferentially amplified (<32% of total reads) in one library, but not in the other two libraries. There were no significant differences in term of read numbers and sequence length distribution between the second and the third libraries (data not shown). Thus only the data of the second library were further analyzed. In total, 34,244,779 reads were generated by Solexa sequencing of the egg small RNA library (Table S1). In the library, <30.4 million reads were high-quality clean reads, which could be merged into 3,053,121 unique tags (Table S1 and S2). The redundancy of the library was 89.9 (Redundancy = 1002[Total Unique Clean Reads/Total high-quality Clean Reads6100]), which was higher than that of our previous small RNA libraries [36,38]. The match rate of the library was more than 70% (Table S2), which is dramatically higher than that of our previous libraries SjE30 and SjE45 (both at <1%), which were constructed with total small RNA isolated from egg-trapped liver tissues [38]. Therefore, after eliminating contaminated reads derived from host tissue, the data presented here should reflect the authentic repertoire of sncRNA in S. japonicum eggs. Since the egg small RNA library has the same order of magnitude of reads number as other available S. japonicum small RNA libraries, it is possible to compare the profile and expression level of sncRNAs in the egg stage with that of other developmental stages.
Length distribution of small RNA reads in different S. japonicum libraries Here, we investigated the length distribution of small RNA reads from the egg small RNA library along with libraries of other developmental stages that perfectly matched the draft genome sequences of S. japonicum. As shown in Figure 1, the dominant species of small RNAs within S. japonicum are between 18 and 23 nt, which ruled out the piRNA (normally 29-31 nt) pathway in this species. The length distribution of the reads in SjE presented a pattern similar to that of SjC and SjH, at both total and unique levels, which featured as that the percentage of 22-nt and 23-nt reads, particularly the latter, was significantly higher at total level compare to that of unique level ( Figure 1A, B, and D).

Classification of sncRNAs in different small RNA libraries of S. japonicum
The sncRNA transcripts in different small RNA libraries of S. japonicum were systematically analyzed using a more rigorous bioinformatic pipeline than previously described [38]. As shown in Figure 2, the percentage of miRNAs in the egg libraries accounted for <24%, which was slightly higher than that in the SjF library (9.65%), but lower than that in other three libraries (SjC, SjL, and SjH). The percentage of reads mapping to the 28S, 18S, and 5.8S rRNA genes, as well as the intergenic spacer sequences among these genes in the egg libraries was also modest (<25%) when compared with that in the SjL, SjH, SjM, and SjF libraries. Unexpectedly, the reads derived from tRNAs were dramatically expanded in the egg libraries, compared with that of the other libraries. The tRNA-derived small RNA reads accounted for 23.0% of the RNA population in SjE, suggesting that there may be a specific processing mechanism of tRNA transcripts in the egg stage of the parasite. In addition, the TE-derived siRNAs in the egg libraries primarily originated from two types of retroelements, long terminal retrotransposons (LTRs) and long interspersed nucleotide elements (LINEs). There was no significant variation in the transcription of TE-derived siRNAs among different developmental stages of the parasite.

MiRNAs expressed in S. japonicum eggs
The clean reads from the six libraries were aligned to the 55 S. japonicum miRNA precursors in the Sanger miRBase [43,44] (Release 18). In the egg libraries, 75 out of 78 known S. japonicum mature miRNAs were detected, significantly more than that of previous egg libraries (18 and 25 known mature miRNAs were found in the SjE30 and SjE45 libraries, respectively), suggesting that those miRNAs with low expression in eggs were detected in this study (Table S3). Of the miRNAs identified, sja-miR-71b-5p, sja-miR-71, sja-miR-1, sja-miR-36-3p, and sja-124-3p were the most abundant members at the egg stage ( Figure 3A). These five miRNAs accounted for approximately 86% of all known miRNAs in the SjE library, which further supports our earlier finding that there is a strongly biased expression of particular miRNA families in each particular developmental stage of the parasite [38] ( Figure 3). A similar phenomenon was also observed in other species, such as Clonorchis sinesis, in which members from the miR-71 family accounted for one third of the reads in the adult stage [37]. In the parasitic nematode Trichinella spiralis, members derived from the miR-1 and let-7 families were predominantly expressed in larvae [45]. Combining the TPM value of miRNAs (Table S3) and Northern blot analysis (Figure 4), we found that members of the sja-miR-71 family were the most highly expressed ones in the egg stage, implying that these miRNAs may play important regulatory functions during this stage.
The miR-36 family has so far been observed only in helminthes [46]. A conserved ortholog of sja-miR-36-3p was identified in S. mansoni by computational prediction [29]. Recently, Liu et al. also detected a putative miR-36 family member (tsp-Novel-21) in T. spiralis, which was mainly expressed in the adult worms [45]. The alignment of the Sj-miR-36-3 sequence with orthologs from other organisms is shown in Figure S2. All the orthologs shared a conserved seed sequence ''CACCGGG'' except bma-miR-36a and bma-miR-36b. In C. elegans, miR-36 was one of the eight functionally redundant members of the cel-miR-35 family (cel-miR-35,42). Previously, Alvarez-Saavedra and colleagues comprehensively analyzed the function of the miR-35 family members in C. elegans, and found that mutation of seven members of this family led to developmental suppression, including embryonic and larval lethality [47]. More recently, this family was found to regulate the G1/S transition of intestinal cells and promote proliferation of germ cells in C. elegans [48]. Further, in the parasitic nematode of Ascaris suum, each asu-miR-36 family member was expressed in a narrow window from early to middle embryogenesis, implying that each member from this family may finely regulate the development of the parasite [41]. Considering the conserved seed sequence and its high expression in eggs compared to adult worms, it is postulated that sja-miR-36-3p may plays a similar role during embryonic development of S. japonicum. However, sja-miR-36-3p was also abundantly expressed in cercariae and lung-stage schistosomula (Table S3), suggesting that subtle regulatory mechanisms may be exerted at different developmental stages by one schistosomal miRNA.
tRNA-derived small RNA in S. japonicum Previously, tRNA-derived RNA fragments (tRFs) precisely processed from mature or precursor tRNAs were detected in prostate cancer cell lines by ultra-high-throughput sequencing [49]. More recently, tRFs were identified in the halophilic   archaeon Haloferax volcanii [50] and in the plant pathogenic fungus, Magnaporthe oryzae [51], indicating the existence of tRFs in various organisms. However, because the tRFs were found to be processed from different mature tRNAs by specific endonucleases (e.g., ELAC2 and Dicer) under stress responses and possibly on other occasions [52], tRFs may be suppressed at other developmental stages of S. japonicum, which may have prevented the identification of this group of small RNAs in our previous studies [36,38]. Here, after mapping to the predicted S. japonicum tRNA sequences (sja.trna.bed), we unexpectedly found that small RNAs derived from tRNAs were abundantly present in the SjE library. The percentage of tRNA-derived small RNAs was significantly upregulated in the egg stage compared with that of other stages ( Figure 2). For the first time, we defined a panel of highly expressed Sj-tRFs processed from mature tRNAs of S. japonicum ( Table 1). The secondary structures of the tRNAs for generating these Sj-tRFs are shown in Figure 5A-K.
Previous studies indicated that tRFs can be derived from different positions within tRNAs [49]. Lee et al. categorized tRFs into three series, of which tRF-5 and tRF-3 were aligned to the 59 and 39 ends of mature tRNA, respectively, whereas the tRF-1 series are located within pre-tRNA, and their 59 ends start precisely after the 39 ends of the mature tRNA sequence [49]. Here, we found that the tRFs in S. japonicum were preferentially processed from the 59 end of mature tRNAs (Table 1 and Figure 5A-K), which represent an extensive terminal and asymmetric processing of tRNA, as recently reported [53]. Furthermore, like that observed in miRNAs, isoforms from particular tRNA(s), such as Sj-tRF-001-1,Sj-tRF-001-4, were also commonly identified ( Table 1).
The tRF, Sj-tRF-001-1, which could be processed from mature tRNA-Ala (AGC) and tRNA-Ala (UGC) ( Figure 5A), represented the most abundant read deposited in the SjE library. However, the hybridization signal corresponding to 24-nt small RNAs was stronger than that of 19-nt ( Figure 5J, left panel), which could be due to the cross-hybridization to some other tRFs homologous to Sj-tRF-001-2. In addition, Sj-tRF-002 and Sj-tRF-006, which are respectively derived from tRNA-Ala (CGC) and tRNA-Cys (GCA), are highly homologous to Sj-tRF-001-1, with only one base mismatch. Therefore, the signal at the position of 19-nt may  reflect the expression of all three tRFs ( Figure 5J, left panel). Sj-tRF-003-1 and Sj-tRF-003-2 exhibited a biased expression in eggs as well as their pre-and mature tRNAs ( Figure 5J, middle panel). Moreover, a relatively low expressed tRF, Sj-tRF-010, was also detected mainly in the egg stage ( Figure 5J, right panel). Although the biogenesis and function of tRFs remains to be further clarified, one exciting finding reported by Lee et al, was that tRF-1001 can regulate cell proliferation in prostate cancer cell lines [49]. As tRFs are predominantly expressed in the egg stage of S. japonicum, we postulate that some of these tRFs may play a role in regulating the development of schistosomal embryos.
Small RNAs mapping to S. japonicum rDNA repeats Deep analysis of the small RNA identified a large proportion of sncRNAs that were derived from rDNA, especially in the hepatic schistosomula and adult worms (Figure 2). A similar result was also observed by Xue et al., who observed an over-representation of fragmented rRNAs in the pool of short RNAs in their S. japonicum library [32]. They suggested that processing of the 28S large rRNA subunit (the phenomenon known as ''nick in vivo'' in flatworms) may be the source of the large proportion of rRNA fragments. After precisely mapping the small RNAs to the ribosomal DNA repeat sequences (5.8S-ITS2-28S-IGS-18S-ITS1) of S. japonicum, we confirm that the majority of these small RNAs were derived from the 28S rRNA ( Figure 6). Notably, the amount of small RNA reads, 59-CUGACCUCGGAUCAGAC-GUGAU(U)-39, derived from the 59-terminus of 28S rRNA was predominant in the SjE library. Since the highly integrity of the total RNA extracted from the eggs ( Figure S3), it is suggested that the high-abundant rRNA fragments derived from specific regions of the 28S rRNA were not due to the random degradation of rRNA transcripts. Whether these small RNAs are functionally processed products of rRNA, as in the case of tRNA-derived sncRNAs, remains to be determined.

Expression of 29 well-defined TE-derived endo-siRNAs in S. japonicum
TE sequences constitute 21.84% of the S. japonicum genome [54], and the mobile genetic elements (MGEs) in the S. japonicum genome, including short interspersed nucleotide elements (SINEs)like retrotransposons, LTR retrotransposons, non-LTR retrotransposons, and penelope-like retrotransposons, have been characterized in several studies [54][55][56]. Some of these elements are still actively mobile within the genome [55], implying that they must be under tight regulation to maintain the stability of the genome [57]. We recently demonstrated that TE-derived siRNAs in S. japonicum were at least partially associated with the Argonaute protein, Ago2 [26]. In this study, the reads from the six libraries were mapped to the 29 well-defined TEs, including 18 LTR retrotransposons (SjCHGCS1-18), 7 non-LTR retrotransposons (SjCHGCS19-22, SjR1, SjR2, and Sjpido), 3 Penelope-like retrotransposons (Sj-penelope1-3), and Gulliver. The expression levels of siRNAs derived from these retrotransposons in the related libraries were presented based on their TPM value ( Table 2). We found that siRNAs derived from these 29 well-defined TEs were stably expressed at different developmental stages, though slightly suppressed in the male adult worms (Figure 7). However, the expression of siRNAs derived from different TEs was relatively diverse; for instance, the siRNAs originating from Gulliver, SjCHGCS18, SjCHGCS21, SjCHGCS22, and Sj-penelope3 were present nearly at undetectable levels, whereas those from SjCHGCS3 were highly expressed. In addition, the expression of both sense and antisense siRNAs from one specific TE showed a symmetrical pattern, with the exception of SjCHGCS17, from which the antisense siRNAs are predominant (Table 2). No major variation in expression level was observed in all libraries, supporting the earlier speculation that the main function of TEderived siRNAs was to maintain the genomic stability of S. japonicum [26].
In summary, we have further dissected the expression characteristics of the small RNAome in the egg stage of S. japonicum. Strong biased expression patterns of certain miRNA family members were observed, of which, the expression of sja-miR-71, sja-miR-71-5p, and sja-miR-36-3p were prominent in this pathologically-related stage. Transfer RNA (tRNA)-derived small RNA fragments, precisely processed mainly from the 59 side of tRNA transcripts, were identified for the first time as a novel class of small RNA in S. japonicum, which exhibited a significant stagebiased expression pattern, indicating their potential regulatory function in this stage. The most highly expressed tRF, Sj-tRF-001-1, has the potential to serve as an egg stage-specific bio-marker. The TE-derived siRNAs, which showed less variation in expression among different stages, also appeared to be an important constituent of the small RNA population, and is likely to protect the integrity of the genome against retroelements. The data in this study provide novel insights into the small RNAome of S. japonicum, which will facilitate a deeper understanding of the biology of this important parasitic pathogen.

Animals and Parasites
S. japonicum-infected Oncomelania hupensis snails were purchased from Jiangxi Provincial Institute of Parasitic Diseases, Nanchang, China. The cercariae were shed by exposing the snails to light conditions. A total of six New Zealand white rabbits were randomly assigned to two groups. Each rabbit was percutaneously infected with,1,200 cercariae. Hepatic schistosomula and mixed sex adult worms were recovered at 2 and 6 weeks post-infection, respectively, by hepatic-portal perfusion from the infected rabbits. Male and female worms were manually separated with the aid of a light microscope, and washed three times with phosphate buffered saline (PBS). Liver tissues were also obtained from the infected rabbits at 6 weeks post-infection. All procedures carried out on animals within this study were conducted following the animal husbandry guidelines of the Chinese Academy of Medical Sciences and with permission from the Experimental Animal Committee of the Chinese Academy of Medical Sciences with the Ethical Clearance Number IPB-2011-6.

Egg isolation
The schistosomal eggs were isolated by an improved sieving and enzymatic method [58]. The egg-trapped liver tissues were chopped with a scalpel blade and homogenated to a smooth consistency in 500 ml ice-cold PBS. The suspension was successively passed through 80, 120, 160, 200, and 260 mesh metal sieves, and finally a 320 mesh nylon screen. After repeated washes with PBS, the eggs on the nylon screen were collected in a 50 ml tube. The eggs were washed three times by discarding the tissue debris-containing suspension after natural sedimentation on ice. The pellet was resuspended in 50 ml PBS containing 10 mg collagenase B, 125 mg trypsin, 10 mg penicillin, and 20 mg streptomycin, then incubated at 37uC for 3 h with gentle shaking. The sample was then centrifuged at 1,500 rpm at 4uC for 5 min, and the supernatant was removed. This washing procedure was repeated twice more. The egg pellet was resuspended in 2 ml PBS and layered on the top of a Percoll column (containing a mixture of 2.4 ml of Percoll and 9.6 ml of 0.25 M sucrose) in a 15 ml tube. The tube was centrifuged at 2,000 rpm at 4uC for 5 min. Liver debris remaining in the supernatant was removed. The eggs were resuspended in 2 ml PBS for two more Percoll separations. The egg pellet was washed for 3 times with PBS then transferred to 1.5 ml tubes. The purity and integrity of the eggs was examined with the aid of a light microscope.

Total RNA preparation
After centrifuging at 12,000 rpm for 1 min, the egg pellet was ground in liquid nitrogen. Total RNA from the eggs was extracted using Trizol reagent (Invitrogen, CA, USA) according to the manufacturer's protocol. Total RNA from hepatic schistosomula, male adult worms, female adult worms, and normal rabbit livers were also isolated using Trizol reagent. RNA quantification and quality were examined with a Nanodrop ND-1000 spectrophotometer (Nanodrop Technologies, Wilmington, DE) and standard agarose gel electrophoresis. All RNA samples were stored at 280uC until use.

Small RNA library construction and deep sequencing
The egg total RNA sample was evaluated with an Agilent 2100 Bioanalyzer before library construction ( Figure S3). Three small RNA libraries were constructed as described previously [36]. Small RNAs between 15-30 nucleotides (nt) were recovered from a 15% TBE-Urea polyacrylamide gel electrophoresis (PAGE), and ligated into Illumina's proprietary 59 and 39 adaptors. The product was converted into single-stranded cDNA using Superscript III reverse transcriptase (Invitrogen, CA, USA). The cDNA was then amplified with Illumina's small RNA primer pair using Phusion high-fidelity DNA polymerase (NEB) in 18 PCR cycles. The purified PCR products were sequenced using Illumina's Genome Analyzer platform at the BGI (Beijing Genomics Institute, Shenzhen, China). Bioinformatics analysis of reads from different small RNAs libraries of S. japonicum Raw datasets were produced by deep sequencing of three libraries. After primary analysis, only the data generated from one library, designated as SjE, were further analyzed. The data were simultaneously analyzed with the previous small RNA datasets from cercariae, lung-stage schistosomula, hepatic schistosomula, male adult worms and female adult worms. First, the low quality reads, adaptor null reads, insert null reads, 59 adaptor contaminants, and reads with poly(A) tail were filtered. Adapter sequences were then trimmed from both ends of clean reads. Clean reads were obtained after all identical sequences were counted and merged as unique sequences These unique sequences were mapped onto the S. japonicum genome of SGST (http:// lifecenter.sgst.cn) using the program SOAP version 2.20 [59]. We investigated the length distribution of the perfectly matched small RNA reads in the six libraries [38]. Further, these small RNA reads were categorized using an optimal bioinformatic pipeline. In our previous study [38], analysis was focused mainly on miRNAs, whereas here, the unique sequences originating from snoRNAs (small nucleolar RNA) were in the first step filtered out [60], and sequences of rRNAs and tRNAs were investigated separately. In detail, 28S, 18S, and 5.8S rRNA sequences and rRNA intergenic spacer sequences (GenBank Accession Number: Z46504.4 [48], AY157226.1 [61], FJ852569.1 [62], and EU835685.1 [63]) of S. japonicum were retrieved from the NCBI GenBank database [64]. The putative tRNA gene sequences of S. japonicum were downloaded from http://www.bioinf.uni-leipzig. de/publications/supplements/08-014 [28]. Reads from different small RNA libraries mapped to those rRNA and tRNA sequences, other than to the data of Rfam as in the previous study [38], were respectively defined as rRNA-derived and tRNA-derived small RNAs. The secondary structures of mature tRNAs were predicted using an on line algorithm, ARAGORN [65]. The remaining perfectly matched reads were then BLAST-searched against the 77 known mature miRNAs of S. japonicum deposited in the Sanger miRBase [43,44] (Release 18) using the program Patscan [66], and were further BLAST-searched against the conserved and novel S. japonicum miRNAs reported in our previous study [38]. Next, the reads were matched to the transposable elements in the S. japonicum genome predicted by using REPET software (http://urgi. versailles.inra.fr/index.php/urgi/Tools/REPET), in the order of LINE (Long Interspersed Elements), SINE (Short Interspersed Elements), LTR (Transposable elements with Long Terminal Repeats), TIR (Terminal inverted repeat), MITE (Miniature inverted-repeat transposable elements), and unknown TEs. The remaining small RNAs were aligned to S. japonicum predicted mRNA sequences (sjr_mRNA.fasta) downloaded from SDSPB using SOAP 2.20 aligner, and perfectly matched reads were retained as mRNA-related siRNA. Finally, the remaining reads were labeled as unknown small RNAs. We employed IDEG6 [67] to identify miRNAs or tRFs showing statistically significant difference in relative abundance (as reflected by TPM, transcripts per million) between any two small RNA libraries. The general Chi-square test was applied to determine whether one particular miRNA or tRF was significantly differentially expressed between any two samples (P value, = 0.01) (Table S3 and S4).
To further characterize the small RNAome, full length sequences of 29 classes of retrotransposons [54,56,68,69] were retrieved from the NCBI GenBank database [64]. The small RNA reads from the libraries SjC, SL, SjH, SjF, SjM, and SjE were mapped to the sequences of rDNA repeat and the above mentioned retrotransposons. The abundance of rDNA-derived small RNAs or retrotransposon-derived siRNAs was reflected based on their expression values (TPM). A set of graphs depicting the distribution and abundance of these small RNAs was constructed as previously described [26]. All sequence data of the six small RNA libraries have been submitted to NIH Short Read Archive with the Accession numbers of SRR786675 (for SjE), SRR786666 (for SjC), SRR786671 (for SjL), SRR786672 (for SjH), SRR786673 (for SjM), and SRR786674 (for SjF).

Confirmation of sncRNA expression by Northern blot
The 59-DIG-labeled miRCURY LNA probes were synthesized by Exiqon (Vedbaek, Denmark) (Http://www.exiqon.com) (Table  S5). Northern blot analysis was performed as described previously [70]. Total RNAs (10 mg each) from different S. japonicum stages were resolved by 15% denaturing PAGE (7 M urea). The samples were then transferred to neutral nylon membranes (Hybond-NX, GE) by capillary with 206SSC, and cross-linked to the membrane using an EDC (1-ethyl-3-(3-dimethylaminopropyl) carbodiimide) method [71]. The blots were rinsed thoroughly with double distilled water and pre-hybridized at 37uC for 3 h in DIG Easy Granule (Roche). Hybridization was carried out in fresh DIG Easy Granule containing 1 nM DIG-labeled LNA probe at the recommended temperature (RNA Tm230uC) overnight. Blots were then washed sequentially in a low stringency buffer (26SSC, 0.1% w/v SDS) and a high stringency buffer (0.16SSC, 0.1% w/v SDS) at the hybridization temperature. After briefly rinsing in washing buffer, the blots were incubated in blocking buffer at room temperature for at least 2 h (DIG washing and blocking buffer Set, Roche). Subsequently, the blots were incubated with a 10,000-fold dilution of anti-DIG-AP Fab fragment (Roche) in blocking buffer at room temperature for 30 min then washed 5 times for 15 min each in washing buffer. Blots were then rinsed in detection buffer for 5 min. Anti-DIG-AP was detected using CDPstar chemiluminescent substrate for alkaline phosphatase (Roche). Blots were stripped by boiling for 1 min at 100uC in 10 mM Tris-HCl, 5 mM EDTA, and 0.1% SDS and probed up to three times.