Oncomelania hupensis retains its ability to transmit Schistosoma japonicum 13 years after migration from permissive to non-permissive areas

Background The East Route Project (ERP) of the South-to-North Water Diversion Project (SNWDP) stretches across schistosomiasis endemic and non-endemic areas in China, which may lead to the dispersal of Oncomelania hupensis, the intermediate host of Schistosoma japonicum, from permissive areas along the Yangtze River Basin to non-permissive areas in northern China. A previous survey demonstrated that O. hupensis could survive and breed for 13 years (12 generations) after being transferred to a non-permissive area, and could be infected by S. japonicum. However, it is not clear if the migrated snails will change their ability to transmit S. japonicum. Methods We infected mice with the cercariae released from the infected transferred snails bred in Jining city of Shandong Province (non-permissive areas) for 13 years. The mice in the control group were infected with cercariae derived from the snails collected in their original habitat (Jiangdu county of Jiangsu Province, permissive areas). Then, we explored the pathogenicity to mice including worm burden, liver egg count and pathology. Additionally, the gene expression profiles of the adult male and female worms recovered from the infected mice were analyzed by RNA sequencing. Results The worm burden, liver egg count and pathology of the mice infected with cercariae released from transferred snails bred in non-permissive areas for 13 years showed no significant differences, when compared with the control cercariae. Slight changes occurred at the transcription level between adult male and female worms recovered from mice infected with cercariae derived from snails bred in permissive and non-permissive areas. Only fourteen genes were significantly differentially expressed in the comparison of adult female worms, and no significantly differentially expressed gene was found in the comparison of adult male worms. Conclusions Our findings strongly suggest that transferred snails did not change their schistosomiasis transmission ability and the worms derived from them retained the original pathogenicity, even after migrating from permissive to non-permissive areas for 13 years. Therefore, a long-term surveillance system of snails along the SNWDP is urgently needed to prevent the diffusion of O. hupensis and reduce the risk of transmission of schistosomiasis.


