Differential expression of microRNA between normally developed and underdeveloped female worms of Schistosoma japonicum

Eggs produced by bisexual infected mature female worms (MF) of Schistosoma japonicum are important in the transmission of the parasite and responsible for the pathogenesis of schistosomiasis. The single-sex infected female worms (SF) cannot mature and do not produce normal eggs; also they do not induce severe damage to the host. In this study, the microRNA (miRNA) expression profiles of 25d MF and 25d SF were investigated through Solexa deep-sequencing technology to explore the developmental mechanisms of schistosome female worms. There were 36 differentially expressed miRNA, 20 up-regulated and 16 down-regulated found in MF/SF worms, including some development related miRNA such as bantam (ban), let-7, miR-124, miR-8, miR-1, miR-7. There were 166 target genes of up-regulated miRNA and 201 target genes of down-regulated miRNA after comparing the target gene prediction software results with RNA-Seq transcriptome results. Analysis of the target genes shows that different ones are involved in MF and SF worms in Gene Ontology terms, with a similar situation in KEGG. This observation indicates that different genes regulated by differentially expressed miRNA take part in MF and SF and lead to differential sexual status. This means that the sexual status of female worms is regulated by miRNA.


Introduction
Schistosomiasis is an important helminth infection that mainly occurs in developing countries. It is endemic in 74 countries, with roughly 120 million individuals being severely affected [1]. No successful vaccine is currently available for this disease. Praziquantel is the only and most extensively used drug for treatment, but it is ineffective against young worms [2] and there are concerns that drug resistance to praziquantel may be developing [3]. Consequently, there is an urgent need to identify novel drug targets and/or develop alternative therapeutic strategies.
In schistosomiasis, the severely pathological injuries to the final host are caused by the amounts of eggs produced by mature normally developed female worms (MF). These eggs are also responsible for the spread of the disease [4], but the underdeveloped female worms (SF) which come from single-sex infection cannot develop normally or produce normal eggs. So these single-sex infected worms (female or male) will not seriously damage the host, and will not induce the prevalence of schistosomiasis. Based on the important role of eggs produced by mature females, the developmental mechanisms of female worms is the research target for many schistosome researchers.

Open Access
*Correspondence: yameijin@shvri.ac.cn 1  MicroRNA (miRNA) are a class of small non-coding RNA (18-25 nucleotides in length) that are generated from endogenous transcripts [5]. They function as guide molecules for post-transcriptional gene regulation by base-pairing to their messenger RNA (mRNA) targets [6], they bind to the target and induce mRNA translation repression and/or exonucleolytic decay. Through these actions, miRNA regulate gene expression during development, differentiation, proliferation, death, and metabolism in many organisms [7]. To date, several groups have carried out studies to identify miRNA in schistosomes to provide some insight into the role of miRNA in schistosome development [8][9][10].
In this study, the expression profiles of miRNA of 25d MF and SF Schistosoma japonicum (S. japonicum) were investigated. Based on the analysis of the differentially expressed miRNA and their predicted target genes, the correlation between miRNA and the mature status of female worms was studied.

Ethics statement
All animal experiments were conducted in accordance with the guidelines of the Committee for the Care and Use of Laboratory Animals of Shanghai Veterinary Research Institute, Chinese Academy of Agricultural Sciences (permit no. SYXK-2016-0010). The study protocol was approved by the Ethics and Animal Welfare Committee of the Shanghai Veterinary Research Institute, Chinese Academy of Agricultural Sciences.

Host animals and parasites
BALB/C male mice, 4-to 6-week-old, were used for worm collection. S. japonicum (Chinese strain) worms were maintained in mice and Oncomelania hupensis hupensis snails at the Shanghai Veterinary Research Institute, Chinese Academy of Agricultural Sciences (SHVRI, CAAS). Cercariae of S. japonicum were shed from infected snails after exposure to artificial light for 5 h at 25-28 °C. The 25d MF were collected from bi-sexual infected mice at 25 days post-infection by perfusion. Worms were washed five times in phosphate buffered saline (PBS, pH 7.4) at 4 °C to remove excreted and secreted proteins and other interfering substances and were then incubated at 4 °C for 30 min. The paired worms spontaneously separated into females and males, then the female worms were collected. The 25d SF were obtained after perfusion from mice infected by female cercariae. The sexes of cercariae shed from individual snails were identified by polymerase chain reaction (PCR) based on the patent: CN 101597645. Then the female worms were washed five times in PBS at 4 °C to remove excreted and secreted proteins and other interfering substances, then collected.

