Comprehensive Transcriptome Analysis of Sex-Biased Expressed Genes Reveals Discrete Biological and Physiological Features of Male and Female Schistosoma japonicum

Schistosomiasis is a chronic and debilitating disease caused by blood flukes (digenetic trematodes) of the genus Schistosoma. Schistosomes are sexually dimorphic and exhibit dramatic morphological changes during a complex lifecycle which requires subtle gene regulatory mechanisms to fulfil these complex biological processes. In the current study, a 41,982 features custom DNA microarray, which represents the most comprehensive probe coverage for any schistosome transcriptome study, was designed based on public domain and local databases to explore differential gene expression in S. japonicum. We found that approximately 1/10 of the total annotated genes in the S. japonicum genome are differentially expressed between adult males and females. In general, genes associated with the cytoskeleton, and motor and neuronal activities were readily expressed in male adult worms, whereas genes involved in amino acid metabolism, nucleotide biosynthesis, gluconeogenesis, glycosylation, cell cycle processes, DNA synthesis and genome fidelity and stability were enriched in females. Further, miRNAs target sites within these gene sets were predicted, which provides a scenario whereby the miRNAs potentially regulate these sex-biased expressed genes. The study significantly expands the expressional and regulatory characteristics of gender-biased expressed genes in schistosomes with high accuracy. The data provide a better appreciation of the biological and physiological features of male and female schistosome parasites, which may lead to novel vaccine targets and the development of new therapeutic interventions.


