Correlates of severity of disease in Macaca mulatta infected with Plasmodium cynomolgi

TL;DR

This study investigates the immune and transcriptomic responses in Macaca mulatta infected with Plasmodium cynomolgi, identifying molecular signatures associated with severe malaria to better understand disease severity.

Contribution

It provides an integrative analysis linking gene expression, cytokine profiles, and immune responses to disease severity in primate malaria models.

Findings

Severe infections show up-regulation of PAMP and cytokine response genes.

Distinct transcriptional profiles are associated with severe versus mild malaria.

Molecular signatures identified may inform understanding of severe malaria mechanisms.

Abstract

Characterization of host responses associated with severe malaria through an integrative approach is necessary to understand the dynamics of a \textit{Plasmodium cynomolgi} infection. In this study, we conducted temporal immune profiling, cytokine profiling and transcriptomic analysis of five \textit{Macaca mulatta} infected with \textit{P. cynomolgi}. This experiment resulted in two severe infections, and two mild infections. Our analysis reveals that differential transcriptional up-regulation of genes linked with response to pathogen-associated molecular pattern (PAMP) and pro-inflammatory cytokines is characteristic of hosts experiencing severe malaria. Furthermore, our analysis discovered associations of transcriptional differential regulation unique to severe hosts with specific cellular and cytokine responses. The combined data provide a molecular and cellular basis for the development of severe malaria during \textit{P. cynomolgi} infection.

Figures

Full text

1 Introduction

Malaria remains a public health challenge, responsible for approximately 400,000 (236,000–635,000) death in 2015 [1]. Out of the five human Plasmodium species capable of causing malaria, Plasmodium falciparum and Plasmodium vivax account for the majority of human malaria infections. P. falciparum is responsible for the majority of malaria-related mortality and is most prevalent in sub-Saharan Africa [1]. Although P. vivax infection does not typically result in mortality, this parasite causes substantial morbidity outside of sub-Saharan Africa and exhibits the propensity to cause severe disease. The mechanisms that underly vivax malaria pathogenesis, however, remain poorly understood partially due to experimental constraints such as the lack of an in vitro culture system and a rodent model of P. vivax infection [2, 3, 4].

Plasmodium cynomolgi is a non-human primate parasite that infects old world monkeys and is capable of recapitulating clinical and histopathological findings of vivax malaria patients[5, 6, 7]. It is both genetically and physiologically similar to P. vivax [7, 8]. For instance, both parasites exhibit 48-h erythrocytic cycle during blood stage infection [9], preferential infection of reticulocytes [10] and form hypnozoites, which are dormant, liver stage forms that can activate and cause relapse infections [11]. Due to the difficulty of studying P. vivax pathogenesis, P. cynomolgi infection of rhesus monkeys (Macaca mulatta ) has being used to better understand hypnozoite-caused relapse [12].

Host clearance of malaria parasites without complication requires a concerted effort between inflammatory and anti-inflammatory cytokines; their balance and timing are critical in determining clinical outcomes [13]. The association between cytokine, transcriptomic and immune response, and clinical outcomes has been extensively studied in P. falciparum infections [14, 15, 16, 17, 18], in comparison, much less is known for P. vivax infections.

With the aim to better characterizing P. vivax infection in humans, a time series experiment where four Macaca mulatta were infected with P. cynomolgi was conducted as a part of the Malaria Pathogen Host Interaction Center (MaHPIC) project [19]. This experiment captured host transcriptomic, cellular and cytokine response to P. cynomolgi infection. The subjects within this study responded to the infection in different manners and resulted in two cases of severe malaria, and two cases of mild malaria. The presence of both severe malaria and mild malaria provided us the opportunity to characterize host responses associated with clinical outcomes.

Here we present the combined analysis of cellular, transcriptome and cytokine data collected during this experiment with an emphasis on comparing host responses during severe malaria and mild malaria. Our analysis shows associations of severe malaria with differential up-regulation of pro-inflammatory signals independent of adaptive immune response, whereas moderate malaria is associated with transcriptomic up-regulation of complement and heme metabolism related genes. These distinct characteristics demonstrate that P. cynomolgi induced severe malaria is strongly associated with up-regulation of inflammatory signaling through transcriptional regulation.

2 Results

2.1 Study Design

