Comparative Transcriptome Analyses of Schistosoma japonicum Derived From SCID Mice and BALB/c Mice: Clues to the Abnormality in Parasite Growth and Development

Schistosomiasis, caused by the parasitic flatworms called schistosomes, remains one of the most prevailing parasitic diseases in the world. The prodigious oviposition of female worms after maturity is the main driver of pathology due to infection, yet our understanding about the regulation of development and reproduction of schistosomes is limited. Here, we comparatively profiled the transcriptome of Schistosoma japonicum recovered from SCID and BALB/c mice, which were collected 35 days post-infection, when prominent morphological abnormalities could be observed in schistosomes from SCID mice, by performing RNA-seq analysis. Of the 11,183 identified genes, 62 differentially expressed genes (DEGs) with 39 upregulated and 23 downregulated messenger RNAs (mRNAs) were found in male worms from SCID mice (S_M) vs. male worms from BALB/c mice (B_M), and 240 DEGs with 152 upregulated and 88 downregulated mRNAs were found in female worms from SCID mice (S_F) vs. female worms from BALB/c mice (B_F). We also tested nine DEGs with a relatively higher expression abundance in the gonads of the worms (ovary, vitellaria, or testis), suggesting their potential biological significance in the development and reproduction of the parasites. Gene ontology (GO) enrichment analysis revealed that GO terms such as “microtubule-based process,” “multicellular organismal development,” and “Rho protein signal transduction” were significantly enriched in the DEGs in S_F vs. B_F, whereas GO terms such as “oxidation–reduction process,” “response to stress,” and “response to DNA damage stimulus” were significantly enriched in the DEGs in S_M vs. B_M. These results revealed that the differential expression of some important genes might contribute to the morphological abnormalities of worms in SCID mice. Furthermore, we selected one DEG, the mitochondrial prohibitin complex protein 1 (Phb1), to perform double-stranded RNA (dsRNA)-mediated RNA interference (RNAi) in vivo targeting the worms in BALB/c mice, and we found that it was essential for the growth and reproductive development of both male and female S. japonicum worms. Taken together, these results provided a wealth of information on the differential gene expression profiles of schistosomes from SCID mice when compared with those from BALB/c mice, which were potentially involved in regulating the growth and development of schistosomes. These findings contributed to an understanding of parasite biology and provided a rich resource for the exploitation of antischistosomal intervention targets.