Introduction
Schistosomiasis, caused by infection with blood flukes (digenetic trematodes) of the genus Schistosoma, remains one of the most serious parasitic diseases worldwide, afflicting more than 200 million people, with close to 800 million at risk [1,2]. This debilitating disease causes an annual number of disability-adjusted life years (DALYs) lost of up to 3.3 million in 2010, ranking it as third in the list of global neglected diseases [3]. Three main species, Schistosoma mansoni, S. haematobium and S. japonicum are of clinical relevance. Currently, no practical anti-schistosome vaccine is available and mass chemotherapy with a single effective drug, praziquantel, combined with morbidity management, are the primary strategies adopted for the treatment and control of schistosomiasis [4,5].
Schistosomes have a complex lifecycle involving an aquatic snail as an intermediate host and a mammalian definitive host [2]. In contrast to other trematode species, these parasites are unique in that they exhibit sexual dimorphism, and they thus represent a valuable model for invertebrate conjugal biology research. The availability of schistosome transcriptome [6,7] and genome sequences [8][9][10] for the three major schistosome species, provides a wealth of resources to allow the dissection of gene profiles during development and between the sexes. In this respect, a variety of high-throughput techniques have been widely employed in the study of schistosomes, including the use of microarrays [11][12][13][14][15][16], serial analysis of gene expression (SAGE) [17][18][19], digital gene expression (DGE) [20], and RNAseq [21,22] with each method presenting distinct advantages and disadvantages. These pioneering studies revealed expression patterns and features of developmental-, gender-dependent, tissue-specific, strain-specific and host-associated gene expression within schistosome parasites [11,13,[23][24][25], which have contributed substantially to our understanding of their biology. However, SAGE and DGE suffer from quantifying low abundance mRNA transcripts [20], both of which may omit genes responsible for vital functions present within tissue-specific expression. This is a particular concern given that schistosomes are multi-cellular organisms. The interpretation of DNA microarray results depends on the quality of genetic information contained within the DNA sequences used for probe design. The first generation DNA microarrays used for schistosome studies were designed based on EST transcripts and the data obtained from these chips only provide a compromise interpretation due to the poor annotation of these sequences [11][12][13]. To provide more comprehensive gene profiles during development or between the sexes of schistosomes, a second generation DNA microarray with a well-curated design of probes based on transcriptomic and genomic sequences is required.
MicroRNAs (miRNAs) are small non-coding RNA molecules, which exert important gene regulatory functions at the post transcriptional level. The identification of schistosome miR-NAs has been carried out over the past five years using techniques including cloning methods to deep-sequencing. Comprehensive miRNA expression profiles within several discrete developmental stages of schistosomes, as well as between different sexes have been presented with high accuracy and coverage [26][27][28][29][30][31]. In humans, it has been estimated that miRNAs may regulate the expression of~60% of the total coding mRNA transcripts [32]. Plausible regulatory roles in the development and sexual maturation of schistosomes have been suggested for several miRNAs [26,27]. However, comprehensive functional annotation for individual schistosome miRNAs is still unavailable. In silico prediction represents a high-through approach to achieve this objective, but different prediction tools with different scoring criteria, have led to differing outputs with variable false positive and false negative rates. Recently, a novel experimental approach, using high-throughput sequencing of RNA isolated by cross-linking and immunoprecipitation (HITS-CLIP), has been used to identify miRNA targets in adult S. japonicum worms, although the results from this study were inconclusive [33].
We have constructed a second generation DNA microarray for transcriptomic study of S. japonicum based on S. japonicum and S. mansoni genomic and transcriptomic sequences with multiple probes designed against each target sequence (both forward and reverse) [34][35][36]. By employing this powerful microarray printed with the most comprehensive coverage of probes, we focused on the identification of sex-biased expressed genes and predicted potential miRNA targets against these genes. The study presents a global view of the expressional and regulatory features of gender-associated genes in S. japonicum, and provides novel insights on schistosome conjugal biology.

Ethical statement
All procedures performed on animals within this study were conducted following animal husbandry guidelines of the Chinese Academy of Medical Sciences and with permission from the Experimental Animal Committee (Institute of Pathogen Biology, CAMS) with Ethical Clearance Number IPB-2011-6.

Parasite materials
S. japonicum-infected Oncomelania hupensis were provided by Hunan Institute of Parasitic Diseases, Yueyang, China. Cercariae were shed from these snails under light stimulation and used to percutaneously infect female New Zealand rabbits. Mixed adult worms were also obtained from infected rabbits by perfusion at 6 weeks post-infection (p.i.). Male and female worms were separated manually with the aid of stereomicroscope [27]. All parasite samples were soaked in RNAlater (Ambion, CA, USA), and stored at -80°C until total RNA extraction.

Total RNA isolation
Total RNAs were isolated from male and female S. japonicum using RNeasy Mini kits (QIA-GEN, GmbH, Hilden, Germany) according to the manufacturer's instructions. Potential contaminating genomic DNA was removed from RNA samples using Turbo DNA-free kit (Ambion, CA, USA). The quantity and quality of the RNA samples were assessed by a Nano-DropND-1000 spectrophotometer (NanoDrop Technologies, Wilmington, DE) and denaturing agarose gel electrophoresis.

Microarray construction and hybridization and subsequent data analysis
A schistosome genome-wide microarray was employed for analysing the gene expression profiles of male and female S. japonicum with three biological replicates. The details regarding the design and construction of the microarray, the hybridization method, and feature extraction have been reported previously [34][35][36][37][38]. Briefly, a total of 20,194 S. japonicum target sequences collected for creating an array. For each target sequences, 3 or 4 pairs of 60-mer complementary oligonucleotide probes (forward and reverse probes) were designed. Probes with random sequences were printed as negative controls (background signal) and eight spike-RNA probes from the intergenic sequence of yeast were used as hybridization controls. Microarrays were printed in a 12×135 K feature format (Roche NimbleGen) with a total of 145,000 probes representing 41,982 features. cDNA was labelled with a fluorescent dye (Cy3-dCTP) using a cRNA Amplification and Labelling Kit (CapitalBio, Beijing, China) [39]. Hybridization was performed using three biological replicates for all samples (CapitalBio, Beijing, China). Procedures of array hybridization, washing, scanning, and data acquisition were carried out according to the NimbleGen Arrays User's Guide. The arrays were scanned using a MS200 scanner (Nim-bleGen Systems) at 2-μm resolution, and NimbleScan software (NimbleGen) was used to extract fluorescent intensity raw data from the scanned images. Normalized gene expression data were generated using the Robust Multichip Average (RMA) algorithm [40,41]. Outlier probes were identified and their contribution was reduced at the reported gene expression level, a process which has been shown to improve the sensitivity and reproducibility of microarray results [41]. Then, the expression value of a gene is a weighted average of all forward or reverse probe sets when both background correction and quantile normalization are performed. Raw data and the normalized data have been deposited at the public domain Gene Expression Omnibus under the accession number for the platform GPL18617, and series GSE57143.

Bioinformatics analysis on sex-differentially expressed genes
Potential gender-biased expressed genes of S. japonicum were initially retrieved from the NCBI database (http://www.ncbi.nlm.nih.gov/sites/batchentrez) based on fold-changes (FC) of the mean of the weighted intensity value of forward or reverse probe set between genders (FC 2, three biological replicates). Genes were further considered differentially expressed by FC values from both forward and reverse probe sets 2 between genders (p<0.05, Student's t-test [35,42], without p-value adjustment for multiple testing); only those genes with a mean of signal intensity >100 at least in one gender were included for further investigation. Heat maps were created based on the signal intensities of forward gene or EST sequences using HemI 1.0 software [43]. Gene sets were then functionally annotated using Blast2GO [44]. The gene collection was reannotated using the BLAST program based on the annotation of their homologous sequences from S. mansoni, S. haematobium, Clonorchis sinensis and Echinococcus granulosus, deposited in NCBI database. For hypothetical proteins, conserved protein domains were further searched against the NCBI CDD database (v3.14) [45] for possible improved annotation.

Quantitative real-time PCR
A total of 50 gender-associated and 10 non-gender-associated genes were selected for validation using qRT-PCR as described previously [34]. One microgram male or female total RNA were reverse transcribed into first-strand cDNA using a SuperScript III Reverse Transcriptase Kit (Invitrogen) with oligo (dT) 15 primer according to the manufacturer's instructions. The cDNA products were diluted 20-fold with nuclease-free water before undertaking the qPCR. Each 25 μl PCR reaction contained 12.5 μl of 2×Brilliant II SYBR Green QPCR Master Mix (Agilent, USA), 1 μl cDNA, 1 μl of the forward and reverse primer pair (S1 Table), and 10.5 μl of sterile water. PCR cycling conditions were as follows: 95°C for 10 min, followed by 40 cycles of 30 s denaturation at 95°C and 1 min annealing and extension at 60°C. A dissociation step (95°C for 15 s, 60°C for 1 min, 95°C for 15 s, and 60°C for 15 s) was performed to confirm the amplification specificity for each gene. 26S proteasome non-ATPase regulatory subunit 4 (PSMD4), a reliable reference gene for transcriptomic analysis of S. japonicum [34,46], was employed as a control gene in the assays. The PCR primers were designed using Primer Express 3.0 software (Applied Biosystems, Foster City, USA). PCR reactions were performed in technical triplicates on the 7300 Real-Time PCR system (Applied Biosystems). The relative expression level of each gene was analysed using SDS 1.4 software (Applied Biosystems). Melt curves for the genes tested are shown in S1 Fig Correlations between the microarray and qPCR results for 50 gender-associated genes were checked with the Spearman's correlation coefficient.

miRNA target prediction
The miRNA target sites were predicted using PITA [47], and RNAhybrid [48]. Gene sequences were downloaded from the NCBI website; mature miRNAs were downloaded from the miR-Base (release 19.0) (TPM (transcripts per million) >10 in adult male and female worm libraries [26]). Target sites were first predicted by PITA with the following cutoffs: 1) minimum seed size: 6, and single G:U wobble allowed for seed size of 7 and 8; 2) sites with microRNA-target hybridization energy ΔG duplex score -15 kcal/mol and combined interaction energy ΔΔG score -10 kcal/mol; which were further filtered by RNAhybrid: minimum free energy (mfe) -20 kcal/mol. Target site location within mRNA transcripts (5'-UTR, CDS, or 3'-UTR) was further determined by the annotation available in NCBI database.