The experimental design is shown in Figure 1. Transcriptomic, cellular and cytokine data were measured from peripheral blood samples collected during six time points. Time point 1 (TP1) corresponds to the pre-infection baseline for each subject. Time point 2 (TP2) corresponds to the acute primary infection when parasite counts peaked. At time point 2, two subjects, S1 and S2 developed signs of severe malaria including anemia and thrombocytopenia, which required each animal to be sub-curatively treated to prevent possible severe complications. The other two subjects M1 and M2 experienced mild disease manifestations and did not receive sub-curative treatment. All animals received curative blood-stage treatment with artemether on Day 28 to ensure that subsequent parasitemias were due to relapses and not recrudescing parasitemias. A detailed description of the experimental set-up can be found in the publication by Joyner et al. [19].

2.2 Differential Gene Regulation during Primary Infection

To characterize the differences of transcriptomic responses induced by P. cynomolgi infection between the subjects experiencing severe malaria (severe hosts) and the ones experiencing mild malaria (mild hosts), we conducted RNA-Seq analysis for each clinical group between TP1 and TP2 using the DESeq algorithm [20]. Furthermore, between-group transcriptome analysis was also carried out for TP2. The transcriptome profiles differ the most from the healthy baseline (TP1) at TP2 and TP3 shown by principal component analysis (Figure 2). Additionally, the transcriptome profiles of the severe hosts were more similar to each other than to the mild hosts at TP2 as shown by both principal component analysis and hierarchical clustering using the top 50$\%$ most variably expressed genes (Figure 2).

Mild hosts showed similar degrees of transcriptome perturbation during the onset of primary infection (698 DEGs using an absolute fold-change threshold of $>$ 1.5 and q-value $<$ 0.05) relative to severe hosts (511 DEGs). 148 genes were identified by both analyses. At TP2, 1337 genes were differentially expressed between the mild hosts and severe hosts. 380 of these DEGs were temporally differentially expressed in severe hosts between TP1 and TP2. 220 DEGs were temporally differentially expressed in mild hosts between TP1 and TP2.

To identify differentially regulated transcriptome responses to P. cynomolgi between severe hosts and mild hosts, DEGs were classified into four groups: differentially up-regulated genes in severe hosts (DUGs-S), differentially up-regulated genes in mild hosts (DUGs-M), differentially down-regulated genes in severe hosts (DDGs-S) and differentially down-regulated genes in mild hosts (DDGs-M). Genes were classified into each of these groups based on their between group differential expression at TP2 and their with-in group differential expression between TP1 and TP2. The specific classification scheme is shown in Figure 3 and Figure 4.

2.3 Biological Pathways Perturbed During Primary Infection

To identify the biological pathways associated with each of the four groups of differentially regulated genes, we used MySigDB to conduct enrichment analysis [21]. Using a q-value cut-off of 0.05, DUGs-S are enriched in 205 pathways with the most significantly enriched pathway being cytokine signaling in the immune system. DDGs-S are enriched in only 24 pathways (Figure 3) with the most significantly enriched pathway being immunoregulatory interactions between lymphocyte cells (LCS) and non-lymphocyte cells (Non - LCS). DUGs-M is enriched in 51 pathways with the most significantly enriched pathway being heme metabolism, and DDGs-M are enriched in 6 pathways. The number of enriched pathways within each gene set is proportional to the gene set size. Both DUGs-S and DUGs-M show enrichment in innate immune-related response and interferon gamma related response. However, complement and coagulation-related genes are differentially up-regulated in the mild hosts and differentially down-regulated in the severe hosts.

2.4 Differentially Regulated Genes Demonstrate Distinct Correlation Pattern with Cell Populations and Cytokine Abundance

To identify the similarities among temporal profiles of DUGs-S, DUGs-M, cytokine abundances and cell populations, we constructed two correlation matrices and identified 5 clusters (Figure 5). DUGs-S and DUGs-M exhibited distinct correlation patterns with respect to cytokine abundances and cell populations. Two clusters were identified for DUGs-S, all but 11 genes were clustered into C1 and the rest into C2. C1 contained the majority of cytokine abundances along with a few cell populations such as the percentage of CD14 CD16 monocytes out of total monocyte, the percentage of granulocytes out of total lymphocytes and parasite concentration. Furthermore, we performed enrichment analysis of transcriptional factors for C1 and identified a list of enriched transcriptional factor targets. The top 3 transcriptional factor targets enriched are ISRE, IRF7, and STAT5. Incidentally, IRF7 was also observed to be transcriptionally differentially up-regulated in severe hosts. Unlike C1, C2 contained most of the adaptive and innate immune cell populations but very few genes. In fact, the majority ( $>60\%$) of the correlations between the genes in C1 and the adaptive immune cell populations in C2 are negative.

