Impact of Sickle Cell Trait Hemoglobin on the Intraerythrocytic Transcriptional Program of Plasmodium falciparum

Abstract

ABSTRACT

Sickle-trait hemoglobin (HbAS) confers nearly complete protection from severe, life-threatening falciparum malaria in African children. Despite this clear protection, the molecular mechanisms by which HbAS confers these protective phenotypes remain incompletely understood. As a forward genetic screen for aberrant parasite transcriptional responses associated with parasite neutralization in HbAS red blood cells (RBCs), we performed comparative transcriptomic analyses of Plasmodium falciparum in normal (HbAA) and HbAS erythrocytes during both in vitro cultivation of reference parasite strains and naturally occurring P. falciparum infections in Malian children with HbAA or HbAS. During in vitro cultivation, parasites matured normally in HbAS RBCs, and the temporal expression was largely unperturbed of the highly ordered transcriptional program that underlies the parasite’s maturation throughout the intraerythrocytic development cycle (IDC). However, differential expression analysis identified hundreds of transcripts aberrantly expressed in HbAS, largely occurring late in the IDC. Surprisingly, transcripts encoding members of the Maurer’s clefts were overexpressed in HbAS despite impaired parasite protein export in these RBCs, while parasites in HbAS RBCs underexpressed transcripts associated with the endoplasmic reticulum and those encoding serine repeat antigen proteases that promote parasite egress. Analyses of P. falciparum transcriptomes from 32 children with uncomplicated malaria identified stage-specific differential expression: among infections composed of ring-stage parasites, only cyclophilin 19B was underexpressed in children with HbAS, while trophozoite-stage infections identified a range of differentially expressed transcripts, including downregulation in HbAS of several transcripts associated with severe malaria in collateral studies. Collectively, our comparative transcriptomic screen in vitro and in vivo indicates that P. falciparum adapts to HbAS by altering its protein chaperone and folding machinery, oxidative stress response, and protein export machinery. Because HbAS consistently protects from severe P. falciparum, modulation of these responses may offer avenues by which to neutralize P. falciparum parasites.

IMPORTANCE Sickle-trait hemoglobin (HbAS) confers nearly complete protection from severe, life-threatening malaria, yet the molecular mechanisms that underlie HbAS protection from severe malaria remain incompletely understood. Here, we used transcriptome sequencing (RNA-seq) to measure the impact of HbAS on the blood-stage transcriptome of Plasmodium falciparum in in vitro time series experiments and in vivo samples from natural infections. Our in vitro time series data reveal that, during its blood stage, P. falciparum’s gene expression in HbAS is impacted primarily through alterations in the abundance of gene products as opposed to variations in the timing of gene expression. Collectively, our in vitro and in vivo data indicate that P. falciparum adapts to HbAS by altering its protein chaperone and folding machinery, oxidative stress response, and protein export machinery. Due to the persistent association of HbAS and protection from severe disease, these processes that are modified in HbAS may offer strategies to neutralize P. falciparum.

INTRODUCTION

In regions across Africa where malaria is endemic, red blood cell (RBC) variants are highly prevalent and can confer significant protection against severe malaria (1). Heterozygosity for sickle hemoglobin (Hb) provides one striking example of the protection afforded by an RBC variant: sickle-cell trait (HbAS) reduces the risk of severe, life-threatening malaria by over 90% (2). The functional underpinnings of this extensive protection afforded by HbAS against severe disease hold potential for uncovering mechanisms of parasite pathogenesis that can be exploited by future interventions.

Candidate protective mechanisms reported thus far are diverse and provide evidence that HbAS disrupts mediators of pathogenesis and the parasite life cycle in infected RBCs (iRBCs). Parasite maturation may be inhibited by either low oxygen tension (3–6) or host microRNAs (miRNAs) that are enriched in HbSS and HbAS RBCs (7), though others have reported equivalent growth in HbAA and HbAS iRBCs under the same conditions (8, 9). Multiple lines of evidence support a mechanism in which cytoadherence, a central property of Plasmodium falciparum’s pathogenesis, is reduced in HbAS iRBCs. HbAS reduces both the display on the iRBC surface of P. falciparum erythrocyte membrane protein 1 (PfEMP-1) and knob density (10). This is supported by the observation that HbAS impairs export of proteins, including PfEMP-1, to the iRBC surface via Maurer’s clefts, which sort and traffic parasite proteins (11). The attenuated export of surface proteins and concomitant reduction in cytoadherence would compromise the capacity of HbAS iRBCs to sequester in the microvasculature endothelium, thereby diminishing endothelial activation—as has been observed in vitro (12)—as well as subsequent events in the microvasculature that can lead to severe disease in vivo (13). These downstream phenotypes of HbAS iRBCs are the result of incompletely understood upstream parasite responses to the altered RBC environment.

The natural protection against severe falciparum malaria conferred by HbAS provides an opportunity to identify parasite mechanisms that are neutralized and therefore critical to pathogenesis. In this study, we used transcriptome sequencing (RNA-seq) to perform an unbiased comparative transcriptomic survey of the P. falciparum transcriptional changes that occur in HbAS iRBCs. Owing to the reduced cytoadherence observed in HbAS iRBCs, we hypothesized that gene products functioning in cytoadherence would be downregulated in HbAS. Therefore, we investigated parasite transcription over the intraerythrocytic development cycle (IDC) in vitro in tightly synchronized parasite populations in HbAA and HbAS iRBCs as well as in freshly collected parasites from Malian children with uncomplicated malaria with HbAA or HbAS.

RESULTS

Parasite maturation and transcriptome characteristics.

We interrogated two reference Plasmodium falciparum strains that originated where the sickle allele of β-globin is common (14): 3D7 from West Africa (15) and FUP from Uganda (16). Based on morphology assessed by microscopy, parasites matured similarly in HbAA and HbAS iRBCs for both 3D7 and FUP (Fig. 1A). Parasites in HbAA and HbAS iRBCs transitioned from rings to trophozoites at 27 h postinvasion (hpi) (3D7) and 21 hpi (FUP).

