Survey of transcripts expressed by the invasive juvenile stage of the liver fluke Fasciola hepatica

Construction of a newly excysted juvenile F. hepatica cDNA library

In order to identify the genes expressed during the invasion process of the platyhelminth F. hepatica, we constructed a full length enriched cDNA library using a modified protocol based on selective amplification of capped polyadenylated RNA species. Since the starting parasite material was limiting, a modified size fractioning step of the products was introduced in order to improve the yield [28] (Additional File 1). More than four thousand reads were produced and analyzed using the Partigene pipeline [29]. Quality and vector trimming drastically reduced the starting 4319 ESTs to 1684 high quality sequences, mainly due to the presence of multimers of the adapters used in the generation of the libraries (see methods). This setback could be expected considering the minimal amount of starting material, and might be corrected using 5' blocked adapters in lower concentrations.

Comparison and validation of the FhNEJ ESTs with other databases

In order to establish if the obtained contig sequences correspond to validated transcripts, we compared them to different available databases, including ESTs from the adult F. hepatica stage, predicted coding sequences from selected organisms with complete genomes, and transcriptomes of other eukaryotes representing the main lineages in the metazoan diversity (Additional File 2). To compare the data obtained from the juvenile stage to the adult sequences, we retrieved and analyzed using the Partigene pipeline more than 10,000 EST reads from F. hepatica adult worms available at the Wellcome Trust Sanger Institute, obtaining 4089 contigs (1879 clusters and 2210 singletons), 58% of them showing significant blast hits (E value < 1e^-05) with publicly available databases (Additional File 3). These results are very similar to a recently reported analysis of the same dataset performed using a different pipeline [13].

More than half of the juvenile contigs (55.3%) were also found in adult ESTs (Figure 1). A set of 91 juvenile contigs (17.6%), also present in adults, showed no homology to sequences in other databases, suggesting that they might represent Fasciola specific transcripts expressed in diverse stages of the parasite life cycle. On the other hand, there are several juvenile contigs that are absent from the adult database, although represented in other organisms suggesting that they might represent stage specific transcripts (Figure 1). A set of 114 juvenile contigs (22,1%) were common to all other organisms searched indicating core eukaryotic functions such as ribosomal proteins and common enzymes. The absence of some of them from the adult dataset might suggest that the representation of the adult libraries is still partial. Interestingly, 64 contigs (12.4%) are shared only within flatworms, corresponding to conserved uncharacterized transcripts that might be relevant to parasitism. Also 56 contigs (10.9%) are shared only within metazoans and absent in the non metazoan choanoflagellate Monosiga brevicollis, suggesting that they represent metazoan innovations.

To further characterize conservation patterns between different metazoan lineages, we analyzed the distribution of tblastx hits by three-way comparisons using the Simitri program [31]. As expected, the F. hepatica predicted genes are more similar to homologues from other trematodes rather than cestodes and turbellaria, and to all flatworms rather than other protostomes, supporting the monophylectic origin of flatworms (Figure 2A, C). Consistent with the reports from the schistosomes genomes, we detected slightly more shared genes (being them also more similar) with the complete genomes of vertebrates than with insects and nematodes [18, 19]. These results further support the idea that ancient genomes were gene rich, and that lineage specific gene gain and loss events were frequent during metazoan evolution, particularly within the ecdysozoans [32]. While the relevance of genes shared between trematodes and their hosts has been highlighted, since they may be crucial for parasite adaptation to the host [33], the inverse situation (genes present in the parasite but absent in their hosts), might provide relevant candidates for anti-parasitic intervention.

Additionally, since we included in the study partial genomes from other lophotrochozoans (annelids and mollusks) we were able to compare the Fasciola dataset to these organisms and other phyla. This is relevant since flatworm position in modern phylogeny is still debated, being placed either within or as sister group of the lophotrochozoa [34–36]. The conserved set of liver fluke genes is almost equally distant from cnidarians, mollusks, and annelids, but slightly closer to the two lophotrochozoans than the model ecdysozoans or vertebrates (Figure 2B, D, and Additional File 4). The trend in this (and all other comparisons performed) were maintained when including the 4089 F. hepatica adult contigs suggesting that the effects observed might not be due to sampling bias (data not shown). The comparisons here presented are consistent with the placement of flatworms basal to the lophotrochozoans.