RNA extraction, small RNA library construction and sequencing
Total RNA was extracted using TRIzol reagent (Invitrogen Life Technologies, Carlsbad, CA, USA) according to the manufacturer's instructions. The quality of RNA was measured using a NanoDrop 1000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA). The RNA samples were stored at -80 °C.
The construction of small RNA libraries was carried out using the TruSeq Small RNA Sample Prep Kit (Illumina, San Diego, CA, USA) according to the manufacturer's protocol. Briefly, total RNA was subjected to polyacrylamide gel electrophoresis (PAGE) and then the 18-30 nt size range of RNA was enriched for purification. The purified small RNA was subjected to sequencing with an Illumina Genome Analyzer (Illumina) according to the manufacturer's instructions. Briefly, proprietary adapters were ligated to the 5ʹ and 3ʹ termini of the small RNA, and the ligated small RNA were then used as templates for cDNA synthesis, and further converted into 62-75 nt single-stranded cDNA. The cDNA were amplified using illumina's 3′ adaptor reverse primer and 5′ adaptor forward primers. The purified PCR products were sequenced using Solexa's proprietary sequencingby-synthesis method. DNA sequencing was performed with an Illumina Genome Analyzer at the BGI (Beijing Genomics Institute, Shenzhen, China).

Data processing and mapping sequence reads to the reference genome
All raw datasets produced by deep sequencing from each library (25d MF, 25d SF) were subjected to further analysis. Clean reads were obtained after removal of low quality reads, insert null reads, adaptor null reads, 5′ adaptor contaminants, reads with poly (A) tails and reads shorter than 18 nt. Adapter sequences were then trimmed from both ends of the clean reads. The small RNA with sizes ranging from 18-30 nt were used for further analyses. The clean reads in each library were mapped onto the draft S.japonicum genome sequences (ftp://down:lsbi@ lifec enter .sgst.cn:2121/subje ctDat a/schis tosom a/) and (ftp://ftp.sange r.ac.uk/pub/patho gens/Schis tosom a/ manso ni/genom e/smav5 .2.chr.fa) using the Short Oligonucleotide Alignment Program (SOAP). The small RNA that matched with known rRNA, tRNA, small nuclear RNA (snRNA), and small nucleolar RNA (snoRNA) deposited in the Rfam database (ftp://selab .janel ia.org/ pub/Rfam/) and NCBI GenBank (https ://www.ncbi. nlm.nih.gov/GenBa nk) were excluded. Then, we aligned the predicted miRNA to known miRNA or homologous miRNA in the miRBase (version 21, https ://micro rna. sange r.ac.uk) to identify conserved miRNA Other small RNA were further analyzed to predict new miRNA by exploring the secondary structure, the Dicer cleavage site and minimum free energy of the unannotated small RNA tags that could be mapped to the genome using MIREAP (https ://sourc eforg e.net/proje cts/mirea p).

Bioinformatics analysis of S. japonicum small RNA
Small RNA tags were aligned to the miRNA precursor/mature miRNA of corresponding species in miR-Base20. Detailed information of the alignment, including the structure of known miRNA precursors, and the length and count of tags from the sample, among others, were collected. To make every unique small RNA mapped to only one annotation, we followed the following priority Ruler RNA etc.: rRNA etc. (in which Genbank > Rfam) > known miRNA > repeat > exon > intron 3. The total rRNA proportion should be less than 40%. The expression of the miRNA in two samples were visualized by plotting the Log2-ratio figure and the scatter plot. The procedures adopted were as follows: (1) the expression of miRNA in two samples (25d MF and 25d SF) were normalized to get the expression in transcripts per million (TPM) with the normalization formula (normalized expression = actual miRNA count/total count of clean reads × 1 000 000); and (2) the fold-change and P-value were calculated from the normalized expression, and log2-ratio and scatter plots were generated using the fold-change formula (fold change = log2 (25d MF/25d SF)).

MiRNA quantification by quantitative RT-PCR analysis
Differentially expressed miRNA were examined using qPCR with SYBR green. The first cDNA strand was synthesized from 0.2 μg total RNA using PrimeScript ® RT reagent Kit (Perfect Real Time, TaKaRa, Otsu, Japan, code: DRR037A). Stem-loop qRT-PCR was performed to quantify the sex-biased expressed miRNA. Stem-loop RT primers (Table 1) were used to reverse-transcribe mature miRNA to cDNA. The RNA templates for the qPCR were performed on the same samples used for Solexa deep sequencing. The 20 μL reaction RT system contained 2 μL of total RNA (0.2 μg), 0.4 μL (10 μM) of each individual stem-loop RT primer, 0.4 μL ROX Reference Dye II (50 ×), 10 μL SYBR Premix Ex Taq (2 ×), and 6.8 μL Easy Dilution. The RT primers of the miRNA, the forward primers of miRNA and the common reverse primer were designed as shown in Table 1. The cDNA were synthesized by incubation for 30 min at 16 °C, 45 min at 37 °C, and 15 s at 85 °C. The products were amplified using SYBR ® Premix ExTaq (TliRNaseH Plus, TaKaRa, code: RR420A) in an ABI Prism 7500 sequence detection system (Applied Biosystems, Waltham, MA, USA) according to the manufacturer's instruction. Our system contained 2 μL of cDNA from the RT reaction product, 5 μL of 2 × SYBR Premix ExTaq (TaKaRa), and 0.4 μL of 10 μM forward and reverse primers. The NADPH gene was used in the qRT-PCR reaction as an endogenous reference [8]. Relative levels of gene expression were calculated using the 2 − △△Ct method. All assays were performed in triplicate.

Prediction and bioinformatics analysis of differentially expressed miRNA target genes
The regulatory effect of differentially expressed miRNA in 25d MF and 25d SF were predicted by analyzing the functions of their predicted target genes and by comparing the proportion of their predicted target genes in differentially expressed genes. In this study, RNA hybrid (https ://140.109.42.4/cgi-bin/RNAhy brid/ RNAhy brid.cgi) [11] were used to predict the target genes of differentially expressed miRNA. The predicted target genes were filtered with minimum free energy (MFE < = -20 kcal/mol). These target genes were refiltered by comparison with the RNA-Seq (Quantification) transcriptome results of SF and MF. The worms were the same sets used in miRNA deep sequencing (Data not published); the genes that were not consistent in the gene expression profile and miRNA target gene profile were excluded. Then those corresponding to target genes were mapped to the terms in the Gene Ontology (GO) database, to identify significantly enriched GO terms in DEG (differentially expressed genes) compared with the genome background. Pathway enrichment analysis was used to identify significantly enriched metabolic pathways or signal transduction pathways in DEG compared with the whole genome background. Detailed pathway information was determined with the KEGG database.

Host animals and parasites
Cercariae were collected from individual snails and their sex was identified by PCR. There is a 153 bp specific band from female S. japonicum cercariae, and no band from male S. japonicum cercariae (Figure 1). For those snails that contain two sex cercariae, the SF female worms were guaranteed by sacrificing the mice that were infected by cercariae from the same snail when the worms were collected.

Preliminary analysis of raw reads of miRNA in S. japonicum 25d MF and 25d SF
The cDNA libraries derived from 18-30 nt long RNA isolated from 25d MF and 25d SF were constructed and sequenced using an Illumina (Solexa) DNA sequencer, yielding a total of 9 358 376 and 9 314 275 clean reads corresponding to 97.74% and 97.29% of high-quality sequence for S. japonicum 25d MF and 25d SF, respectively ( Table 2). The length distribution of the miRNA from the two libraries mainly ranged from 20-24 nt (Figure 2).

Predictions and bioinformatics analysis of the predicted target genes for differentially expressed miRNA
To analyze the effect of the differential expression of miRNA on 25d MF and 25d SF, the target genes of the differentially expressed miRNA between the two samples were predicted and analyzed. The target genes of up-or down-regulated miRNA were predicted and analyzed separately. Based on the target gene prediction software, 4728 and 3906 target genes of up-and downregulated miRNA were acquired respectively. After comparison with the RNA-Seq transcriptome results, there were 166 target genes of up-regulated miRNA and 201 target genes of down-regulated miRNA left respectively. The bioinformatics analysis for target genes of either up-or down-regulated miRNA was applied separately. Not all the target genes can be applied to GO terms, for example, only 41 out of 166 target genes of up-regulated miRNA and 76 out of 201 target genes of down-regulated miRNA were logged in the GO terms of molecular function (Additional file 1). The target genes of up-or downregulated miRNA were mainly concentrated on the same GO terms (Additional file 2 and Additional file 3). For molecular-function, they were both concentrated on binding, catalytic activity, molecular transducer activity, transporter activity ( Figure 5). But the specific target genes of each GO term were different for up-or downregulated miRNA and there were also target genes of up-and down-regulated miRNA concentrated on different GO terms ( Figure 5). For biological process and cellular component GO terms, the same situation was found (Additional file 4 and Additional file 5).
Some target genes of up-or down-regulated miRNA were selected to analyze their functions. The results show that target genes were not only regulated by only one miRNA, such as the target gene Sjc_0005980, but that a target gene can be involved in more than one GO term (Additional file 6 and Additional file 7). This means they participate in multiple important biological processes.
Pathway enrichment analysis also testified to this situation. KEGG results show that some target genes of up-and down-regulated miRNA take part in the same pathway, but the individual target genes in up-and down-regulated miRNA are different. In addition also some target genes of up-or down-regulated miRNA are involved in various pathways, however, only target genes of up-regulated miRNA are involved in these pathways, such as starch and sucrose metabolism, carbohydrate digestion and absorption, TGF-beta signaling pathway, chemokine signaling pathway, and so on. For the pathways, protein processing in the endoplasmic reticulum, N-glycan biosynthesis, ribosome biogenesis in eukaryotes, etc., only target genes of down-regulated miRNA take part (Table 4).

Discussion
The MF worms of S. japonicum can develop to sexual maturation and maintain this mature state by pairing with male worms [14]. The SF worms were considerably different in size and structure of inner tissues or systems from those MF worms [15]. Although the long term infected SF can lay very few eggs, these non-viable eggs do not contain developing or live miracidia so there are no soluble miracidial antigens that can induce the host's immune response, therefore the pathological reaction in infected final host and the schistosomiasis spreading in the environment will be relieved or stopped [16]. During the development of S. japonicum, males and females begin to pair at 15 or 16 d post-infection, and females begin to lay eggs from 24 d post-infection [14]. So S. japonicum will develop totally mature worms at 25 d post-infection. In this study, the miRNA profiles of 25d MF and 25d SF were compared, some differentially expressed miRNA in the two samples were found, and some clues about the correlation between differentially expressed miRNA and the mature status of schistosome female worms appeared by analyzing the differentially expressed miRNA and their target genes.
The adult female and male worms have obvious gender dimorphism, which can be identified by naked eyes, but the sex of miracidium, sporocyst and cercariae cannot be identified in this way, as they do not have obvious gender dimorphism. One snail can be infected by several miracidia, but in the majority of snails (65.7%), only one miracidium can develop into a mother sporocyst [17]. The proportion of single-or mixed-sex infections of parasites within snail individuals in the field has been reported for S. japonicum, most cercariae shed from an individual snail were single-sex [18]. Based on this situation, in this study, the sex of cercariae overflow from an individual snail was confirmed by PCR and infected mice were sacrificed, to ensure that SF worms were recovered.
There are significant morphological differences between normal developed and underdeveloped female S. japonicum. In our previous research, the underdeveloped female S. japonicum was significantly shorter and thinner than normal developed, and the mature ootid in the indistinct ovary of the underdeveloped female could not be observed [19]. These distinct morphological changes may be associated with miRNA regulation, in this study; there are 20 up-regulated and 16 downregulated miRNA in 25d MF. The differentially expressed  [20].
Bantam is a miRNA that is functionally well characterized in Drosophila. In Drosophila, it is a developmentally regulated miRNA that controls cell proliferation and regulates the proapoptotic gene hid [21]. Bantam was also identified with genes promoting organ growth and systemic body growth by connecting insulin signaling, and ecdysone production [22]. Boulan et al. found that lower ban activity induced growth cessation of late ecdysoneproducing cells that led to the cells entry into a metamorphosis stage, whereas ban's high activity promotes the larvae's systemic growth in Drosophila [22]. In our study, ban was also up-regulated in developed females and down-regulated in underdeveloped females (Table 3). These phenomena are also observed in Caenorhabditis elegans, the deletion of ban orthologs in Caenorhabditis elegans cause severe growth defects and prevent entry into the dauer state [23].
In this study, two miRNA, miR-7 and miR-1, which have been characterized as tumor suppressors in many kinds of human cancers, were up-regulated in 25d SF. They were found to play a key role in cell proliferation and survival either by being involved in the insulin-like growth factor 1 (IGF1) receptor signaling pathway or epidermal growth factor (EGF) receptor signaling pathway. EGF is an early growth response gene that is induced by mitogenic stimuli, including serum and purified factors such as IGF, EGF, and TGF-beta [24]. Let-7 and miR-124 were also up-regulated in 25d SF in our study. The Let-7 miRNA family has an important function in cell differentiation and cancer [25]. Worringer et al. suggested that let-7 based pathways are an inhibitory influence on the  reprogramming process through a regulatory pathway involving prodifferentiation factors, including EGF receptor 1 [26]. MiR-124 also plays a key role in cell differentiation, especially in the nervous system [27]; it was also thought to be a growth suppressor [28]. MiR-8, recognized as a critical mediator for several well characterized developmental regulatory networks such as Wnt, Notch, and TGF-beta in Drosophila or mammalian cells [29,30], is up-regulated in 25d SF too. In addition, the Wnt, Notch, and TGF-beta pathways are shown to be involved in the parasite development, oogenesis and embryogenesis in schistosoma mansoni [31][32][33][34]. The up-regulation of these differentiation-associated miRNA in 25d SF may be related to the underdevelopmental status of SF. Sexual immaturity may mean more potential for differentiation when a single female pairs with her mate.
The bioinformatics analysis of the predicted target genes for differentially expressed miRNA shows that the target genes distributed in different GO terms, which means these differential functional genes, take part in many important biological processes. The KEGG results show that more target genes in 25d MF are involved in the regulation of actin cytoskeleton, focal adhesion, protein digestion and absorption, tight junction, and calcium signaling pathway metabolic processes whereas more target genes in 25d SF are involved in lysosomes, oocyte meiosis, long-term potentiation and salivary secretion pathways. Only target genes of 25d MF were detected in starch and sucrose metabolism, galactose metabolism, carbohydrate digestion and absorption, axon guidance, TGF-beta signaling pathway, and chemokine signaling pathways whereas only target genes in 25d SF were found in pathways involved in protein processing in the endoplasmic reticulum, N-glycan biosynthesis, ribosome biogenesis eukaryotes, gap junction, dopaminergic synapse, spliceosome, RNA transport, GnRH signaling pathway, and long-term depression. All these are associated with the developmental stages of female worms.
Some receptors in the Wnt, Notch, IGF, EGF, TGFbeta signaling pathways which interact with miRNA have been characterized in both S. mansoni and S. japonicum [35,36], these signaling pathways have also been implicated in critical roles in schistosome development by the S. japonicum genome study [37], but direct evidence for miRNA regulated signaling pathways remains to be experimentally clarified. This study shows that some miRNA could regulate the sexual development of female worms and that it analyzed target genes. This will shed some clues on the understanding of the development mechanism of schistosome female worms.  Table 4 Comparison of predicted target genes participating in regulating metabolic processes between up-regulated miRNAs and down-regulated miRNAs in 25d MF.