FIG 1

P. falciparum maturation and the overall transcriptional program are conserved in HbAS iRBCs. (A) The relative abundance of ring, trophozoite, and schizont stage parasites from light microscopy readings in HbAA (solid line) and HbAS (dotted line) iRBCs. 3D7 and FUP parasites mature at similar rates in HbAA and HbAS iRBCs. Counts are based on 500 parasites from each sample per time point and were assessed by a reader masked to time point and RBC type. (B) Principal-component analysis of transcript expression following variance-stabilizing transformation based on the 500 transcripts with the highest variance across all samples. Samples cluster according to parasite strain and developmental stage. (C) Expression of 5,499 transcripts passing thresholds in the ASM276v2 transcriptome in 3D7 and FUP isolates. In both 3D7 and FUP, the ordering of expression is largely conserved between HbAA and HbAS samples. Transcripts are ordered vertically by the peak expression time in HbAA RBCs for each parasite strain. The columns along the x axis depict the normalized expression values of each transcript per time point. Expression values were averaged across replicates and are normalized around the mean for each transcript over the time series and depicted as a z score of standard deviations from the mean. Values range between −1.5 (cyan) and 1.5 (yellow). Scripts for RNA-seq read preparation and quantification are available on our GitHub page.

Principal-component analysis (PCA) based on transcript counts clustered samples by parasite line and developmental stage (Fig. 1B). Notably, the principal components that explain the vast majority of the transcriptional variance between samples did not separate parasite transcriptomes in HbAS iRBCs from those in HbAA iRBCs. Similarly, visualization of normalized expression of 5,499 parasite transcripts suggests that the progression of P. falciparum’s IDC transcriptional program remains largely intact in HbAS iRBCs (Fig. 1C). These data indicate both that the parasite populations in each sample were highly synchronous and that HbAS does not broadly disrupt P. falciparum’s transcriptional program during the IDC.

Transcriptional synchrony of parasites in HbAS RBCs.

P. falciparum maturation in HbAS iRBCs was described previously to be delayed compared to that in HbAA iRBCs (6). To assess parasite maturation at a finer scale, we computed, for the subset of transcripts that had a single transcription peak, the temporal shift of each transcript’s expression peak in HbAS iRBCs compared to HbAA (Fig. 2A). The distribution of peak shifts between HbAA and HbAS over the entire time course was centered around zero for both 3D7 and FUP (Fig. 2B), with slight but nonsignificant mean peak shift delays in HbAS RBCs for both 3D7 (−1.62 h) and FUP (−0.91 h) parasites.

FIG 2

Shifts in peak expression time between HbAA and HbAS. (A) Peak expression time of the subsets of mono-peaking parasite transcripts in HbAA and HbAS iRBCs for 3D7 (left pair) or FUP (right pair) parasites. For each pair, transcripts are ordered vertically by the peak expression time in HbAA iRBCs. In both 3D7 and FUP, the ordering of expression is largely conserved between HbAA and HbAS samples. The columns along the x axis depict the normalized expression values of each transcript per time point. Expression values were averaged across replicates and are normalized around the mean for each transcript over the time series and depicted as a z score of standard deviations from the mean. Values range between −1.5 (cyan) and 1.5 (yellow). (B) Temporal shifts for peak expression time point for the subset of mono-peaking parasite transcripts in HbAA compared with HbAS iRBCs for 3D7 (left) or FUP (right) parasites. For each transcript, the time between its expression peak in parasites in HbAA and HbAS iRBCs was calculated and averaged across β-globin replicates for each strain. Each x-axis category indicates the number of transcripts that normally peak at this time point in HbAA iRBCs, and the y values indicate the distributions of the shifts of the peaks of these transcripts in HbAS iRBCs. Circles are sized relative to the number of transcripts with that y value. The top histogram indicates the number of parasite transcripts that peak in HbAA iRBCs at each x-value time point, and the right histogram indicates the number of transcripts with peak shifts of each y value. The distribution of peak shifts (displayed along the right y axis) is centered around zero, indicating that most transcripts peak at the same time in parasites growing in HbAA and HbAS iRBCs. See Jupyter Notebook “expressionPeakChanges.ipynb” on our GitHub for the full workflow.

Differentially expressed transcripts identified in HbAS in vitro.

HbAS impacts P. falciparum’s transcriptional program in the IDC at various time points in 3D7 and FUP isolates (Fig. 3A and B), though in both strains we observed an absence of differential expression through much of the ring stage, with transcriptional divergence apparent in the mid- to late-trophozoite stage that accelerated as parasites entered the schizont stage (Fig. 3B).

FIG 3

Differential expression analysis reveals parasite transcript expression in HbAS iRBCs diverges from HbAA parasites primarily in the schizont stage. (A) Scatterplots at 3, 12, 24, 36, 42, and 48 h postinvasion of the log₂ fold change of individual transcripts for HbAS parasites compared to HbAA parasites (hpi) in 3D7 (x axes) and FUP (y axes) parasites. Points in the upper right and lower left quadrants signify transcripts that are differentially expressed in the same direction between parasite strains, i.e., those in the upper right quadrant are upregulated in both 3D7 and FUP in HbAS iRBCs, and those in the lower left quadrant are downregulated in both 3D7 and FUP in HbAS iRBCs. Point color indicates the statistical significance of DESeq2 differential expression analysis for each transcript: white points have an adjusted P value of >0.05 in both 3D7 and FUP, red points have an adjusted P value of <0.05 in 3D7 but not FUP, blue points have an adjusted P value of <0.05 in FUP but not 3D7, and purple points have an adjusted P value of <0.05 in both 3D7 and FUP. (B) Black bars indicate counts of transcripts in 3D7 (top) and FUP (bottom) classified as downregulated (left column) or upregulated (right column) at a significance level of an adjusted P value of <0.05 at each time point, from 3 hpi to 48 hpi. Each plot is colored according to the predominant stage of the parasites at each time point: ring (red), trophozoite (blue), or schizont (green). There were minimal shared differentially expressed transcripts until 39 hpi, after which 966 transcripts were differentially expressed in both 3D7 and FUP (493 downregulated and 460 upregulated). (C) Scatterplot of transcript mean fold changes within gene ontology (GO) terms that are enriched for differentially expressed transcripts in HbAS iRBCs in both 3D7 (x axis) and FUP (y axis) parasites. Color indicates GO term, shading indicates hpi (darker = later hpi), and size of point indicates the proportion of transcripts within the GO term that are differentially expressed at an adjusted P value of <0.05. Significant expression changes of GO terms occurred at later time points during the schizont stage, including upregulation of Maurer’s cleft transcripts in HbAS iRBCs and concurrent downregulation of transcripts involved in the endoplasmic reticulum, the tricarboxylic acid cycle, cysteine-type peptidase activity, and nucleosomes.

