De Novo Transcriptome Analysis of Oncomelania hupensis after Molluscicide Treatment by Next-Generation Sequencing: Implications for Biology and Future Snail Interventions

The freshwater snail Oncomelania hupensis is the only intermediate host of Schistosoma japonicum, which causes schistosomiasis. This disease is endemic in the Far East, especially in mainland China. Because niclosamide is the only molluscicide recommended by the World Health Organization, 50% wettable powder of niclosamide ethanolamine salt (WPN), the only chemical molluscicide available in China, has been widely used as the main snail control method for over two decades. Recently, a novel molluscicide derived from niclosamide, the salt of quinoid-2',5-dichloro-4'-nitro-salicylanilide (Liu Dai Shui Yang An, LDS), has been developed and proven to have the same molluscicidal effect as WPN, with lower cost and significantly lower toxicity to fish than WPN. The mechanism by which these molluscicides cause snail death is not known. Here, we report the next-generation transcriptome sequencing of O. hupensis; 145,008,667 clean reads were generated and assembled into 254,286 unigenes. Using GO and KEGG databases, 14,860 unigenes were assigned GO annotations and 4,686 unigenes were mapped to 250 KEGG pathways. Many sequences involved in key processes associated with biological regulation and innate immunity have been identified. After the snails were exposed to LDS and WPN, 254 unigenes showed significant differential expression. These genes were shown to be involved in cell structure defects and the inhibition of neurohumoral transmission and energy metabolism, which may cause snail death. Gene expression patterns differed after exposure to LDS and WPN, and these differences must be elucidated by the identification and annotation of these unknown unigenes. We believe that this first large-scale transcriptome dataset for O. hupensis will provide an opportunity for the in-depth analysis of this biomedically important freshwater snail at the molecular level and accelerate studies of the O. hupensis genome. The data elucidating the molluscicidal mechanism will be of great benefit in future snail control efforts.