Results and Discussion
Global view of gender-biased expressed genes in S. japonicum Based on the results generated from the use of a microarray with the most comprehensive and informative probe design to date, signal intensities from 4,303 and 6,224 sequences were upregulated (FC 2) in male and female adult worms, respectively, which enabled us to retrieve 2,459 (1,344 and 1,115 male and female-biased expressed genes, respectively) potential gender differentially expressed genes from NCBI database (S2 Table). These gene sets further underwent screening with stringent criteria (See Materials and Methods and Table 1). Since alternative splicing and bidirectional transcription are frequent events in the S. japonicum transcriptome [20,46], these criteria can significantly increase the accuracy but may sacrifice sensitivity of the screening procedure to some degree. The procedure finally led us to identify 685 and 430 mRNA transcripts, and 130 and 86 expressed sequence tags (ESTs), exhibiting gender biased expression in male and female worms, respectively (Table 1 and S3-S6 Tables). In comparison with other reports [10,13,21,49], a relatively small number of RNA transcripts (~12%, 85 male-biased and 51 female-biased) were identified as previously "reported" gender-associated genes (S3 and S4 Tables), thereby allowing us to further explore these novel gender-associated genes in S. japonicum.
In general, the number of male-associated genes was higher than female-associated genes, yet more transcripts in female-associated genes exhibit a stronger biased (greater fold change) expression (Fig 1, Fig 2A and 2B). For example, 7% female-associated genes show a strong biased expression (fold change >100) compared with male worms, while none of male-associated genes presented a fold change >100. The percentage of genes showing a high fluorescence intensity (>10,000) in male-and female-biased expressed genes were 13.7%, and 17.5% respectively ( Fig 2C), although this parameter may be affected by the GC composition of the 60-mer probes. A similar result was obtained when analysing the gender-biased EST sequences (S2 Fig), but with a low number of genes (about 1/5 compared to the mRNA data). Furthermore, most of these EST sequences were annotated as either unknown or as hypothetical protein (81 (62.30%) and 47 (54.65%) in the adult male-and female-biased ESTs, respectively) (S5 and S6 Tables). This may reflect the fact they are short sequences from the 5'-untranslational regions (UTR) and 3' UTR of mRNA transcripts [6]. More importantly, a comparison of the mRNA and EST data highlights the power of our second generation S. japonicum DNA microarray in profiling gene expression, since the design of the first generation of S. japonicum chip was based on EST data only. A subset of these ESTs overlaps with those gender-biased mRNA transcripts; i.e., Aromatic-L-amino-acid decarboxylase, 22.6 kDa tegumental membrane-associated antigen, putative wnt inhibitor frzb2, semaphorin-5B, 16 kDa calcium-binding protein, ancient ubiquitous protein 1, myosin heavy chain, paramyosin, calponin-3 and Annexin A3, listed in both the male-biased mRNA and EST transcripts, whereas TES (Trematode Eggshell Synthesis) domain containing protein, UV excision repair protein RAD23, alanine aminotransferase 2 and DNA replication licensing factor mcm7-A were listed in both the female-biased mRNA and EST transcripts. We then focused on analysing the gender-biased mRNA data further.