We identified the transcripts that were differentially expressed in the same direction in both 3D7 and FUP at each time point (Fig. 3A), which indicated a consistent late-stage effect of HbAS on transcriptional phenotypes across parasite strains (Fig. 3B). In a gene ontology (GO) enrichment analysis of this shared set of transcripts (17) (Fig. 3C), the only GO term that was enriched for upregulated transcripts in both 3D7 and FUP in HbAS iRBCs at any time point in our time series experiments was that for Maurer’s clefts (Fig. 4A), which are parasite-derived membranous structures that sort and traffic proteins and in HbAS iRBCs are both dysmorphic and dysfunctional (8). Several critical components of Maurer’s clefts, including the membrane-associated histidine-rich protein 1 (MAHRP1; PF3D7_1370300), and skeleton-binding protein 1 (SBP1; PF3D7_0501300), are encoded by transcripts that were among the most significantly upregulated in HbAS RBCs.

FIG 4

Differential expression of Maurer’s cleft, histone, and SERA transcripts in HbAS RBCs. (A) Heat maps of log₂ fold change between HbAS and HbAA of Maurer’s cleft transcripts in 3D7 and FUP parasites. Transcript IDs were extracted from the Maurer’s cleft GO term (GO:0020036). (B) Heat maps of log₂ fold change between HbAS and HbAA of histone transcripts in 3D7 and FUP parasites. Transcript IDs were extracted from the nucleosome GO term (GO:0000786). (C) Heat maps of log₂ fold change between HbAS and HbAA of SERA transcripts in 3D7 and FUP parasites. The color bar is restricted to log₂ fold change considered significant in this study [|log₂(HbAS/HbAA)| > log₂(1.5)].

The most downregulated GO term in both 3D7 and FUP in HbAS iRBCs was that of the nucleosome (GO:0000786), which comprises eight transcripts that encode the core histones and histone variants in P. falciparum. In both 3D7 and FUP growing in HbAS RBCs, all histone transcripts were downregulated at 48 hpi (Fig. 4B). Similarly, downregulated transcripts were enriched in the tricarboxylic acid (TCA) cycle (GO:0006099) and the endoplasmic reticulum (ER) (GO:0005783) at late time points in 3D7 and FUP in HbAS iRBCs. The ER GO term encompasses a wide range of proteins, including heat shock proteins, phosphatases, peptidyl-prolyl cis-trans-isomerases, and plasmepsins. Of 53 transcripts in the ER GO term, 30 transcripts were downregulated in both strains in HbAS iRBCs between 45 and 48 hpi.

Among molecular function GO terms, downregulated transcripts were enriched only in that of cysteine-type peptidase activity (GO:0008234), which primarily comprises the serine repeat antigen (SERA) proteases that play a pivotal role enabling parasite egress from the iRBC. SERA2 to -7 were significantly downregulated in 3D7 and FUP parasites in HbAS iRBCs at late time points (Fig. 4C).

Temporal dysregulation of transcript expression in HbAS iRBCs.

To identify the transcripts that are temporally dysregulated, we analyzed our time series data using dynamic time warping (DTW) (Fig. 5A), which determines between two series the optimal alignment of a transcript’s expression and assigns a score as a measure of similarity (low score) (Fig. 5B, top row) or difference (high score) (Fig. 5B, bottom row). Overall, DTW scores comparing transcript temporal expression between HbAA and HbAS RBCs were close to normally distributed with a positive skew in 3D7 and in FUP (Fig. 5A). Between strains, mean transcript DTW scores comparing HbAA versus HbAS time series were positively linearly correlated when computed for 3D7 or FUP (r = 0.524; P < 0.0001), indicating a conserved impact on transcript timing across strains.

FIG 5

