Genome-wide transcriptome profiling reveals functional networks involving the Plasmodium falciparum drug resistance transporters PfCRT and PfMDR1

Abstract

Background: The acquisition of multidrug resistance by Plasmodium falciparum underscores the need to understand the underlying molecular mechanisms so as to counter their impact on malaria control. For the many antimalarials whose mode of action relates to inhibition of heme detoxification inside infected erythrocytes, the digestive vacuole transporters PfCRT and PfMDR1 constitute primary resistance determinants.

Results: Using gene expression microarrays over the course of the parasite intra-erythrocytic developmental cycle, we compared the transcriptomic profiles between P. falciparum strains displaying mutant or wild-type pfcrt or varying in pfcrt or pfmdr1 expression levels. To account for differences in the time of sampling, we developed a computational method termed Hypergeometric Analysis of Time Series, which combines Fast Fourier Transform with a modified Gene Set Enrichment Analysis. Our analysis revealed coordinated changes in genes involved in protein catabolism, translation initiation and DNA/RNA metabolism. We also observed differential expression of genes with a role in transport or coding for components of the digestive vacuole. Interestingly, a global comparison of all profiled transcriptomes uncovered a tight correlation between the transcript levels of pfcrt and pfmdr1, extending to dozens of other genes, suggesting an intricate regulatory balance in order to maintain optimal physiological processes.

Conclusions: This study provides insight into the mechanisms by which P. falciparum adjusts to the acquisition of mutations or gene amplification in key transporter loci that mediate drug resistance. Our results implicate several biological pathways that may be differentially regulated to compensate for impaired transporter function and alterations in parasite vacuole physiology.

Fig. 1

Alignment of all time series to the 3D7 reference transcriptome reveal differences in growth rates between the profiled parasite lines. a Plot of all time courses before temporal alignment of the transcriptome data. For each transcriptome data set, the original 8 time points of each 48-hour time course (corresponding to a single intra-erythrocytic developmental cycle (IDC)) were aligned to the transcriptome of 3D7 that earlier had been elucidated at a one-hour resolution [46]. b Phaseograms of all parasite lines in the FCB set after dynamic time warping. Genes (on the Y-axis) were ordered according to their phase of expression computed using Fast Fourier Transform (FFT), which interpolated expression levels at unassayed time points. The resulting output generated a smooth sine curve with a single peak of maximum expression (shown in red, spanning a heatpmap gradient to the minimum transcript levels shown in blue) for each gene across the entire 48-hour IDC. c Phaseograms of the three parasite lines in the 7G8 set after dynamic time warping. d Principal Component Analysis (PCA) confirming the proper alignment of all transcriptomes. The ~3,000 to ~3,500 genes represented in each transcriptome were used to generate a matrix of values that was reduced to a tri-dimensional projection using PCA. The first two principal components accurately recover the progression of the lines through the 48-hour lifecycle while the third principal component focuses on differences between the lines as symbolized by the size of the circles. Circle colors refer to the parasite lines depicted in (a)

Fig. 2

Differential gene expression detected by chromosome mapping. a Chromosome mapping of gene expression fold differences between FCB and FCB^{pfmdr1_3′KD} shows differential transcription of a large fragment on chromosome 5 corresponding to the pfmdr1-containing 100-kb amplicon region. For each gene located on chromosome 5, the expression fold difference between FCB and FCB^{pfmdr1_3′KD} was mapped using the chromosomal coordinates. Interestingly, PF3D7_0507500 (a subtilisin-like protease) and PF3D7_0501200 (a parasite-infected erythrocyte surface protein) show higher expression in FCB in comparison to both FCB^{pfmdr1_3′KD} and 106/1 (see panel B). The inset shows the genes within the pfmdr1-containing amplicon region. b Chromosome mapping of gene expression fold differences between FCB and 106/1 for each gene on chromosome 5 with the inset showing genes within the pfmdr1-containing amplicon region. c Chromosome mapping of gene expression fold differences between 7G8 and 7G8^pfcrt_T76K shows differential gene expression of a large number of genes located on chromosome 2 (inset shows a cluster of genes coding for proteins involved in parasite invasion). d Chromosome mapping of gene expression fold differences between 7G8 and 7G8^pfcrt_CTL for each gene on chromosome 2. In each pairwise comparison, the fold difference in gene expression was calculated as the difference between the areas under the curve computed on the full 384 time point data

Fig. 3

