Insufficient compensation for anaemia during acute cynomolgi malaria in rhesus macaques
The parasitological and haematological data used in this analysis were published previously and made publicly available in Joyner et al. [9]. Briefly, five rhesus macaques were infected with P. cynomolgi M/B strain sporozoites and followed for up to 100 days. One macaque developed severe disease and irreversible kidney failure that required euthanasia; thus, this animal was removed from the current analysis [9]. Seven BM samples were collected longitudinally from each macaque, and these samples were classified into the following infection points for analysis in this manuscript: prior to inoculation (pre-infection), during acute primary infection, after the peak of parasitaemia (i.e. post-peak), before and between relapses (i.e. inter-relapses), and during relapses (Fig. 1a).
The rhesus macaques infected with P. cynomolgi M/B strain developed varying levels of anaemia during the primary infection (Fig. 1b; [9]). Despite the significant decrease in haemoglobin levels during acute infections, the circulating reticulocytes did not increase significantly until after the peak of parasitaemia or after sub-curative blood-stage treatment had been administered (Fig. 1b). These results suggested that the BM was not appropriately responding to the anaemia during acute infections but was restored after the peak of parasitaemia. Interestingly, relapses did not result in changes in haemoglobin levels or peripheral reticulocyte counts (Fig. 1b).
To assess whether each animal was appropriately compensating for anaemia, the reticulocyte production index (RPI) was calculated as previously described [9]. The RPI adjusts for alterations in reticulocyte maturation time in an anemic individual and provides a metric to assess if an appropriate compensatory response is being made by the BM to combat anaemia. An RPI value of under 2 when haemoglobin levels are decreasing suggests that the BM is not mounting an appropriate compensatory response for anaemia whereas an RPI of greater than 3 when haemoglobin is decreasing indicates an appropriate compensation. In support of the haemoglobin and reticulocyte data, there was not an appropriate compensatory response by the BM during acute infections (Fig. 1b). However, an appropriate response was mounted after the peak of parasitaemia and treatment, though haemoglobin levels continued to decrease (Fig. 1b).
Erythropoietin (EPO) is critical for inducing erythropoiesis during normal RBC turnover and pathological conditions. Thus, the lack of response by the BM during acute infection could have been due to dysregulation of erythropoietin production. In contrast to this hypothesis, EPO levels were significantly increased during acute infection (Fig. 1c). Collectively, these analyses demonstrated that there was not an appropriate compensatory response by the BM in the face of the developing anaemia early during acute infections.
Acute malaria, but not relapses, leads to substantial changes in the bone marrow transcriptome
The next goals of this study were to analyse the transcriptional changes in the BM that may account for the lack of compensatory response by the BM during acute infection and to determine if relapses induced changes in the BM since anaemia was not observed during relapses. RNA-Seq was performed using BM MNCs from the four rhesus macaques infected with P. cynomolgi at the time points indicated in Fig. 1a and classified for downstream analysis as described above.
Hierarchical clustering was performed with the R mixOmics package [9] to determine the similarity and/or differences in the BM transcriptome among the animals for each infection point. Using this approach, the acute infection measurements clustered together, but no other infection points clustered with a discernable pattern, suggesting that factors other than the BM transcriptome were influencing the grouping of these stages (Additional file 1A). After determining that 26.7% of the gene expression variance was attributable to each individual, the variance associated with each individual was removed using a supervised normalization of microarray (SNM) transformation (Additional file 1B).
After transformation, 42.3% of the variance of the dataset was now attributable to the infection point (Additional file 1C). Hierarchical clustering of the data after accounting for the variance attributable to individual animals maintained clustering of the acute infection stage (Fig. 2a). However, this also led to a more robust clustering of the other infection stages, although the clusters were not completely correlated with infection points (Fig. 2a). This suggested that acute malaria caused significant and coherent changes in the BM transcriptome, but there did not appear to be unique changes during relapses, post-peak, or inter-relapse infection periods.
Since hierarchical clustering indicated that transcriptional profiles at post-peak and inter-relapse infection points were generally similar to pre-infection profiles, similarities and differences were formally evaluated between each infection point. The number of differentially expressed genes (DEGs) between infection points was determined using analysis of variance (ANOVA). Compared to pre-infection, 1266 genes were differentially expressed during acute infection, 315 after the peak of parasitaemia, and only 201 genes during relapses (Fig. 2b). For relapses, 75 of these genes were in common with acute infection (Fig. 2c). During acute infection, 1148 of the DEGs were upregulated whereas 118 were downregulated compared to pre-infection. Relapses showed much less variance in gene expression even for genes with large absolute upregulation and high significance during acute infection. Thus, while acute infection caused significant perturbations in the BM transcriptome, relapses had a small effect.
Pathways and processes altered in the bone marrow during acute malaria
Gene ontology analysis was employed to understand the pathways and processes altered during acute infection. Approximately 37 biological processes were significantly changed in the BM during acute infection (Additional file 2). The identified pathways could largely be consolidated into biological processes related to cellular stress (e.g. response to endoplasmic reticulum stress), increases in transcription (e.g. RNA processing), and ongoing immune responses in the BM (e.g. Type I IFN, IFNγ, cell response to cytokine signaling, etc.). The pathways and processes of the other infection points were also analysed but none were significantly enriched; this is likely due to the low number of differentially expressed genes during these infection stages.
Pathway analysis using MetaCore (Thomson Reuters) and Molecular Signatures Database (MSigDB) [29, 30] identified additional molecular pathways that were increased during acute infection (Fig. 3). The pathways identified were similar between both databases and were skewed towards upregulation (Fig. 3a–d). This is likely due to the low number of down-regulated DEGs during acute infection (Fig. 2b). Full lists of pathways identified in each database are provided in Additional files 3, 4, 5 and 6.
The enriched pathways and gene sets from both databases highlighted the ongoing immune response in the BM during acute infection; in particular, pathways involved in cytokine signaling were overrepresented (Fig. 3a, b). The most predominant cytokine signatures to emerge from the analysis of acute infection in both databases were Type I (i.e. IFN-α/β) and Type II interferon (i.e. IFN-γ), which were identified by both databases (Fig. 3a, b). A heatmap of DEGs identified in these pathways shows upregulation during acute infection but not other infection points (Fig. 3e, f). Although the related pathway was identified, changes in expression of the ifnα and ifnβ1 genes themselves were not captured by RNA-Seq during acute infection because the genes did not pass the FPKM threshold. ifnγ expression, on the other hand, was detectable and significantly increased (p < 0.001) during acute infection and after the peak of parasitaemia (p < 0.05) (Fig. 3g). To ensure that gene expression was indicative of the protein levels, the concentrations of IFNγ and IFNα in the plasma were measured via multiplex assay. Consistent with the transcriptional data, IFNγ levels were increased (p = 0.008; Fig. 3i) and IFNα levels did not change during acute infection although a significant decrease was detected after the peak of parasitaemia (Fig. 3h).
Other cytokine signatures identified as transcriptionally upregulated in the BM during acute infection included IL-10 and IL-27, which are known to have effects on the BM during acute malaria, and genes involved in the IL-3 signalling pathway (Fig. 3b) [31, 32]. Interestingly, the IL-5 signalling pathway was the only cytokine signalling pathway identified as being downregulated during acute infection, which could be due to the skewing towards a Th1 response (Fig. 3d).
This analysis also revealed the enrichment of pathways associated with pathogen recognition receptors (PRRs) that are important for sensing intracellular and extracellular microbial components as well as host danger molecules. Specifically, the analysis highlighted the importance of engagement of Toll-like Receptors, NOD-Like receptors and RIG-I/MDA5 (Fig. 3a, b).
Intermediate and non-classical monocytes may negatively impact the erythroid lineage during acute malaria
Bone marrow is made up of a heterogeneous population of cells that together yield the transcriptional profiles measured here. Weighted Gene Set Co-expression Network Analysis (WGCNA) was employed to investigate which cells types in the BM may potentially be responsible for the observed changes in the BM transcriptome during acute malaria [33]. WGCNA was performed using DEGs identified during the acute infection and flow cytometry data collected at each infection point. The gating strategy for the cellular subsets used in WGCNA analysis is shown in Additional file 7. The populations included in the analysis were as follows: classical monocytes, intermediate monocytes, non-classical monocytes, B-cells, T-cells, granulocytes and a mixed population of erythroid progenitors (Additional file 8). Samples from other infection points were not included in this analysis to allow the acute infection to drive the formation of the coexpression modules.
The dominant modes of variance of each module were used to evaluate the association of the modules with cellular subsets. WGCNA identified four gene coexpression modules (Fig. 4a). The modules coloured turquoise and blue were composed of 490 and 335 DEGs, respectively, and were positively correlated with the number of intermediate and non-classical monocytes (p < 0.05) in the BM. The turquoise module was negatively correlated with B cell numbers (r = −0.52, p = 0.01). However, unlike monocytes, this was a negative correlation and is likely related to the disruption of the BM environment where B-cells develop. Other immune cell types identified by flow cytometry were not correlated with any of these co-expression modules. These results suggest that changes in the BM transcriptome during acute malaria may have been primarily related to monocytes.
To determine whether there was a link between the identified co-expression modules and suppression of the erythroid progenitor population, the relationship of the co-expression modules with the erythroid progenitor population was examined. All four co-expression modules were correlated with a drop in erythroid progenitors (p < 0.05 in all cases). Since only the turquoise and blue modules were associated with immune cells (monocytes), this suggested that the correlation of the yellow and brown co-expression modules represented non-immunological or diffuse immunological pathways.
Intermediate and non-classical monocytes are associated with pathways upregulated during acute malaria in the bone marrow
To test the hypothesis that monocytes were a major driving force of the transcriptome changes in the BM, pathway enrichment analysis was first performed on the genes in each module [29, 30]. The MSigDB database was used for this analysis; results obtained with this database were comparable to those obtained using MetaCore (Fig. 4b–e). The turquoise module was enriched with pathways involved in the immune response (and interferon signalling pathways), general cellular response pathways (e.g. unfolded protein response), and RNA metabolism (Fig. 4b). The blue module was enriched in pathways related to the immune response, including cytokine signalling by IFNγ, activation of NF-κB, and the IL-2 pathway (Fig. 4c). Each of these pathways are known to be involved in monocyte activation [34]. This is consistent with the hypothesis that the activation status of the intermediate and non-classical monocyte populations in the BM during acute infection may have influenced the BM transcriptome. Indeed, the intermediate and non-classical monocyte populations were expanded in the marrow during acute infection, but these differences did not reach statistical significance likely due to the sample size (Additional file 9).
The blue module was also enriched for the PRR pathways activated through NOD like receptors (NLRs) and RIG-I like receptors (RLRs) (Fig. 4c). This suggests that these pathways, which were originally identified by DEG analysis, may be active in intermediate and non-classical monocytes in the BM in response to Plasmodium infection.
To validate the transcriptomic signatures related to cytokine signalling and identify the cytokines that may negatively impact the erythroid progenitors, WGCNA analysis was repeated using the flow cytometry, transcriptomics and cytokine data. Seventeen of 45 cytokine concentrations were found to be positively correlated with the turquoise and blue modules. Transitively, these cytokines were thus also negatively associated with the erythroid progenitor population (Fig. 5). Although negatively correlated with the erythroid population in the BM, these cytokines were positively associated with the frequency of non-classical and intermediate monocytes (Fig. 5). IFNγ was correlated with the blue, brown, and turquoise modules, and IFN signalling was enriched in the turquoise and blue modules, which provided confirmatory support for the transcriptome analysis and the validity of the WGCNA approach (Fig. 4a, c).
The brown module did not correlate with particular immune cell types, and this is likely due to the diverse set of pathways represented in this module (Fig. 4d). The yellow module highlighted IFN signalling pathways (Fig. 4e) but was not correlated with protein levels of IFN. Interestingly, a signature driven by erythropoietin (EPO) was also enriched in the yellow module (Additional file 10). Although an EPO signature is to be expected given that EPO plays an important role to increase BM output of RBCs during anaemia, this signalling was inversely correlated in this study with the erythroid progenitors, suggesting that despite the elevated EPO levels in the plasma and signalling in the BM, an appropriate compensatory response could not be made. Indeed, these data are consistent with previous literature showing that malaria induces an EPO-independent anaemia [35].
Transcriptional networks related to erythropoiesis are disrupted in the bone marrow during acute malaria
Given insufficient compensation for anaemia during acute malaria, the transcriptome data was mined to determine if transcriptional networks related to erythropoiesis were disrupted during acute P. cynomolgi infections. Correlation analysis (CA) was performed using the erythroid population frequencies determined by flow cytometry and BM MNC gene expression to identify genes that may be related to the decrease in erythroid progenitors observed in the BM (Fig. 6b). Samples from pre-infection, acute primary, and post-peak were utilized for this analysis to determine changes that occurred before, during, and after the onset of anaemia.
547 genes from the same RNA-Seq dataset utilized in the previous analyses were identified as positively correlated (Spearman’s correlation coefficient >0.6; p < 0.05) at FDR < 0.2 with the frequency of the erythroid progenitor population in the BM (Additional file 11). After identifying the related genes, the iRegulon pipeline was utilized to explore what transcription factors may be influencing the genes that were positively correlated with the change in erythroid progenitor cell frequencies during acute infection. iRegulon identified multiple transcription factors including CCAAT/Enhancer Binding Protein (CEBPD), GATA Binding Proteins (GATA1/GATA2) and TEA Domain Transcription Factor 4 (TEAD4) as having their targets enriched in the genes that correlated with erythroid progenitor cell frequencies (Fig. 6a).
This analysis focused on GATA1 and GATA2 because these proteins are master regulators of erythropoiesis, had the largest number of target genes identified in the dataset, and had two of the highest NES scores, which are a metric of statistical significance in iRegulon (Fig. 6a) [36, 37]. First, the expression of GATA1 and GATA2 was examined directly to understand if the gene expression of these proteins was changed during infection, which could directly impact the erythropoietic response during infection. GATA2 gene expression did not change throughout the initial infection (Fig. 6d). In contrast, GATA1 gene expression was not changed during acute infection but was significantly increased after the peak of parasitaemia when erythropoiesis was restored (Fig. 6c). Thus, transcription of GATA1/GATA2 did not appear to completely or directly explain the correlation between erythroid progenitors and hundreds of genes that are their targets, though the increases in GATA1 and progenitor levels at post-peak suggest GATA1 transcription could still be involved.
Next, the genes regulated by GATA1 and GATA2 were examined to identify whether the genes regulated by these proteins were differentially regulated even though transcription of GATA1/GATA2 was not altered significantly during acute infection. The target genes of GATA1 and GATA2 correlated with the small but significant decrease (p = 0.045) in the frequency of the erythroid progenitor population in the BM during acute infection (Fig. 6b). Notably, the expansion of other cell types in the marrow during acute infection could drive this decrease in percentage of erythroid progenitors. Importantly, the genes regulated by GATA1 and GATA2 were downregulated during acute infection, but upregulated post-peak when the erythroid progenitor population in the BM was restored to pre-infection levels and an increase in peripheral reticulocytes was observed, indicating restored erythropoiesis (Fig. 6e, f). Therefore, direct suppression of GATA1/GATA2 transcription does not appear to be responsible for the lack of compensatory erythropoiesis in the periphery, and instead, it seems that these proteins’ functions may be compromised at the post-translational level during acute infection.