Three clusters were identified for DUGs-M. C3 contained 37 genes and the majority of the pro-inflammatory cytokines along with reticulocyte population and parasite concentration. Enrichment analysis of C3 reveals that it is enriched in heme metabolism related genes. C4, on the other hand, contained the rest of DUGs-M along with caspase 3 CD 8 and CD 4 T cell populations, monocyte and granulocyte populations. C4 is enriched in complement and coagulation-related genes. C5 is composed entirely of adaptive immune cell types such as B cell and traditional T cell populations. Transcription factor analysis of C3 and C4 did not reveal any enrichment of transcription factor targets. Despite the fact that both heme metabolism and complement-related genes are differentially up-regulated at TP2 in mild hosts, their expression patterns are associated with different cell populations and cytokine abundances. Furthermore, their temporal expressions differ over the entire experiment, in contrast to the uniform temporal expression patterns of DUGs-S.

2.5 Cytokine Perturbation during Primary Infection

To characterize the differences of cytokine responses induced by P. cynomolgi infection between the clinical groups, two-way ANOVA was used to identify significantly changed cytokines. 4 cytokines out of all 44 cytokines measured displayed significantly different abundances between the two clinical groups. All 4 cytokines have a higher abundance in severe hosts at TP2 (Figure 6).

To determine whether the differences in cytokine abundances can be attributed to transcription levels, Pearson’s correlation coefficients between the time series of each of the cytokines and their respective transcripts were calculated (Figure 6, 7). The correlations of IL-6, MIG, MIP-1 and SICAM1 (soluble ICAM1) to their transcripts are 0.94, 0.77, 0.83 and 0.48 respectively. The correlation between SICAM1 and its transcript within the severe hosts is 0.59, which is significantly higher than that of the mild hosts (correlation of 0.17). Additionally, MIG, SICAM1, and MIP-1 only differ in cytokine abundances but not transcriptional abundances. The transcription regulators of IL-6, MIG, MIP-1 and SICAM1 were identified using the TRRUST database [22] and analyzed for differential up-regulation in severe hosts. One transcription factor of IL-6, JUN, and one transcription factor of ICAM1, TWIST2, were identified as differentially up-regulated in severe hosts at TP2.

Interestingly, 13 cytokines exhibited negative correlation between their protein abundances and transcriptional abundances. One of these 13 cytokines, IL-2 is transcriptionally differentially up-regulated in mild hosts at TP2. Additionally, all the genes belonging to DUGs-S showed positive correlations to their protein abundances. MCP1 and IL1RA exhibit the highest correlation between transcript abundances and protein abundances. Both IL1-beta and IL1RA showed differential up-regulation in severe hosts at the transcriptional level, however, their transcript-protein correlation differs. IL-1-beta showed a 5 fold transcriptional change which is associated with a 4 fold increase protein abundance. IL-1RA showed 7 fold increase in transcripts but $\sim 8000$ fold change in protein abundance. The ratio of IL1 to IL1RA is lower in severe hosts (0.0045) in comparison to mild hosts (0.04).

3 Discussion

In this study, we analyzed temporal transcriptomic, cellular and cytokine profiles of four M. mulatta infected with P. cynomolgi to characterize the host response to infection. We leveraged the fact that two of the subjects experienced severe malaria during the primary infection and the other two subjects experienced mild malaria to understand the association of host responses with P. cynomolgi disease severity. The diverse data types allow us to cross-validate transcriptional response with actual protein abundance in relation to cellular populations.