Introduction
The freshwater snail Oncomelania hupensis is the only intermediate host of Schistosoma japonicum, the agent of the most virulent form of schistosomiasis [1,2], which is endemic mainly in China and, to a lesser extent, in Indonesia and the Philippines. In China, continuous control efforts have achieved a remarkable reduction in the prevalence of the infection and the burden of disease among humans, making China one of the most successful countries in terms of schistosomiasis control in the world [3][4][5]; however, the intermediate host O. hupensis is still present in numerous areas covering a total of 3.69 billion m 2 [6], mainly in the lake and marshland regions along the Yangtze River. In addition, as the ecology of central and southern China are changing drastically as a result of the construction of the Three Gorges Super Dam, the South-to-North water-transfer project and climate change, new habitats have arisen, and the distribution of O. hupensis, and thus the distribution and transmission of schistosomiasis, have changed [7][8][9][10]. As a consequence, the infection risk to humans and livestock is severe and unpredictable. As schistosomiasis is a neglected tropical disease mainly governed by the geographical distribution of O. hupensis, snail control, as a major component of the schistosomiasis control program, is still key in attempts to interrupt the transmission of the disease [11,12] and may be an important component of an integrated strategy for the elimination of schistosomiasis [7,13,14].
Although many approaches, including environmental modification, molluscicide treatment, biological control and social intervention, have been used to control the intermediate host snails in China, the best, most cost-effective approach is still not clear [7,14]. Molluscicide treatment remains the most widely used method due to its wide coverage, simple procedure and fast action. As niclosamide is the only molluscicide recommended by the World Health Organization [15][16][17] due to its high efficiency, low toxicity and comparatively low environmental contamination, 50% wettable powder of niclosamide ethanolamine salt (WPN) has been extensively used as the only chemical molluscicide in fields in China with endemic schistosomiasis for more than 2 decades. Although WPN has been repeatedly used for long periods, the molluscicidal mechanism is still not well known. Some histochemical studies have shown that WPN immersion can decrease the activity of enzymes in snails, such as phenol oxidase [18], cytochrome C oxidase (CCO), lactate dehydrogenase (LDH), succinate dehydrogenase (SDH), cholinesterase (CHE) and nitric oxide synthase (NOS) [19], leading to hepatopancreas and cerebral ganglion cell damage [20]. These findings suggest that WPN might kill snails by blocking nervous system transmission and metabolic pathways [21]. All of these studies are limited to several enzymatic activities and organ-structure changes in snails after molluscicide treatment. The molecular basis of molluscicidal action is not yet known, hindering identification of valuable target genes and pathways. Park et al. found that niclosamide can cause mitochondrial fragmentation and may contribute to the apoptotic and autophagic cell death of HeLa cells [22], which has yet not been confirmed in snails.
In addition, studies have shown that WPN application has physical complications that include rapid precipitation and difficult dispersion, as well as complications such as toxicity to aquatic animals and expense. To overcome the precipitation and dispersal problems without reducing molluscicidal activity, larger amounts of WPN must be used in the field, leading to increased cost and environmental toxicity [23,24]. Dai et al. reported a new molluscicidal formulation of niclosamide, a 25% suspension concentrate of niclosamide mixture (SCN) containing niclosamide, alkyl polyglycosides, glycerol sodium carboxymethyl cellulose and sodium benzoate [25]. SCN has improved molluscicidal efficiency, with better dispersion, lower cost and less acute toxicity to fish than WPN; however, its biodegradability [26], effect on chromosomal aberrations and genotoxicity [27] remain to be evaluated. Because niclosamideis currently the only chemical available for snail control, the possibility of the emergence of resistance to this compound in the intermediate host snails has received considerable attention. Although resistance has not been demonstrated [12,28], the response of O. hupensis snails to niclosamide has exhibited regional variation in China [12,28] which make the resistance development possible. In our previous study, a new molluscicide derived from niclosamide, the salt quinoid-2', 5dichloro-4'-nitro-salicylanilide (Liu Dai Shui Yang An, LDS), was developed and shown to have the same molluscicidal effects as WPN, with lower cost, significantly lower toxicity to fish and a simpler production process [29]. LDS was then recommended as a new choice of molluscicide for use in a variety of habitats. Further studies are required on the molluscicidal mechanism, toxicity and differential effects on protein expression of LDS and WPN to maximize the efficacy of molluscicidal treatment and reduce environmental pollution and the possibility of the emergence of niclosamide resistance.
The objective of this study was to present the first known large-scale transcriptomic dataset generated from O. hupensis using Illumina technology, with emphasis on the gene-expression profiles of O. hupensis after treatment with the molluscicides LDS and WPN. The assembled and annotated sequences and gene-expression profiles will allow us to identify most of the genes expressed in O. hupensis without prior genome information and discover valuable resources for the analysis of the genes and metabolic pathways involved in the molluscicidal mechanism.

Snail Collection and Molluscicide Assay
O. hupensis specimens were collected in September 2012 in a location (30°08 0 N, 112°22 0 E) in Jiangling County, Hubei Province, China. This zone is a marshland beside the Yangtze River and is not a national park or other protected area or private land. Snail control in this zone is supervised by the Institute of Schistosomiasis Control of Jinzhou County, Hubei Province, China. Our collection was permitted by and carried out with the assistance of the institute. No other specific permission was required for us to collect and use these snails. O. hupensis is not on the list of endangered or protected species. A total of 3,000 snails were brought back to the laboratory. To ensure compatibility, after they were rinsed in dechlorinated tap water and kept for 24 h in the lab, snails with higher vitality ranging from 6 to 8 mm in shell length were selected and then maintained on wet filter paper in a Petri dish. The Petri dish was covered with plastic mesh and placed in an incubator at 25°C that was supplied with fresh humidified air, and the filter paper was changed every three days.
The molluscicidal effects of LDS and WPN (LDS and WPN were obtained from Hubei Provincial Center for Diseases Control and Prevention) were assessed by immersion of the snails in these substances at room temperature (25 ± 1°C). In previous molluscicide experiments, LDS and WPN with active concentrations from 0.05 mg/L to 1 mg/L have been used to estimate molluscicidal activity. More than 50% of the snails died after 12 h treatment, and 100% of the snails died after 24 h treatment at a concentration of up to 0.2 mg/L. To ensure that enough samples could be collected for next-generation sequencing under effective molluscicide treatment, active concentrations of 0.1 mg/L LDS and WPN were used. Twenty active, mature O. hupensis were exposed to 80 mL in individual 100 mL flasks. In the control group, the same number of snails was kept in dechlorinated tap water under the same conditions. After immersion in LDS, WPN and H 2 O for 6 h, 12 h, 24 h and 48 h, separately, the snails were washed with dechlorinated water, mortality was assessed, and the snails were then crushed for the determination of the infections status [30]. All treatments were performed in triplicate and survival analysis (Kaplan-Meier/Log Rank (Mantel-Cox)) in IBM SPSS Statistics 20 was used to estimate the survival time of snails in immersion test.
The living snails from the molluscicide assay were further washed in diethylpyrocarbonate (DEPC)-treated water. The soft bodies of the snails were dissected out individually on glass slides and examined microscopically to verify infection. The snails that were confirmed to have no helminthic infection were used in this study. The whole soft tissue of each individual was carefully washed 3 times in ice-cold DEPC-treated water containing 0.3% NaCl and then pooled into each group before storage in liquid nitrogen until further use.

