Anti-schistosomal intervention targets identified by lifecycle transcriptomic analyses

Abstract

Background: Novel methods to identify anthelmintic drug and vaccine targets are urgently needed, especially for those parasite species currently being controlled by singular, often limited strategies. A clearer understanding of the transcriptional components underpinning helminth development will enable identification of exploitable molecules essential for successful parasite/host interactions. Towards this end, we present a combinatorial, bioinformatics-led approach, employing both statistical and network analyses of transcriptomic data, for identifying new immunoprophylactic and therapeutic lead targets to combat schistosomiasis.

Methodology/principal findings: Utilisation of a Schistosoma mansoni oligonucleotide DNA microarray consisting of 37,632 elements enabled gene expression profiling from 15 distinct parasite lifecycle stages, spanning three unique ecological niches. Statistical approaches of data analysis revealed differential expression of 973 gene products that minimally describe the three major characteristics of schistosome development: asexual processes within intermediate snail hosts, sexual maturation within definitive vertebrate hosts and sexual dimorphism amongst adult male and female worms. Furthermore, we identified a group of 338 constitutively expressed schistosome gene products (including 41 transcripts sharing no sequence similarity outside the Platyhelminthes), which are likely to be essential for schistosome lifecycle progression. While highly informative, statistics-led bioinformatics mining of the transcriptional dataset has limitations, including the inability to identify higher order relationships between differentially expressed transcripts and lifecycle stages. Network analysis, coupled to Gene Ontology enrichment investigations, facilitated a re-examination of the dataset and identified 387 clusters (containing 12,132 gene products) displaying novel examples of developmentally regulated classes (including 294 schistosomula and/or adult transcripts with no known sequence similarity outside the Platyhelminthes), which were undetectable by the statistical comparisons.

Conclusions/significance: Collectively, statistical and network-based exploratory analyses of transcriptomic datasets have led to a thorough characterisation of schistosome development. Information obtained from these experiments highlighted key transcriptional programs associated with lifecycle progression and identified numerous anti-schistosomal candidate molecules including G-protein coupled receptors, tetraspanins, Dyp-type peroxidases, fucosyltransferases, leishmanolysins and the netrin/netrin receptor complex.

Figure 1. Egg and daughter sporocyst parasite stages display the most divergent transcriptional profiles within the Schistosoma lifecycle.

Correlations of M (log₂ (Alexafluor647/Alexafluor555)) ratio values between intra- (lifecycle replicate) and inter- (life-stage comparison) parasite samples were calculated. A heatmap of the correlations, calculated on data normalised within each DNA microarray, is illustrated. Red represents a high correlation and blue depicts a low or lower correlation. Un-scaled data were clustered by Euclidean distance and shown as dendrograms above and to the left of the heatmap. * Indicates female sample replicate manually removed.

Figure 2. Clusters associated with the schistosome digenetic and dioecious states are the most statistically significant transcriptional patterns identified.

Normalised and p value corrected data (see Methods ) were used to generate statistically significant gene lists (Table 1, Datasets S2 & S3). A) A heatmap of differentially expressed genes significant at an adjusted p value of <0.05 in at least one bimodal life stage comparison is illustrated (red – highly expressed, blue – weakly expressed). Log₂ intensity values were clustered by Euclidean distance and shown as a dendrogram on top of the heatmap. Three major clusters are identified (I – expressed predominantly in the definitive host, II – expressed predominantly in the intermediate host and III – expressed differentially between the adult sexes). B) Specific examples of some transcripts found in each of the three major clusters additionally showing >32-fold differential expression (DE) in at least one bimodal life-stage comparison. A full list of transcripts showing >32 fold DE in at least one bimodal life-stage comparison is found in Dataset S3. Bar graphs represent average normalised AF647 intensity values (Dataset S4) for the chosen transcripts. Array ID represents unique 50-mer oligonucleotide identifier on the DNA microarray (list of parent sequences where 50-mer oligonucleotide was designed is found in Dataset S1); BLASTx ID represents name of closest match to parent oligonucleotide sequence in NCBI database; SchistoGeneDB ID represents gene model (version 4.0) identifier of the mapped 50-mer oligonucleotide.

Figure 3. A cohort of constitutively expressed gene products is identified from DNA microarray analysis of the Schistosoma lifecycle.

Normalised log₂ AF647 intensity data were used to identify gene products constitutively expressed throughout the schistosome lifecycle (Dataset S5). A) Heatmap of transcripts exhibiting the lowest levels of variation (1% variation across lifecycle). Log₂ intensity values were clustered by Euclidean distance and shown as a dendrogram to the left of the heatmap (red – highly expressed, blue – weakly expressed). Specific examples of gene products constitutively expressed are illustrated (full list included in Dataset S5). B–D) Oligonucleotide sequences (Array ID) mapping to specific S. mansoni gene models (version 4.0) in SchistoGeneDB were analysed for Gene Ontology enrichment (p<0.01) using a hypergeometric test. B) Biological process enriched terms; C) Molecular function enriched terms and D) Cellular component enriched terms.