Genes with significant expression fold changes organized by categories for pairwise comparisons of the transcriptomes in the pfmdr1 data set. Gene expression fold changes were calculated as the difference between the areas under the curve generated on the full 384 time point data and were then normalized using a baseline created by a background pool of 11 transcriptomes. For each gene, all 55 possible pairwise combinations of transcriptomes created a distribution of expression fold changes, which was used as a normalization factor. a Genes showing higher expression in the parental line FCB in comparison to pfmdr1-knock-down strain FCB^{pfmdr1_3′KD}. To identify significant differences in gene expression, we applied a threshold of a fold change >1.5 and a normalized fold change at least 2 standard deviations (SD) above the mean. b Genes showing higher expression in FCB in comparison to 106/1. c Genes showing lower expression in the parental line FCB compared to the knock-down strain FCB^{pfmdr1_3′KD}. We applied a threshold of a fold change <0.6 and a normalized fold change 3 SD below the mean. d Genes showing lower expression in FCB in comparison to 106/1. e Venn diagrams depicting the overlap of genes with higher (fold change >1.5, left) or lower (fold change <0.6, right) expression in FCB for the FCB vs. FCB^{pfmdr1_3′KD} and FCB vs. 106/1 pairwise comparisons

Fig. 4

Hypergeometric Analysis of Time Series (HATS) for pairwise comparisons of the transcriptomes in the pfmdr1 data set. HATS was performed by combining FFT to align the time series with a hypergeometric distribution approach based on random permutations to compute an enrichment score for each gene set. Normalization of the enrichment scores was performed using an expression baseline constituted by a background pool of 11 transcriptome data sets (see Methods). We applied a threshold of normalized gene set rank >0.9 for gene sets we considered positively enriched and normalized gene set rank <0.1 for gene sets considered negatively enriched. Gene sets were built using GO, KEGG and the Malaria Parasite Metabolic Pathway (listed in Additional file 26: Table S17). a Gene sets that are enriched in FCB in comparison to FCB^{pfmdr1_3′KD} (mean ranking values >0.9). b Gene sets that are enriched in FCB in comparison to 106/1 (mean ranking values >0.9). c Gene sets that are depleted in FCB in comparison to FCB^{pfmdr1_3′KD} (mean ranking values <0.1). d Gene sets that are depleted in FCB in comparison to 106/1 (mean ranking values <0.1). e Stage-associated enrichment analysis in FCB in comparison to FCB^{pfmdr1_3′KD}. Gene sets with significant enrichment (orange) or depletion (blue) for three windows of the parasite IDC correspond to rings (R, 12–16 h post-invasion), trophozoites (T, 24–28 h post-invasion) and schizonts (S, 32–36 h post-invasion). f Stage-associated enrichment analysis in FCB in comparison to 106/1. Gene sets with significant enrichment (orange) or depletion (blue) during the same three stages as in (E). g Venn diagrams depicting overlaps of gene sets significantly enriched in FCB (mean rank of enrichment score >0.9) for both FCB vs. FCB^{pfmdr1_3′KD} and FCB vs. 106/1 pairwise comparisons, as identified from our HATS analysis. h Venn diagrams depicting overlaps of gene sets significantly depleted in FCB (mean rank of enrichment score <0.1) for both FCB vs. FCB^{pfmdr1_3′KD} and FCB vs. 106/1 pairwise comparisons, as identified using HATS

Fig. 5

Genes with significant expression fold changes organized by categories for pairwise comparisons of the transcriptomes in the pfcrt data set. Gene expression fold changes were calculated as the difference between the areas under the curve generated on the full 384 time point data and were then normalized using a baseline created by the background pool of 11 transcriptomes. For each gene, all 55 possible pairwise combinations of transcriptomes created a distribution of expression fold changes, which was used as a normalization factor. To identify significant differences, we applied a threshold of a fold change >1.5 and a normalized fold change at least 2 standard deviations (SD) above the mean for higher gene expression, and a threshold of a Fold change <0.6 and a normalized fold change 3 SD below the mean for lower expression. a Genes with higher expression in the parental line 7G8 in comparison to the recombinant control strain 7G8^pfcrt_CTL. b Genes showing higher expression in 7G8 compared to the recombinant back-mutant strain 7G8^pfcrt_T76K. c Venn diagrams depicting the overlap of genes with higher (fold change >1.5, left) or lower (fold change <0.6, right) expression in 7G8 for both 7G8 vs. 7G8^pfcrt_CTL and 7G8 vs. 7G8^pfcrt_T76K pairwise comparisons. d Genes with lower expression in comparison with the recombinant control strain 7G8^pfcrt_CTL. e Genes with lower expression in the parental line 7G8 in comparison with the recombinant back-mutant strain 7G8^pfcrt_T76K