RNA Extraction, mRNA-seq Library Construction and Illumina Sequencing
Total RNA was extracted from 20 living snails for each group separately using TRIzol Reagent (Ambion, Life Technologies) according to the manufacturer's instructions. Any trace genomic DNA was removed by treating RNA samples with DNase I (Fermentas) prior to cDNA synthesis. The integrity and size distribution of RNA were then verified using an Agilent 2100 Bioanalyzer. RNA samples with RNA integrity numbers (RIN) 7.5 were used for cDNA library preparation. The concentration of RNA in each sample was determined using a NanoDrop 2000 spectrophotometer.
The cDNA library was prepared for each group according to the Illumina protocol. In brief, poly (A) mRNA was purified from the total RNA using oligo (dT) magnetic beads and chopped into short fragments using divalent cations under elevated temperature. The cleaved RNA fragments were reverse transcribed into first-strand cDNA using random primers and then synthesized into double-stranded cDNA. From the cDNA, a paired-end library was synthesized for each group using the mRNA-Seq Sample Preparation Kit (Illumina) according to the manufacturer's instructions. Short fragments were purified with the QIA quick PCR Purification Kit (Qiagen), which was also used for continued end repair and 'A' base addition. These DNA fragments were then ligated into adapters and purified through gel separation. Finally, the adaptorligated libraries were amplified by PCR for sequencing. Illumina sequencing was performed at Beijing Berry Genomics Co., Ltd., using the HiSeq 2000 platform. The transcriptome datasets are available from the NCBI Sequence Read Archive (SRA) under accession number SRP041729.

Sequence Analysis and Functional Annotation
Raw data generated by Illumina sequencing were transformed by base calling and preprocessed to remove adaptor fragments and low-quality reads to yield clean reads for analysis. These short clean reads were assembled into contigs, transcripts and unigenes using the programs Velvet and Oases. RPKM was calculated by RNA-Seq by Expectation-Maximization (RSEM) program and used to normalize the abundances of the unigenes. The P and Q value (< 0.05) of expression fold change was calculated in Empirical analysis of digital gene expression data in R (edgeR) and used to identify the unigenes that were differentially expressed between molluscicide-treated groups.
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 of < 10 −5 . Gene names were assigned to each protein sequence based on the best BLAST hit. To obtain the Gene Ontology (GO) annotation of unigenes, the Blast2GO program was used. WEGO software was used to perform GO functional classifications for all unigenes and to explore the macro-distribution of gene functions. The metabolic pathways of the unigenes were predicted by Kyoto Encyclopedia of Genes and Genomes (KEGG) mapping.