HbAS impacts the temporal expression of transcripts governing gene regulation. (A) Scatterplot of dynamic time warping (DTW) scores in time series data of individual transcripts in HbAA compared to HbAS samples in 3D7 (x axis) and FUP (y axis) parasites. For each strain, scores were computed pairwise between HbAA and HbAS samples and averaged. Transcripts with dynamic time warping scores greater than two standard deviations from the overall mean (points to the right of the red line for 3D7 and points above the blue line for FUP) were considered temporally dysregulated. Histograms indicate overall distribution of transcript DTW scores in either 3D7 (top) or FUP (right). (B) The Euclidean distance scores and warping path of selected transcripts that correspond to the highlighted circles in panel A. Mean-normalized transcript expression is plotted along the x and y axes: that for HbAA iRBCs is plotted from top to bottom along the left y axis, and that for HbAS iRBCs is plotted from left to right along the upper x axis. The heat map matrix is colored according to the Euclidean distance between the mean-normalized expression values between HbAA and HbAS at the corresponding time points. Dynamic time warping scores are cumulatively computed along the warping path (white line). Transcripts for the ring-infected erythrocyte surface antigen (RESA, PF3D7_010220) and the putative transcription factor TFIIB (PF3D7_0110800) (top) had very low DTW scores between HbAA and HbAS samples in both strains, while those for the Plasmodium exported protein (hyp2; PF3D7_0220500) and a glycerol kinase (PF3D7_1351600) had high DTW scores, indicating dissimilar transcriptional profiles between HbAA and HbAS RBCs (bottom). (C) Heat maps of mean-normalized expression of the 155 (3D7; top) and 162 (FUP; bottom) transcripts with DTW scores exceeding two standard deviations from the mean score in HbAA (left) and HbAS (right) iRBCs (points beyond the dashed lines in panel A). Transcripts for each parasite strain are order vertically by their peak expression time in HbAA iRBCs. The ordering of peak expression in HbAS samples is not conserved in temporally dysregulated transcripts identified by DTW. Transcripts are ordered by peak expression time in HbAA samples and displayed as a z score of standard deviations from the mean. (D) Scatterplot of mean DTW scores within GO terms that are enriched for temporally dysregulated transcripts in HbAS iRBCs in both 3D7 (x axis) and FUP (y axis) parasites. Color indicates GO term, shading indicates hpi (darker = later hpi), and size of point indicates the proportion of transcripts within the GO term that are temporally dysregulated as defined by a DTW score greater than 2 standard deviations from the overall mean of DTW scores (points beyond the dashed lines in panel A). Transcripts involved in gene regulation processes demonstrate higher DTW scores. See Jupyter Notebook “DTW.ipynb” on our GitHub for full workflow.

We observed in HbAS iRBCs 155 temporally dysregulated transcripts in 3D7 and 162 in FUP (Fig. 5C). In contrast to the generally unperturbed expression patterns of the overall transcriptome (Fig. 1C), these subsets of transcripts demonstrate aberrant expression profiles in HbAS iRBCs (Fig. 5C). Transcripts were grouped by GO terms and, for each term, the mean DTW score was calculated between parasites in HbAA and HbAS. We observed that transcripts associated with antigenic variation and cytoadhesion (i.e., the variant surface antigens) were among the most conserved in the timing of their expression (Fig. 5D). In contrast, the highest DTW scores were recorded for transcripts involved in transcription by RNA polymerase II (GO:0006366), RNA splicing (GO:0008380), regulation of translational initiation (GO:0006446), and ribosome biogenesis (GO:0042254), which collectively suggest that HbAS may have posttranscriptional consequences beyond what can be observed in our data.

We focused more closely on ring-stage expression by screening with the temporal alignment Kendall-Tau (TAKT) algorithm, which measures the overall similarity and temporal shift between two curves (18), and then computing for these candidate transcripts DTW scores. Using this approach in ring-stage parasites, we identified 9 transcripts in 3D7 and 11 in FUP that were dysregulated between parasites in HbAA and HbAS, though none of these were shared between the two strains.

In vivo transcriptomes of P. falciparum in children with HbAS.

We next compared parasite transcriptomes between 16 uncomplicated P. falciparum episodes in HbAS children matched to 16 episodes in HbAA children in Kéniéroba, Mali. We first used stage-specific ring (n = 7 transcripts), trophozoite (n = 4), and schizont (n = 6) transcripts to classify 24 infections as ring stage and 12 as trophozoite stage (Fig. 6A). PCA and hierarchical clustering supported these classifications (Fig. 6B). Late-stage trophozoites are rare in circulation (19), so we re-estimated the stage of each infection using an independent published mixture model, which suggested that the trophozoite classifications were early trophozoites (20) and otherwise confirmed initial findings.

FIG 6

Differential expression profiling of parasite transcripts in freshly collected infections in Malian children with uncomplicated malaria with HbAA or HbAS. (A) Proportion of parasite stage-specific transcripts expressed in each patient’s infection (rows) in HbAA (left) and HbAS (right) children. Expression of stage-restricted transcripts for ring (red), trophozoite (blue), and schizont (green) stages were measured as TPM, summed within each stage, and expressed as a proportion of stage for each infection. We observed a nearly complete absence of schizont-specific transcript expression, and most infections were predominantly either rings or trophozoites with little admixture. See Jupyter Notebook “Genes_by_stage.ipynb” on our GitHub for the full workflow. (B) Principal-component analysis of transcript expression from in vivo isolates following variance-stabilizing transformation based on the 500 transcripts with the highest variance across all samples. Clusters separate according to parasite stage identified from stage-specific transcript expression. (C and D) Volcano plots of differential transcript expression conditioned by parasite stage. The x axis shows the log₂ fold change of HbAS versus HbAA children, with positive values indicating increased expression in HbAS. The y axis shows −log₁₀ P values of each transcript, and horizontal dashed lines indicate P value thresholds that capture the transcripts with adjusted P values of <0.05. (E) Scatterplots of the log₂ fold change of individual transcripts for severe-malaria parasites compared to uncomplicated-malaria parasites in Gambian children (y axes) and HbAS compared to HbAA parasites (x axes). Points in the upper right and lower left quadrants signify transcripts that are differentially expressed in the same direction between studies; i.e., those in the upper right quadrant are upregulated in both severe malaria and HbAS, and those in the lower left quadrant are downregulated in both severe malaria and HbAS. All transcripts included were statistically significant in severe malaria versus uncomplicated malaria (26). Red points have an adjusted P value of <0.05 in trophozoite parasites from HbAS versus HbAA Malian children.

Differential expression in vivo.

The only differentially expressed transcript in ring-stage infections was cyclophilin 19B (cyp19B; PF3D7_1115600), a peptidyl-prolyl cis-trans-isomerase (PPIase) (Fig. 6C): compared to children with HbAA, cyp19B in those with HbAS was downregulated 3-fold (adjusted P = 0.023010). Quantitative PCR (qPCR) confirmed that cyp19B is reduced in ring-stage parasite transcriptomes in HbAS. Cyclophilins function in protein folding and trafficking, in RNA processing, and generally as molecular chaperones (21) and are most highly expressed in schizont-stage parasites in vitro (22).