qPCR validation of DNA microarray data
In order to validate the microarray results, a subset (50, 4.48%) of the gender-differentially expressed genes was selected for validation using qPCR. Generally, the fold changes obtained with the qPCR assay were higher than these obtained by the microarray signals, especially for these extremely biased (high fold change) genes (Fig 3A), which is a phenomenon common in microarray validation experiments [14,50]. Strong correlations were observed between the two methods (for male-biased expressed genes, r = 0.9419, p<0.0001; for female-biased expressed genes, r = 0.9041, p<0.0001) ( Fig 3B). Further, the expression of 10 non-gender-associated genes was also validated by qPCR, which showed good consistency with the DNA microarray data (S3 Fig).

Putative function prediction by GO analysis
To predict the potential function of these sex-biased genes, different functional categories were further determined by Gene Ontology [51] (Fig 4, S7 and S8 Tables). Of the biological process categories, genes involved metabolic and biosynthetic processes are more active in female worms compared to male worms, indicating that the nutritional acquisition is more crucial for female worms, probably reflective of its status of oviposition which requires abundant nutrition for the laying of thousands of eggs per day per worm pair. This finding was consistent with that of a similar transcriptomic study of S. mansoni [21], showing that cellular protein modification process, DNA metabolic process and catalytic activity were the top three enriched categories in females. Of the molecular function categories, more genes associated with protein, ion, small molecule and carbohydrate derivative binding; transmembrane and substrate-specific transporter activity were more highly expressed in male S. japonicum worms, indicating more active host-schistosome interplay (i.e., host ligand-receptor interaction) and energy metabolism than in females. In comparison, assembly, calcium ion binding, protein binding, receptor activity, potassium ion transport and regulation of transmembrane transport were found to be significantly enriched GO categories in S. mansoni males [21]. In the cellular component categories, gene products localised to membrane regions are more abundant in adult S. japonicum males, while gene products localised to membrane-bounded organelles more enriched in adult females.

