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

Abstract

Background: 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.

Results: NGS of the whole transcriptome resulted in 19,534,766 sequence reads, yielding a total of 40,260 transcripts (N₅₀ = 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).

Conclusions: 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.

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⁻⁵) 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⁻⁵), 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.

Figure 2

Nonsynonymous/synonymous substitution rate of the liver fluke transcripts. Ka/Ks ratio distribution (%) of a total of 12,394 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).

Figure 3

Relative abundances of the liver fluke transcripts in various biological functions based on manual inspection of protein motifs. The relative abundances (%) of all the functionally categorised transcripts within whole transcriptome (a), VR group and its PSR and PDR subgroups (b) in various functions are demonstrated. Nucleic acid binding/transcription and signaling mechanisms were highly pronounced functions of the whole transcriptome, with relative transcript abundances of 20.91% and 16.66%, respectively (Figure 3a). 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.

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).

Figure 5

Relative abundances of the virulence and immunomodulation related transcripts determined by the transcript subgroups. The Venn diagram shows percentage distribution of the virulence and immunomodulation related transcripts (VIRs, n = 71), determined by the exclusive homology with the pathogen databases (PDR subgroup), observation of the signs of positive selection (PSR subgroup) and/or cytokine signaling relation (CSR subgroup) among the immunomodulation categorised transcripts. This diagram indicates that VIRs showing the signs of positive selection (Ka/Ks > 1; in PSR subgroup) are evolved towards cytokine signaling. The Venn diagram was drawn using eulerAPE v3 (http://www.eulerdiagrams.org/eulerAPE/) [118].

Figure 6

Distribution of putative virulence and immunomodulation-related orthologous genes by the immunomodulation property. The number of the virulence and immunomodulation-related transcripts and the corresponding putative orthologous genes based on the predicted immunomodulation property are demonstrated. The number of the putative orthologous genes showing MHC I homology (n = 14) was the highest, which was followed by the others showing sequence homologies with TGF-β receptor (n = 11), T-cell or Toll/interleukin-1 receptors (n = 5 for both receptors), and the relation with TGF-β stimulation (n = 4).