Compositional characteristics of F. hepatica predicted proteins

The average G+C content of the F. hepatica ESTs (both juvenile or adult) was 45%, a value substantially higher than in S. mansoni and S. japonicum (34%) [37]. Since variation in GC content can result in skewed codon usage [38], we analyzed the frequency of codons and amino acids of the predicted protein coding sequences in all F. hepatica available assemble ESTs (NEJ and adult stage), and compared it to those observed in other trematodes. As indicated in Figure 3A, there is a detectable difference in codon frequency, between the schistosomes and the other trematodes (including F. hepatica). Schistosomes prefer the most AU rich codon of each synonymous family, and are also strongly biased against C or G in the third codon position confirming early predictions obtained with limited gene sets [39]. More striking is the fact that significant differences were also found at the amino acid level, where schistosomes uses less Arg, Ala and Gly, and are enriched in Asn, Ile and Ser (Figure 3B). In a recent paper the tRNA complement of S. mansoni and S. japonicum is analyzed, but no significant correlation between tRNA copy number with the overall codon usage were found in any of the species [40]. The biological and evolutionary significance of the differences here observed is not clear, and deserves further consideration. In any case, these results raise the question that schistosomes might represent a more divergent than expected model for other trematodes.

Gene Ontology classification and functional annotation

Gene Ontology (GO) provides a useful way of classifying and annotating sequence information. Our analysis of the F. hepatica juvenile dataset showed up to 179 NEJ contigs with GO assignment. The molecular function classification showed a predominance of the binding category overlapping with almost all other categorizations, followed by enzymes (catalytic activity) and structural components. The discrimination within the binding class showed three main divisions of similar relevance, two overlapping with enzymes and ribosomal proteins and a set identified as protein and DNA binding associated with regulatory functions (Figure 4A). The more represented biological process categories were linked with metabolism, regulation and development (Figure 4B), showing a consistent assignment of GO cellular components (data not shown). Functional annotation of predicted proteins showed a general representation of the diverse biological functions. Proteases and antioxidant enzymes should be highlighted since they have long been under scrutiny for their putative involvement in invasion and immune evasion processes [9, 10, 41–49]. Novel proteins included ribosomal proteins (Additional File 5) several factors associated with protein and gene expression, cell signaling and apoptosis, as well as orthologues of candidate antigens that induce protection against other helminthiasis. They include tetraspanin-like protein [50], a membrane spanning protein located at the tegument of S. mansoni, Sm22.6 tegument antigen [51], and venom allergen-like (VAL) proteins, a candidate vaccine antigen against Necator americanus and Ancylostoma caninum[52–55].

Relevant molecules for parasitism

Secreted cathepsins were among the more represented transcripts in juvenile ESTs, and also in the adult dataset (Table 2, Additional File 6). A more detailed analysis of these transcripts showed that different isoforms are are being expressed by the invading and adult stage. While cathepsins L3, L4 and L6 are detected in the juvenile ESTs, they are absent from the much larger adult dataset (with the exception of cathepsins L4). Proteomic analysis have shown that cathepsins L1 and L2 are clearly predominant in adults, in agreement with the relative abundance of their transcripts in the adult EST database [13] (Additional File 6), and it has been proposed that the repertoire of cathepsin Ls gradually change from those expressed in juveniles to a different set characteristic of the adults worms [13–15]. Interestingly, it has recently been reported that the juvenile predominant cathepsin L3 has a strong collagenase activity, that might result essential for the invasion process [61], while the "adult" cathepsin L1 is involved in hemoglobin degradation [14].

We found evidence that within the less characterized cathepsin B gene family a similar phenomenon might be taking place. The cathepsin B forms that appear as frequent in juveniles are quite distinct to the cathepsin B transcripts found in the adult stage dataset (Additional Files 7 and 8), suggesting that they might also be functionally distinct; cathepsin B1 functions as a digestive enzyme in the juvenile gut [62].