Analysis of the top 40 gender-biased expressed genes in S. japonicum
We then further analysed the highest, by fold change, gender-biased expressed genes for each sex (Top 40, Tables 2 and 3). Of the male-specific genes, the majority were previously uncharacterized; thus we inferred their putative functions based on their homology to other species. An example being gamma-crystallin related domain-containing protein (FN317557, ranked 3) which is a Ca 2+ binding protein, whose biological functions are not fully understood. However, it has been suggested that crystallin in mammals may function in protecting retinal neurons from damage caused by environmental and/or metabolic stress [52]. Aromatic-L-amino-acid decarboxylase (AY812557, Top 5) participates in dopamine and serotonin (5-HT) neurotransmitter synthesis [53]. It has been shown that the expression of its ortholog in the male worms of S. mansoni was pairing-dependent [54]. Neuroglian (AY811209, Top 10) is a homologue of a L1-type neural adhere molecule in vertebrates. In Drosophila, it has been shown that neuroglian is expressed in a variety of embryonic cells, and one of its primary functions is as a central coordinator of synaptic growth [55]. With semaphorin-5B (AY808492, Top 11), its vertebrate equivalents act as positive axonal growth guidance cues as well as functioning in the regulation of synapse morphology [56]. Similarly, spondin-1 (AY812421, Top 22) is an extracellular matrix protein, with its homologue in C. elegans also acting in axon guidance and fasciculation in motoneurons [57]. The over-expression of these genes in adult males of S. japonicum collectively indicates that the neuronal activities in this sex are more robust, or play a more prominent biological function, when compared to female schistosomes.  In addition, the Top 4 gene (AY811231) encodes a kelch-like protein 10, which may be related to spermiogenesis, since its vertebrate equivalents are testis-specific, and it has been suggested to participate in protein ubiquitination and subsequent proteasomal degradation during spermiogenesis [58]. Another interesting gene highly expressed in males is the zinc finger transcription factor Cubitus interruptus (AY810451, Top 16), whose counterpart in Drosophila can act as a mediator of hedgehog signal transduction to facilitate normal osteoblast differentiation via regulating the expression of the decapentaplegic gene [59], the ortholog of mammalian bone morphogenetic protein (BMP) 2. A homologue of BMP has been characterized in S. mansoni, which is also more abundant in male parasites, but its localization remains elusive [60]. Furthermore, nebulette (AY223389, Top 31) is a nebulin-like protein, its vertebrate counterparts interacting with a number of skeletal proteins such as actin, tropomyosin, alpha-actinin, Xin, and XIRP2 [61], and function in regulating the assembly and lengths of thin filaments in skeletal muscle [62], indicating this protein is related to the motility of the parasite. Previous research in vertebrates has indicated that nephrin is a structural component of the slit diaphragm [63]. However, a recent study showed that the nephrin homologue in planarians is probably expressed in neoblasts, but not in flame cells and neurons [64]. It would be of value to determine the molecular localization of nephrin in schistosomes which may provide further clues to its precise cellular function.
A number of trematode eggshell synthesis (TES) domain containing proteins, an asparagine-rich antigen Pfa35-2, a histidine-rich glycoprotein, an extracellular superoxide dismutase (Ex-SOD), a female-specific 800 protein, two distinct of tyrosinase homologues are listed in the top 40 female-associated genes (Table 3). A finding consistent with previous studies on schistosomes, with the potential molecular functions for some of these genes having been suggested [10,11,15,16,21,49,65]. The data presented here thus show the consistency of high-throughput gene profiling technologies in confirming these highly differentially expressed genes. Intriguingly, it has been shown that DNA vaccination with one of these genes, Cu-Zn SOD, induced a high reduction of worm burden in mice infected with S. mansoni [66]. Further, it is noteworthy to consider some of the novel genes listed in Table 3, which exhibit a variety of molecular functions. For example, the top 10 gene encodes a UV excision repair protein, RAD23, which is involved in nucleotide excision repair (NER). This may be linked to the extensive DNA damage caused by antioxidants released during the detoxification of hemoglobin byproducts in females, which is discussed further below. The gene encoding a large neutral amino acids transporter small subunit 2 (Top 17) was expressed specifically in females, which likely reflects their special physiological status requiring a considerable amount of large neutral amino acids [67]. Furthermore, poly(rC)-binding protein 3 (Alpha-CP3) (Top 25) is a member of the KH-domain containing protein, which binds to RNA transcripts via a C-rich pyrimidine region, and further mediates post-transcriptional regulatory activities. Previously, miRNAs, another post-transcriptional regulatory factor, have also been shown to exhibit sex-biased expression in schistosomes [26,30]. This observation indicates that both Alpha-CP3 and miRNAs are responsible for post-transcriptional gene regulation in S. japonicum by targeting different gene sites. Genes encoding cytoskeleton and motor proteins are more up-regulated in male adult worms In general, previous studies on schistosomes have shown consistently that genes encoding cytoskeleton and motor proteins (i.e., actin, titin, alpha-actinin, dynein light chain, myosin heavy chain, paramyosin, tropomyosin, fimbrin, and troponins) are up-regulated in male worms [10,13,15,16,20,49]. This was also reflected by enriched GO analysis of S. mansoni adult worms showing that microtubule, microtubule-based process and troponins complex were significantly enriched GO categories in males [21]. These data thus indicate a role for males in the physical support of females, which enables the latter to save energy and to focus on nutrient acquisition and egg-laying. More cytoskeleton genes, such as microtubule-associated protein 2, actin related protein 2/3 complex, tensin, laminin were identified in this study as a result of using our comprehensive probe design. In addition, other genes encoding cytoskeleton related proteins, such as PDZ-, multiple PDZ-or PDZ and LIM domain-containing proteins (FN317962, AY810295, AY811780, AY812842, AY815664, AY808539, AY812903 and AY814003) (S3 Table), which contribute to the formation and maintenance of cell complex scaffolding were highly expressed in adult males. The majority of these genes were 2-4 fold more highly expressed in male worms compared with females, which supports the notion that male worms render physical support to females to facilitate their migration against the blood flow from the portal liver sites to the smaller mesenteric circulation where they lay their eggs [11]. Further, members of the PDZ domain-containing protein, such as GIPC3 and Scribble, have been suggested as potential drug targets based on non-canonical protein interaction, though they do not display a sex-biased expression [68,69].