1 School of Basic Medical Sciences, Wuhan University, Wuhan, China, 2 Laboratory Animal Center, School of Medicine, Wuhan University, Wuhan, China Schistosomiasis, caused by the parasitic flatworms called schistosomes, remains one of the most prevailing parasitic diseases in the world. The prodigious oviposition of female worms after maturity is the main driver of pathology due to infection, yet our understanding about the regulation of development and reproduction of schistosomes is limited. Here, we comparatively profiled the transcriptome of Schistosoma japonicum recovered from SCID and BALB/c mice, which were collected 35 days post-infection, when prominent morphological abnormalities could be observed in schistosomes from SCID mice, by performing RNA-seq analysis. Of the 11,183 identified genes, 62 differentially expressed genes (DEGs) with 39 upregulated and 23 downregulated messenger RNAs (mRNAs) were found in male worms from SCID mice (S_M) vs. male worms from BALB/c mice (B_M), and 240 DEGs with 152 upregulated and 88 downregulated mRNAs were found in female worms from SCID mice (S_F) vs. female worms from BALB/c mice (B_F). We also tested nine DEGs with a relatively higher expression abundance in the gonads of the worms (ovary, vitellaria, or testis), suggesting their potential biological significance in the development and reproduction of the parasites. Gene ontology (GO) enrichment analysis revealed that GO terms such as "microtubule-based process," "multicellular organismal development," and "Rho protein signal transduction" were significantly enriched in the DEGs in S_F vs. B_F, whereas GO terms such as "oxidation-reduction process," "response to stress," and "response to DNA damage stimulus" were significantly enriched in the DEGs in S_M vs. B_M. These results revealed that the differential expression of some important genes might contribute to the morphological abnormalities of worms in SCID mice. Furthermore, we selected one DEG, the mitochondrial prohibitin complex protein 1 (Phb1), to perform double-stranded RNA (dsRNA)-mediated RNA interference (RNAi) in vivo targeting the worms in BALB/c mice, and we found that it was essential for the growth and reproductive development of both male and female S. japonicum worms. Taken together, these results provided a wealth of information on the differential gene INTRODUCTION Schistosomiasis, caused by infection with the blood-dwelling endoparasites of the genus Schistosoma (the only digenetic trematodes), remains one of the most serious parasitic diseases worldwide, afflicting more than 200 million people and generating annual losses of 1.7-4.5 million disability-adjusted life years (DALYs) of humans, with about 800 million at risk (King et al., 2005;Steinmann et al., 2006;Gray et al., 2010;Ross et al., 2013;WHO, 2016). Schistosomes have a complex life cycle involving both a snail intermediate and a vertebrate definitive host and have coevolved a long-lived and intricate relationship with their human hosts (Harrison and Doenhoff, 1983;Amiri et al., 1992;Nakazawa et al., 1997;Cheever et al., 1999;Hernandez et al., 2004;Blank et al., 2006;Liu et al., 2006;deWalick et al., 2012), as well as an indispensable interplay between the adult male and female parasites during their development and maturation (Den Hollander and Erasmus, 1985;Cai et al., 2016;. As the only members of the trematodes, schistosomes have evolved to form separate sexes-dioecious (Platt and Brooks, 1997;Kunz, 2001). Furthermore, mating is required to induce sexual development and maturation of the female worms as a prerequisite for egg production (Kunz, 2001;Ross et al., 2002;Hu et al., 2004;Quack et al., 2006;Beckmann et al., 2010;Lu et al., 2016). The latter finally leads to inflammatory processes in the affected organs of their host (such as the gut and liver of the hosts infected with Schistosoma japonicum or Schistosoma mansoni parasites) due to egg deposition, culminating in the formation of granulomas and fibrosis (Ross et al., 2002). The clinical signs of schistosomiasis are dependent on the maturation stage of the parasites and their eggs. In humans, acute infection is characterized by a debilitating febrile illness (Katayama fever) that usually occurs before the appearance of eggs in the stool, with a peak at 6-8 weeks post-infection. In chronic diseases, eggs trapped in various tissues induce the formation of granulomatous inflammation, which subsequently results in fibrosis that causes many of the pathological conditions. It is believed that parasites definitely cause damage to their host, and the host resists against the invading parasites by its immune response (Han et al., 2009;The Schistosoma japonicum Genome Sequencing and Functional Analysis Consortium, 2009;deWalick et al., 2012;Duan et al., 2014;Morais et al., 2018;Roderfeld et al., 2018). However, several studies have already observed that host factors such as the immune factors can modulate parasite development and maturation (Amiri et al., 1992;Hernandez et al., 2004;LoVerde et al., 2009;Beckmann et al., 2010;deWalick et al., 2012). Moreover, several studies have discovered that schistosomes showed retarded growth, development, and reproduction in immunodeficient mammalian hosts, leading to attenuated pathogenesis to hosts with decreased ovulation, such as attenuated hepatic granulomas and fibrosis formation (Amiri et al., 1992;Davies et al., 2004;Cheng et al., 2008;Tang et al., 2013). We then ask how the immunodeficiency of SCID mice could lead to the retarded growth and development of the schistosomes residing in them and what the significant changes in the gene expression of worms under their aberrant morphological changes are. However, few studies are available investigating the molecular regulation of schistosomes that is involved in their abnormal growth and development in immunodeficient hosts, although some studies focusing on the host's factors are available (Davies et al., 2001;Blank et al., 2006;Lamb et al., 2010).
Numerous researchers had endeavored to explore the genes participating or involved in the regulation of the development and reproductive maturation of schistosomes by comparative transcriptome, proteome, and metabolome research and gained great discoveries in the growth and reproductive biology of schistosomes (Gupta and Basch, 1987;Haseeb et al., 1989;Siegel and Tracy, 1989;Haseeb and Eveland, 1991;Grevelding et al., 1997;Haseeb, 1998;Liu et al., 2006;Cogswell et al., 2012;Leutner et al., 2013;Protasio et al., 2013;Sun et al., 2014;Han et al., 2015a,b;Cao et al., 2016;Picard et al., 2016;Zhai et al., 2018;Li et al., 2019;Mao et al., 2019). Recently, we have also investigated and reported the differential metabolomic profiles of S. japonicum parasites from SCID and BALB/c mice in order to explore the molecular events involved in the aberrant morphologies in the growth and reproduction of schistosomes from a metabolomic level. Some important metabolic clues to what is involved in the morphological abnormalities in the growth and reproductive development of schistosomes from SCID mice have been discovered, such as the potential influence of the meiotic process due to a disturbed retinol metabolism, which is involved in meiosis initiation . A comprehensive understanding of the gene expression alterations would also help us find more clues about the molecular basis associated with the morphological alterations of S. japonicum parasites from SCID mice. Synchronously, here, we further explored and report the differentially expressed transcriptomes of S. japonicum parasites from SCID and BALB/c mice by RNA sequencing (RNA-seq) . We have identified some differentially expressed genes with higher expression levels in the gonads of the worms, i.e., ovary and vitellaria of female worms, and testes of male worms, suggesting their potentially great biological significance in the development and reproduction of schistosomes. Moreover, we also report the double-stranded RNA (dsRNA)-mediated RNA interference (RNAi) in vivo of the mitochondrial prohibitin complex protein 1 (Phb1) of the worms. S. japonicum Phb1 was found highly expressed in the vitellaria of female worms and was found essential for the growth and reproductive development of both male and female S. japonicum worms.
These results, to our knowledge, provide a wealth of information on the differential gene expression profiles of schistosomes from SCID mice when compared with those from BALB/c mice. These findings also provide a rich and valuable resource for schistosomes and schistosomiasis research, such as exploring potential anti-schistosome drug and vaccine targets.

Parasite Collection for RNA-seq
Schistosoma japonicum-infected O. hupensis were purchased from the Institute of Parasitic Disease Control and Prevention, Hunan Province, China. Six-to 7-week-old female BALB/c and SCID mice in SPF grade were purchased from Beijing Hua Fu Kang BioScience Co., Ltd. through WUCAE. The mice were maintained by WUCAE and were allowed to stay at a room temperature of 25 • C with free access to food and water for 12 days in order to adapt to the new environment. The release of S. japonicum cercariae, artificially infecting mice, and the collection of worms on the 35th day post-infection were performed according to the methods described in our previously published article . After being washed with phosphate-buffered saline (PBS), the female and male worms were collected separately and equal amounts of parasites (with 20 worms in each sample) were divided into several sterile cryopreservation tubes and stored in TRIzol reagent (Invitrogen, United States). The worm samples were labeled as "S_F" for the female worms from SCID mice, "S_M" for the male worms from SCID mice, "B_F" for the female worms from BALB/c mice, and "B_M" for the male worms from BALB/c mice. The worms were then stored provisionally at −80 • C until further use.

Total RNA Extraction, Library Construction, and Sequencing
Four schistosome samples (labeled as S_F, S_M, B_F, and B_M and containing 20 worms in each sample) were separately homogenized in glass homogenizers, and their total RNAs were isolated using TRIzol reagent according to the manufacturer's instructions (Invitrogen, United States).
The potential contaminating genomic DNA was removed from the RNA samples using recombinant DNase I (RNase-free) according to the manufacturer's protocol (Takara, China). The prepared RNA samples were packed in dry ice and delivered to the Beijing Compass Biotechnology Co., Ltd. (Beijing, China) for commercial RNA-seq analysis. Before library construction, the integrity of the RNA samples was assessed qualitatively by denatured agarose gel electrophoresis and quantitatively by an Agilent 2100 Bioanalyzer, which reported no significant degradation of the total RNA samples. The purity and the quantity of the total RNA were measured by using a Nanodrop 2000 spectrophotometer (NanoDrop Technologies, Wilmington, DE, United States). RNA quality assessment was performed by a professional technician of the Beijing Compass Company.
For transcriptome library construction, 3 µg of the total RNA from each sample was used as the input material for the subsequent RNA sample preparations. Ribosomal RNA (rRNA) was depleted using the Epicentre R Ribo-zero rRNA Removal Kit. Sequencing libraries were generated and RNA sequencing was performed according to the manufacturer's protocol using the NEBNext R Ultra Directional RNA Library Prep Kit for Illumina (Illumina R , NEB, United States). Briefly, first-strand complementary DNA (cDNA) was synthesized using a random hexamer primer and M-MuLV Reverse Transcriptase. Secondstrand cDNA synthesis was subsequently performed using DNA polymerase I and RNase H. In the reaction buffer, dNTPs with dTTP were replaced by dUTP. After end repair, adenylation of the 3 end of the DNA fragments, and adapter ligation using the NEBNext Adaptor with a hairpin-loop structure, each library was enriched by 15 cycles of PCR, after which the library fragments were purified with the AMPure XP system (Beckman Coulter, Beverly, MA, United States) in order to obtain the strandedspecific library with select cDNA fragments of the preferred 150to 200-bp length. The quality of the libraries was assessed on an Agilent Bioanalyzer 2100 system. The libraries were quantified with a Qubit 2.0 Fluorometer (Life Technologies, Grand Island, NY, United States) and sequencing was performed on an Illumina Hiseq TM 2500 platform (Illumina).
Library construction and sequencing was performed by professional technicians of the Beijing Compass Biotechnology Co., Ltd 1 .

Data Processing
Firstly, the raw reads of the transcriptomic sequencing were separated based on the barcodes and processed to remove adapters, poly-N-containing reads (N% > 10%), and lowquality reads (quality score ≤ 20) using custom perl scripts. Because the schistosomes were obtained from experimentally infected mice, the reads were aligned to the mouse genome using Bowtie 2 (v0.12.7) and Tophat (v2.0.0) (Trapnell et al., 2009;Langmead and Salzberg, 2012;Kim et al., 2013). Non-aligned reads were considered to be S. japonicumspecific. Clean reads were matched to the transcriptomes of S. japonicum (PRJEA34885) in WormBase ParaSite 2 and assembled using Trinity, with default parameters . Transcripts ≤ 150 bp in length were discarded. Finally, the reconstructed unigenes were annotated by blastx to public databases [the Swiss-Prot protein database, Caenorhabditis elegans transcriptome, Clonorchis sinensis transcriptome, the UniProt protein database, S. mansoni transcriptome, Schistosoma haematobium transcriptome, S. japonicum V4.0, and the NCBI non-redundant (nr) database] with a threshold of 10 −10 . The number of clean reads for each transcript was calculated and then normalized to fragments per kilobase of exon model per million reads (FPKM), which associates the read numbers with the gene expression levels. Corset was applied to hierarchically cluster short transcripts into long genes for the downstream analysis (Davidson and Oshlack, 2014). The assembled transcriptome was annotated by seven public databases, including Nr, Nt, Pfam 3 (Bateman et al., 2002), KOG/COG 4 , Swiss-Prot 5 , Kyoto Encyclopedia of Genes and Genomes (KEGG 6 ) (Kanehisa et al., 2008), and Gene ontology (GO 7 ). The Coding-Non-coding Index (CNCI), Coding Potential Calculator (CPC), and Pfam were used to distinguish mRNAs from long non-coding RNAs (lncRNAs) (Bateman et al., 2002;Kong et al., 2007;Sun et al., 2013). CNCI profiles can effectively distinguish protein-coding and non-coding sequences independent of the known annotations by adjoining nucleotide triplets (Sun et al., 2013). CPC searches the sequences with the known protein sequence database to clarify the coding and non-coding transcripts mainly through assessing the extent and quality of the open reading frame (ORF) in a transcript (Kong et al., 2007). Each transcript can be translated in all three possible frames and the Pfam Scan (v1.3) used to identify the occurrence of any of the known protein family domains documented in the Pfam database (Bateman et al., 2002). Transcripts without coding potential were a candidate set of lncRNAs, which will be reported and discussed elsewhere.

Profiling of Differentially Expressed Genes and Gene Annotation
The quantification of differentially expressed genes (DEGs) in each sample was calculated by Cuffdiff (v2.1.1) (Trapnell et al., 2012), and transcripts with a P-adjust < 0.05 were considered differentially expressed. The gene sequences were functionally annotated using Blast2GO at www.blast2go.org (Gotz et al., 2008;Liu et al., 2015), and the output provided combined graphics for the three categories of gene ontology (GO) terms: biological processes, molecular functions, and cellular components. The KEGG automated annotation server 8 was used to assign pathwaybased functional orthology to the DEGs (Moriya et al., 2007;Kanehisa et al., 2008), in which a P-value less than 0.05 was considered to be statistically significant. Signal peptides were predicted using the program SignalP 4.1 server 9 employing both the neural network and hidden Markov model (Nielsen, 2017), and transmembrane helices were predicted using TMHMM 2.0 10 (Chen et al., 2003).

Validation of Differential Expression by qRT-PCR
A subset of the DEGs with significant biological significance from the list of differential transcriptomes was selected for further validation by quantitative reverse transcriptase PCR (qRT-PCR). The total RNA of the above four schistosome samples collected in the same experiment was isolated using TRIzol reagent (Invitrogen, United States) and treated with recombinant RNasefree DNase I (Takara, China) according to the manufacturer's instructions. For each sample, 1 µg of the total RNA was used to synthesize the first-strand cDNA using a Reverse Transcriptase Kit (TaKaRa, Dalian, China) with oligo(dT) 18 primers in a final volume of 20 µl, followed by dilution of the transcribed products by adding another 20 µl ddH 2 O. qRT-PCR was performed in an optical 96-well plate on StepOne Plus Real-Time PCR System (Applied Biosystems, Thermo Fisher Scientific, United States) using SYBR R Green PCR Master Mix (TaKaRa, Dalian, China) according to the manufacturer's instructions. Each real-time PCR reaction (in a final volume of 20 µl) contained 10 µl of 2 × SYBR R Green Real-Time PCR Master Mix, 0.25 µl of each primer (10 µM, the forward and reverse primers), 1 µl of cDNA, 0.4 µl of ROX Reference dye (50×), and 8.1 µl of sterile distilled water. In parallel, for each sample, 1 µl of sterile distilled water, as the blank template, was included as the negative control. The cycling conditions included an initial denaturation and activation at 95 • C for 3 min, followed by 45 cycles at 95 • C for 10 s and 60 • C for 20 s. All amplifications were followed by a dissociation curve analysis of the amplified products by a dissociation step (95 • C for 15 s, 65 • C for 10 s, and 95 • C for 10 s) to confirm the amplicon specificity for each gene. Specific primers of the candidate genes were designed using the NCBI/Primer-BLAST 11 , with specific parameters set as PCR amplicon length of 100-200 bp, melting temperature (T m ) of approximately 60 • C, and primer pair specificity checking against the Refseq mRNA (database) of Schistosoma (taxid:6181) (Organism), and were commercially synthesized by Sangon Biotech (Shanghai, China) Co., Ltd. The gene expression levels were normalized with the 26S proteasome non-ATPase regulatory subunit 4 (PSMD4) (Liu et al., 2012), which was validated as a reliable reference gene in the transcriptomic analysis of S. japonicum elsewhere (Liu et al., 2012). The relative expression of the candidate genes was calculated using the 2 − Ct method (Livak and Schmittgen, 2001). All reactions were run in duplicate. All primers used are listed in Supplementary Table S5.
Statistical significance in the gene expression comparison was analyzed using the IBM SPSS Statistics 20 software. The relative data are expressed as the mean ± SD, and differences between groups were determined for statistical significance using Student's t-test. A P-value of ≤0.05 was considered to be of statistical significance.

Relative Expression of Some Differentially Expressed Genes of Important Biological Significance in Life Cycle Stages and Gonads
Another batch of S. japonicum-infected O. hupensis, also purchased from the Institute of Parasitic Disease Control and Prevention of Hunan Province, was used to release cercariae in order to infect twenty 6-to 7-week-old SPF-grade female BALB/c mice, which were purchased from Beijing Hua Fu Kang BioScience Co., Ltd. through WUCAE. The cercariae were collected to infect BALB/c mice, with 50 ± 2 cercariae per mouse, via percutaneous exposure according to the methods described above. Schistosome worms at 21, 28, 35, and 49 days post-infection were collected by hepato-portal perfusion of five mice each time, and the eggs at 35 and 49 days post-infection were harvested from the livers of mice by comminution and layered filtration with a nylon mesh. The worm gonads, the testes of male worms, and the ovaries and vitellaria of female worms were separated by cutting under a dissecting microscope (Supplementary Figure S5).
For total RNA isolation and examination of the expression of some DEGs of important biological significance through life cycle stages in definitive hosts and the gonads of worms, real-time PCR was performed according to the method described above.

dsRNA-Mediated Knockdown of SjPhb1 in vivo
In accordance with the standard protocols for in vivo RNAi in adult schistosomes (Li et al., 2018), specific SjPhb1 dsRNA of approximately 580 bp, which was designed with the online RNAi Design Tool 12 , was generated with DNA templates (Sjp_0046680) produced by PCR with specific primers that contained a T7 RNA polymerase promoter sequence at the 5 end (SjPhb1-RNAiU: 5 -GCGGATCCTAATACGACTCACTATAGGGACCTACAAACT GTGAATATAACG-3 ; SjPhb1-RNAiL: 5 -CCGCTCGAG TAATACGACTCACTATAGGGGCAGGAAGGTTAAGAAGTG -3 ). This was done by in vitro transcription using the T7 RiboMAX Express RNAi system (Promega, United States) according to the manufacturer's instructions. The irrelevant enhanced green fluorescent protein (EGFP) dsRNA of approximately 640 bp was also generated using pEGFP-N1 plasmids as DNA template amplification and used as the negative control . The concentration of dsRNA was measured in (Biofuture, United Kingdom). The size and integrity of the resulting dsRNA was confirmed by gel electrophoresis in 1.2% agarose MOPS, which was conducted to confirm that single RNA bands were of the correct size.
Fifteen infected BALB/c mice (50 ± 2 cercariae per mouse) were designed to be randomly allocated into three groups, with five mice in each group, as follows: SjPhb1 dsRNA injection group; EGFP dsRNA injection group set as the negative group; and 0.7% NaCl injection as the blank control group. SjPhb1 dsRNA or EGFP dsRNA (10 µg) dissolved in 300 µl of sterilized 0.7% NaCl solution was injected into the corresponding groups of infected mice immediately after infection (day 1), which were then given additional injections every 5 days until day 26 (days 6, 11, 16, 21, and 26) (Supplementary Figure S6). The blank control group of infected mice received injection with 300 µl of 0.7% NaCl at each time point mentioned above. But during the experiment process, one infected mouse assigned to the EGFP dsRNA injection group was injected with SjPhb1 dsRNA at the first dose by mistake, so that mouse was later assigned to the SjPhb1 dsRNA injection group and received SjPhb1 dsRNA injection; thus, the sample size is six mice for the SjPhb1 dsRNA injection group and four mice for the EGFP dsRNA injection group. In addition, one mouse in the 0.7% NaCl injection group suffered accidental death when immobilized for injection, so the sample size for this group is four mice.
On the 42nd day post-infection, the three groups of mice were sacrificed and worm samples were collected for recording worm count, percent of paired worms, and size measurements (worm length), and the efficiency of SjPhb1 dsRNA-mediated specific silencing of SjPhb1 was confirmed by qRT-PCR. The livers of mice were also collected for pathological evaluation of egg count estimation, granuloma, and fibrosis. Specifically, to examine whether SjPhb1 (RNAi) affected the overall growth of the worms, the parasites were immobilized on ice and imaged on a LY-WN-HPCCD, 5M Pixels High-Speed HPCCD-5 Color Microscope Camera connected to a VistaVision trinocular dissecting microscope at × 10 magnification. Parasite length was measured in digital images using the Image-Pro Plus 5.0 software (Media Cybernetics, United States). Quantitative analysis of parasite length was conducted using GraphPad Prism 6.0 software. For egg counting, 0.5 g of the left lateral lobe of the liver of each mouse was pulverized and digested in 10 ml of 5% KOH solution overnight at 37 • C. Eggs in 10 µl of digest in triplicate for every sample were counted under the microscope, and the number of eggs deposited in every gram of liver per couple of worms was finally calculated . To examine whether SjPhb1 knockdown in worms affected the pathogenicity of the worms, hematoxylin and eosin (H&E) and Masson staining were performed on the paraffin sections of mouse liver to evaluate the egg granuloma formation and fibrosis (Tao et al., 2018;Zhu et al., 2018).
The statistical significance in the RNAi effects comparison was analyzed using the IBM SPSS Statistics 20 software. The relative data are expressed as the mean ± SD, and differences between groups were determined for statistical significance using Student's t-test. A P-value of ≤0.05 was considered to be of statistical significance.

Transcriptome Profiles of Schistosomes From SCID and BALB/c Mice
The four total RNA samples, including the female worms recovered from BALB/c mice (B_F), male worms recovered from BALB/c mice (B_M), female worms recovered from SCID mice (S_F), and male worms recovered from SCID mice (S_M), have prominent 18S and 5S ribosomal bands on agarose gels (Figure 1A), while the 28S rRNA band is not present due to a known gap region within the molecule (in vivo nick of the L-rRNA of S. japonicum worms total RNA) (van Keulen et al., 1991;Chen and Qiu, 1993;Li and Shi, 1997;Lu et al., 2017). According to the Bioanalyzer 2100 analysis, the 28S rRNA peak is also absent in the samples, the 18S rRNA peak is present as the main peak, and the second peak is the weaker 5S peak ( Figure 1B). These indicated that all RNA samples possessed high integrity and purity and could be used for further experiments (Figure 1 and Supplementary Table S1).
A separate sequencing was performed on the double-stranded cDNA (dscDNA) libraries generated using the RNA transcripts of the 5-week-old male and female S. japonicum worms recovered from SCID and BALB/c mice. The Q30 was >90% for each sample ( Table 1). The RNA-seq data of the transcriptomes obtained based on the total RNA of the schistosomes with high quality revealed that a total of 103.6, 142.6, 90.3, and 132.3 million raw reads were obtained in B_F, B_M, S_F, and S_M samples, respectively. And 97.6, 134. 9, 84.5, and 124.6 million clear reads were obtained after quality filtration using TopHat2 (Kim et al., 2013) in B_F, B_M, S_F, and S_M samples, respectively (Table 1), in which reads matching transposon, mitochondrial, and ribosomal genes were filtered out. The mapping rates of clean reads to the reference S. japonicum genomic data (PRJEA34885) available at the WormBase ParaSite database 13 (Howe et al., 2016;Lee et al., 2018)  46.77% of which were uniquely mapped to the reference genome ( Table 1). Of the above total, 3.51, 5.52, 3.79, and 4.54% of clean reads belonged to putative alternatively spliced fragments in B_F, B_M, S_F, and S_M samples, respectively, and the remaining clean reads corresponded to unique transcript fragments with no evidence of alternative splicing (Table 1). Approximately 20.36,29.65,24.13,and 25.66% of the mapped clean reads in the above four samples were mapped to the annotated genes (protein-coding genes) ( Table 1), respectively.

Sample Correlation Analysis
Pearson's correlation analyses of the global expression profiles of schistosome samples revealed a stronger similarity between schistosome worms of the same sex from SCID and BALB/c mice (R 2 > 0.90) than that between worms of different sexes from the same host (R 2 < 0.80) (Figure 2A and Supplementary Figure S1). Similar results were obtained in the hierarchical clustering analysis of the global expression profiles of schistosome samples between the compared groups ( Figure 2B). As expected, the distance of correlation between worms of the same sex from different hosts was found to be shorter than that between the worms of opposite sex from the two kinds of hosts.

Differential Gene Expression Profiles of Schistosomes
In total, 11,183 gene transcripts were detected using FPKM > 0.01 as the threshold value by cross-referencing the S. japonicum gene coordinates in the WormBase ParaSite database with the coordinates of the clean reads that matched those S. japonicum genes ( Table 2 and Supplementary Table S2). Of these, 9,993 genes were expressed in female worms from SCID FIGURE 1 | Qualitative and quantitative integrity identification of the total RNA of the Schistosoma japonicum worm samples. (A) Qualitative integrity identification of the total RNA of the S. japonicum worm samples by denatured agarose gel electrophoresis. Of the total RNA from each sample, 3 µl was loaded on 1.2% denaturing agarose gel with ethidium bromide (EB). The voltage of the electrophoresis was 65 V. The gel was observed and photographed with an ultraviolet transmission analyzer after the electrophoresis was finished. The 28S rRNA bands are absent, the 18S rRNA bands are present as the main bands with appropriate peaks, and the second bands are the weaker 5S bands. Therefore, the in vivo nick of the L-rRNA of S. japonicum worms total RNA was observed and confirmed here, which indicated good integrity of the total RNA of the S. japonicum worm samples in this study. M: RNA ladder. 1: S_F; 2: S_M; 3: B_F; 4: B_M.
(B) Quantitative integrity identification of the total RNA of the S. japonicum worm samples by an Agilent 2100 Bioanalyzer. The peaks of 28S rRNA are also absent and the peaks of 18S rRNA are present as the main peaks for all four samples. Therefore, the in vivo nick of the L-rRNA in S. japonicum worms total RNA was also observed and confirmed here, which further indicated the good integrity and purity of the total RNA of the S. japonicum worm samples used here. B_F female worms from BALB/c mice, B_M male worms from BALB/c mice, S_F female worms from SCID mice, S_M male worms from SCID mice.  mice, 10,697 were expressed in male worms from SCID mice, 9, 760 were expressed in female worms from BALB/c mice, and 10,801 genes were expressed in male worms from BALB/c mice ( Table 2 and Supplementary Table S2). Moreover, 9,243 genes were commonly expressed in all four parasite forms ( Table 2 and  Supplementary Table S2).
Two independent pairwise comparisons between schistosomes of the same sex recovered from SCID and BALB/c mice were performed to identify their differentially expressed genes. In detail, significant differential expression between female worms from SCID mice and female worms from BALB/c mice (S_F vs. B_F), and between male worms from SCID mice and male worms from BALB/c mice (S_M vs. B_M), was computed with a false discovery rate (FDR) of 5% and a q value (P-adjusted value) < 0.05. In brief, 298 DEGs were obtained from a pairwise comparison of the schistosome samples (i.e., S_F vs. B_F and S_M vs. B_M) between SCID mice and BALB/c mice 5 weeks post-infection. The distributions of the up-and downregulated genes between the paired comparisons are displayed as scatter plots and heatmaps (Figures 3A-D). Specifically, 62 DEGs with 39 upregulated and 23 downregulated mRNAs were found in S_M vs. B_M ( Figure 3A and Supplementary Table S3), and 240 DEGs with 152 upregulated and 88 downregulated mRNAs were found in S_F vs. B_F (Figure 3B and Supplementary Table S4). It is reasonable that female worms from SCID and BALB/c mice had more differentially expressed genes than the male worms from the two hosts, which was also evidenced by the sample correlation analysis results wherein a weaker correlation was found between S_F and B_F than that between S_M and B_M (0.915 vs. 0.94). Of these DEGs between S_M vs. B_M, the top five DEGs downregulated to undetectable levels in S_M worms are Sjp_0125510 (cytochrome c oxidase subunit 2), Sjp_0130060 (cytochrome oxidase subunit I), Sjp_0117620 (26S proteasome non-ATPase regulatory subunit), Sjp_0005030 (SJCHGC08493 protein), and Sjp_0043150 (SJCHGC01077 protein After integrating the two subsets of differentially expressed genes, moreover, a Venn diagram revealed that four DEGs are present in both comparisons ( Figure 3E): they are Sjp_0117620 (26S proteasome non-ATPase regulatory subunit 10), Sjp_0127870 (mitochondrial proton/calcium exchanger protein; leucine zipper-EF-hand-containing transmembrane protein 1), Sjp_0130040 (dynein heavy chain, axonemal), and Sjp_0094540 (hypothetical protein containing proteinase inhibitor I25, cystatin domain) ( Table 3). Two DEGs, Sjp_0117620 and Sjp_0094540, showed the same change direction in male and female worms from SCID mice when compared with those in worms from BALB/c mice, with the former gene downregulated and the latter gene upregulated in worms from SCID mice ( Table 3). The other two DEGs, Sjp_0127870 and Sjp_0130040, showed opposite change directions in male and female worms from SCID mice, with both genes downregulated in male worms but upregulated in female worms from SCID mice ( Table 3).

Validation of RNA-seq Data by Quantitative RT-PCR
Using 26S proteasome non-ATPase regulatory subunit 4 (PSMD4) as a housekeeping gene (Liu et al., 2012), a subset of 35 DEGs from both S_F vs. B_F and S_M vs. B_M comparisons were assayed by qRT-PCR to validate the data of the RNAseq profiles (Supplementary Table S5 and Supplementary Figure S2). These selected genes were identified as DEGs in only S_F vs. B_F, in only S_M vs. B_M, or in the overlapping of both, and all of them were randomly selected among those that were annotated to ensure their biological relevance. In detail, we selected 19 DEGs just in S_F vs. B_F (Supplementary Figure S2A). Another two DEGs from both S_F vs. B_F and S_M vs. B_M were also selected for consideration. The RNAseq data for the selected DEGs correlated moderately with the qPCR data in S_F vs. B_F, with a concordance rate of 66.67% ( Figure 4C). We selected 14 DEGs just in S_M vs. B_M (Supplementary Figure S2B). The other two DEGs from both S_F vs. B_F and S_M vs. B_M, i.e., 26S proteasome subunit P28related (Sjp_0117620) and leucine zipper-EF-hand-containing transmembrane protein 1, mitochondrial (Sjp_0127870), were also selected for consideration. The RNA-seq data for the selected DEGs correlated moderately with the qPCR data in S_M vs. B_M, with a concordance rate of 56.25% (Figure 4D).

GO and KEGG Pathway Analyses of Differentially Expressed Genes
To identify significantly enriched gene categories and obtain a better understanding of the enriched functions of the differentially expressed genes, we performed GO term enrichment analyses (Ashburner et al., 2000). Significantly enriched GO categories (P < 0.05) probably identified gene groups related to parasite development and the morphological abnormalities of schistosomes from SCID mice.
Firstly, 113 differentially expressed genes between S_F vs. B_F were annotated with GO terms in three independent categories (Supplementary Tables S4, S7). Among the biological process ontology, the top three enriched GO     Tables S4, S7). Secondly, 37 differentially expressed gene sequences between S_M vs. B_M were annotated with GO terms in three independent categories (Supplementary Tables S3, S8). Among the biological process ontology, the top three enriched GO categories were oxidation-reduction process (GO:0055114), response to stress (GO:0006950), and response to DNA damage stimulus (GO:0006974) (Figure 5B, Supplementary  Figure S4, and Supplementary Table S8). In the case of the cellular components ontology, the top three enriched GO categories were extracellular region (GO:0005576), extracellular region part (GO:0044421), and nucleoplasmic THO complex (GO:0000446) (Figure 5B and Supplementary  Table S8). On the basis of the molecular function ontology, the top three enriched categories were metal ion binding (GO:0046872), cation binding (GO:0043169), and substratespecific transmembrane transporter activity (GO:0022891) ( Figure 5B and Supplementary Table S8). GO analysis also revealed an enrichment of DEGs between S_M vs. B_M that are mainly involved in important biological processes and cellular components, such as genes enriched in oxidation-reduction process, the decreased cytochrome c oxidase subunits 1 and 2 (cox1 and cox2) (Sjp_0130260 and Sjp_0130060), increased superoxide dismutase [Cu-Zn] (Sjp_0042990), decreased ATPase subunit 6 (Sjp_0122570), and genes enriched in response to stress or DNA damage stimulus, the increased UV excision repair protein RAD23 (Sjp_0118730) and histone H3.3 (Sjp_0121450) (Supplementary Tables S3, S8).
KEGG pathway and enrichment analysis indicated that the differentially expressed genes between S_F vs. B_F were enriched in purine metabolism, RNA polymerase, fatty acid biosynthesis, fatty acid degradation, pyrimidine metabolism, pyruvate metabolism, fatty acid metabolism, glycolysis/gluconeogenesis, and oxidative phosphorylation ( Figure 6A and Supplementary Table S9). KEGG pathway and enrichment analysis also revealed that the differentially expressed genes between S_M vs. B_M were enriched in oxidative phosphorylation, riboflavin metabolism, tyrosine metabolism, cysteine and methionine metabolism, propanoate metabolism, pyruvate metabolism, glycolysis/gluconeogenesis, biosynthesis of secondary metabolites, microbial metabolism in diverse environments, etc. (Figure 6B and Supplementary Table S10). Most of the enriched pathways had a P-value greater than 0.05 due largely to the limited number of DEG hits in the pathways.

Expression Profiles of DEGs With Important Biological Significance in Schistosome Life Cycle Stages and Gonads
In order to predict the potential important biological significance of the identified DEGs of schistosomes in their development during the life cycle stages in the definitive hosts, nine of the differentially expressed genes of important biological significance were selected for profiling their expression patterns through their life cycle in BALB/c mice (Figure 7) and in the gonads (Figure 8) by qRT-PCR. The worms, through their life cycle in mouse, were collected at 21, 28, 35, and 49 days post-infection (dpi), and the gonads (testes of the male worms and ovaries and vitelline glands of the female worms) were separated and collected from worms at 28, 35, and 49 dpi (Supplementary Figure S5). The results revealed that, firstly, mitochondrial prohibitin complex protein 1 (Phb1, Sjp_0046680), von Willebrand factor type A domain-containing protein (Sjp_0073890), dynein heavy chain (Sjp_0014800), and retrotransposon-like protein 1 (Sjp_0125250) were highly expressed in mature female worms at 35 and/or 49 dpi (Figure 8). Secondly, bone morphogenetic protein antagonist noggin (Sjp_0105360), glutamate receptor ionotropic (iGluR), NMDA 2B (Sjp_0082750), and transient receptor potential channel (Sjp_0091900) were highly expressed in male worms. Lastly, dead end protein homolog 1 (Dnd1, Sjp_0122230), von Willebrand factor type A domain-containing protein (Sjp_0073890), dynein heavy chain (Sjp_0014800), and transient receptor potential channel (Sjp_0091900) were highly expressed in eggs at 35 and/or 49 dpi. Additionally, Dnd1 (Sjp_0122230), noggin (Sjp_0105360), aspartyl protease protein (Sjp_0049310), dynein heavy chain (Sjp_0014800), and retrotransposon-like protein 1 (Sjp_0125250) were highly expressed in the testes of male worms and the vitellaria of female worms. von Willebrand factor type A domain-containing protein (Sjp_0073890) and transient receptor potential channel (Sjp_0091900) were specifically expressed in female gonads-ovaries and/or vitellaria. Glutamate receptor ionotropic (iGluR), NMDA 2B (Sjp_0082750) was highly expressed in the testes of male worms.

dsRNA-Mediated RNAi in vivo Revealed SjPhb1 Is Essential to Maintain the Growth and Development of Schistosomes
In order to investigate and confirm the biological function of the identified DEGs between S_F vs. B_F and S_M vs. B_M involved in the regulation of the growth and development of schistosomes, we silenced the expression of S. japonicum Phb1 gene (SjPhb1, Sjp_0046680), the function of which was predicted to be important in the development and reproduction of schistosomes, as above, using dsRNA-mediated RNAi in vivo through a mouse tail vein injection immediately after infection Li et al., 2018). Prohibitin proteins are multifunctional proteins located mainly at the inner membrane of the mitochondria and expressed in universal species. Prohibitin 1 (PHB1) and prohibitin 2 (PHB2) are the two highly homologous subunits of the eukaryotic mitochondrial PHB complex. They are interdependent on the protein level, and loss of one simultaneously leads to the loss of the other. The homologous gene of prohibitin in other animals was reported to play a vital role in the mitochondria's function, cell proteolysis, senescence, and apoptosis and is essential for embryonic viability and germline function-spermatogenesis and folliculogenesis (Artal-Sanz et al., 2003;He et al., 2008;Tyc et al., 2010;Peng et al., 2015;Chowdhury et al., 2016;Jin et al., 2016;Nguyen et al., 2016;Chai et al., 2017;Xu et al., 2017;Wang D. et al., 2017).
In this experiment, BALB/c mice infected with S. japonicum worms were injected with specific SjPhb1 dsRNA via tail vein immediately after infection, followed by an additional injection every 5 days for up to 26 dpi. The mice were finally sacrificed on day 42 after infection for effects assessment, which was in contrast to the negative (EGFP dsRNA treatment) and blank controls (0.7% NaCl treatment) (Supplementary Figure S5). Figure S7A and Supplementary Table S11) and pairing rate (Supplementary Figure S7B and Supplementary Table S11) were observed among the three groups, although a decreasing trend in worm burden could be detectable in the SjPhb1 dsRNA treatment group (Supplementary Figure S7A and Supplementary Table S11). To our excitement, dsRNA-mediated RNAi in vivo led to obvious morphologic changes in worms of the specific SjPhb1 dsRNA treatment group, wherein a significant decrease in worm length was detected in both female and male worms when compared with those of the two control groups (Figures 9a-g and Supplementary Table S11). qRT-PCR found that dsRNA treatment still led to a decrease in SjPhb1 transcript levels by approximately 12% in female worms (P > 0.05) and 29% in male worms (P < 0.05) on the 42nd day post-infection, after a final injection on the 26th day (Figure 9h), suggesting a sustained knockdown role of SjPhb1 dsRNA, although the dsRNA would inevitably decrease gradually since the last injection. Pathological examination of the egg granuloma and fibrosis in the livers of mice by H&E (Figures 10a-c) and Masson staining (Figures 10df) found a decreased size of egg granuloma and attenuated fibrosis formation in the livers of mice of the SjPhb1 dsRNA treatment group, which could be mainly attributed to the Transcriptome Analysis of S. japonicum From Mice decreased egg deposition in livers by 42.3% (Figure 10g and Supplementary Table S11).