The identification of differentially up-regulated cytokines and transcripts in the severe hosts combined with the downstream pathway analysis demonstrate that severe malaria within our experiment is associated with elevated pro-inflammatory cytokine levels (MIP, IL6, MIG and SICAM1) and the up-regulation of innate immunity and anti-viral responses related gene sets (TLR Pathway, NOD Pathway, and RIG-I Pathway). Additionally, the positive correlations between the abundance of MIP, IL6, MIG and SICAM1 and their respective transcripts suggest that the differential up-regulation of these four cytokines are associated with transcriptional up-regulation. On the other hand, only two transcription factors that activate the production of these cytokines have been identified to be differentially up-regulated. The lack of detection of other transcription factors associated with these cytokines could be attributed to the lack of statistical power due to small sample size, or that the up-regulation of other transcription factors happened prior to our sampling. Additionally, most of the transcriptional regulators that control the expression of the acute phase inflammatory cytokines identified here are controlled by post-translational modifications, and thus, an increase in transcription is not alway observed even though signaling may be occurring. Our characterization of severe malaria caused by malaria infection through elevated IL-6 levels and up-regulation of innate and anti-viral related response have also been observed in studies of P. falciparum infection severity [15, 23, 14, 24], suggesting likely shared mechanism of recognition of foreign dsRNA and DNA material by the innate immune cells. Correlation analysis of DUGs-S, cell population, and cytokine abundance demonstrates that differential up-regulation of innate immunity and anti-viral responses related gene sets are accompanied by the increase in pro-inflammatory cytokine abundances. The high correlation between cytokine abundance and DUGs-s suggest that the increase in pro-inflammatory cytokines are not independent of transcriptional up-regulation. The negative correlations between DUGs-S and adaptive immune cell types suggest possible innate immune suppression by the adaptive immune system [25, 26] or cross-talk between innate immune response and adaptive immune response.

On the other hand, the mild hosts are characterized by the differential up-regulation of the complement, coagulation and heme metabolism related genes. Interestingly, the complement and coagulation pathways are differentially down-regulated in the severe hosts. Complement activation has been previously linked to children experiencing severe malaria during P. falciparum infection [27, 28], however, our result indicate that within the context of P. cynomolgi infected M. mulatta, elevated levels of complement related transcripts are associated with mild malaria. Correlation analysis of DUGs-M, cell population, and cytokine abundance demonstrates that heme metabolism and complement-related genes are associated with different cell populations and cytokine signals. Complement-related genes are associated with both innate and adaptive immune cells and heme metabolism related genes are associated with parasitemia and inflammation. Interestingly, SICAM1 which is differentially up-regulated in severe hosts at time point 2 shows a significantly higher transcript-cytokine correlation within the severe hosts than that of the mild hosts, which suggest that the abundance of SICAM1 can be further attributed to post-transcriptional activities within the mild hosts such as less cellular damage, inflammation or endothelial activation. Furthermore, IL-1B was observed to be transcriptionally differentially up-regulated in severe hosts, but the large difference of IL-1/IL1RA ratio between the two clinical groups suggest that IL-1 signaling might actually have been hindered in severe hosts.

In summary, we provided transcriptomic, cellular and cytokine evidence associated with severe disease outcome during P. cynomolgi infection in M. mulatta. We observed that severe malaria is associated with differential up-regulation of innate viral response related genes and pro-inflammatory cytokine, specifically elevated level of IL-6. Mild malaria is associated with up-regulation of the complement pathway strongly associated with monocyte and neutrophil populations. Our analysis provided a molecular and cellular basis for the development of severe malaria during P. cynomolgi infection. Larger and more frequent sampling studies are needed to validate our findings and specifically determine the importance of innate immune activation and the role of complement in controlling disease severity.

4 Material and Methods

4.1 Experimental Setup and Data Collection

A detailed description of the experimental set up and the generation of cytokine and immune profiles can be found in the publication by Joyner et al. [19]. The procedure for RNA-seq data collection and the library size normalized data can be found on Gene Omnibus under series number GSE99486. The daily clinical data of the experiment is stored in PlasmoDB using the identifier E04MalariaClinical. The cytokine and immune profiling data are deposited on ImmPort under the identifier SDY1015.

4.2 Bioinformatics Analysis

Differential expression analysis was carried out using MATLAB’s implementation of DESeq algorithm on library size normalized read counts. Enrichment analysis of identified DEG was conducted using MySigDB web service. The gene sets tested for enrichment included canonical pathways, Hallmark pathways, BioCarta pathways and KEGG pathways. A q-value cut-off of 0.05 was used. Clustering of transcriptomic, cellular population and cytokine abundances was done using k-means clustering. The number of clusters used was determined by the silhouette value of different numbers of clusters (1-6).