Comparative transcriptome profiling approach to glean virulence and immunomodulation-related genes of Fasciola hepatica

Fasciola hepatica causes chronic liver disease, fasciolosis, leading to significant losses in the livestock economy and concerns for human health in many countries. The identification of F. hepatica genes involved in the parasite’s virulence through modulation of host immune system is utmost important to comprehend evasion mechanisms of the parasite and develop more effective strategies against fasciolosis. In this study, to identify the parasite’s putative virulence genes which are associated with host immunomodulation, we explored whole transcriptome of an adult F. hepatica using current transcriptome profiling approaches integrated with detailed in silico analyses. In brief, the comparison of the parasite transcripts with the specialised public databases containing sequence data of non-parasitic organisms (Dugesiidae species and Caenorhabditis elegans) or of numerous pathogens and investigation of the sequences in terms of nucleotide evolution (directional selection) and cytokine signaling relation were conducted. NGS of the whole transcriptome resulted in 19,534,766 sequence reads, yielding a total of 40,260 transcripts (N50 = 522 bp). A number of the parasite transcripts (n = 1,671) were predicted to be virulence-related on the basis of the exclusive homology with the pathogen-associated data, positive selection or relationship with cytokine signaling. Of these, a group of the virulence-related genes (n = 62), not previously described, were found likely to be associated with immunomodulation based on in silico functional categorisation, showing significant sequence similarities with various immune receptors (i.e. MHC I class, TGF-β receptor, toll/interleukin-1 receptor, T-cell receptor, TNF receptor, and IL-18 receptor accessory protein), cytokines (i.e. TGF-β, interleukin-4/interleukin-13 and TNF-α), cluster of differentiations (e.g. CD48 and CD147) or molecules associated with other immunomodulatory mechanisms (such as regulation of macrophage activation). Some of the genes (n = 5) appeared to be under positive selection (Ka/Ks > 1), imitating proteins associated with cytokine signaling (through sequence homologies with thrombospondin type 1, toll/interleukin-1 receptor, TGF-β receptor and CD147). With a comparative transcriptome profiling approach, we have identified a number of potential immunomodulator genes of F. hepatica (n = 62), which are firstly described here, could be employed for the development of better strategies (including RNAi) in the battle against both zoonotically and economically important disease, fasciolosis.