Validation of Differential Gene Expression
Eleven unigenes that were differentially expressed and two unigenes that were not differentially expressed between the H 2 O-treated group and the LDS-treated group and had significant BLAST hits with the homologs were chosen for validation using real-time PCR with gene-specific primers. Before real-time PCR, semi-quantitative PCR was carried out to confirm the presence or absence of fragments in the original extracted fractions. Total RNA samples were the same as the samples prepared for cDNA library construction. The first-strand cDNAs were random-hexamer primed and synthesized using the RevertAid First Strand cDNA Synthesis Kit (Thermo). Oligonucleotide primers were synthesized based on the sequences determined for differentially expressed candidate genes, and 0.4 μM of each primer was used in the following PCR protocol: 94°C for 3 min followed by 39 cycles of 94°C for 15 s, optimal annealing Tm for 15 s, and 72°C for 20 s. The primers used to amplify the O. hupensis 18S (F: 5'-CGTCC TTTTGGTGACTCTGG-3'; R: 5'-TGGATGTGGTAGCCGTTTCTC-3') were designed based on the O. hupensis 18S ribosomal RNA gene partial sequence (AF367667). O. hupensis 18S was also used in parallel amplification to monitor the transcription of this internal control for normalization. Real-time PCR for each transcript of interest (for primers, see S1 Table) was performed in triplicate in a 20 μL reaction for each reaction, and each reaction was repeated three times, using SYBR Green Real-time PCR master mix (TOYOBO) and a Bio-Rad CFX96 thermal cycler following the manufacturers' protocols. Each of the PCR products was verified by gel electrophoresis as having a clear single band. The mRNA levels were normalized to O. hupensis 18S, the fold changes were calculated by comparing the mRNA of LDS-and WPN-treated snails to that of H 2 O-treated snails using the comparative delta-delta-Ct method [31], and the results were statistically analyzed using One-way ANOVA and then Post Hoc Multiple Comparisons (Dunnett) in SPSS 13.0.

Molluscicide Assay
After the O. hupensis were immersed in LDS or WPN at an active concentration of 0.1 mg/L or H 2 O for 48 h, snail mortality increased with exposure time both in the LDS-and WPN-treated groups ( Table 1, S1 Fig.). The survival time of snails were significant difference (P = 0.0000) among three groups by Kaplan-Meier/Log Rank (Mantel-Cox) of survival analysis. By pairwise comparisons (Log Rank (Mantel-Cox)), statistically significant differences were showed in the survival time of LDS-H 2 O paired group (P = 0.0000) and WPN-H 2 O paired group (P = 0.0000), but no significant difference (P = 0.230) in LDS-WPN paired group (Table 1). To sample enough snails in the subsequent experiment, snails that were still living after molluscicide treatment for 12 h were selected.

Functional Annotation
Of the all obtained unigenes, 47,906 (18.84%) had significant BLAST hits to homologs in one or more databases (Nt, Nr or Swiss-Prot), with size ranging from 201 bp to 22,915bp, and an average length of 1445 bp. For most of unigenes which a blast hit is unavailable, the length of unigenes ranged from 203 bp to 10,124 bp, with an average length of 490 bp, but all of them contained protein-coding regions with the analysis of GETORF in EMBOSS 3.0. GO analysis was carried out on the putative proteins blasted against the Swiss-Prot database. Of the 35,198 unigenes that could be annotated as putative proteins in Swiss-Prot, 14,860 were assigned one or more GO terms. For the control (H 2 O-treated) group, of 21,620 unigenes matching putative proteins, 8,927 unigenes were assigned one or more GO terms. Of 26,996 and 34,331 unigenes matching putative proteins in the LDS-and WPN-treated groups, 11,173 and 14,545 unigenes were assigned one or more GO terms, respectively ( Table 2). All GO assignments fell into broad categories for all three major GO functional domains (Cellular Components, Molecular Function and Biological Process), as presented in Fig. 1. The subcategories in each group are organized in similar patterns. Most of the unigenes in the Cellular Components category were classed as "cell part", "cell" and "organelle"; "binding" and "catalytic activity" dominated in the Molecular Function category. The Biological Process category was composed primarily of "cellular process," followed by "metabolic process," "biological regulation" and "response to stimulus".
KEGG analysis demonstrated the biological pathways in which unigenes are involved. A total of 4,686 unigenes of the 35,198 putative proteins were mapped to 250 KEGG pathways (S2 Table) in six categories and 41 subcategories (Fig. 2). Of the 250 KEGG pathways, there were 89 pathways involved in metabolism, 22 pathways involved in genetic information processing, 56 pathways involved in organismal systems, 15 pathways involved in cellular processes, 18 pathways involved in environmental information processing and 50 pathways involved in human diseases. Among these categories, "metabolic pathways" (1,311 unigenes; 27.98% of KEGG-annotated unigenes), "purine metabolism" (357 unigenes; 7.62% of KEGG-annotated unigenes), "pyrimidine metabolism" (279 unigenes; 5.95% of KEGG-annotated unigenes) and "oxidative phosphorylation" (200 unigenes; 4.27% of KEGG-annotated unigenes) were dominant (S2 Table). The gene catalog provided a comprehensive understanding of the gene-transcription profiles of O. hupensis and a valuable foundation for screening differentially expressed genes after treatment with molluscicides.