Genes involved in neuronal activities are more up-regulated in male adult worms
Motor activity in schistosomes is closely controlled by the neuronal system. In addition to the male-biased transcripts related to neurotransmitter synthesis, synapse growth and axon development, discussed earlier, additional genes involved in neuronal activities were also found enriched in male worms, including ionotropic glutamate receptor (AY815670), neuron navigator 3 (AY808520), neurogenic locus notch protein-like protein (AY809231), excitatory amino acid transporter (AY810837), neuronal calcium sensor 2 (FN317645), and synaptic vesicle membrane protein VAT-1-like protein (AY811071) (S3 Table). The data reinforce the active neuronal activities in male parasites, which could be linked to the fact that the adult male worms are directly exposed to the cardiovascular system and need to monitor and respond to environmental cues from the host [20]. One should be aware that schistosome neuronal system is not only responsible for motor activity, but also plays an essential role in a wide variety of biological processes, such as cercarial penetration, blood feeding and digestion, waste disposal, reproductive activities, and egg excretion [70], highlighting the significance of the neuronal signaling pathways for parasite survival. Several neuronal receptors, such as the glutamate receptor (SmGluR) [71], serotonin receptor (Sm5HTR) and G protein-coupled acetylcholine receptor, have been suggested as potential targets for novel drug development against S. mansoni [72,73]. Further, it has been shown that two glutamate receptor genes (GRIN1 and NMDA receptor) were up-regulated in paired males treated with PZQ [74]. The expressed products of the neuronal pathway-related genes identified here may represent alternative targets for drug development against the schistosome parasites.
Genes involved in amino acid metabolism, nucleotide biosynthesis and gluconeogenesis are more up-regulated in female adult worms As indicated by the GO analysis, metabolic and biosynthetic processes are more vigorous in female worms (Fig 4A). This is supported by the fact that genes encoding enzymes participating in a variety of metabolic pathways were actively transcribed within this sex. Several members of the venom allergen-like (VAL) family (i.e., VAL 27 (FN318592) and 28 (AY815621)) were also significantly more highly expressed in female parasites. Previously, esophageal secreted proteins encoded by micro exon gene (MEG) 4.1, 4.2, and 14 and VAL-7 in S. mansoni have been shown to play a vital role in erythrocyte lysis and tethering and killing of leucocytes [75]. Recently, SjMEGs 4.1, 8.2, 9, 11 and VAL-7 have been suggested as potential targets of the self-cure process based on observations with the Rhesus Macaques animal model [76]. It would be interesting to determine the localization and function of VAL-27 and 28, which may play a specific role in the biology of female parasites. Increased amino acid metabolism could be expected in females based on the observation that genes involved in amino acid transport (i.e., large neutral amino acids transporter (FN327074) and L-amino acid transporter (FN313722)) and amino acid metabolism (i.e., putative L-asparaginase (AY814032), gamma-glutamylcyclotransferase (AY814775), alanine aminotransferase (AY915267) and aspartate-ammonia ligase (FN326707)) were up-regulated in this sex (S4 Table). In addition, genes related to nucleotide biosynthesis (i.e., adenylosuccinate synthetise (AY816019), ribonucleoside-diphosphate reductase subunit M1 (FN330781) and hypoxanthine-guanine phosphoribosyltransferase (AY915002)) were enriched in adult females, emphasising the increased DNA synthesis is evident during vitellocyte differentiation in female parasites. In respect to energy metabolism, it has been shown that, in schistosomes, two typical facilitated diffusion glucose transporter proteins 1 (GTP1) and 4 (GTP4) are responsible for transporting glucose from the exterior to the inside of worms [77,78]. However, no sexually biased expression was observed for these genes in the current study, in contrast to genes encoding enzymes involved in gluconeogenesis, phosphoenolpyruvate carboxykinase (AY813371) and fructose-1,6-bisphosphatase 1 (FN318294), showing a 2.2 and 2.8 fold up-regulation, respectively, in female than in male worms. These observations indicate that the female parasite may have a relatively increased ability to generate glucose from non-carbohydrate carbon substrates, and this characteristic may need to be considered when targeting glucose metabolism for potential vaccine candidates.