Background
Fasciola hepatica, the liver fluke, is a digenetic trematode helminth, causing highly damaging hepatobiliary disease (fasciolosis) in mammalians including economically important ruminants (such as cattle and sheep) and humans [1]. Fasciolosis results in a significant economic loss in livestock industry worldwide and more cases have been reported in humans in different countries [1][2][3][4]. The disease is currently treated with anti-helminthics (such as triclabendazole), but the observed anti-helminthic drug resistance [5] necessitates more effective strategies for the treatment and/or the prevention of fasciolosis.
For the elucidation of infection mechanisms of the liver fluke and the development of better strategies in dealing with fasciolosis, the most critical step is the identification and characterisation of the genes which are important for the establishment of parasitism. The transcriptome of F. hepatica was reported in a previous publication where a previous sequencing platform (454 NGS, Roche) was employed for the purpose of describing general biological characteristics of the parasite [6]. However, the sequence data from that study is likely to be still further from encapturing the entirety of the transcriptome profile of F. hepatica and the virulence factors of the parasite, particularly those related to host immunomodulation, are worthwhile for additional investigation with current NGS platforms (such as HiSeq 2000, Illumina). The current NGS technologies produce greater transcriptomic data in comparison with the previous sequencing systems and increase the chance of detection of parasite transcripts which are expressed at lower levels (relative to housekeeping transcripts) but with significant importance in immune evasion.
One of the most promising approaches for determining the virulence factors of parasitic organisms is in silico comparison of parasites' transcripts with publicly available data [7][8][9][10]. A recent data resource, the helminth secretome database (HSD) [11,12], contains a broad repertoire of excretory/secretory (ES) protein sequences (n > 18,000) of various parasitic helminths including trematodes, cestodes and nematodes. ES proteins of endoparasites are, in general, thought to play vital roles in the establishment of infection [11][12][13][14]. Additionally, Vaccine Investigation and Online Information Network (Violin; http://www.violinet. org), provides gene and protein sequences that are affiliated with infection mechanisms of various micro (virus, bacteria, protozoon) and macro (helminth) infectious organisms [15]. Compared to other resources providing whole genome data of numerous pathogens, both HSD and Violin contain filtered data better suiting to glean pathogen-related molecules of infectious organisms. However, the major caveat of the specialised databases is the presence of insufficiently refined data (such as house keeping genes/proteins in particularly HSD database) that are related with regular physiological events in both parasitic and non-parasitic organisms, but not indeed linked to virulence.
Parasitism genes which are not directly involved in virulence, but rather associated with regular physiological mechanisms, could be uncovered by sequence homology with taxonomically similar, free-living (non-parasitic) organisms [9]. Very recently, a large nucleotide sequence collection from free-living/non-infectious trematodes (taxonomically close to F. hepatica) of Dugesiidae family (including Dugesia sp. and Schmidtea sp.) has become publicly available at DNA Data Bank of Japan (DDBJ; www.ddbj.nig.ac.jp). Additionally, a comprehensive data for a well studied free-living model nematode, Caenorhabditis elegans, is freely accessible from a regularly updated resource, WormBase (http://www.wormbase.org).
The data of non-parasitic organisms from current resources have been useful tool to investigate the genes that are under directional selection through the comparative analysis of nucleotide diversity by assessing nonsynonymous/synonymous (Ka/Ks) substitution rates [16][17][18]. It is a well known fact that parasitism related genes important for the evasion of defensive systems of host and the corresponding genes in the host are under constant selective pressure favoring nucleotide subsitutions [18]. For example, the strong influence of directional selection in the evolution of avirulence genes to gain immunomodulatory properties has been clearly reported in multiple studies [19][20][21]. Furthermore, lineage-specific sub-or neofunctionalisation of genes which are vital for the establishment and maintenance of parasitism could be identified through comparative genomics [21][22][23].
Apart from the analysis of sequence homology with infectious and non-infectious organisms, the virulence genes associated with modulation of host immune responses (such as cytokine signaling) can be identified by the manual inspection of encoded protein motifs in public databases [24,25].
To date, in vivo and in vitro studies have shown that the liver fluke modulates host immune responses for enhancing its virulence [17,18,26]; however, which genes of the parasite imitate the components of host immune system have not yet been elucidated in detail.
The main purpose of this study was to glean virulence and immunomodulatory F. hepatica genes through comparative transcriptome profiling with the transcriptomes of non-parasitic related organisms by focusing on the genes which are evolved in lineage-specific manner, under positive selection and show similar motifs of host immune system genes involved in cytokine signaling.

Results
Transcriptome profile of F. hepatica Workflow illustrating the experimental steps of the study is demonstrated in Figure 1. From a total of 19,534,766 sequence reads, generated by the sequencing instrument (HiSeq 2000, Illumina) with paired-end 2X 100 bp reading, 81,090 contigs (contig N 50 = 377) were de novo assembled, of which, 40,260 transcripts were annotated with blast searches (blastx and blastn/tblastx) as described in Additional file 1. The obtained base number in this study was approximately 12.5 times higher than that reported in the previous transcriptome study of F. hepatica [6]. The transcript N 50 was 522 bp and the length of a total of 7,861 transcripts was equal or greater than the observed N 50 length in this work. In the present study, F. hepatica G + C content was 48.01%, which was similar to that reported in the previous related studies of F. hepatica (44,5%, 47.0 ± 14.1%) [6,27]. The identified transcript sequences in our study corresponded to a total of 28,142 unique accession numbers [n = 24,243 (NCBI related), n = 3,899 (GeneDB, CHGC, and SchistoDB related)].
The identity of species was genetically confirmed by the presence of a F. hepatica transcript (showing significant similarity with previously known species-specific heat shock protein 70 of F. hepatica [28], #ABS52704.1; E value = 8.00 −21 ), in addition to the morphological identification of the isolated parasite. The amino acid sequence for this protein at the correct frame of the transcript sequence was valine (V-599) as in F. hepatica, but not leucine (L-599) as in F. gigantica [28], confirming the species specificity of the isolated parasite in this study. In terms of drug resistance, the isolated parasite in the present study was found to be susceptible to albendazole (an anthelmintic benzimidazole drug). This was extrapolated by the comparative analysis of the translated aminoacid sequences of a F. hepatica transcript (annotated with tubulin beta-2 of F. hepatica; #CAP72050.1; E value = 0) with the drug susceptibility associated amino acid residues (N-165; F-167; E-198; F-200; R-241) of tubulin beta-2 of F. hepatica [29].
Nonsynonymous/synonymous substitution rate of F. hepatica transcripts A total of 16,832 orthologous pairs (E value < 10 −3 ) could be subjected to nonsynonymous/synonymous substitution rate analysis (13,288 F. hepatica transcripts showing Figure 1 Experimental design. The transcriptome study is consisted of wet lab procedures and data analysis processes. The parasite transcript sequences were analysed in terms of positive selection, cytokine signaling and pathogen database homology. The parasite sequences showing the signs of positive selection (Ka/Ks > 1) were termed PSRs (positive selection related transcripts). The cytokine signaling relation for the inferred parasite proteins were detected by the protein motif information and the related transcripts were termed CSRs (cytokine signaling related transcripts). Non-virulence-related transcripts (NVTs) were determined by the sequence homology of the parasite transcripts with the data of the non-pathogen associated databases (DDBJ, WormBase; E value < 10 −5 ) by excluding the non-parasite organism homologous PSRs and CSRs. The parasite transcripts showing the exclusive sequence homology with the pathogen-associated databases (HSD and Violin) (E value < 10 −5 ), but not being NVTs, were termed PRDs (pathogen database related transcripts). The parasite transcripts including PSRs, CSRs and PDRs were considered virulence-related transcripts (VRs). Immunomodulation categorisable VRs (PSRs/PDRs and all CSRs) were termed virulence and immunomodulation-related transcripts (VIRs). Asterisk (*) indicates that data from other resources (NCBI, UniProt, GeneDB, Gene Ontology and referred publications) were used for functional categorisation of some PSRs and PDRs.
homology with the sequences of Dugesiidae species and the remaining parasite transcripts showing homology with the sequences of C. elegans). Ka/Ks ratio was calculable for a total of 12,394 transcript pairs (Additional file 2). Ka/Ks analysis revealed that majority of the analysed F. hepatica transcripts (89.63%) have Ka/Ks < 0.5, hinting a purifiying selection against nonsynonymous changes as expected, a minority of the transcripts (3.37%) have Ka/Ks > 1, and a tiny portion of the transcripts (0.24%) have Ka/Ks = 1 (Figure 2). F. hepatica transcripts with Ka/ Ks > 1 were hereafter called transcripts under positive selection (termed PSRs within PSR subgroup). Only the small percentage of the orthologous transcripts were under positive selection, which confirms the hypothesis that some of the genes are diversing from the former versions to possess virulence capabilities as previously suggested [18].
The level of homology of F. hepatica transcripts in cytokine signaling The detailed analysis of InterProScan descriptions for all functionally categorised transcripts revealed that some of the protein motifs (family, domain or functional site), inferred from a group of transcript sequences (n = 35), suggesting possible involvement of them in cytokine signaling (named CSRs under CSR subgroup).
Sequence homology of F. hepatica transcripts with the specialised secondary databases and virulence-related transcripts of the parasite Approximately half of the total liver fluke transcript number (i.e. 51.87%) showed sequence homology with nucleotide or protein sequences of non-parasitic organisms including Dugesiidae species and C. elegans (E value < 10 −5 ) (Table 1a). The level of sequential homology of F. hepatica transcripts with Dugesia sp. and Schmidtea sp. was slightly higher (~1.2%) than that with the transcripts of C. elegans (37.77 %), but this increased to 46.36% when considering the parasite transcripts sequentially homologous to any of Dugesiidae species (E value < 10 −5 ). The parasite transcripts which were similar to the sequences of Dugesiidae species or C. elegans (E value < 10 −5 ) but not showed the signs of positive selection (Ka/Ks > 1) and/or cytokine signaling relation, a total of 20,483 liver fluke sequences (50.88%), were named non-virulence-related transcripts (NVTs). Based on the transcript number, approximately a quarter of F. hepatica transcriptome (26.56%) showed homology with the sequences from helminth secretome database (HSD ; Table 1a), but only a minor part of these transcripts (3.1%; n = 1,251) were out of the category of NVTs that these transcripts were named HTs (E value < 10 −5 ; Table 1b). A smaller percentage of F. hepatica transcriptome showed homology with the data of Vaccine Investigation and Online Information Network (Violin), but only a minority of these transcripts (0.29%; n = 117) were not observed in the category of NVTs and these transcripts were termed VTs (E value < 10 −5 ; Table 1b). A number of non-NVTs (n = 23) were found sequentially homologous to both HSD and Violin databases (HTs/VTs; E value < 10 −5 ). All HTs, VTs and HTs/VTs were called PDR(s) [pathogen database related transcript(s)] under PDR subgroup.
Overall, 4.15 % of all identified F. hepatica transcripts was assumed to be virulence-related transcripts (ascribed as VRs under VR group) on the basis of the degree of homology with sequences from publicly available databases; i) pathogen-associated data (PDR subgroup; n = 1,391; E value < 10 −5 ), ii) observation of the signs of positive selection (PSR subgroup; n = 418; Ka/Ks > 1) and iii) cytokine signaling relation (CSR subgroup; n = 35; manual inspection), yielding a total of 1,671 transcripts (excluding transcripts commonly determined by the different analyses) (Additional file 3). A number of the transcripts with Ka/Ks > 1 (n = 169) were also observed in PDR subgroup (n = 141 for HTs, n = 5 for VTs, and n = 23 for HTs/VTs). Of the cytokine signaling related transcripts (CSRs), one transcript showed the signs of positive selection and sequence homology with the HSD database (#20661), three transcripts (#6733, #23314 and #64440) showed the signs of positive selection without sequential similarity of any pathogen related databases, and the others showed sequence homology with the non-parasitic databases (n = 21). A number of CSRs (n = 10) showed non-similarity with any specialised databases used in this study. Dugesiidae or C. elegans orthologous F. hepatica transcripts is shown. The majority of the transcripts [89.63% (87.90% at P < 0.05)] were found to be associated with negative selection (Ka/Ks < 1). A group of transcripts [6.75% (1.88% at P < 0.05)] could be under directional selection (0.5 < Ka/Ks < 1). Another group of transcripts (0.24%) were related to neutral selection. Only a small portion of the transcripts (3.37%, n = 418) showed the signs of positive selection (Ka/Ks > 1) and, of this, a minority (0.64%, n = 79) reached statistical significance (P < 0.05).

Protein function profile for whole transcriptome and VR group
The InterProScan search for whole transcriptome of the parasite revealed a total of 5,089 unique accession number indicating protein families, domains or functional sites. Protein sequences inferred from a total of 20,160 transcripts (50.07% of all the annotated transcripts) were categorised in various functional groups by considering the biological functions of the parasite in previous publications [24,25,27]. The majority of the protein sequences (93.23%) were functionally categorisable with InterProScan information and the rest (particularly virulencerelated transcripts) were classifiable on the basis of the information from the other resources such as Gene Ontology, UniProt, NCBI or referred publications (Additional file 4). The analysis indicated that the abundant transcripts were mostly involved in nucleic acid binding/ transcription (20.91%) and signaling (16.66%), and only Numerical values for the F. hepatica transcripts which show similarity with nucleotide/protein sequences of the pathogen-associated and non-pathogen-associated databases (a) and the numbers of the virulence-related F. hepatica transcripts (VRs) (b) are listed. A total of 18,663 liver fluke transcripts (46.36%) and an additional 2,218 transcripts (5.51%) showed homology with sequences of Dugesiidae species and C. elegans, respectively, yielding a total of 20,881 transcripts (51.87%) (E value < 10 −5 ) (a). Overall, 20,483 liver fluke transcripts (50.88%) (orthologous to the non-parasitic organisms, E value < 10 −5 ) were determined to be non-virulencerelated (NVTs) after excluding the transcripts showing the signs of positive selection (Ka/Ks > 1) (n = 377) and the relationship with cytokine signaling (n = 21). Approximately a quarter of the total liver fluke transcripts (n = 10,694) showed sequence homology with the HSD data and 4.30% of the total liver fluke transcripts showed sequential similarity with the Violin data. A total of 1,251 transcripts (3.11%) and another set of 117 transcripts (0.29%) were exclusively homologous to the HSD (HTs) and Violin data (VTs), respectively but not identified in the category of NVTs (b). A small number of VRs (n = 23) were common for both HTs and VTs (HTs/VTs). In total, 1,391 VRs (HTs, VTs and HTs/VTs) were included in PDR subgroup (containing pathogen database related transcripts). A total of 246 VRs showing the signs of positive selection (Ka/Ks > 1) were observed in PSR group alone. Some of VRs (n = 31) were only identified by CSR subgroup. A number of PSRs (n = 168) were common for PDR subgroup. Three of VRs were detected by both PSR and CSR subgroups. Only one of VRs was determined by all the subgroups (PDR, PSR and CSR). Overall, a total of 1,671 transcripts, identified at least by one of the subgroups, were predicted to be virulence-related. The percentage values indicate the ratio of transcript number to total transcript number (n = 40,260). NVTs: Non-virulence-related transcripts.
2.08% of all the categorised transcripts (i.e. 419 transcripts) was associated with immunomodulation ( Figure 3a). Biological functions of the virulence-related transcripts under VR group, and PDR and PSR subgroups were found to be mostly related to nucleic acid binding/transcription, signaling and unknown mechanisms, respectively ( Figure 3b). Interestingly, the relative abundance of the transcripts with unknown mechanisms in PSR subgroup (26.79 %) was higher than those in PDR subgroup (19.12%), suggesting possible novel functions enhancing the virulence of the parasite (Figure 3b). The relative quantity of immunomodulation related transcripts within PDR subgroup (2.66 %) was higher more than twice, compared to those within PSR subgroup (1.20 %) (Figure 3b), likely hinting the importance of high diversity at lineage/genus level for gaining immunomodulation capability. The protein function of all the transcripts under CSR group was inherently related to immunomodulation.

Subcellular localisation profile for whole transcriptome and VR group
Subcellular localisation signals of the parasite protein sequences were mostly related to cytosol (21.68%), extracellular space (20.57%) and nucleus (19.15%) (Figure 4a). Similarly, most of the detected subcellular localisation signals of VRs were associated with the cytosol, extracellular and nucleus parts, in order ( Figure 4b). The extracellular localisation signal was more oftenly observed within VR group (20.42%) and PDR subgroup (20.80%) and CSR subgroup (21.80%), in comparison to PSR subgroup (17.81%). Further details about the subcellular localisation for the transcripts are shown individually in Additional file 5.
Virulence and immunomodulation-related transcripts and genes of F. hepatica Of VRs (n = 1,671), the immunomodulation categorised PDRs, PSRs, and all CSRs, a total of seventy-one transcripts, corresponding to 64 putative genes, were named virulence and immunomodulation-related transcripts (VIRs) under VIR set ( Table 2). Further details about the sequence characteristics including available InterProScan accesion number for VIRs are shown in Additional file 6. The majority of VIRs were specifically detectable through the level of homology with the pathogen related databases (PDR subgroup) (49.3%) and cytokine signaling relation (CSR subgroup) (43.7%) and the minor part of VIRs, observable in PDR subgroup (1.4%), CSR subgroup (4.7%) or both (1.4%), showed the signs of positive selection (PSR subgroup) ( Figure 5). ). This pattern was similar within VR group except the relative transcript abundance of unknown mechanisms was higher within this group (20.47%), PSR subgroup (26.79%) and PDR subgroup (19.12 %). The transcript proportions of amino acid/protein metabolism, catalytic activity, energy, neuron system and signaling within PSR subgroup were higher (around 1.5-5 times) than those within PDR subgroup. The relative transcript ratios for immunomodulation, ion/haem binding and proteolysis/protein degradation functions were higher (about 2-5 times) in PDR subgroup in comparison with PSR subgroup.
A number of VIRs (n = 15) showed sequential similarities with MHC I and an important part of VIRs indicated the relationship with TGF-β signaling based on sequence homologies with TGF-β, TGF-β receptor or other proteins associated with stimulation or inhibition of TGF-β ( Figure 6). The other sequence homologies were associated with various immunomodulatory molecules including T-cell receptor, toll/interleukin-1 receptor, TNF receptor, cluster of differentiations (i.e. CD48, CD59, CD 147), IL-18 receptor accessory protein, interleukin-4/interleukin-13, TNF-α, modulators of T-cell function, suppressors of cytokine signaling and of IKBKE 1, or molecules involved in other immunomodulation-related mechanisms (including IL-10 stimulation, leukocyte mediated cytotoxicity, proline binding or macrophage migration inhibition, toll like receptor 4 regulation; Figure 6). The majority of VIRs (n = 64) individually were potentially located at extracellular space (n = 63) or localised within plasma membrane (n = 1) ( Table 2), possibly indicating direct interactions of the parasite proteins with host immune system.
Gene ontology categories related to whole transcriptome, VR group/subgroups and VIR set Gene ontology categories (i.e. biological process, molecular function and cellular component) for whole transcriptome of the liver fluke, in comparison with VR group, its subgroups (CSR, PDR, PSR) and VIR set are shown in Table 3. The relative transcript abundances of the cellular (GO:0009987; 31.27%) and metabolic processes (GO:0008152; 26.68 %) were found to be higher, in comparison to other biological processes in the total transcriptome, and this pattern was similar for VR group and its PDR and PSR subgroups. However, the proportions of the transcripts in response to stimulus (GO:0050896; 21.82%) and immune system process (GO:0002376; 9.09%) were the highest within VIR set than those within VR group and its subgroups. Furthermore, VIR set has higher sequence proportions in molecular transducer activity (GO:0060089; 20%) and receptor activity (GO:0004872; 6.15%), compared to VR group and its subgroups except CSR subgroup. The relative abundance in the GO cellular compartmant category within VIR set was skewed to membrane (GO:0016020; 68.29%).
Biological pathways related to whole transcriptome, VR group/subgroups and VIR set A total of 96 different KEGG biological pathways (195 enzyme types) were determined in whole transcriptome of the parasite, in which the purine (map00230; 16.62%) and pyrimidine metabolisms (map00240; 7.93%) were the biological pathways with higher abundant transcript numbers, where the relative abundances of the other pathways, Figure 4 Subcellular localisation signals of the liver fluke transcripts. The relative amount (%) of the detected subcellular localisation signals within whole transcriptome (a), VR group and its subgroups including CSR subgroup, PDR subgroup and PSR subgroup (b) are demonstrated. The relative abundances for the cytosol, extracellular part and nucleus signals were found high levels within the total transcriptome, 21.68%, 20.57%, 19.15%, respectively. This profile was similarly found within VR. The relative percentages of the extracellular signal in CSR subgroup, PDR subgroup and VR group are higher (~2%) than those in PSR subgroup. The abundances of the signals of mitochondria&peroxisome, endoplasmic reticulum and plasma membrane were proportionally slightly higher (about 0.5-1.5%) in PSR subgroup than those within VR group and other subgroups (CSR and PDR subgroups).        Gene ontology categories (parental 2) for the whole transcriptome (WT) and VR group with its subgroups (CSR, PDR, PSR) and VIR set are shown. The relative transcript abundances of response to stimulus (GO:0050896; 21.82%) and immune system process (GO:0002376; 9.09%) were predominant in VIR set while the proportional transcript abundances for the cellular process (GO:0009987) and metabolic process (GO:0008152) were remarkably higher within WT, VR group, PDR and PSR subgroups (around 19-32%). The relative transcript abundances of molecular transducer (GO:0060089) and receptor activities (GO:0004872) appeared to be much higher within VIR set (20%, 6.15%, respectively) and CSR subgroup (23.21%, 7.14%, respectively), compared to the other subgroups, VR group and WT. The proportional transcript abundance for the membrane (GO:0016020) was predominant within VIR set (68.29%) and CSR subgroup (93.75%), but the cellular part (GO:0005623) (around 34.02-36.02%) was with the highest transcript ratio within VR group and its PDR and PSR subgroups and WT. Numerical values and brackets indicate the total number of transcript sequences and the relative transcript abundance (%) within each transcript classification, respectively.
in transcript number, were less than 5% (Additional file 7). The pathway with the highest transcript number was the purine (map00230; around 12-16%) metabolism within PSR and PDR subgroups; however, the relative transcript abundance for the pyrimidine metabolism (map00240) within PDR subgroup was approximately twice than that within PSR subgroup (Table 4). Aminobenzoate degradation (map00627), beta-Alanine metabolism (map00410), glycine, serine and threonine metabolism (map00260) were uniquely identified in PDR subgroup while butanoate metabolism (map00650) and pentose phosphate pathway (map00030) were found specific to PSR subgroup among the transcript subgroups (Table 4).