DISCUSSION
In our previous publication , we have reported that mice lacking functional T and B cells (SCID mice) negatively influenced the growth and development of the schistosomes, which led to the downregulation of the granuloma formation in SCID mice when compared with that in BALB/c mice. Furthermore, in order to explore the molecular mechanisms involved in the regulation of growth and development of schistosomes, we also tested the metabolic profiles of the male and female S. japonicum worms recovered from SCID and BALB/c mice at 5 weeks post-infection using an untargeted liquid chromatography-tandem mass spectrometry (LC-MS/MS)-based high-resolution metabolomic investigation and screened out a series of perturbed metabolites and metabolic pathways potentially involved in the regulation of the growth and development of S. japonicum worms . In this study, synchronously, a comparative RNA-seq investigation was also performed to identify the distinct biosignatures in the transcriptomic profiles of male and female S. japonicum worms recovered from SCID and BALB/c mice, respectively.
In the results, Pearson's correlation analyses and the hierarchical clustering analysis of the schistosome samples found a relatively weaker correlation in S_F vs. B_F than that in S_M vs. B_M. This indicates that the female schistosome worms were affected more severely than did the male worms in SCID mice when compared with those in BALB/c mice. This was also verified by the subsequent finding that more differentially expressed genes were obtained in S_F vs. B_F than in S_M vs.
B_M. This could be expectable and reasonable as the growth and development of female worms were affected by male worms as well as the host's factors, i.e., the sexual maturation of female worms depends on the pairing with male worms. A similar finding was also obtained in our previous metabolomics research on the differential metabolites in S_F vs. B_F than in S_M vs. B_M . The overall results of the qRT-PCR verification of the RNAseq data showed that more than half (23/37, 62.16%) of the tested DEGs had consistent expression patterns between the RNA-seq and qRT-PCR data (Figure 4 and Supplementary Table S6), suggesting a moderate consistency rate between the qRT-PCR and RNA-seq data. The concordance rates between the RNAseq data and the qRT-PCR results were in the range 50-70%, which we consider to be a generally successful validation as it is already known that all transcriptome techniques, including microarray, RNA-seq, and qPCR, have their inherent pitfalls that affect accurate quantification and cannot be fully controlled (Marioni et al., 2008;Mortazavi et al., 2008;Malone and Oliver, 2011;Shi and He, 2014;Wang C. et al., 2014).
It was reported that DLC-1 of C. elegans and Plasmodium falciparum played an important role in growth and germ cell proliferation and gametogenesis in vivo (Daher et al., 2010;Ellenbecker et al., 2019). DLC-1 (Sjp_0057280), which is involved in microtubule-based process/microtubule-associated complex/microtubule cytoskeleton, was also identified as differentially expressed in female schistosomes from SCID mice when compared with those from BALB/c mice. This indicates a disturbed process of growth and germ cell proliferation and gametogenesis in female worms from SCID mice, which needs further research to test and confirm in the future. In addition, the enriched Rho protein signal transduction and regulation of Rho protein signal transduction (GO:0007266 and GO:0035023), with increased ephexin-1 (Sjp_0110200 and Sjp_0096350, which function as guanine nucleotide exchange factors for the important Rho-type GTPases), were also reported to be involved in a number of cell signaling pathways regulating actin cytoskeleton organization, gene transcription, cell cycle progression, and membrane trafficking (Hotchin and Hall, 1996;Tapon and Hall, 1997;Hall, 1998;Aspenstrom, 1999;Hall and Nobes, 2000;Santos et al., 2002;Vermeire et al., 2003;Quack et al., 2009;Sit and Manser, 2011;Singh et al., 2017). 5-HT affects the feeding behavior and obesity in the central nervous system. On the other hand, peripheral 5-HT may also play an important role in regulating glucose and lipid metabolism (Watanabe et al., 2014(Watanabe et al., , 2016. A decrease in serotonin receptor 5HTR (Sjp_0082160) was detected in female worms from SCID mice, which suggests a weakened glucose and lipid metabolism in these worms. The protein HSP70 is an important part of the cell's machinery for protein folding and helps protect cells from stress (Neumann et al., 1993;Maresca and Kobayashi, 1994;Jayasena et al., 1999;Brochu et al., 2004;Hatherley et al., 2014;Charnaud et al., 2017). The homolog HSP70 (Sjp_0026380) transcript was detected as increased in female worms from SCID mice, which indicates an unfavorable environment of SCID for schistosomes. Wnt genes, encoding signaling molecules, play important roles in regulating patterning and morphogenesis during vertebrate gastrulation (Rauch et al., 1997;Kilian et al., 2003;McIntyre et al., 2013;Martyn et al., 2018). An increase of Wnt5 in female worms (Sjp_0003340) from SCID mice indicates an increased need of Wnt5 in the development of female worms from SCID mice, as SjWnt5 was predicted to play a role in the regulation of parasite muscle development and the development of the reproductive organs of both sexes based on its prominent expression in the subtegumental musculature and acetabulum musculature of schistosomula and adult worms, the testes of the male and the ovary as well as the vitellaria of the female (Ta et al., 2015). The trematode eggshell synthesis protein (Sjp_0070860) is important to the eggshell structure (Bobek et al., 1988;Ebersberger et al., 2005), the disturbance of which would result in influenced egg formation and release by female worms.
In the identified DEGs in S_M vs. B_M, meanwhile, cytochrome c oxidase, which is a large transformation protein complex found in the mitochondria and catalyzes the electron transport from cytochrome c to oxygen coupled to proton translocation, was reported to have important roles in energy metabolism of cell respiration and stress-induced apoptosis, and cytochrome c oxidase deficiency usually leads to muscle weakness and movement problems (Kadenbach et al., 2004;Pang et al., 2011;Douiev and Saada, 2018;Ramzan et al., 2019). Therefore, the decrease in cytochrome c oxidase in male worms from SCID mice probably results in attenuated energy metabolism and supply and increased apoptosis and finally manifests as developmental retardation or defects (Baden et al., 2007;Hatakeyama and Goto, 2017). In addition, the increase in UV excision repair protein RAD23, which is a nucleotide excision repair protein and plays an important role in nucleotide excision DNA repair through the ubiquitin-proteasome pathway when the genetic material is damaged by UV irradiation or other genotoxic stresses (Mueller and Smerdon, 1996;Schauber et al., 1998;Xie et al., 2004;Gong et al., 2006;Dantuma et al., 2009;Wade et al., 2009;Zhou et al., 2015;Lahari et al., 2017;Verma et al., 2019), also indicates potential genotoxic stresses (or unfavorable living environments) exerted by the hemato-microenvironment of SCID mice on the worms.
In summary, therefore, these differential or disturbed gene expressions associated with structural components and biological/molecular function could be inferred as the potential causes of growth and development retardation of female schistosomes from SCID mice (Anderson et al., 2015;Han et al., 2015a,b). Interestingly, furthermore, the majority of DEGs between S_F vs. B_F were categorized into the cellular components ontology, while most of the DEGs between S_M vs. B_M were categorized into the molecular function ontology (Figure 5). This phenomenon suggests that more influence or effects from the host on the female worms mainly focus on the structural composition, while more effects from the host on the male worms mainly focus on the molecular function. This could be addressed as the growth and development of the female worms being dependent on interactions with the male worms as well as their hosts (Shaw et al., 1977;Den Hollander and Erasmus, 1985;Gupta and Basch, 1987;Haseeb et al., 1989;Siegel and Tracy, 1989;Haseeb and Eveland, 1991;Kunz et al., 1995;Kunz, 2001;Hernandez et al., 2004;Han et al., 2009;deWalick et al., 2012).
In order to test and verify the biological function of the identified DEGs in the regulation of the growth and development of schistosomes, we silenced the expression of SjPhb1 (Sjp_0046680) in vivo through dsRNA-mediated RNAi and confirmed that Sjphb1 is essential in governing the parasite growth and development, although further mechanical research remains needed.
To conclude, for the first time, our current research pointed to the relevant biological processes enriched in the S. japonicum female and male worms from SCID mice based on their DEGs compared with those from BALB/c mice, providing a set of genes involved in the biological processes associated with the morphological abnormalities of worms in SCID mice. This research also provided a list of potential antischistosomal targets for drug or vaccine development.

DATA AVAILABILITY STATEMENT
The raw data and normalized gene-level data from the RNAseq discussed in this article have been completely deposited in NCBI's Gene Expression Omnibus and are accessible through GEO Series accession number GSE122317 (https:// www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE122317). All relevant data are within the paper and its Supporting Information files.

ETHICS STATEMENT
The animal study was reviewed and approved by the Institutional Animal Care and Use Committee (IACUC) at Wuhan University Center for Animal Experiments (WUCAE).

AUTHOR CONTRIBUTIONS
RL conceived and designed the experiments. RL, W-JC, FY, Y-DZ, Q-PZ, and H-BT performed the experiments. RL, H-FD, and HJ contributed reagents, materials, and analysis tools. RL analyzed the data and wrote the manuscript. RL and HJ critically revised the manuscript. All authors read and approved the final version of the manuscript.

FUNDING
This work was supported by the Special Project for Schistosomiasis Control of Health and Family Planning Commission of Hubei Province (grant no. WJ2017X002), the Natural Science Foundation of Hubei Province (grant no. 2017CFB239), and the "Fundamental Research Funds for the Central Universities" of China (grant no. 2042017kf0033) to RL. The funders had no role in the study design, data collection and analysis, decision to publish, or preparation of the manuscript.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmicb.2020. 00274/full#supplementary-material FIGURE S1 | Sample relationships revealed by transcriptome profiles between worm samples. (A-F) Scatter plots of Pearson's correlation analysis based on the transcriptome profiles between worm samples. The figures in the matrix R 2 are the square of the correlation coefficient (r) between two samples. FIGURE S2 | Validation of a subset of DEGs by qRT-PCR. Thirty-five differentially expressed genes in female or male worms between SCID mice and BALB/c mice identified by RNA-Seq were selected for validation by qRT-PCR analysis. (A) 21 DEGs between S_F vs. B_F were selected for validation. (B) 16 DEGs between S_M vs. B_M were selected for validation. Bars represent standard deviation of the mean for the duplicate biological replicates. Sjp_0117620 and Sjp_0127870 were differentially expressed both S_F vs. B_F and S_M vs. B_M identified by RNA-Seq. The PSMD4 (FN320595) was used for internal normalization. The Student's t-test was used to calculate the statistical significance of the expression differences among the compared parasite groups (difference with a P-value < 0.05 was considered to be statistically significant).  FIGURE S6 | Schematic overview of the experimental design of RNAi in vivo. All biological replicates of BALB/c mice infected with S. japonicum cercaria (n = 5 for each group -SjPhb1 dsRNA injection group, EGFP dsRNA injection group, and 0.7% NaCl injection group.) were maintained in the same conditions, except that they were received SjPhb1 dsRNA injection, EGFP dsRNA injection, and 0.7% NaCl injection through tail vein immediately after infection then given an additional injection every 5 days for up to 26 days, respectively. TABLE S1 | Purity and quantity of total RNA of the S. japonicum worms sample.        TABLE S11 | (a) dsRNA-mediated RNAi effects assessment (worm burden). (b) dsRNA-mediated RNAi effects assessment (pairing rate). (c) dsRNA-mediated RNAi effects assessment (egg count). (d) dsRNA-mediated RNAi effects assessment (egg count).