We observed in trophozoite-stage patient samples differential expression of 74 transcripts (Fig. 6D), among which we did not observe any enrichment for GO terms. The most upregulated transcript in HbAS trophozoites was a putative protein phosphatase (PF3D7_1127000; log₂ fold change = 5.00; adjusted P = 1.08 × 10⁻¹⁸). Also among the 10 most upregulated transcripts in HbAS trophozoites were the gene for the 130-kDa glycophorin binding protein (GBP130; PF3D7_1016300) and the two glycophorin binding protein homologs (GBPH, PF3D7_1401000; GBPH2, PF3D7_1301200), each upregulated more than 8-fold in HbAS trophozoites. Additionally, one of the most significantly upregulated transcripts in HbAS trophozoites (log₂ fold change = 1.55; adjusted P = 0.0077) encodes the heat shock protein 70-x (HSP70-x; PF3D7_08317000), which is thought to chaperone export of PfEMP-1 via J-dots and Maurer’s clefts (23), similar to the 2-fold overexpression of the critical Maurer’s cleft components MAHRP1 (PF3D7_1370300; adjusted P = 0.041) and MAHRP2 (PF3D7_1353200, adjusted P = 0.034).

Among the 41 downregulated transcripts in trophozoites in HbAS children, the most highly significant was that encoding a putative acid phosphatase (PF3D7_1430300; log₂ fold change = −2.37; adjusted P = 1.77 × 10⁻⁶). We also observed among the top 10 downregulated transcripts the genes for putative autophagy protein 5 (ATG5) (PF3D7_1430400; adjusted P = 0.0076), which associates with autophagosome-like structures (24), and for the DNA repair protein RAD5 (PF3D7_1343400; adjusted P = 0.0091), which has been linked to delayed clearance following artemisinin treatment (25).

Because HbAS confers protection specifically from severe malaria, we compared our data with those from a prior study comparing parasite gene expression in severe and uncomplicated malaria (26), which identified 236 differentially expressed transcripts. Seven of these transcripts were also differentially expressed between trophozoite-stage parasites in HbAS and HbAA children (Fig. 6E). Of these seven, six were transcripts overexpressed in severe malaria. One of these, gbp130 (PF3D7_1016300), was upregulated in children with HbAS. Five of the six transcripts that were upregulated in severe malaria were downregulated in children with HbAS (PF3D7_0932100, PF3D7_1203500, PF3D7_0524100, PF3D7_0919800, and PF3D7_1430300). Among the transcripts with statistically significant upregulation in severe malaria but downregulation in HbAS RBCs was that of a putative acid phosphatase (PF3D7_1430300), which was the most statistically significantly downregulated transcript in our HbAS trophozoite samples compared to HbAA trophozoites.

Intersection of aberrant transcripts in HbAS across analytic approaches.

Across analyses of both differential expression and temporal dysregulation in our time series experiments, only 15 candidate transcripts were aberrantly expressed in HbAS by 3D7 and FUP parasites (Fig. 7A). These include PTEX88, which was upregulated in both 3D7 and FUP, which is noteworthy owing to its critical role in exporting proteins from the parasitophorous vacuole membrane (PVM) that remodel the iRBC. Of these candidates, 10 overlapped candidates identified from in vivo infections in the same direction. The sole transcript differentially expressed in in vivo ring-stage parasites, cyp19B, was also downregulated in both 3D7 and FUP in HbAS (Fig. 7B to D).

FIG 7

Intersection of aberrant transcripts across analytical and experimental approaches. (A) Venn diagram of transcripts that were differentially expressed and temporally dysregulated in 3D7 and FUP in in vitro time series experiments. The 15 transcripts that were differentially expressed and temporally dysregulated in 3D7 and FUP are visualized in a heat map of log₂ fold change between HbAS and HbAA samples over 48 h. This list includes seven unannotated transcripts (PF3D7_0921900, PF3D7_1423700, PF3D7_1117300, PF3D7_1414200, PF3D7_0925900, PF3D7_1120000, and PF3D7_1310400), an exported protein (PF3D7_0220500), two rRNA processing proteins (PF3D7_1250000 and PF3D7_0818400), a translation initiation factor (PF3D7_0310600.1), a nuclear preribosomal assembly protein (PF3D7_1124800), a nucleoporin (PF3D7_0905100), a dimethyladenosine transferase (PF3D7_1415800), and the translocon component PTEX88 (PF3D7_1105600). (B) Venn diagram of transcripts that are up- and downregulated in both in vitro time series and in vivo samples from Malian children. The 10 transcripts that are shared between both experiments are visualized in a heat map, constrained to the conditions specified for panel A. (C) Log₂ fold change of transcripts over the in vitro time series of transcripts of biological significance that were included in panels A and B. (D) Normalized expression of transcripts of transcripts of biological significance in in vivo samples from panel B.

Among the genes that were significantly upregulated in vitro and in vivo was gbp130, whose expression was increased more than 30-fold in HbAS trophozoites compared to HbAA. In our in vitro time series, we found that this gene is significantly upregulated in HbAS FUP at 48 hpi, and one HbAS 3D7 replicate also demonstrated upregulation at 48 hpi. In addition, the glycophorin-binding protein homologs, gbph and gbph2, were significantly upregulated in vivo in trophozoite-stage parasite isolates and in vitro in both HbAS 3D7 and HbAS FUP. In vitro, hsp70-x (PF3D7_0831700) demonstrated a pattern of elevated transcript abundance at late time points in 3D7 and FUP, although statistical significance was achieved only in FUP (Fig. 7C and D). Additionally, mahrp1 and mahrp2 were upregulated in HbAS in vitro and in vivo (Fig. 7C and D). Among downregulated transcripts in vivo and in vitro, we observed a palmitoyltransferase (PF3D7_0609800), a thiamine-phosphate pyrophosphorylase (PF3D7_0614000), an exported protein (hyp11, PF3D7_1102900) and two unannotated proteins of unknown function (PF3D7_0316200 and PF3D7_1449200). The protein phosphatase (PF3D7_1127000) that was the most highly upregulated transcript in HbAS trophozoites was highly expressed in HbAS and was also highly expressed and slightly upregulated in HbAS FUP, although it did not pass the threshold for statistical significance. In 3D7, however, this transcript was expressed at low levels in all samples.