Sequences encoding detoxifying enzymes like thireodoxin peroxidase (TPx), superoxide dismutase (SOD), thioredoxin 2, glutation S-transferases and a novel glutathione peroxidase not previously reported in F. hepatica were also found in juveniles, stressing their relevance for immune evasion [72]. In flatworms thioredoxin and glutathione peroxidases are the main enzymes involved in detoxifying reactive oxygen species produced by host immune effector cells [73, 74].

Secreted and surface proteins that may modulate host interactions are considered as relevant targets for vaccine or anti-parasitic drug design [75]. SignalP analysis identified putative signal peptides in 60 NEJ predicted proteins, while 52 had an N-terminal signal anchor peptide. Several putative secreted proteins were novel (with no significant hits) or conserved only in trematodes but not detected in other taxa. Some of these transcripts were among the more represented ESTs in juveniles (Table 2). The repeated detection of these transcripts in partial datasets from diverse trematodes support the notion that they are truly highly expressed genes in trematodes, and may be important mediators for parasitism. We selected Contig FHC00023, a predicted secreted protein of unknown function that is the most frequent in the juvenile ESTs with no homologies outside F. hepatica for further analysis. By real time PCR we found that this transcript is predominantly expressed in the invasive stage confirming the in silico observation (Additional File 9 panel C). The putative ORF is characterized by repeated Ser and Thr residues predicted to be glycosilated, and in further analysis showed faint homology with mucins. Parasite-specific proteins (with no counterparts in vertebrates) like these are ideal targets for development of therapeutic agents since they would have no cross-reactivity with host molecules. The elucidation of the function of these proteins is an important task. The growing availability of functional genomics tools like RNA interference in F. hepatica and model trematodes [76–78] offers some hope this can be accomplished.

Parasites

