Abstract

Background

   The freshwater snail Oncomelania hupensis is the obligate intermediate
   host for Schistosoma japonicum in China. Transcriptomic examination of
   snail–schistosome interactions can provide valuable information of host
   response at physiological and immune levels.

Methods

   To investigate S. japonicum-induced changes in O. hupensis gene
   expression, we utilized high-throughput sequencing to identify
   transcripts that were differentially expressed between infected snails
   and their uninfected controls at two key time-point, Day 7 and Day 30
   after challenge. Time-series transcriptomic profiles were analyzed
   using R package DESeq 2, followed by GO, KEGG and (weighted gene
   correlation network analysis) WGCNA analysis to elucidate and identify
   important molecular mechanism, and subsequently understand
   host–parasite relationship. The identified unigenes was verified by
   bioinformatics and real-time PCR. Possible adaptation molecular
   mechanisms of O. hupensis to S. japonicum challenge were proposed.

Results

   Transcriptomic analyses of O. hupensis by S. japonicum invasion yielded
   billion reads including 92,144 annotated transcripts. Over 5000
   differentially expressed genes (DEGs) were identified by pairwise
   comparisons of infected libraries from two time points to uninfected
   libraries in O. hupensis. In total, 6032 gene ontology terms and 149
   KEGG pathways were enriched. After the snails were infected with S.
   japonicum on Day 7 and Day 30, DEGs were shown to be involved in many
   key processes associated with biological regulation and innate immunity
   pathways. Gene expression patterns differed after exposure to S.
   japonicum. Using WGCNA, 16 modules were identified. Module-trait
   analysis identified that a module involved in RNA binding, ribosome,
   translation, mRNA processing, and structural constituent of ribosome
   were strongly associated with S. japonicum invasion. Many of the genes
   from enriched KEGG pathways were involved in lysosome, spliceosome and
   ribosome, indicating that S. japonicum invasion may activate the
   regulation of ribosomes and immune response to infection in O.
   hupensis.

Conclusions

   Our analysis provided a temporally dynamic gene expression pattern of
   O. hupensis by S. japonicum invasion. The identification of gene
   candidates serves as a foundation for future investigations of S.
   japonicum infection. Additionally, major DEGs expression patterns and
   putative key regulatory pathways would provide useful information to
   construct gene regulatory networks between host-parasite crosstalk.

   Keywords: Oncomelania hupensis, Schistosoma japonicum, Invasion,
   Transcriptome

Introduction

   Schistosomiasis is a zoonotic parasitic disease estimated to affect
   over 250 million people worldwide [[39]1]. It is transmitted by snails
   infected with worms of the trematode genus Schistosoma. Oncomelania
   hupensis (Gastropoda: Rissooidea), mainly distributed in China’s
   Yangtze river basin and its south, is the only intermediate host of
   Schistosoma japonicum, which plays a key role in the transmission of
   schistosomiasis in China [[40]2, [41]3]. The number of Schistosomiasis
   patients and the burden of the disease have dramatically reduced since
   the establishment of People’s Republic of China (P.R.C), together with
   the national schistosomiasis control strategy shift from the morbidity
   control strategy to an integrated strategy in 2004 [[42]4–[43]6].
   However, the disease remains endemic in some regions and contributes to
   a major public health concern in China. Moreover, there are still O.
   hupensis snails distributed in the ecologically complicated
   environments, such as lakeshore, marshland and mountainous areas. The
   breeding and spread of O. hupensis is the hotbed for the transmission
   and retransmission of Schistosomiasis, even after elimination.

   The life cycle of schistosomiasis is complicated. Its life cycle
   necessitates the presence of an intermediate host-a fresh water mollusk
   [[44]7]. After released from the egg, the miracidium, the first larval
   stage in the life cycle, begin to search for its specific intermediate
   host. It penetrates the snail, subsequently multiply asexually into
   multicellular sporocysts and later into cercarial larvae [[45]8,
   [46]9]. The S. japonicum miracidium develop into sporocysts larva in
   the O. hupensis within a week. The glands, ganglia, eye spots and
   apical papilla get degenerated in this stage, except the ciliated
   epidermis replaced by neodermis (tegument). After about 30 days
   infected snails release free swimming cercariae in response to
   sunlight. These can penetrate the skin of the mammalian host within
   12–24 h and finally invade the definite host [[47]10, [48]11]. In
   addition to the parasite itself, it is well-known that the infection in
   O. hupensis depends on many factors, such as parasite-snail
   compatibility, environmental aspects, individual snail defense capacity
   and specificity, and so on. Thus, the schistosome-snail interaction is
   a key area of research of biomedical significance.

   With the help of tremendous breakthroughs in molecular research, the
   application of next-generation sequencing (NGS) technology in genomics,
   transcriptomics, epigenomics and metagenomics have been performed to
   drive forward our understanding of many processes and potential
   regulatory mechanisms in schistosomiasis involved both the parasite and
   host. The genome of Biomphalaria glabrata published in May 2017
   [[49]12], which has dramatically improved the efficiency of gene
   discovery in snail. In a recent study Wang et al. tried to determine
   modification of neuropeptides in the snail B. glabrata ganglia nervous
   system before and after infection with Schistosoma mansoni miracidia
   [[50]13]. They found that the neuropeptides and precursor proteins,
   particularly those involved in snail reproduction were heavily down
   regulated and less abundant in prepatent snails compared to
   non-infected snails. For O. hupensis, transcriptomic studies by RNA-seq
   have been performed after molluscicide treatment. Many genes involved
   in key processes associated with biological regulation and innate
   immunity have been identified, which would benefit for elucidating the
   molluscicidal mechanism [[51]14]. Another study by comparing the
   differences in gene expression between the hilly and marshland snails
   revealed that there is a potential relationship between expression
   profiling and ecological adaptation of the snail that may have
   implications for schistosomiasis control in future [[52]15]. Our team
   reported preliminary analysis of raw data of O. hupensis before and
   after S. japonicum invasion, and 178,436,865 clean reads with a Q20
   percentage of 87.90% were obtained [[53]16]. However, we failed to
   conduct in-depth analysis due to sequence contamination of S.
   japonicum. Nevertheless, the key molecular mechanisms involved in the
   parasite–snail interaction, the subsequent larva development, and the
   main molecular events occurring during infection establishment in the
   intermediate host have not been investigated so far.

   Here, we provide and further analyze global view of the transcriptomes
   at two key time point to determine which physiological changes and
   immune responses occur in the snail after infections with S. japonicum
   miracidia. Understanding the pathogen–host interactions and critical
   molecular events on these time nodes during infection establishment in
   prepatent snails will provide new insights for the development of novel
   markers of detecting infected snails and anti-schistosome strategies.