Finally, we examined the set of candidate aberrant transcripts identified in vivo that were similarly expressed between Hb genotypes in vitro, a set which may indicate parasite responses that are specific to physiologic conditions. These 17 in vivo-specific transcripts include upregulation of atg5 (PF3D7_1430400) and rad5 (PF3D7_1343400), suggesting that autophagy and DNA repair, as well as the differential expression of enzymes whose specific functions are currently not well understood, including the putative protein phosphatase PF3D7_112700 and the putative acid phosphatase PF3D7_1430300, are impacted by HbAS in vivo.

DISCUSSION

HbAS consistently protects against severe P. falciparum malaria by mechanisms that remain obscure. In this study, we investigated the transcriptional changes that occur in P. falciparum in response to HbAS, and overall, the parasite’s transcriptional program remains largely unperturbed in HbAS iRBCs. However, hundreds of transcripts demonstrated differential expression in HbAS, with the bulk of these changes occurring late in the IDC. Combining our in vitro experiments with parasite transcript expression in vivo, we observed impacts of HbAS on transcripts encoding the parasite’s protein chaperone and folding machinery, oxidative stress response, and protein export machinery.

Transcripts that were most clearly temporally dysregulated in HbAS iRBCs were those that are components of gene regulation: the GO biological processes with the highest DTW scores were those for transcription by RNA polymerase II, RNA splicing, regulation of translation initiation, and ribosome biogenesis (Fig. 5D). Mistimed expression of these transcripts could result in a cascade of downstream expression changes. In particular, temporal dysregulation of the subunits of RNA polymerase II, the complex responsible for the production of all parasite mRNA, could result in altered transcript abundance. Furthermore, the downstream impact of temporal dysregulation of transcripts involved in RNA splicing, ribosome biogenesis, and translation initiation could lead to significant alterations in the HbAS parasite proteome. These potential posttranscriptional effects of HbAS on the parasite would be consonant with our parallel study of the effect of HbAS on the transcription of var genes and the presentation on the RBC surface of P. falciparum erythrocyte membrane protein 1 (PfEMP1) (27). In those experiments, HbAS allows var transcription but attenuates both the presentation of PfEMP1 on the RBC surface and the adhesion of infected RBCs to extracellular ligands, suggesting that posttranscriptional impairment of var translation or transport is an important component of attenuated pathogenesis.

Differential expression in HbAS emerged primarily during the schizont stage in our in vitro time series. The molecular functions and cellular components that were most downregulated in HbAS include cysteine-type peptidase activity, including the SERA proteins and the core histones and histone variants (Fig. 3C). The SERA proteins play a role in erythrocyte rupture due to their localization to the parasitophorous vacuole (PV) in late stages of the IDC and to a papain-like cysteine peptidase domain that is common to all SERA proteins (28). The highly expressed SERA5 (29) and SERA6 are essential to blood-stage parasite growth (30, 31), and SERA5 inactivation disrupts schizont rupture and merozoite release (32). Disruption of individual SERA expression can result in compensatory overexpression of other SERAs (33), though our observation of consistent and significant downregulation of SERA2, -3, -4, -5, -6, and -7 (Fig. 4C) precludes this in HbAS and suggests that SERA-mediated egress is impaired in HbAS. IDC propagation also relies on a highly regulated epigenome (34), which is normally regulated by histones transcribed late in the asexual blood stage of development for nucleosome packaging in merozoites (35, 36). Transcripts encoding P. falciparum histones were consistently underexpressed in HbAS iRBCs (Fig. 4B), which could produce an altered chromatin landscape in the new merozoites and alter the blood stage transcriptional program after subsequent invasions. Though we did not observe large transcriptional changes in in vivo parasites from HbAS children, such an effect in vivo could be masked if such a defective transcriptional program inhibits these merozoites from successfully parasitizing an RBC. Taken together, these observations could collectively reduce the efficiency of propagation. With confirmation during consecutive in vitro cycles in HbAS RBCs of the abundance of SERA5 transcripts and protein, this could explain why HbAS is associated with lower parasite densities in vivo (37–39).

HbAS results in dysmorphic and dysfunctional Maurer’s clefts (8, 40, 41). Surprisingly, we found that transcripts related to Maurer’s clefts were among the most upregulated transcripts both in vitro and in vivo. P. falciparum exports proteins, including PfEMP-1, to the iRBC surface in knob-like protrusions that mediate its cytoadherence. Maurer’s clefts serve as the hub through which these proteins are trafficked (42). In particular, sbp1, mahrp1, and mahrp2 were upregulated in both of our in vitro strains, with the latter two also upregulated in vivo in trophozoite-stage parasites (Fig. 7D). In the asexual blood stage, sbp1 and mahrp1 are nonessential, but individual disruption of either gene leads to dysmorphic Maurer’s clefts and the absence of PfEMP-1 on the iRBC surface (43, 44). In the IDC, mahrp2 is essential, and the MAHRP2 protein localizes to electron-dense tubular structures known as tethers that connect Maurer’s clefts to the parasitophorous vacuole membrane (PVM) and the RBC membrane (45, 46). These data suggest the possibility that gene expression of Maurer’s clefts is under transcription-translation feedback control that is contingent upon the proper downstream assembly and function of its protein complexes. In this scenario, parasites in HbAS iRBCs would continue to express Maurer’s cleft genes to compensate for the dysmorphic and dysfunctional Maurer’s clefts.