Genes involved in cell cycle processes are more up-regulated in female adult worms
A female worm must pair with a male to become completely sexually mature, when the reproductive organs, mainly the ovary and the vitelline glands, undergo terminal differentiation. The vitellarium contributes two thirds of the mature female body volume, and a select number of cells within this structure undergo stage 2 and 3 differentiation, and further terminal differentiation [19,79]. This transcriptional basis of this phenomenon was further supported in our study by observation that genes associated with cell differentiation were highly expressed in female worms (S4 Table). Examples of such genes include G2/mitotic-specific cyclin-B3 (AY809873), Polo-like kinase (FN317236), abnormal spindle-like microcephaly-associated protein (AY812148), cyclindependent kinase 1 (AY815214), cell division cycle 20 (fizzy)-related protein (AY223249), Cyclin-T2 (FN317410) and regulator of chromosome condensation (AY810273). Further examples include those associated with cell cycle arrest in response to DNA damage or spindle abnormalities, including cell cycle checkpoint control protein RAD9B (AY812096), mitotic spindle assembly checkpoint protein MAD2A (AY814258), checkpoint protein HUS1 (AY813370), serine/threonine-protein kinase chk2 (FN313971), mitotic checkpoint serine/threonine-protein kinase BUB1 beta (AY808857). It has been shown that apoptosis is an important cellular process in schistosomes [80]. The activities of caspase-3 and -7, both central proteolytic enzymes involved in this process, were active across different developmental stages of S. japonicum, with a peak expression in the schistosomula 14 days p.i. [80]. Here, the expression of apoptosis-related genes, caspase 7 (AY813428), programmed cell death protein 2 (AY814013) and 4 (AY814519), as well as serine/threonine-protein kinase pim-1 (FN317924), was more extensive in female worms, which may represent an instinct response to get rid of damaged cells.
Genes involved in DNA synthesis and genome fidelity and stability are more up-regulated in female adult worms Use of an autoradiographic method has shown that the pairing of male and female parasites impacts on the DNA synthesis in females, but not in males [38]. DNA synthesis may present an important cellular process as a consequence of vitellocyte differentiation. Here, we found a set of genes involved in DNA replication processes, were preferentially expressed in female worms (S4 Table). These included DNA replication licensing factor mcm2 (AY815400), mcm4 (AY914892), and mcm7-A (AY815974), DNA polymerase alpha subunit B (AY811018), Origin recognition complex subunit 4 (AY812655), DNA replication factor Cdt1 (FN313910) and DNA replication complex GINS protein PSF2 (AY815313), were all up-regulated in females. In addition, transcripts for a number of DNA damage repair related genes (i.e., UV excision repair protein RAD23 (AY811322), DNA repair protein RAD51 (AY812723), DNA polymerase epsilon subunit 2 (AY815035), and DNA mismatch repair protein msh2 (AY814227) and msh6 (AY811972)) as well as genes encoding a chromosome transmission fidelity protein (AY810243) and a mini-chromosome maintenance complex-binding protein (AY810626) were also enriched in female parasites. These observations potentially reflect the need for repairing DNA damage caused by oxygen radicals released during the process of hemoglobin digestion, and the guarantee of chromosomal fidelity during vitellocyte mitosis and/or egg embryonic development. Furthermore, the gene encoding Argonaute 2, a protein which binds small interfering RNA, was found over-expressed in adult females, which is consistent with the results of a previous study [46]. The expression of the Ago2 ortholog in S. mansoni was observed predominantly in the gonads (particularly in the posterior ovary) [81], and SjAgo2 has been shown to play a vital role in germline cell maintenance via suppression of the activity of transposable elements (TEs) [46].