Methods

Sample collection and preparation

   Oncomelania hupensis specimens were originally collected from Guichi
   district, Chizhou city, Anhui province during April 2008 to February
   2012 (Fig. [54]1a). The collected samples were reared in the insectary
   for a week, and confirmed by cercaria shedding test. After excluding a
   natural infection, all snails were used for subsequent investigations.
   The freshly hatched miracidia were harvested from the liver of New
   Zealand white rabbits as describe in [[55]17], then pooled and placed
   in a container with de-chlorinated water. Subsequently, snails were
   separately placed in a 24-well culture plate of de-chlorinated water,
   and then added 10 miracidia to each well, resulting in a snail to
   miracidia ratio of ~1:10. Snails were illuminated under light for 6 –8
   h at 25 − 30 °C conditions. Snails were then removed and placed into
   20 × 30 cm trays and cultured. De-chlorinated water was added daily to
   the culture trays to keep the moisture at a relative humidity of 85%.
   Throughout the experiment, snails were checked daily and dead ones were
   removed with forceps and counted. Finally, 10 snails were collected at
   days 7 and 30 post-infection. Infection was assessed by compression of
   tissues between two glass plates after removing of shell and examined
   under microscope to confirm the presence of cercariae or sporocysts
   (Fig. [56]1b). The snail soft tissues were pooled (at least three
   replicates), and then immediately flash frozen and stored for further
   analysis.

Fig. 1.

   [57]Fig. 1
   [58]Open in a new tab

   a Geographical location of the sampling site in Guichi district in
   Anhui province; b overall workflow for collection of non-infected and
   infected O. hupensis. Control group is Group N(naïve stage), while
   infection group include Group A (initial infection stage, 7 dpi) and
   Group C (late infection stage, 30 dpi)

Transcriptome preprocessing and submission

   Briefly, total RNAs were extracted separately from three different
   groups, designated as N (N1–N3, naïve stage), A (A1–A3, initial
   infection stage, 7 dpi) and C (C1–C3, late infection stage, 30 dpi),
   respectively, using TRIzol reagent (Invitrogen, USA) according to the
   manufacturer’s protocol. Genomic DNA was removed by treating RNA
   samples with DNase I (Fermentas) prior to cDNA synthesis. RNA
   quantification and quality were measured by the Nanodrop ND-1000
   spectrophotometer (Nanodrop Technologies, Wilmington, DE), combining
   with Agilent 2100 Bioanalyzer and denaturing gel electrophoresis.

   Total RNA (3 µg) from each group was prepared according to the Illumina
   protocol as described previously. In brief, Magnetic beads with
   oligo-dT were used to combine the poly-A of the mRNA for purifying the
   mRNA from the total RNA. Subsequently, poly (A) mRNA was chopped into
   short fragments using divalent cations under elevated temperature. The
   fragments were used to synthesize first-strand cDNA with random
   primers, and first-strand cDNA was transformed into double-strand cDNA
   by using RNase H and DNA polymerase I according to the manufacturer’s
   instructions. Short fragments of desirable lengths (200–300 bp) were
   purified using the QIAquick PCR Extraction Kit (Qiagen, Valencia, CA,
   USA), which was also used for continued end repair and ‘A’ base
   addition. The end-repaired DNA fragments were ligated with sequencing
   adapters through A and T complementary base pairing. AMPure XP beads
   (Beckman Coulter, Shanghai, China) were used to remove unsuitable
   fragments. The multiplexed cDNA libraries were checked using PicoGreen
   (Quantifluor™-STfluorometerE6090, Promega, CA, USA) and
   fluorospectrophotometry (Quant-iT59PicoGreen dsDNA Assay Kit;
   Invitrogen, P7589) and quantified with Agilent 2100 (Agilent 2100
   Bioanalyzer, Agilent, 2100; Agilent High Sensitivity DNA Kit, Agilent,
   5067–4626). The sequencing library was gradually diluted and quantified
   to 4–5 pM and sequenced using the Illumina NextSeq™ 500 platform at
   Shanghai Personal Biotechnology Co., Ltd. The transcriptome datasets
   are available from the NCBI Sequence Read Archive (SRA) under accession
   number SRR9598637–SRR9598645.

De novo transcriptome assembly and functional annotation

   Clean reads were obtained from raw data by removing adaptor sequences,
   low quality reads (Q < 20) and reads containing adapter or poly-N with
   the help of cutadapt [[59]18] and fqtrim [[60]19]. De novo assembly of
   clean data was assembled based on clean reads using Trinity [[61]20] to
   generate transcripts. For further analysis, the assembled sequences
   were searched against the NCBI non-redundant nucleotide database (Nt),
   non-redundant protein database (Nr) and Swiss-Prot database with an
   E-value < 10^−5. For contigs annotated as the same references, we