Of particular interest among the six transcripts that were significantly downregulated in HbAS both in vitro and in vivo was cyp19B, which in our in vivo samples was approximately 3-fold downregulated. CYP19B localizes to the parasite cytosol and exhibits PPIase activity, participating in protein folding and potentially acting as a molecular chaperone via protein trafficking and regulation of multiprotein complexes (22). CYP19B is bound by cyclosporine (CsA), which inhibits in vitro both the PPIase activity of CYP19B (47) as well as the growth of P. falciparum (22). Interestingly, cyp19B was identified as the most upregulated gene in artemisinin-resistant (Art-R) parasites in a Southeast Asian population transcriptomic survey, in which its expression was the most positively correlated single transcript with increased parasite clearance half-life (48), where its upregulation correlated with a robust parasite stress response to the free-radical damage initiated by artemisinin. In contrast, in our experiments cyp19B was underexpressed in HbAS RBCs, despite the presumptively oxidizing and harsh environment of the HbAS RBC. CYP19B is also a member of the Plasmodium reactive oxidative stress complex (PROSC), and we observed a downregulation of a majority of the members of PROSC in HbAS in vitro at late time points in the IDC in both 3D7 and FUP, including cyp19B, hsp70-2, grp94, pdi8, and pferc, though not in vivo. Additionally, the relatively high levels of cyp19B expression in HbAA in vivo ring-stage samples were unexpected, given that both transcript expression (in our in vitro time series) and protein expression (49) are highest in the late trophozoite and schizont stages.

Among our in vivo parasites, parasites were either predominantly rings or trophozoites with no minor population of >20%, suggesting that circulating parasites in in vivo infections are broadly synchronized. This divergence accounted for most of the variance in transcriptomes, underscoring the need for stage-matched comparisons of parasite transcription. In vitro, most of the transcriptional changes we observed occurred in schizont-stage parasites. However, given that schizont-stage parasites are likely sequestered in the microvasculature and thus not in circulation, we could not confirm that similar changes occur in vivo. Finally, though we observed relatively few differentially expressed transcripts in HbAS children—particularly among the ring stages that predominated—the intersection of statistically significant differential expression between our trophozoite-stage data compared with transcriptomic data from severe malaria indicates that HbAS reduces the expression of a subset of transcripts that are associated with severe malaria. The functions of these transcripts are largely obscure, but our analyses collectively suggest their functional role in severe malaria.

Despite the expectation that HbAS would reduce microvascular cytoadherence and thereby produce more circulating trophozoites, we observed similar proportions of trophozoite-predominant infections in HbAS (7/16) and HbAA children (5/16) (Fig. 6A and B) and no evidence of late trophozoite or schizont gene expression in our HbAS patient samples (Fig. 6A). This may have been influenced by our sampling scheme, and since parasites are thought to be somewhat synchronous in natural infections in prior studies (50) and in our data, we would expect to find only ring-stage and early trophozoite parasites in our patients.

It is surprising that we observed conservation of the overall order of the transcriptional program in the IDC despite the presence of substantial temporal dysregulation in HbAS of gene regulatory GO terms such as RNA polymerase II and translation initiation. Interestingly, our in vitro time series data in HbAS RBCs mirrored key aspects of the transcriptional changes observed in P. falciparum lines with mutations of the drug resistance transporter genes pfmdr-1 and pfcrt (51), with an overall transcriptomic concordance followed by late-stage differential expression. This suggests the possibility that P. falciparum’s transcriptional program in the IDC is sufficiently robust to adapt to fitness challenges which primarily manifest by modulating the expression levels of particular genes while still producing the necessary transcripts on time to complete the IDC. The maintenance of the parasite’s transcriptional program relies on a robust architecture that is only partially understood, highlighted by recent reports that Plasmodium spp. possess intrinsic oscillators that govern the periodicity of their IDC program (52, 53). The alterations we observed in transcript abundances may be a consequence of the temporal dysregulation of gene regulatory products we discovered with DTW analysis. Such subtle alterations could, with successive rounds of infection, be amplified into fitness-mediating phenotypic changes. A need for successive rounds of infection to manifest transcriptional divergence in HbAS RBCs could also explain the paucity of overlap between the lists of differentially expressed transcripts in vitro and in vivo. Alternate explanations are the altered fitness landscape in the bloodstream of partially immune hosts, the sequestration in vivo of the most transcriptionally divergent schizont stages, or an interaction in vivo between HbAS and parasite genotype.

Our study had several limitations. Our in vitro time series experiments were conducted for one round of infection over 48 h in order to ensure tight synchrony of parasites, and therefore we could not observe impacts of perturbations in histones and chromatin modifications on subsequent infections. Our in vitro experiments were performed in single gas mixture at 1% O₂, while in vivo, P. falciparum encounters in circulation a range of oxygen tension that can differentially impact its growth in HbAS iRBCs in vitro (3–6). Though prior studies have reported a growth defect in HbAS RBCs at 1% O₂, we did not observe such a defect over one IDC.

Our data encompass P. falciparum’s transcriptional response to infection in HbAS RBCs, highlighting multiple avenues by which HbAS may neutralize P. falciparum. The potential neutralizing effect of these aberrations will require testing in systems of pathogenic phenotypes, including extracellular adhesion, endothelial activation, and immune evasion. These transcriptional perturbations induced in HbAS provide new targets to neutralize parasite mediators of disease and suggest that downstream changes in the proteome and epigenome are candidates for exploration.

MATERIALS AND METHODS

Subjects and samples.