Fig. 6

Hypergeometric Analysis of Time Series (HATS) for pairwise comparisons of the transcriptomes in the pfcrt data set. As for the pfmdr1 dataset, we performed HATS analysis by combining FFT to align the time series with a hypergeometric distribution approach based on random permutations to compute an enrichment score for each gene set. Normalization of the enrichment scores was performed using the expression baseline constituted by the background pool of 11 transcriptome data sets. As before, we applied a threshold of normalized gene set rank >0.9 for gene sets we considered positively enriched and normalized gene set rank <0.1 for gene sets considered negatively enriched. The gene sets were built using GO, KEGG, and the Malaria Parasite Metabolic Pathway. a Gene sets that are enriched (left, mean ranking values >0.9) or depleted (right, mean ranking values <0.1) in 7G8 in comparison with 7G8^pfcrt_CTL. b Gene sets that are enriched (left, mean ranking values >0.9) or depleted (right, mean ranking values <0.1) in 7G8 in comparison to 7G8^pfcrt_T76K. c Venn diagrams depicting overlaps of gene sets significantly enriched (left) or depleted (right) in 7G8 for both 7G8 vs. 7G8^pfcrt_CTL and 7G8 vs. 7G8^pfcrt_T76K pairwise comparisons

Fig. 7

Analysis of the transcriptional differences between 7G8^pfcrt_CTL and 7G8^pfcrt_T76K. a Genes showing higher expression (fold change >1.5 and normalized fold change ≥2 SD above the mean) in the recombinant control 7G8^pfcrt_CTL in comparison to the recombinant back-mutant strain 7G8^pfcrt_T76K. Gene expression fold changes were calculated as the difference between the areas under the curve generated on the full 384 time point data and were then normalized using a baseline created by the background pool of 11 transcriptomes. b Genes showing a lower expression (fold change <0.6 and normalized fold change ≤3 SD below the mean) in the recombinant control 7G8^pfcrt_CTL in comparison to the recombinant back-mutant strain 7G8^pfcrt_T76K . c HATS analysis was used to identify gene sets that are enriched in 7G8^pfcrt_CTL in comparison to 7G8^pfcrt_T76K (mean ranking values >0.9). Enrichment scores were normalized using the expression baseline constituted by the background pool of 11 transcriptome data sets. d Gene sets that are depleted in 7G8^pfcrt_CTL in comparison to 7G8^pfcrt_T76K (mean ranking values <0.1), as identified from HATS analysis. e Stage-associated enrichment analysis. Gene sets with significant enrichment (orange) or depletion (blue) in 7G8^pfcrt_CTL in comparison to 7G8^pfcrt_T76K for three windows of the parasite IDC corresponding to rings (R, 12–16 h post-invasion), trophozoites (T, 24–28 h post-invasion) and schizonts (S, 32–36 h post-invasion) using stage-enrichment analysis

Fig. 8

Investigation of pfcrt and pfmdr1 expression networks: Gene expression correlations between all P. falciparum genes and pfcrt/pfmdr1. a Cluster heatmap of expression data for P. falciparum genes. The hierarchical clustering was generated using Pearson correlation coefficients (PCC) calculated using log₂-transformed and normalized expression values of 2,600 genes across 55 pairwise comparisons of 11 parasite transcriptome data sets. Heatmaps of the hierarchical clustering show several domains of high correlation between gene pairs, including a large cluster corresponding to strong interactions between plasmepsin X (PF3D7_0808200) and multiple genes involved in cell motility such as myosin (PF3D7_1246400), kinesin-like protein (PF3D7_0724900) and invasion like rhoptry (PF3D7_0414900) or IMC1-related (PF3D7_0304100) proteins. Two smaller clusters (boxes 1 and 2) are shown including one containing both pfcrt and pfmdr1. b Venn diagrams showing the number of genes whose expression highly correlates (PCC >0.7) (top) or anti-correlates (PCC <−0.7) (bottom) with that of both pfcrt and pfmdr1. PCC values were calculated using log₂-transformed and normalized expression values of 2600 genes across 110 pairwise comparisons of 11 parasite transcriptome data sets