Differentially Expressed Unigenes Involved in the Response to Molluscicide Exposure
To investigate changes in gene expression and identify the critical genes in the O. hupensis response to molluscicide treatment, the expression levels of unigenes in each library were calculated using the RPKM algorithm. The RPKM value can be used directly to compare the differences in gene expression among samples. The frequency of the RPKM distribution in three samples has a similar pattern (S2 Fig.); most of unigenes belong to RPKM interval (0, 10), and the RPKM interval (90, 100) is the minimum frequency for each sample. A Q value (justified P value) ⩽ 0.05 was used as a threshold to evaluate the significance of differences in  unigene expression. According to quantification analysis (S3-S5 Tables), a total of 254 unigenes showed significantly differential expression between any pair, whereas no unigene showed a significant difference among all three groups. The amounts of differentially expressed unigenes with each molluscicide treatment are presented in Fig. 3A. Among these unigenes, compared to the H 2 O-treated group, 33 and 12 unigenes were up-regulated and down-regulated at 12 h after LDS exposure, respectively, and 83 and 119 unigenes were up-regulated and down-regulated at 12 h after WPN exposure, respectively. Forty-two unigenes overlapped among differentially expressed unigenes of paired groups (Fig. 3B) To understand the functions of differentially expressed unigenes, GO functional and KEGG pathway analyses were performed for each paired group, although the number of unigenes that could be annotated was limited. For 45 unigenes differentially expressed between the H 2 Otreated and LDS-treated groups, 12 unigenes had significant BLAST hits with homologs in one or more databases (Nt, Nr or Swiss-Prot). Among these unigenes, 10 were annotated as putative proteins in Swiss-Prot, two were assigned one or more GO terms, and three were mapped to KEGG pathways. Of the 202 differentially expressed unigenes in the H 2 O-WPN paired group, 105 unigenes had significant BLAST hits with homologs in one or more databases and 88 were annotated as putative proteins in Swiss-Prot. Among these 88 unigenes, 27 were assigned GO terms and 13 unigenes were mapped to KEGG pathways. Thirty-eight out of 49 unigenes differentially expressed between the LDS-and WPN-treated groups had significant BLAST hits to homologs in known databases. Six of these 38 unigenes were assigned GO terms, and seven were mapped to KEGG pathways (S6 Table). When comparing all GO terms of differentially expressed genes in the pairwise groups, "cellular process" and "metabolic process" were the most frequently enriched terms involved in Biological Process for any pairwise group, namely, H 2 O-LDS, H 2 O-WPN and LDS-WPN. "Cell part" and "cell" were the most frequently enriched terms involved in Cellular Components for any pairwise group. For GO terms classified as Molecular Function, the most frequently enriched terms were associated with "binding" and "catalytic activity", although no unigene involved in "catalytic activity" was enriched in LDS-WPN (Fig. 4).