Genes involved in glycosylation are more up-regulated in female adult worms
Glycosylation in schistosomes is a complex process which plays an essential role in host-pathogen interplay, particularly in terms of immune evasion and modulation [82]. A comprehensive glycomic analysis has revealed that the dominant N-glycans structure dynamically changes during the development of S mansoni. For example, N-glycans with Galβ1-4 (Fucα1-3) GlcNAc (LeX) and core-xylose motifs are abundant in cercariae, but are lost rapidly after entry to the mammalian host, while GalNAcβ1-4GlcNAc (LDN)-motifs gradually became predominant during the transition of schistosomula to adult worms. Further, fucosylated motifenriched N-glycans are presented during egg development [83]. Also, it has been shown that tri-antennary type glycans are predominant in adult females compared with adult males [84].

Hypothetical genes in male and female adult worms
Within those gender-associated genes, a wide array of genes was annotated as hypothetical protein (155 (22.6%) and 95 (22.1%) in adult male-and female-biased expressed genes, respectively). These hypothetical genes may encode schistosome-specific proteins that lack homologous domains with other species, but limited attention has been paid on this gene set. Further research on these genes and their expressed products may further the discovery of new vaccine candidates and drug targets.
miRNA target prediction against gender-biased expressed genes miRNA profiles have been established across the different developmental stages and different sexes of S. japonicum and S. mansoni, and the potential function for some miRNAs have been suggested based on the profiling data, but their precise roles, such as how they regulate potential targets, remain elusive. Target prediction is an important pipeline in order to learn about the function of miRNAs. Previous miRNA target prediction has been carried out on schistosomes mainly within the 3' UTR of mRNAs [29,85]. However, some studies have shown that the target sites are not limited to the 3' UTR and can be located within the CDS and even the 5' UTR [86,87]. In addition to canonical target sites, non-canonical sites (i.e., "non-seed" sites [88,89]) have widespread biological functions, which undoubtedly increases prediction complexity. Here, we have focused on the mechanism whereby miRNAs can potentially regulate the expression of gender-associated genes.
Combining the algorithms of PITA [47] and RNAhybrid [48], putative miRNA target sites were predicted against the full length mRNA transcripts that exhibited gender differential expression. Half of these were predicted to contain miRNA target sites (Table 1, S9 and S10 Tables). On average, 1.6 putative miRNA target sites were predicted per individual gene, with most sites located within the CDS (~70%) of both male and female biased expressed genes ( Fig 5A). Work with HEK293 cells has shown that within the exonic crosslink-centered regions, 50% of sites correspond to the CDS compared with 46% to 3' UTRs [90]. In another study, it was shown that 41% and 40% miRNA binding sites were located in the 3' UTRs and CDS, respectively, in the human brain [91]. The over-concentration of binding sites in the CDS reported here may have been caused by the fact that the 3' UTR may be fractured in a group of S. japonicum mRNA transcripts (only about 10% gender-associated mRNA transcripts have a complete 3' UTR based on the presence of a poly(A) tail).The seed type plays has an important impact on miRNA regulatory function [92]. Here in our analysis, no mismatch was allowed in the seed site, and a single G:U wobble was only allowed for seed sizes of 7 and 8. The percentage of target sites was gradually decreased in size type of 6:0:0, 7:0:0 to 8:0:0, while the percentage of target sites was similar for seed types of 7:0:1 and 8:0:1 (Fig 5B). Individually, sja-let-7, sja-miR-1 sja-miR-7-5p, sja-miR-3479-5p sja-miR-190-5p, sja-miR-71 and sja-miR-71b-5p have the most putative sites within the sex-biased expressed genes (Fig 5C) of which sja-let-7, sja-miR-1 sja-miR-7-5p are male-biased miRNAs, while sja-miR-71b-5p is female- Sex-Biased Expressed Genes of Schistosoma japonicum biased [26]. In contrast, few target sites have been predicted for sja-miR-125b and sja-bantam, two miRNAs abundantly expressed in male and female worms, respectively, indicating that they may regulate non-gender-associated genes.

Conclusions
In this study, we present the most comprehensive transcriptomic perspective on sex-expressed genes in S. japonicum, which sheds further light on key biological and physiological features of the male and female parasites. Furthermore, we present a global view on how miRNAs potentially modulate the expression of gender-associated genes via miRNA target sites prediction. This study provides novel insights on schistosome conjugal biology, which may help in the discovery of new anti-fecundity vaccine candidates and drug targets against this persistent pathogen.