HbAA and HbAS blood samples were collected under institutional review board (IRB) number Pro00007816 for in vitro experiments. Genotypes were confirmed by Sanger sequencing. Field samples were collected in an observational study of the efficacy of artemether-lumefantrine for the treatment of uncomplicated malaria in children in Kéniéroba, Mali (ClinicalTrials.gov identifier NCT02645604), where HbAS reduces the risk of uncomplicated malaria by 34% (54). Inclusion criteria were age between 2 and 17 years and febrile uncomplicated falciparum malaria as assessed by light microscopy. From all participants, prior to treatment, we collected venous blood, passed this through cellulose columns (55), and stored up to 2 ml of the flowthrough in RNAprotect (Qiagen). Parasite density was estimated from light microscopy, hemoglobin genotype was assessed using high-performance liquid chromatography (HPLC), and ABO blood group was assessed by agglutination assay. For transcriptional analysis, we selected all available samples from children with HbAS and matched each of these to a sample from a child with HbAA with regard to age (within 1 year), month of episode, parasite density, ethnic background, and, if possible, ABO blood type. The field study was approved by the Institutional Review Board of the University of Sciences, Techniques, and Technologies of Bamako (IRB00001983).

Parasite culture.

P. falciparum parasites of the 3D7 and FUP strains were cultured in human RBCs in RPMI 1640 supplemented with Albumax II, l-glutamine, hypoxanthine, and d-glucose at a 2% hematocrit in 1% O₂–5% CO₂ in N₂ at 37°C (56). iRBCs were kept in filtered-lid cell culture flasks inside a modular hypoxia incubator chamber and mixed on a nutator. P. falciparum parasites from either strain were synchronized using Percoll centrifugation, and the schizont fraction was then split equally into separate flasks of HbAS or HbAA. Following a 3-h incubation, these cultures were treated with sorbitol and allowed to progress through the intraerythrocytic cycle in parallel, with sampling every 3 h for light microscopy and RNA preservation. Parasite density and maturation were assessed following Giemsa staining by readers masked to hemoglobin type and hours postinvasion.

Library preparation.

Total RNA was isolated using TRIzol and the RNeasy minikit (Qiagen) after DNase treatment and quantified by a Qubit high-sensitivity RNA assay (Thermo Fisher Scientific) before storage at −80°C. Libraries were prepared from total RNA with the Kapa stranded mRNA-seq library prep kit. Libraries from each of the two time series experiments (n = 64/experiment) were sequenced on a full flow cell on the NovaSeq 6000 S2 platform with 150-bp paired-end reads. Field study samples (n = 32) were sequenced on a full flow cell of the NovaSeq 6000 S1 platform with 50-bp paired-end reads.

Transcriptome quantification.

Reads were assessed for quality with FastQC, trimmed and quality-filtered with Trimmomatic (57), depleted of reads mapping to the human genome (GRCh38.p13) with STAR (58) and mapped and quantified with Salmon (59) using the Plasmodium falciparum 3D7 transcriptome (ASM276v2). Quantification files of parasite gene expression were summarized with the tximport package in R (60). Principal-component analysis data were generated in R with the plotPCA function from DESeq2 and then plotted with ggplot2. For heat map visualization, transcripts per million (TPM) were normalized around the mean for each transcript over the time series and depicted as a z score of standard deviations from the mean.

Peak shift analysis.

Analyses of the shift of transcript peak were conducted by first identifying for each transcript in each time series the time point at which the transcript achieved its maximum expression and then computing for each transcript the difference between time series of the transcript’s peak. We repeated this approach using only transcripts with a single peak over the time series, which were identified using the find_peaks function from the signal processing submodule in SciPy. We compared the distributions of peak shift values between samples using Kruskal-Wallis one-way analysis of variance.

Differential expression analysis.

All differential expression analyses were performed with the DESeq2 (61) R package using unnormalized read counts per transcript. Transcripts were considered differentially expressed if, as determined by DESeq2 with a false-discovery rate set to 0.05, their adjusted P value was <0.05 and |log₂(HbAS/HbAA)| was >log₂(1.5) in the time series and, with a false-discovery rate set to 0.1, the adjusted P value was <0.05 and |log₂(HbAS/HbAA)| was >log₂(1.5) in the field samples.

Gene ontology enrichment analysis.

Gene ontology enrichment analysis was performed with GOATOOLS (17). Gene IDs and associated GO IDs were downloaded from NCBI (ftp://ftp.ncbi.nlm.nih.gov/gene/DATA/gene2go.gz), and P. falciparum’s annotations were extracted with the taxonomic identifier 36329. Using the transcripts identified as differentially expressed by DESeq2 according to the conditions described above as input, gene ontology terms were accepted for Benjamini-Hochberg-adjusted P values of <0.05.

Analyses of transcript temporal regulation.

Dynamic time warping (DTW) of time series transcript expression data was computed with the tslearn toolkit (62). We considered transcripts with mean DTW scores across β-globin replicates that exceeded two standard deviations above the mean DTW scores of all transcripts in HbAA versus HbAS comparisons temporally dysregulated.

We used the temporal alignment Kendall-Tau algorithm (TAKT) (18) to identify transcripts within each parasite that had similar dynamics but were temporally shifted between β-globin types during the ring stage (18). Measures of similarity and associated significance scores were produced for a gene between similar β-globin type replicates and a P value of ≤0.05 was used to identify transcripts with similar dynamics across all β-globin types within each parasite.

Stage classification of in vivo parasites.

We identified the set of transcripts shared between 3D7 and FUP from our in vitro time series that peaked within rings, trophozoites, or schizonts and then selected transcripts that were highly expressed at their peak, discarding those that did not exceed 1,000 TPM. We also estimated parasite stage with a previously published mixture model (20). Our count data were converted to reads per kilobase per million (RPKM) and log₂ transformed to match the staged expression data (63) used for stage estimation.

Data availability.

The RNA-seq data are available in GEO (accession no. GSE163144; BioProject no. PRJNA685106). All scripts, protocols, analyses, and supplemental output are available on GitHub (64). See Supplemental Text S1 for more experimental, processing, and analytic detail.