Validation of Illumina Expression Patterns by qRT-PCR
To confirm the reliability of the sequencing analysis and to verify the relationship of the unigenes after molluscicide treatment, 13 unigenes in the H 2 O-LDS pairwise group were selected and detected by qRT-PCR, including 11 differentially expressed unigenes and 2 unigenes that showed no significant difference in sequencing. Before qRT-PCR, all candidates were detected in Oncomelania cDNA and showed a band of the expected size through agarose-gel electrophoresis by semi-quantitative PCR. The profile of their expression (Fig. 5) was consistent with the results of the sequencing. The expression of two unigenes with no different expression in sequencing also showed no significant difference after LDS and WPN exposure in real-time PCR, with the exception of unigene 351626, which were up-regulated after WPN exposure, and the expression level was significantly different, consistent with the sequencing result. Five of those 11 differentially expressed unigenes were down-regulated, and six of those 11 unigenes were up-regulated after LDS or WPN exposure, and the expression levels were both significantly different compared to the H 2 O-treated group, with the exception of unigene 148608, which showed no significant difference in the H 2 O-WPN pairwise group, consistent with the sequencing result. Although unigene 313329 showed no significant difference (Q = 0.078) in expression level between H 2 O and WPN by sequencing, the difference was confirmed to be significant by real-time PCR (P = 0.007) and transcriptome sequencing (P = 4.53×10 -5 ).