Fasciola hepatica metacercariae were obtained in our laboratory from experimentally infected Lymnaea viatrix snails and maintained encysted on 0.4% carboxymethyl cellulose until use. Excystment of metacercariae was performed as described previously [15]. Briefly, metacercariae were placed in a 100-μm filter and incubated 5 min with 1% sodium hypochlorite, washed exhaustively with PBS and incubated at 39°C for up to 3 hours in a solution prepared by mixing equal volumes of A (0.4% sodium taurocholate, 120 mM NaHCO₃, 140 mM NaCl, pH 8.0 and B (50 mM HCl, 33 mM L-cysteine). The emerging NEJs were collected in a 20 μm-filter with RPMI-1640 medium and used for RNA extraction.

RNA extraction, ligation of RNA adaptors and cDNA synthesis

Total RNA from 1200 NEJs was prepared using the Micro to Midi RNA Extraction Kit (Invitrogen), according to the manufacturer's protocol. Two hundred nanograms of total RNA were used for cDNA synthesis using the protocol described [28]. Briefly, the non capped RNAs were dephosphorylated, and the complete mRNAs were decapped by a pyrophosphatase treatment, and later ligated to the GeneRacer RNA oligo to introduce a 5' priming site in complete mRNAs. After this treatment, first strand synthesis was performed with the reverse transcriptase Superscript III (Invitrogen) using GeneRacer oligo-dT primer (5'GCTGTCAACGATACGCTAC GTAACGGCATGACAGTG(T)₁₈3').

Preparation of NEJ cDNA libraries

Amplification of full-length cDNAs was performed by PCR using universal forward (GeneRacer 5'Nested: 5'GGACACTGACATGGACTGAAGGAGTA3') and reverse primers (GenerRacer 3'Nested: 5'CGCTACGTAACGGCATGACAGTG3') provided by the GeneRacer kit. PCR was carried out for 30 cycles (94°C, 45 sec; 68°C 45 sec; 72°C 5 min) using Hot Start Taq DNA polymerase (Fermentas). PCR products were size fractionated in three subpopulations (300-800 bp, 800-2000 bp and >2000 bp) by excision from 1% TBE agarose gels, purified with QIAquick Gel Extraction kit (QIAGEN), ligated to pCR4-TA cloning vectors (Invitrogen), electroporated into One Shot TOP 10 Electrocompetent E. coli (Invitrogen), and plated on LB Amp/X-Gal. Recombinant clones from the libraries were randomly picked, grown in Circle Growth medium and stocked at -80°C in 96-well plates in 30% glycerol.

DNA sequencing and Bioinformatics analysis

Clones were cultured in 96 well plates with Circle Growth media and plasmid DNA was purified by alkaline lysis in 96 well plates. DNA was sequenced with M13 reverse primers using the Dyenamic ET Dye Terminator cycle sequencing kit for MegaBace DNA Analysis Systems (GE Healthcare Life Sciences) according to the manufacturer's instructions.

The sequence reads obtained were processed and analyzed using the Partigene pipeline [29]. Briefly, Trace2dbest [79] processed the chromatograms removing low quality (Phred <15, <150 bp) and vector sequences, and the resulting preprocessed ESTs were assembled in a two-step process carried out by CLOBB [80] and Phrap programs; the resulting contigs and singletons were compared to a set of databases maintained locally (listed in Supplementary Table 1) using tBLASTX and BLASTX. Functional categories were analyzed using annot8r [81]. Signal sequence prediction was performed using SignalP3.0 program [82]. Prediction of trans-membrane domains were conducted using TMHMM software [83]. Blast results comparisons were performed with Simitri [31]. All the available ESTs reads from adult stage of F. hepatica available at the Wellcome Trust Sanger Institute http://www.sanger.ac.uk/Projects/Helminths were downloaded and processed with the same pipeline. The juvenile sequences here obtained were deposited at the dbEST with the accessions GT740211 to GT741887.

Codon usage and amino acid frequencies

For F. hepatica adult and juvenile stages, Echinostoma paraensei, Opistorchis viverrini and Clonorchis sinensis, open reading frames were predicted from the assembled EST data using EMBOSS bioinformatics suite [84]. The longest ORF from each EST was retrieved and the predicted protein sequence was blasted against the NCBI nr databank. ORFs with significant hits were kept for further analysis. For S. japonicum and S. mansoni predicted coding regions obtained through the respective genome projects were analyzed. Codon and amino acid usage was calculated using GCUA: General Codon Usage Analysis tool [85]

Real time PCR

Real time PCR experiments were carried out in an Applied Biosystems 7500 Real time PCR System. Ten microliters of different dilutions of cDNA of NEJ and adult parasites were amplified using 0.2 μM each specific primers, 1.5 mM MgCl₂, 25 uM dNTPs, 0.25 U Platinun Taq DNA polymerase (Invitrogen), 1× SYBR Green, 1× PCR buffer in a 20 uL volume reaction. Primers sequences are β-actin (Forward 5'-GTGTTGGATTCTGGTGATGGTGTC-3' and Reverse 5'-CAATTTCTCCTTGAT GTCTCG-3'), FHC00023 (Forward 5'-ATGGTGCGAACGCTAAG-3' and Reverse 5'-GAAGAACGCAACGCCGAAGA-3'), Legumain 1 (Forward 5'-CAAGGATGTTTATGAAGGG-3'and Reverse 5'-TGCTTTGTTCATGCTGGC-3') Legumain 3 (Forward 5'-AGCAGACAAAACCCTTATCGT-3'and Reverse 5'-GGAATAATAGTAGGCGACGTG-3'). Reactions were performed in triplicate using the following PCR amplification conditions, 1 cycle (94°C, 5 min), 40 cycles (94°C, 15 seg; 60°C, 10 seg; 72°C, 15 seg). All results were analyzed using the 2^-ΔΔCt method and β-actin as internal control group [86].

Note added in proof

Recently a separate study describing the generation of more than 500,000 sequences from an adult cDNA library using 454 sequencing was published [87]. However, at the time of writing these sequences have not been made publically accessible and hence a comparative analysis of this dataset was not possible.