Figure 4. Network analysis of the schistosome lifecycle identifies egg and daughter-sporocyst enriched transcripts.

Network analysis was performed as described in Methods using the BioLayout Express^3D package . A) One network graph component containing two major clusters is illustrated (I; egg enriched transcripts and II; daughter sporocyst enriched transcripts). Different Markov clusters (MCL) within the network are represented by different coloured nodes (genes). The network analysis was performed using a Pearson correlation of r = 0.94 and Markov cluster stringency of 1.7 with accompanying details indicated. Nodes included in the identified 387 MCLs are included in Dataset S6, and form the basis of all subsequent descriptions. B) MCL1 includes transcripts highly enriched in the egg lifecycle stage (I). Some representative MCL gene examples are provided in addition to the GO terms significantly enriched in this cluster. C) MCL2 contains transcripts highly enriched in the daughter sporocyst lifecycle stage (II). Some representative MCL gene examples are provided. Bar graphs represent average normalised AF647 intensity values (Dataset S4) for all transcripts in that cluster, with replicate values obtained from independent life-stage hybridizations indicated as identical colours.

Figure 5. Markov clustering expands the repertoire of differentially expressed transcripts associated with the schistosome dioecious and digenetic state.

Network analysis was performed as described in Methods (Pearson correlation, r = 0.94 and Markov cluster stringency  = 1.7). Selected MCLs are included to illustrate transcripts enriched in the intermediate host (A, MCL9), definitive host (B, MCL69) or in adult male and female schistosomes (C, MCL7 - Male > Female; MCL12 - Female > Male). In each case, specific examples of transcripts included in the selected MCLs are listed (full list of transcripts included in each cluster is found in Dataset S6). Enriched Gene Ontology Molecular Function categories are also included where identified. Bar graphs represent average normalised AF647 intensity values (Dataset S4) for all transcripts in that cluster, with replicate values obtained from independent life-stage hybridizations indicated as identical colours.

Figure 6. The minichromosome maintenance (MCM) heterohexamer is coordinately expressed in both mother sporocysts and 3-week adult worms.

A) MCL190 identified 9 gene products sharing a peak of gene expression in the mother sporocyst and 3-wk adult worms. Bar graphs represent average normalised AF647 intensity values (Dataset S4) for all transcripts in MCL190, with replicate values obtained from independent life-stage hybridizations indicated as identical colours. B) List of genes contained within MCL190 demonstrates the presence of 4 out of the 6 transcripts constituting the eukaryotic MCM heterohexamer. C) Heatmap (Spearman's rank correlation coefficient, average linkage clustering, red – highly expressed, blue – weakly expressed) of all 6 MCM heterohexamer subunits (Dataset S10, including two MCM2 paralogs) defined by both array and GeneDB identifiers (ID). * indicates three MCM transcripts not clustered within MCL190. D) Scatterplot of MCM heterohexamer normalised AF647 intensity units in mother sporocyst versus 3-wk worm samples. A correlation coefficient of 0.738 suggests coordinated expression of the MCM heterohexamer with individual outliers (Smp_054840 and Smp_037590) indicated.

Figure 7. Analysis of three important schistosome gene families reveals differential transcription of individual members across the parasite's lifecycle.

Using the latest version of the schistosome genome assembly (version 4.0), all putative fucosyltransferase- (PF00852, Dataset S7), tetraspanin- (PF00335, Dataset S8) and GPCR- (PF00001, PF00002 and PF00003, Dataset S9) family members were identified. DNA microarray oligonucleotides were mapped to specific regions of the predicted genes using the S. mansoni genome browser function of S. mansoni GeneDB . For those predicted gene family members where 50-mer oligonucleotides could be mapped, expression data was extracted and analysed. Heatmap representation of all S. mansoni: A) fucosyltransferases (PF00852), B) tetraspanins (PF00335) and C) GPCRs (PF00001, PF00002 and PF00003) differentially expressed across the parasite's lifecycle (all data included in Datasets S7–9. Heatmaps (red – highly expressed, blue – weakly expressed) were created from normalised AF647 intensities using a Spearman's rank correlation coefficient and average linkage clustering ( Methods ). Corresponding real time RT-PCR confirmation of selected family members (*) (D – fucosyltransferases: CD154939/Smp_054300 & AF183577/Smp_148850; E – tetraspanins: CD148755/Smp_176870 & CD069616/Smp_155310; F – GPCRs: CONTIG2914/Smp_128940 & CD098416/Smp_099670) is indicated to the right of heatmap.