Discussion
According to the immersion assay, the molluscicidal activity of LDS against O. hupensis at an active concentration of 0.1 mg/L showed no significant difference compared to that of WPN, confirming previous results of lab and field assays by Yuan et al. [29]. The major difficulties in using molluscicides to control snails are the cost and the toxicity to fish. Previous studies imply that LDS has the least acute toxicity to fish among the molluscicides accepted for use by the Chinese government [25,29]. In addition, based on its application in 23 counties in 8 provinces, LDS has 13.17-36.84% lower costs than WPN (data not shown). LDS is thus strongly recommended as a new choice of molluscicide. It is necessary to understand the molluscicidal mechanism of LDS and determine the difference between LDS and WPN, which will have implications for identifying potential key molecules related with metabolism, apoptosis and cell death after chemical treatment. This will help us to understand metabolic pathways or apoptotic programs which may be vital to snail development and death, to snail control finally.
With continued computational developments, high-throughput mRNA sequencing technology with de novo assembly is especially suitable for gene-expression profiling in those species that do not have closely related reference genomes [32,33]. Here, by profiling the transcriptome of O. hupensis on the Illumina Hiseq 2000 platform, 29 Gb of coverage with 145,008,667 clean reads was obtained and assembled into 254,286 unigenes by de novo assembly. This work provided the largest genomic dataset for O. hupensis, the first extensive genetic resources for this biomedically important freshwater snail, for which very limited genetic or genomic information was available in public databases. Previously available data included 1,432 EST sequences from an O. hupensis hepatopancreas cDNA library [34] and head-foot SSH-library [35], partial and complete mitochondrial sequences [36], several internal transcribed spacer sequences [37], ribosomal RNA genes and microsatellites [38,39]. Due to the small proportion of mollusk-specific sequences in current databases, only 18.84% (47,906 unigenes) of 254,286 unigenes were successfully matched to homologs in Nt, Nr or Swiss-Prot. To confirm that the transcriptome dataset presented here corresponds to O. hupensis, 10,000 reads were randomly selected and analyzed, 1,717 showed significant BLAST hits with homologs in the Nt database, and 38.7% (665 reads) of which most closely matched O. hupensis.
Based on general GO analysis, 42.2% (14,860 unigenes) of the 35,198 putative protein sequences could be annotated, suggesting that they have relatively conserved functions. Although these sequences comprise a small part of the transcriptome sequences, this identification will certainly improve our understanding of the development of snail organs or tissues, biological regulation and metabolites, stimulus and immune response. Considering the key role of this species in schistosomiasis transmission and control, we focused on the immune-system-related and death-related genes. Many new genes have been identified, including 3,377 unigenes involved in immune-system processes, 790 unigenes involved in apoptosis and 374 unigenes involved in cell death. Not only the gene or protein names and descriptions but also their putative conserved domains, GO terms, and potential metabolic pathways were assessed. For example, 48 unigenes were identified as toll-like receptors (TLRs) in this study, with high identities to TLR1-4, TLR6, TLR8, TLR13 and TLR22 from oysters to humans (data not shown); most of these TLRs have a conserved structure, including leucine-rich repeats (LRR), transmembrane (TM) regions and a Toll-Interleukin 1 receptor (TIR) domain. In addition, 5 unigenes were identified as MyD88, and several unigenes showed high identities to TRAF6, TRAF3, IRAK4, TBK1, NFκB1 and NFκB2. All of these molecules were confirmed to be involved in the TLR signaling pathway [40], which plays a central role in the initiation of innate immunity, including host-cell recognition and response to pathogens [41]. These findings from transcriptome analysis will be beneficial to understanding the role of TLRs in inducing innate immune responses and in shaping the adaptive immunity of freshwater snails. KEGG predictions allowed us to identify 13.3% (4,686 unigenes) of the 35,198 putative proteins mapped to 250 KEGG pathways. Of these 4,686 unigenes, 59 contributed to apoptosis and 41 contributed to TLR signaling. These transcriptome data will provide new insights and facilitate further studies of O. hupensis genes and their functions.
In sequencing transcripts from groups treated with different molluscicides, 45 unigenes showed significant differential expression after 12 h of exposure to LDS, of which 33 were upregulated and 12 were down-regulated; 202 unigenes showed significant differential expression after 12 h of exposure to WPN, of which 83 were up-regulated and 119 were down-regulated. Thus, the gene-expression pattern differed even opposite after exposure to LDS and WPN, as confirmed by the 49 unigenes that were differentially expressed between LDS and WPN exposure. Because the available matched putative proteins were limited, information from GO and KEGG analysis for the differentially expressed unigenes was not sufficient to identify molecules or pathways critical to inducing snail death after molluscicide treatment. However, the unigenes involved in "cellular process" and "metabolic process" showed the most significant differences after LDS and WPN exposure, and unigenes involved in "binding" and "catalytic activity" also showed significant differences after LDS and WPN exposure. This result was consistent with ultrastructural and metabolic enzyme changes observed in snails in our study. After LDS and WPN exposure, the muscle fibers loosened and became disordered with widened gaps, neural fibers severed, cell membranes collapsed, the number of glycogen granules of the endoplasmic reticulum decreased, mitochondria deformed, vacuolization occurred, nuclear membranes ruptured or disappeared, heterochromatin polarized in the hepatopancreas and muscle cells, and the number of endocrine granules decreased and that of secondary lysosomes increased in hepatopancreas cells (data not shown). In addition, the activity of some enzymes showed changes after molluscicide treatment; for example, alanine aminotransferase (ALT), lactate dehydrogenase (LDH) and nitric oxide synthase (NOS) activity decreased dramatically, whereas the activity of superoxide dismutase (SOD) increased followed by decreasing to levels below those in the control group. This result may imply that molluscicides cause defects in cell structure, inhibit neurohumoral transmission and energy metabolism, and cause snail death. The cell defects in LDS-treated snails were more severe than those in the WPN-treated group. This difference is perhaps related with the 49 differentially expressed unigenes, or the opposite regulation pattern on gene expression that require further identification and investigation.
In summary, the first published large-scale transcriptome analysis of O. hupensis, the only intermediate host of S. japonicum, is presented here. A total of 145,008,667 clean reads were assembled into 254,286 unigenes. The availability of this transcriptome dataset forms the basis for further studies of the growth, development, death and immune response of O. hupensis, which will be beneficial for snail control and interruption of disease transmission through biological methods. In total, 254 unigenes showed significantly differential expression after molluscicide treatment. GO and KEGG analyses, as well as morphological analyses, indicated that cell structure defects and the inhibition of neurohumoral transmission and energy metabolism caused by LDS and WPN lead to snail death. These findings will be more fully elucidated after these unigenes have been identified and annotated completely.