Discussion
In this study, we applied detailed in silico analyses to determine the virulence and immunomodulation-related genes of the liver fluke through a comparative assessment of the transcriptome profile with current NGS technology (HiSeq 2000, Illumina). The observed GO categories of F. hepatica transcripts in terms of biological processes were mostly concordant with the previous classification [6]. However, some GO categories, including signaling (GO:0023052; biological process), receptor activity (GO:0004872; molecular function) and membrane (GO:0016020; cellular component) were firstly described in the present study. This indicates overall more comprehensive coverage of total transcriptome profile of this parasite. Beside GO analysis approach with blast2GO, the manual inspection approach, including analysis of all the detected protein motifs, significantly increased the number of the categorised transcript into various functional terms that were previously used in the related studies [24,25,27].
Comparative transcriptome profile analysis of the parasite transcriptome with the sequences from non-parasitic organisms (i.e. Dugesiidae species and C. elegans) revealed that 51.87% of the obtained parasite transcripts shared significant homology, and the remaining (48.13%) were possibly be lineage-or genus-specific. A small part of the transcriptome (3.46%) showed the sequence homology with the pathogen related databases but not with the non-parasite related databases. The important strategy for picking out candidate transcripts with virulence was based on the identification of the transcripts with Ka/Ks > 1 with the assumption of the possible fast evolutionary pattern of parasitism associated genes [18]. This strategy has been successfully employed to identify virulence and immunomodulation-related genes of other parasitic organisms such as Plasmodium and Theileria species [30][31][32][33]. In our study, as a result of Ka/Ks analysis of 12,394 orthologous transcripts, a total of 418 F. hepatica transcripts were found to be with Ka/Ks > 1, hinting at a faster evolutionary rate because of likely involvement of these genes in the process of parasitism. More detailed analysis of the transcriptome with the motifs of proteins known to be associated with cytokine signaling was useful for the elucidation of other important genes from F. hepatica. The similar evolutionary characteristics of previously described virulence genes including cathepsin protease L1 and L2 [34,35], cathepsin protease B [36], thioredoxin/peroxiredoxin [13], glutathione S-transferase (sigma, mu, omega classes) [37][38][39], protein disulfide-isomerase [40] and 14-3-3 protein [41] boosts the authenticity of parasitism associated genes identified in this study. The presence of more transcripts likely encoding immunomodulatory proteins within the pathogen database related subgroup (PDR subgroup), in comparison to the positive selection related subgroup (PSR subgroup), suggests a possible role for nucleotide diversity at the lineage/ genus for the establishment of parasitism. Because of more frequent extracellular localisation of the transcripts within PDR subgroup and the cytokine signaling related subgroup (CSR subgroup) relative to PSR subgroup, it could be postulated that these gene products could directly interact with host proteins. Pyrimidine metabolism and other aminoacid related metabolisms such as aminobenzoate degradation (map00627), beta-Alanine metabolism (map00410), glycine, serine and threonine metabolism (map00260), specific to PDR subgroup, could hint a possible lineage/genus specific nucleotide diversification.
Through the comparative transcriptome analysis of the liver fluke sequences, a set of virulence-related transcripts (n = 71, corresponding to 64 putative genes) were found likely to possess immunomodulatory functional characteristics. To our knowledge, all the virulence and immunomodulatory genes of F. hepatica identified hitherto have not been reported, with the exception of two putative genes (i.e. #7694 and #32989) which are related to CD59 [42].
The proportion of transcripts with receptor activity (GO:0004872) was higher in VIR set, in comparison with VR group and its PDR and PSR subgroups, possibly denoting that VIRs are coevolved with host proteins due to direct interactions. The skewed proportionality in abundance of membrane (GO:0016020) in VIR set supports this further.
Putative parasite genes (n = 14) with sequence homologies to MHC class I receptor had the highest proportion of all identified parasitism genes. The predicted protein sequences of the related transcripts (i.e. #38312, #53490, #55117, #65009, #22528, #61208, #13771, #73151, #59522, #38573 and #58384) cover the region of alpha 1 and 2 domains, but not of alpha 3 domain of MHC class 1 molecule. We speculate that the parasite peptides may be presented to host immune cells (i.e. CD8 T-cells and NK cells) by alpha 1 and 2 domains, but the cytotoxic effects of immune cells could be potentially blocked because of the absence of alpha 3 domain. This proposed mechanism may result in suppression of subsequent pro-inflammatory responses including related CD4+ T-cell differentiation (possibly Th1) and cytotoxic cell killing mechanisms, known to be harmful to the liver fluke [51][52][53].
Another interesting finding was the presence of a number of putative genes with T-receptor homology for a number of putative genes (i.e. #27939, #34729, #44260, #65009 and #73578). The imitation of T-cell receptor by the parasite may interefere with the presentation of the parasite's antigen to T-cells and subsequent cellular and humoral responses. Taken together with the observed homologies related to T-cell receptor and the proteins of TGF-β signaling accentuates the importance of the stimulation of Th2 type responses for the success of parasitism as suggested by another study conducted for Trichuris suis [54].
The identification of a transcript sharing sequence similarity with CD147 (basigin) (i.e. #23314) was another noteworthy finding. CD147 is a membrane protein of suppressed T-reg cells and associated with negative regulation of T-reg associated cytokine signaling and T-cell activation (via impaired expression of IL-2 receptor α-chain CD25) [56,57]. We suggest that F. hepatica transcripts which coevolved with CD147 could be important for the regulation of T-reg associated responses. The other cluster of differentiation (CD) similarity was found to be associated with CD48 (i.e. #59632), stimulating various regulatory factors in B and T lymphocytes through binding molecules such as CD2 and CD244 (2B4) [58].
Two of F. hepatica transcripts, #10702 and #79120, were found to share high level of sequential similarities with macrophage inhibition factor and ctg4a (protein canopy homolog-related 3), respectively. Both proteins are known to be involved in the activation of macrophages during inflammation [59,60]. In relation to this, three genes, #42935/#76935, #72283, #3866, were associated with binding proline (GYF super family, # cl00072) [61,62] which is secreted by Th2 type immune responserelated alternatively activated macrophages promoting the development of fibrosis [63]. Taken together, the parasite proteins appear to be involved in the regulation of the macrophage activation and control of the development of excessive fibrotic tissue, in concordance with observations in clinical studies [64,65].
There are other identified transcripts which showed sequential similarities with a neutrophil protein (i.e. #5195) [66], molecules known to suppressors of cytokine signaling (i.e. #11419, #47252) [67] and of IKBKE1 [1IKK-epsilon and TBK1, influencing type I IFN production (#38525)] [68] or modulators of T-cell (i.e. #7920 and #70639) [69]. These genes could be other components of the immunomodulatory mechanism induced by the parasite, but which specific immune responses they stimulate can not be predicted with the available data.

Conclusions
By comparative analysis of total transcriptome of F. hepatica with publicly open databases, a number of putative genes (n = 62), which are potentially critical for virulence through immunomodulation or associated mechanisms and firstly described in this study. The majority of these genes appeared to be lineage-or genus-specific, suggesting a modus operandi through the enhancement of sequential diversity for genes encoding proteins which are likely to be at the frontline for the establishment of parasitism. In addition, the nucleotide diversity stemming from positive selectional pressure was found to be associated with cytokine signaling mechanisms by relying on the observed homologies with known genes such as toll/ interleukin-1 receptor, TGF-β receptor, CD147 and a S. mansoni orthologous protein containing thrombospondin (type 1) domain (associated with TGF-β stimulation). A significant percentage of the transcripts have a remarkable level of sequential similarity with host immune receptors and cytokines, which are known to be part of array of immunological responses through CD4+ T-helper differentiation, indicating modulation of host immune system via controlling cellular responses associated with the T-helper heterogeneity (T-reg, Th1, Th2 and Th17 in particular). In conclusion, the blockage of the effects of aforementioned parasite proteins with RNAi or other strategies (e.g. vaccine or drug) would be a good approach in the fight against fasciolosis. This may even be important in promoting the efficacy of other immunoprophylactic molecules which are experimentally tested in the prevention of fasciolosis. Apart from the dealing with fasciolosis, synthetic versions of the identified virulence and immunomodulatory genes reported in this study could be important to control undesired pro-inflammatory responses, by considering immunoregulatory effects of the parasite and current therapeutic approaches in the treatment of autoimmune defects (e.g. helminth therapy) [70][71][72][73]. Studies including the synthesis of the indicated genes in prokaryotic and eukaryotic systems and evaluating their immunological effects in vitro and in vivo are underway. The present approach can be used for other studies with similar purposes.

Parasites
Adult liver flukes were collected from the bile ducts of naturally infected cattle in an abottoir (Tuzla, İstanbul) and immediately placed in a 50 ml tube containing warm PBS (Biochrom) as previously described [74]. Intact liver flukes were washed with warm PBS and placed in a flask containing culture medium (DMEM, Invitrogen) and gentamycin (50 μg as a final concentration) (Invitrogen). The flask containing the parasites was kept in a suitable enviroment (37°C, 5% CO 2 ) for 2 hours for the regurgitation of all contents from the parasites' digestive tracts as previously described [6]. After the flukes were washed with PBS (37°C), each fluke was placed in a seperate cryogenic tube (#1620-2700, Seal-Rite), snap-frozen using liquid nitrogen and kept −80°C until use.

RNA extraction
Total RNA from whole body of each fluke was extracted by using RNeasy Protect Mini Kit (Qiagen) with an oncolumn DNAse step (Macherey-Nagel) as previously described [75]. Purity and quantity of the extracted total RNA were analysed using a spectrophotometer (NANODROP 1000, Thermo Scientific) and a fluorescence based system (Qubit 2.0 Fluorometer using Qubit RNA BR assay kit, Invitrogen), respectively [75]. Quality of RNA for each extraction was analysed by a microfluidics capillary based electrophoretic system (Agilent 2100 Bioanalyzer using Agilent RNA 6000 Nano Kit (Agilent Technologies) according to the manufacturer's recommendations except heat-denaturation which breaks 28S rRNA and prevent determination of RNA integrity number (RIN) [75]. Among several RNA extractions, the sample with the highest RIN number (RIN = 10) was used for sequencing.
Next generation RNA-sequencing (RNA-seq) of whole transcriptome of F. hepatica RNASeq library was prepared from 1.25 μg of total RNA with TruSeq RNA Sample Preparation kit (Illumina) according to the kit's user guide (Part#15026495 Ref. D). Briefly, mRNA was denaturated (65°C for 5 min, 4°C hold), eluted, fragmented and primed (elution 1; 80°C for 2 min, 25°C hold, elution 2-frag-prime; 94°C for 8 min, 4°C hold). Double strand (ds) complementary DNA (cDNA) was synthesised using a reverse transcriptase enzyme (SuperScript II Reverse Transcriptase, Invitrogen) and the other required reagents supplied by the kit. After the end repair, 3′ end adenylation and adapter ligation steps, cDNA fragments were enriched with PCR amplification as described by the user guide. All cDNA clean up steps were performed using Agencourt AMPure XP beads (Beckman Coulter) and a magnetic stand (Agencourt Bioscience Corporation). Quality and quantitative parameters of the library were determined by the Agilent High Sensitivity DNA kit (#5067-4626) and KAPA Library Quantification Kit (#KK4844, KAPA Biosystems) using the Agilent 2100 Bioanalyzer and a quantitative PCR system (iQ5, Biorad) according to the manufacturers' instructions, respectively. After ds DNA fragments were denatured with NaOH (0.05 N as final concentration), a bridge PCR amplification for the cluster generation from single-molecule DNA templates was performed on the inside surface of a flow cell (Illumina) by an automated instrument ( De novo assembly of sequence reads and annotation of the contiguous sequences Sequencing images were generated by HiSeq 2000 and the image analysis step was performed by RTA (Real Time Analysis) software (Illumina). The base calling step was performed with RTA or OLB (Off-Line Basecaller) softwares (Illumina). Cluster intensities and noise estimates were used in the analysis. The base sequence from each cluster, a confidence level for each base, and the filtering parameter (whether the read passes filtering) were given as the output for the base calling. The output files (with bcl extension) were converted to compressed fastq files by analysis software, CASAVA 1.8.2 (Illumina). The reads with quality score less than 33 were eliminated and the first 13 bases of each read were trimmed because of the insufficient quality. Reads with length less than 15 bases were removed before the assembly step. A bioinformatic programme, VELVET 1.2.08 [76], was used for the first step of the assembly and shortPaired run was applied where k-mer length was set to 31, insertion length was set to 400, expected coverage was set to 25, coverage cutoff was set to 2, and minimum contig length was set to 100. The output of this step was used as the input for the second step of the assembly process where OASES 0.2.28 [77] programme was applied. As similar to the previous process, insertion length with the value of 400 and coverage cutoff with the value of 2 were used for OASES 0.2.28 run. Contiguous sequences (contigs) were annotated with blast analyses (E value < 10 −5 ). Nucleotide sequences of the contigs were searched aganist publicly available nonredundant protein and nucleotide databases from NCBI using standalone blastx and blastn programs, respectively [24]. The contigs with E values (E value < 10 −5 ) in the blastn analysis were searched against the same database with standalone tblastx programme to predict the correct frame of the contigs. Nucleotide sequences of the remaining unannotated contigs were searched against Schistosoma mansoni database [n = 11,810 for protein (obtained from Martin Aslett, The Wellcome Trust Sanger Institute, United Kingdom), n = 11,912 for mRNA, downloaded from GeneDB (www.genedb.org [78]) on 20 January 2014], S. japonicum database [n = 12,657; v3, both protein and nucleotide, n = 13,469; v4, both protein and nucleotide, n = 17,401; cDNAs, downloaded from Chinese National Human Genome Center (CHGC) at Shanghai, The Schistosoma japonicum Genome Project (http://www.chgc.sh.cn/ japonicum/Resources.html) on 21 January 2014], S. haematobium and S.mansoni databases [ShaeEgypt; n = 13,073 for protein and nucleotide, SmanPuertoRico; n = 3,897for nucleotide/n = 3,896 for protein, downloaded from SchistoDB (http://SchistoDB.net) [79] on 23 January 2014]. All the blast results were extracted without cutoff parameters using Blast Parser (v1.2.6.14; http://geneproject.altervista.org/) and annotated contigs were termed transcripts. Nucleotide sequences of the identified transcripts were conceptually translated into amino acid sequences using Transeq (http://www.ebi.ac.uk/Tools/st/emboss_transeq/) [80] based on the blast matching frames. The transcript N 50 value was determined as previously described [81,82].

Investigation of sequence homology of F. hepatica transcripts with the specialised databases
To detect non-parasitism homology of the liver fluke transcripts, sequence data of Dugesiidae family and of Caenorhabditis elegans were used, because 1) Species of Dugesiidae family and C. elegans are known free-living (non-parasitic) organisms, 2) both Dugesiidae and Fasciolidae families taxonomically belong to same phylum, platyhelminths (flatworms; http://www.ncbi.nlm.nih.gov/), 3) a large number of nucleotide sequences of Dugesia sp. and of Schmidtea sp. in the Dugesiidae family is publicly available [DNA Data Bank of Japan (DDBJ); www.ddbj.nig.ac.jp], 4) C. elegans is a well studied free-living organism and a comprehensive sequence information for this organism is publicly available (WormBase; http:// www.wormbase.org). Nucleotide sequences of Dugesia sp. (n = 72,225) and of Schmidtea sp. (n = 82,784) were downloaded from the DDBJ resource on 10 and 11 February 2014, respectively. In addition, a small number of nucleotide sequences (n =125) for other organisms belonging to the Dugesiidae family, including Cura sp., Girardia sp., Neppia sp., Romankenkius sp. were downloaded from the same resource (11 February 2014). Protein coding nucleotide sequences (cds) and protein sequences of C. elegans (c_elegans.PRJNA13758.current.* and c_elegans.PRJNA13758.current_development.*, n = 26,769 and n = 26,983, respectively) were obtained from the ftp site of WormBase (http://www.wormbase.org). Nucleotide sequences of all the F. hepatica transcripts were searched against nucleotide sequences of the Dugesiidae species and protein sequences of C. elegans using the standalone tblastx and blastx, respectively (E value < 10 −5 ).
Detection of nonsynonymous/synonymous substitution rate of F. hepatica transcripts Protein coding sequences of F. hepatica transcripts were obtained using GETORF (http://emboss.bioinformatics.nl/cgi-bin/emboss/getorf ), with the consideration of the frame sense and longest sequence length (n > 30), and translated into amino acid sequences with Transeq (Jemboss; v1.5) [89]. The F. hepatica orthologous Dugesiidae sequences were determined with the blastx search based on the parasite's translated amino acid sequences (E value < 10 −3 ). Protein coding sequences of the Dugesiidae sequences were determined by GETORF with the same parameters and the sequences were trimmed from the ends to excise error letters (until the end letter, leaving the minimum sequence length of 18). The sequence alignment for each orthologous transcript was carried out with ClustalW (v2.1) [90] using ParaAT (Parallel Alignment and back-Translation; version 1.0) [91]. The ratio for the number of nonsynonymous substitutions per nonsynonymous site (Ka) to the number of synonymous substitutions per synonymous site (Ks) for the aligned transcripts were estimated using KaKs_Calculator (version 1.2) [92], with the MYN method (a modified version of the Yang-Nielsen algorithm [93,94]. For the remaining F. hepatica transcripts that are not orthologous to Dugesiidae or Ka/Ks ratio was not calculable, the Ka/Ks analysis was performed using the F. hepatica orthologous C. elegans sequences with the same approach. The P value for Ka/ Ks ratio was calculated by the Fisher's exact test by KaKs_Calculator and Ka/Ks ratio with the P value less than 0.05 was accepted statistically significant.

Functional categorisation
Translated amino acid sequences of all the F. hepatica transcripts were analysed with InterProScan 5.0 [95] using blast2GO [96] (version 2.7.0/2.7.1). All available applications in the InterProScan configuration were run, which were; blastProDom (scanning the families in the ProDom database [97]), FPrintScan (scanning the fingerprints in the PRINTS database [98]), HMMPIR (scanning the hidden markow models (HMMs) in the PIR Protein Sequence Database [99]), HMMPFAM (scanning the HMMs in the PFAM protein families database [100,101]), SMART (scanning the HMMs in the SMART domain/ domain families database [102]), HMMTigr (scanning the HMMs in the TIGRFAMs protein families database [103]), ProfileScan (scanning PROSITE profiles [104]), PaternScan (scanning PROSITE profiles with a new version of the PROSITE pattern search software [104,105]), HAMAP (scanning HAMAP profiles [106]), SuperFamily (scanning a library of profile HMMs that represent all proteins of known function [107]), SignalPHMM (predicting signal peptides and/or anchors [108]), TMHMM (predicting transmembrane helices in proteins [109]), HMMPanther (scanning the HMMs in the Panther database [110][111][112]), Gene3D (scanning a large collection of CATH protein domain assignments for ENSEMBL genomes and UniProt sequences [113]), Phobius (predicting combined transmembrane protein topology and signal peptide [114,115]), and Coils (predicting coiled coil regions in proteins [116]). The parasite transcripts were categorised in biological function based on the InterProScan information, considering the order of protein family, domain and functional site (conserved site, active site, binding site or repeat).

Cytokine signaling association
The InterProScan information for the observed protein motifs were manually inspected in terms of the relationship with cytokine signaling at the publicly available database of the European Bioinformatics Institute-InterPro (http://www.ebi.ac.uk/interpro/). The parasite transcripts that are potentially associated with cytokine signaling on the basis of protein family, domain or functional site were reported and termed CSRs (cytokine signaling transcripts).

Subcellular localisation analysis
Subcellular localisations of all the predicted F. hepatica protein sequences were analysed by WoLF PSORT (the value for the 'k used for kNN' was set to 32) [117].

Identification of virulence-and virulence and immunomodulation-related F. hepatica transcripts
The liver fluke transcripts that showed sequential homology with the non-parasite related databases but not showed signs of positive selection and/or cytokine signaling association were termed non-virulence-related transcripts (NVTs). The virulence-related transcripts (VRs) were predicted based on the following criteria; 1) observation of the exclusive homology with the data of the pathogen databases (PDRs), including transcripts sequentially homologous to HSD (termed HTs) and Violin (termed VTs) but not NVTs, 2) demonstration of the signs of positive selection (Ka/Ks > 1; PSRs), 3) detection of the predicted functional protein site which is related to cytokine signaling (proven by protein family, domain or funcitonal site information; CSRs). Some of the virulencerelated transcripts which could not be categorised with the InterProScan data were categorised using the information from the Gene Ontology, UniProt, NCBI and referred publications. The transcripts specific to immunomodulation category was determined on the basis of sequential identity level to known immunomodulatory proteins. The immunodulation categorised PDRs and PSRs and all CSRs were termed virulence and immunomodulationrelated transcript(s) [VIR(s)].