Background
Schistosomiasis japonica, caused by trematode parasite Schistosoma japonicum, is a serious public health issue that leads to human morbidity and mortality, primarily in southern China and large parts of the Philippines [1,2]. Over the past six decades, great achievements have been attained in reducing S. japonicum infections through comprehensive control measures [3,4], but there were still 37,601 schistosomiasis patients in China at the end of 2017, including 1 acute and 29,407 advanced documented cases of schistosomiasis [5]. Over 90% of current schistosomiasis cases occur in the lake and marshland regions of Jiangsu, Jiangxi, Hunan, Hubei and Anhui provinces with vast areas of Oncomelania hupensis habitats [5], which greatly challenge the progress of schistosomiasis transmission control.
The amphibious freshwater snail O. hupensis is the only intermediate host of S. japonicum [6]. The snails become infected when the miracidia hatched from mature eggs penetrate into their bodies and then several larval stages develop, including sporocysts and cercariae [7]. The geographical distribution of snails defines the areas where schistosomiasis is endemic in China, so one of the main approaches is to interrupt the transmission of the disease to control the snails [8][9][10]. In China, O. hupensis is dominated at latitudes below 33° 15′ N [7]. Besides the environmental factors such as temperature, soil type and vegetation [11,12], inevitably some water-conservancy projects affect the dispersal of snails and schistosomiasis epidemic, such as Aswan Reservoir in Egypt, Akosombo Dam in Ghana and Gezira Dam in Sudan [13][14][15][16].
The South-to-North Water Diversion Project (SNWDP), which was put into operation at the end of 2013, is a major strategic project to transfer part of the abundant water resources in the Yangtze River Basin to North and Northwest China. This action changes the situation regarding water shortages in North China, as well as providing flood prevention and drought eluviation in South China [17]. One of the main intakes of the Eastern Route Project (ERP) is in Jiangdu county of Jiangsu Province which is heavily infested with O. hupensis [18]. The route of ERP crosses Baoying county of Jiangsu Province (at 33° 15′ N), the current northern limit zone of O. hupensis distribution in China [19], and then passes northward into Shandong and Hebei provinces [17], nonpermissive areas for O. hupensis survival. Following the construction of the project, studies have been carried out on whether O. hupensis will be brought to the north of China along with the water flow and form new snail habitats across 33° 15′ N, which then can eventually lead to the expansion of S. japonicum endemicity [20][21][22][23][24][25]. However, there is still a lack of long-term longitudinal observations to indicate whether O. hupensis can survive and reproduce in cold areas of northern China and maintain their ability to transmit schistosomiasis.
In our previous study, we transferred O. hupensis from Jiangdu county of Jiangsu Province (permissive areas) to Jining city of Shandong Province (non-permissive areas north of 33° 15′ N), where the ERP of SNWDP passes through. Unlike the results reported from other studies [20,[26][27][28], the snails survived and spawned for 13 years (12 generations) post-migration. During this time period we observed snails infected with miracidia and fed into the natural environment of the lakeshore. Interestingly, the offspring of the transferred snails continued to maintain the susceptibility to S. japonicum infection. In the present study, the cercariae released from the offspring snails were collected and used to infect mice. The worm burden, liver egg count and pathology showed no statistical differences at 50 days post-infection (dpi) compared with the cercariae derived from O. hupensis collected in Jiangdu county of Jiangsu Province which were infected with the same batch of miracidia of S. japonicum and bred under the same conditions. Furthermore, we utilized RNA sequencing to investigate the differences at the transcription level between adult S. japonicum derived from O. hupensis transferred to non-permissive areas and adult parasites derived from snails from their original permissive areas. Unsurprisingly, there were just slight changes at the transcription level. Our findings indicate that snails transferred from a permissive to a non-permissive area retain the original schistosomiasis transmission capacity. Therefore, a long-term surveillance system of a possible schistosomiasis epidemic along the ERP of SNWDP should be constructed to prevent the spread of snails and reduce the risk of schistosomiasis transmission.

Keywords: Schistosoma japonicum, Oncomelania hupensis, Pathogenicity, Transcriptome, Differentially expressed genes, South-to-North Water Diversion Project
The transferred snails were fed in the natural environment of the lakeshore, and survived and spawned for 13 years (12 generations, 2004-2017) (hereinafter referred to as "Shandong snails"). In the marshland of Lake Weishan, we rebuilt a natural environment isolated from adjacent areas as the experimental field, where the Jiangsu snails and Shandong snails were infected with miracidia of S. japonicum (Jiangsu isolate) and bred from July 3rd to November 10th 2017 (130 days). Infected Jiangsu snails and Shandong snails were selected to release cercariae and infect 12 and 13 ICR mice (age, 6-8 weeks; Shanghai SLAC laboratory Animal Co., Ltd., Shanghai, China; 20 cercariae per mouse), respectively. Mice were housed under specific pathogen-free conditions and fed autoclaved food and water ad libitum. Mice were maintained under environmentally controlled conditions (12:12h light:dark photocycle, temperature of 25 °C and relative humidity of 50%). At 50 dpi, mice were sacrificed, and the mature adult worms were recovered by hepatic perfusion. Male and female worms were separated manually under a dissecting microscope, counted for calculating worm burden and then pooled together into four groups: adult male and female worms derived from infected Jiangsu snails and Shandong snails (hereinafter referred to as "JS-M", "JS-F", "SD-M" and "SD-F"). There were three biological replicates in each group and twelve samples of worms were carefully washed and immediately frozen in liquid nitrogen until RNA extraction.

Liver egg count and pathology
Livers of 4 and 6 mice infected by the cercariae were released from Jiangsu and Shandong snails, respectively. They were selected randomly for liver egg count and pathology analysis at 50 dpi. For each mouse, the right quadrate lobe of each liver was removed, weighed, homogenized and digested with 10% NaOH for 10 min at 56 °C. The eggs were counted and the number of eggs per gram of liver tissue per pair worm were calculated. To compare the liver pathology, the left quadrate lobe of each liver from the two groups was fixed in 4% formalin. Paraffin-embedded sections of these samples were prepared and stained with haematoxylin and eosin (H&E).

cDNA library construction and RNA sequencing
Equal amounts of the twelve samples were prepared for cDNA synthesis and RNA sequencing. Before the construction of the cDNA libraries, ribosomal RNAs were removed and poly(A) + mRNAs were isolated with magnetic Oligo-dT beads (Invitrogen), then the libraries were constructed and sequenced by Majorbio Biotech Co., Ltd. (Shanghai, China). Briefly, the libraries were constructed by TruSeq TM RNA Sample Prep Kit (Illumina) using 5 μg of the total RNA for each sample. The synthesized DNA templates were enriched by PCR amplification for 15 cycles. The products were collected and purified via Certified Low Range Ultra Agarose (Bio-Rad, Hercules, USA) gel electrophoresis. Before sequencing, the cDNA libraries were quantified using a TBS-380 micro fluorometer with Picogreen ® reagent (Invitrogen). Clone clusters were produced on an Illumina cBot, using the TruSeq PE Cluster Kit v3-cBotHS, and high-throughput sequencing was carried out on the Illumina HiSeq4000 TruSeq SBS Kit v3-HS (200 cycles).

Bioinformatics analysis of sequencing data
Clean reads were obtained with Trimmomatic v0.38 [29] by removing adapter sequences, empty reads, poly-N containing reads and low-quality sequences with a Q quality score lower than 20 from the raw data. Meanwhile, Q20, Q30 and GC-content of the clean reads were calculated. All of the downstream analysis was based on clean reads with high quality. The improved genome of S. japonicum [30] was used as the reference genome. An index of the reference genome was constructed and paired-end clean reads were aligned to the S. japonicum reference genome using HISAT v2.1.0 [31]. HTSeq v0.11.2 [32] was adopted to count the number of reads mapped to each gene. The genes with mapped reads of all samples > 10 were included for differential expression analysis. Differential expression analysis of the two groups (three replicates of each group) was performed using the DESeq2 R package (1.18.0) [33]. The P-values were adjusted using the Benjamini and Hochberg's approach for controlling the false discovery rate (FDR). Transcripts were determined to be significantly differentially expressed if FDR < 0.05 and > 2-fold change was observed in the level of expression.
Pearson's correlation analysis used for RNA sequencing was performed to examine the correlation of gene expression levels between samples. The correlation of gene expression levels between samples is an important index to test the reliability of experiments and the rationality of sample selection. The closer the correlation coefficient is to 1, the higher the similarity of expression patterns between samples.

qRT-PCR verification of RNA sequencing data
To validate the transcription levels of the genes identified by RNA sequencing, 14 genes (significantly differentially expressed between JS-F and SD-F) were measured using qRT-PCR, based on their involvement in different expression levels in all comparative analyses (JS-M vs SD-M, JS-F vs SD-F). The total RNA was treated with the PrimeScript ™ RT reagent kit with gDNA Eraser (Takara, Dalian, China) to remove genomic DNA and synthesize cDNA template, according to the manufacturer's instructions. Using the RNA sequencing data, primers used for qRT-PCR were designed using Primer Premier 6.0 software (Premier Biosoft International, Palo Alto, CA, USA); these are listed in Table 1. PSMD4 (26S proteasome non-ATPase regulatory subunit 4) was used as an internal reference gene [34]. All qPCR reactions were performed on a LightCycler ® 96 (Roche, Basel, Switzerland) using 2× SYBR green qPCR master mix (Bimake, Houston, USA) according to the manufacturer's instructions. Each 10 μl qPCR reaction mixture comprised a 1 μl (1:4) of cDNA, 5 μl 2× SYBR green master, 0.4 μl (5 μM) of each primer and 3.2 μl ddH 2 O. The qPCR cycle parameters were as follows: 95 °C for 5 min, followed by 40 cycles of 95 °C for 15 s, 60 °C for 30 s; melt curve analysis ranged from 60 °C to 95 °C to ensure that the specific product was amplified in each reaction. The 2 −ΔΔCq method was used to calculate the relative fold change of the differentially expressed transcripts [35]. All RNA samples for qRT-PCR validation were the same as those for the Illumina library synthesis.

Statistical analysis
Data are presented as the mean ± standard deviation (SD). All statistical analyses were performed using a Student's t-test. A probability (P) value of ≤ 0.05 was considered statistically significant. GraphPad Prism software (Version 6, GraphPad Software, La Jolla, CA, USA) was used for all statistical analyses.

Worm burden, liver egg count and pathology
Each mouse in the two groups was infected with 20 cercariae released from Jiangsu and Shandong snails, respectively. At 50 dpi, the adult worms and livers were collected from the two mice groups for comparing the worm burden, liver egg count and pathology. All mice in the two groups were successfully infected with schistosomes. The mean number of adult worms recovered from the two groups were 10.33 ± 1.83 (Jiangsu snails) and 8.77 ± 1.13 (Shandong snails) and the difference was not statistically significant (t = 0.74, df = 1, P = 0.47) (Fig. 2a).
H&E staining showed that there were large numbers of eosinophilic granulomatous nodules in the portal area caused by the eggs excreted from adult worms derived from Jiangsu and Shandong snails. The size of the nodules was different, and some of them had fibrosis contained clusters of Schistosoma eggs. A large quantity of eosinophils and small lymphocytes infiltrated around the eggs. The liver pathology in the two groups of mice was similar (Fig. 2c).

Transcriptomic features of adult S. japonicum derived from O. hupensis bred in permissive and non-permissive areas
The RNA integrity number (RIN) of all twelve RNA samples used for RNA sequencing was > 7. A total of 575,398,346 raw reads were obtained from all twelve cDNA libraries. After eliminating adaptor and low-quality sequences, 570,783,338 clean reads were gained for the subsequent analysis. All clean reads were submitted to the Sequence Read Archive (SRA) database at NCBI (accession No. PRJNA579703). The Q30 percentage (Q30 percentage is the proportion of sequencing bases error rate with quality value < 0.1%) and GC percentage were 95.52% and 43.62%, respectively ( Table 2). The alignment rate to the improved genome of S. japonicum ranged from 83.74% to 98.69%. Over 65% of the clean reads were distributed in exon regions and the rest were distributed in the introns. Pooled clean reads of all samples were mapped to a total of 8578 genes in the S. japonicum genome. Moreover, Pearson's correlation coefficient of gene expression among different samples was greater than 0.9, indicating the high similarity of the gene expression patterns between samples (Fig. 3).

Analysis of differential expression
In the present study, a threshold of FDR < 0.05 and absolute fold change ≥ 2 was used to define significantly differentially expressed genes in pairwise comparisons, including JS-M vs SD-M and JS-F vs SD-F. Of a total of 8578 genes, as shown in the volcano plots (Fig. 4), there was no differentially expressed gene between JS-M and SD-M. Similarly, only 14 genes were differentially expressed between JS-F and SD-F, of which 2 and 12 were significantly upregulated in JS-F and SD-F, respectively. The related information of the 14 differentially expressed genes is listed in Table 3. Interestingly, four of the 12 upregulated genes in SD-F are collagen-related genes (EWB00_010449, EWB00_006199, EWB00_000274 and EWB00_006201) and one of the downregulated genes is lactate dehydrogenase (EWB00_010161).

qRT-PCR validation of RNA sequencing data
To confirm the RNA sequencing data, 14 genes differentially expressed between JS-F and SD-F were further detected by qRT-PCR. According to the RNA sequencing results, the expression levels of EWB00_005036 and EWB00_010161 were upregulated in JS-F, the expression levels of EWB00_002418, EWB00_003584, EWB00_003582, EWB00_008342, EWB00_005617, EWB00_005592, EWB00_010449, EWB00_011245, EWB00_009712, EWB00_006199, EWB00_000274 and EWB00_006201 were upregulated in SD-F, and the expression levels of the 14 genes mentioned above were not significantly different between JS-M and SD-M. After normalizing to the reference gene PSMD4, expression levels determined by qRT-PCR correlated with those obtained by RNA sequencing (r = 0.8400, P < 0.0001), validating the accuracy and reliability of the RNA  sequencing results (Fig. 5). Therefore, the data generated here can be used to investigate the difference in transcription levels between adult S. japonicum derived from O. hupensis bred in both Jiangsu and Shandong Province.

Discussion
The SNWDP stretching across schistosomiasis endemic and non-endemic areas in China, may lead to the dispersal of O. hupensis from permissive areas to non-permissive areas in northern China. Our previous study showed O. hupensis can survive and breed for 13 years (12 generations) after being transferred to a non-permissive area, and maintain the susceptibility to S. japonicum infection (Miao et al., unpublished data) thus highlighting the extreme significance of investigation of the ability transferred snails to transmit S. japonicum.
In the present study, we first evaluated the difference in pathogenicity to mice, using cercariae of S. japonicum derived from the transferred O. hupensis bred in non-permissive areas (Jining City of Shandong Province) and cercariae released from the snails from the original habitat (Jiangsu Province). The results indicated no significant differences between the two groups regarding worm burden, liver egg count and pathology. Furthermore, to explore the differences at the molecular level of the worms developed from cercariae released from Jiangsu and Shandong snails, RNA sequencing technology was used for exploring gene expression profiles of adult S. japonicum (at 50 dpi) derived from the two snail groups. Analysis of differential gene expression showed that no gene was differentially expressed between adult male groups (JS-M vs SD-M) and only 14 genes were differentially expressed between adult female groups (JS-F vs SD-F). This result is not surprising. He et al. [36] found that there was no significant genetic divergence between the snail population in Lake Weishan and its original population in River Yangtze using microsatellite DNA markers, 10 years after migration. These results reconfirm the lack of change in pathogenicity and slight changes in the transcriptional levels of the adult worms  in the present study. On the other hand, the differentially expressed genes in males derived from snails bred in permissive and non-permissive areas were fewer than in females, which may be related to the fact that the gene expression in males is more stable than that in females in the course of development [37]. Of the 14 differentially expressed genes, four of the 12 highly expressed genes in SD-F are collagen-related genes (EWB00_010449, EWB00_006199, EWB00_000274, EWB00_006201). Collagen proteins either become crosslinked into high-strength mature fibers and are distributed in the spaces around the cells or form heterotypic fibrils, distributing in the bones, muscles, tendons, ligaments, skin, cornea and other parts of the organism [38,39]. Besides their structural roles, collagens also have numerous potential functions during growth, repair and morphogenesis, and in pathobiological processes. RNAi silencing of type V collagen in S. japonicum was shown to significantly reduce the hatching rate of eggs and single female spawning rate and also affect the morphology of the worms [40]. One of lowly expressed genes (EWB00_010161) in SD-F is lactate dehydrogenase (LDH), encoding a terminal glycolytic enzyme that catalyzes the interconversion of pyruvate and lactate in the presence of the nicotinamide adenine dinucleotide coenzyme. Schistosome parasitic stages depend mainly on anaerobic energy metabolism and glycolytic enzymes are essential for the survival of parasites [41,42]. Four genes (EWB00_005036, EWB00_003584, EWB00_003582 and EWB00_009712) have no functional domains or unknown domains. The remaining of the highly expressed genes in SD-F are involved in protein modification and transcriptional regulation. In the process of transcriptome sequencing, the changes of the 14 genes were negligible. The materials used for RNA sequencing were from the same group of mice at the same time and the only different factor was the snail origin. Whether the differentially expressed genes are related to the intermediate host, O. hupensis, remains to be further investigated.
Additionally, we performed qRT-PCR to detect the 14 differentially expressed genes for validating the RNA sequencing data and found that the consistency between the two methods was acceptable. To obtain better measurements, there may be a need to further increase the sequencing depth and include more biological replicates for providing greater statistical power.
This study represents a snapshot of comparative transcriptome analysis for adult S. japonicum (at 50 dpi) derived from O. hupensis bred in Jiangsu Province and Shandong Province. In future studies, we could collect sporocysts, freshly released cercariae or schistosomula to carry out a comparative analysis of gene expression.
It is possible that more differences at the molecular level among the comparison groups will be found. In addition, the same effect may be achieved by adding different perspectives of investigation, such as epigenetic analysis or protein modification. On the other hand, considerable variation between some of S. japonicum populations has been reported for genes inferred to be linked to key cellular processes and/or pathogen-host interactions [43]. In our study, we only used one population in the infection experiment which showed no significant differences. However, it is possible to find some differences if using different populations of worms.
Oncomelania hupensis transferred from Jiangdu County (Jiangsu Province) to Lake Weishan (Shandong Province) could survive and spawn for 13 years (12 generations) after migration. The offspring of transferred snails still maintained the susceptibility to S. japonicum infection. In this study, we provide further evidence that the pathogenicity of cercariae released by the snails remained unchanged. Moreover, the adult S. japonicum derived from O. hupensis transferred from permissive to non-permissive areas changed only slightly at the transcriptomic level. Altogether, these results suggest that a long-term surveillance system of schistosomiasis should be established for large-scale water conservancy projects, especially those which pass through schistosomiasis endemic areas and non-endemic areas, to prevent the diffusion of snails and reduce the risk of schistosomiasis transmission.

Conclusions
No difference in pathogenicity to mice was observed between the cercariae of S. japonicum from O. hupensis bred in Jiangsu and Shandong Province, two sites located on the ERP of SNWDP that represent permissive and non-permissive areas for snail breeding, respectively. Unsurprisingly, transcriptome profiles of adult S. japonicum male and female worms (at 50 dpi) derived from the two groups of snails revealed that only 0 and 14 significantly differentially expressed genes were obtained between the comparison groups of male and female worms, respectively. Interestingly, four of the 14 differentially expressed genes between female comparison groups were collagen-related genes and one was lactate dehydrogenase (LDH). However, whether these differentially expressed genes are associated with the intermediate snail hosts, remains to be further investigated. Combined with the results of our previous study, we suggest that O. hupensis transferred from permissive to non-permissive areas retain the original transmission capacity of schistosomiasis. Therefore, a long-term surveillance system of schistosomiasis along the ERP of SNWDP should be constructed to prevent the spread of snails and reduce the risk of schistosomiasis transmission.