Systems analysis of immune responses to attenuated P . falciparum malaria sporozoite vaccination reveals excessive inflammatory signatures correlating with impaired immunity

Abstract

Immunization with radiation-attenuated sporozoites (RAS) can confer sterilizing protection against malaria, although the mechanisms behind this protection are incompletely understood. We performed a systems biology analysis of samples from the Immunization by Mosquito with Radiation Attenuated Sporozoites (IMRAS) trial, which comprised P. falciparum RAS-immunized (PfRAS), malaria-naive participants whose protection from malaria infection was subsequently assessed by controlled human malaria infection (CHMI). Blood samples collected after initial PfRAS immunization were analyzed to compare immune responses between protected and non-protected volunteers leveraging integrative analysis of whole blood RNA-seq, high parameter flow cytometry, and single cell CITEseq of PBMCs. This analysis revealed differences in early innate immune responses indicating divergent paths associated with protection. In particular, elevated levels of inflammatory responses early after the initial immunization were detrimental for the development of protective adaptive immunity. Specifically, non-classical monocytes and early type I interferon responses induced within 1 day of PfRAS vaccination correlated with impaired immunity. Non-protected individuals also showed an increase in Th2 polarized T cell responses whereas we observed a trend towards increased Th1 and T-bet+ CD8 T cell responses in protected individuals. Temporal differences in genes associated with natural killer cells suggest an important role in immune regulation by these cells. These findings give insight into the immune responses that confer protection against malaria and may guide further malaria vaccine development.

Trial registration: ClinicalTrials.gov NCT01994525.

Introduction

Malaria is a devastating disease that results in over 200 million cases and hundreds of thousands of deaths annually. Plasmodium falciparum causes the most serious disease and the most deaths, especially in sub-Saharan Africa and primarily in children [1]. Multi-pronged efforts to eliminate malaria have led to substantial reductions in malaria incidence but the development of drug and insecticide resistance as well as other factors, including the current COVID-19 pandemic, are a challenge to further progress [1, 2]. An effective anti-malarial vaccine has been a long term goal which has proven challenging. Despite exciting recent progress [3], a single approved malaria vaccine exists, the RTS,S subunit vaccine, which elicited 28–33% protection in infants over a 4-year study period [4]. An improved vaccine, especially one that prevents infection, would be a valuable tool in the effort to eliminate this disease. Understanding the immune responses that contribute to vaccine induced immune protection could aid the development of such vaccines.

Sporozoites (SPZs) are the liver-infectious life cycle stage of malaria injected via mosquito bite in natural infections. Many studies in humans and model systems have shown that vaccination with P. falciparum SPZs that have been attenuated by radiation, genetic modification, or drug treatment can result in sterilizing immunity, as determined by subsequent controlled human malaria infection (CHMI) [5–8]. This mode of vaccination aims to elicit immunity against pre-erythrocytic parasite stages, where the biomass of the parasites is low and the infection is asymptomatic. Currently, no universal correlates of protection have been identified and the nature of protective immunity is incompletely understood. Sterilizing immunity is likely to be complex and directed at multiple antigens given that the P. falciparum genome encodes more than 5,300 unique proteins. Available evidence indicates that antibodies against major surface proteins of infecting SPZs, e.g. CSP and TRAP, contribute to protection [9–11]. Animal models have indicated that liver-resident CD8+ T cell responses are important for protection. This is inherently challenging to study in humans, as the human liver is not readily accessible for sampling and a very small fraction of its cells get infected [9, 12, 13]. It is imperative to identify correlates of protection in humans that can aid the improvement of the current vaccines and development of vaccine candidates. To this end, human vaccination and challenge trials with attenuated PfSPZs provide an opportunity to elucidate immune responses that are associated with pre-erythrocytic protection.

In this study, we applied a systems immunology approach to identify correlates of protection that are identifiable up to 28 days after initial vaccination in malaria naïve human trial subjects that participated in the Immunization by Mosquito with Radiation Attenuated Sporozoites (IMRAS) trial [14]. Participants were immunized by PfRAS delivered by mosquito bite with efficacy assessed by CHMI. Five total immunizations were delivered, spaced 4 to 5 weeks apart. The trial was designed with a suboptimal vaccine dose regime to elicit approximately 50% vaccine efficacy to facilitate comparison between protected (P) and non-protected (NP) subjects. Of particular interest in the IMRAS trial is the prime vaccination. IMRAS participants are malaria unexposed, and the initial PfRAS vaccination represents the first time their immune system has been exposed to P. falciparum sporozoites. We hypothesized that the earliest immune responses to PfRAS represent a critical time period determining subsequent development of sterilizing immunity.

Our integrative analysis of whole blood transcriptomics, high parameter flow cytometry and single cell CITE-seq identified numerous vaccine-induced responses including some that correlated with protection. We observed strong negative correlations with protection associated with inflammation, type I interferon (IFN), and signatures related to monocytes and neutrophils, and type 2 polarized T helper cell responses. Differential kinetics in natural killer (NK) cell associated responses and a trend of increased T-helper 1 cells correlated positively with protection. These results suggest that the priming vaccination with radiation attenuated PfSPZs establishes immunological trajectories that can result in protection following additional vaccinations, and show that early inflammatory responses can negatively influence the fate of protective immunity.

Results

The IMRAS cohort analyzed consisted of eleven malaria-naïve adults immunized with five doses of approximately 200 bites from PfRAS NF54 infected mosquitos. All IMRAS participants were between the ages of 18–39, and all but one of the immunized participants were men [14]. The first four doses were given four weeks apart and the final dose was administered five weeks after the fourth. Protection was tested by controlled human malaria infection (CHMI) three weeks after the final vaccination (S1A Fig). Six of the eleven immunized participants were protected, i.e. zero parasitemia after CHMI. Of the five non-protected subjects one developed parasitemia on day 9 after CHMI and four did so on day 13 after CHMI (S1B Fig). There was no significant correlation between the number of PfRAS infectious mosquito bites received and protection status among true-immunized subjects (S1C Fig), and the cumulative number of bites was held at approximately 200 bites across all immunizations for each participant.

PfSPZ vaccination induced broad transcriptome responses

Whole blood transcriptome profiling was performed on all eleven immunized IMRAS participants at 6 timepoints after the initial PfRAS vaccination. We examined transcriptional changes between adjacent timepoints, which we refer to as time “intervals”, namely between days 0–1, 1–3, 3–7, 7–14, and 14–28 after immunization. We conducted linear mixed-effects regression modeling analysis (LMER) of normalized log2-transformed gene expression values to identify significantly responsive genes over all subjects and those that differed between P and NP at each interval (FDR < 0.2 & p < 0.05; see Materials and Methods). P-values calculated using this method were very strongly correlated with p-values calculated using raw RNAseq counts and the R glmmSeq [15] package (S1 Table), indicating consistent results between the two approaches. 90% confidence intervals (CIs) were calculated around model coefficients to label genes as either increased or decreased if the CI was entirely above or below 0, respectively (Fig 1A). This “interval” approach was chosen over an analysis approach relating all later timepoints to day 0/baseline to identify subsequent time points associated with large changes in differentially expressed genes, however, for comparison we have also included the parallel analysis with all timepoints relative to day 0 (S2 Fig). Many significantly responsive genes were observed after the first PfSPZ vaccination; 8170 genes had increased or decreased expression responses over at least at one interval accross all immunized subjects and approximately 10% of these differed significantly between P and NP subjects (Fig 1A).

[Figure omitted. See PDF.]

Fig 1. Vaccine induced gene responses after the first immunization.

A. Barplot and table showing numbers of vaccine-induced genes, including the subset associated with protection with increased (light red, red) or decreased (light green, green) expression over each time interval in all immunized subjects (ALL_UP+PROT_UP, ALL_DN+PROT_DN) (FDR < 0.2, p <0.05, 90% CI > 0 or < 0). Darker colors indicate genes that also differ significantly in expression between protected (P) and non-protected (NP) subjects (PROT_UP, PROT_DN). B. Heatmap of modules significantly enriched for vaccine-induced genes. Each row represents a BTM, each column represents a time interval. Heatmap color shows hypergeometric effect size (ES) of a BTM enriched in genes with increased (red/positive ES) or decreased (blue/negative ES) expression. Non-significant BTMs are shown in white. Assignment of a BTM to a high-level annotation group is indicated by a black square in the corresponding annotation row on the left.

https://doi.org/10.1371/journal.ppat.1010282.g001

The association between vaccine induced gene responses and specific cell populations and immunological processes was determined by testing whether pre-defined coherent blood transcription module sets (BTMs) [16–18] showed enrichment for significantly responsive genes (hypergeometric p < 0.1). We identified 122 BTMs significantly enriched in response genes in at least one time interval. Hierarchical clustering of enriched BTM hypergeometric effect sizes at each interval revealed discrete time-dependent response patterns among the immunized subjects (Fig 1B and S2 Table). BTMs increased during the first day after vaccination were associated broadly with immunity and inflammation, including TLR sensing; antigen processing; interferon and inflammation; monocytes and neutrophils. This was accompanied by decreases in BTMs associated with the cell cycle and T cells. Relatively few enriched BTMs were observed between day 1 and 3. Over subsequent intervals (D3-7, D7-14, D14-28) BTMs associated with monocytes, neutrophils, TLR sensing, inflammation and interferon, decreased sharply. However, we observed an increase in cell cycle-associated responses from D3-7, and an induction of T cell associated BTMs at D14-28. Overall, the priming vaccination resulted in robust and dynamic transcriptional responses in the combined group of P and NP subjects.

Protection associated genes showed distinct response dynamics in P and NP individuals

To identify de-novo co-expressed gene clusters and processes both shared by and specific to specific to P and NP participants, we performed unsupervised clustering of the responses over time of the 1394 genes that were differentially expressed between P and NP (FDR < 1/3 & p < 0.05, Fig 2A; Materials and Methods). Hierarchical clustering was performed separately for P and NP subjects to illustrate the response differences between these two groups and to identify sets of genes with coherent response profiles over time following PfSPZ immunization (Fig 2A and 2C). Immune functions and molecular mechanisms associated with these clusters were identified by the hypergeometric enrichment test using BTMs and by Ingenuity Pathway Analysis (IPA) (see Materials and Methods) (Fig 2B and 2D). We identified four major gene clusters for P subjects (P_1, P_2, P_3, P_4) and five for NP subjects (NP_1, NP_2, NP_3, NP_4, NP_5). The patterns in which these 1394 genes changed over time differed between P and NP, leading to substantial differences in gene composition of most of the major P and NP clusters (Fig 2A, 2C and 2E). A total of 39 signficantly enriched BTMs and 159 significant IPA pathways were identified between P and NP (S3 Table).

[Figure omitted. See PDF.]

Fig 2. Vaccine induced protection associated gene responses.

A,C. Heatmaps showing 1,394 protection-associated genes ordered by hierarchical clustering based on expression in A. P and C. NP subjects for each time interval. Black sidebars indicate genes associated with an IPA pathway or BTM. Colored sidebars indicate P and NP gene clusters. Expression values were z-score transformed in rows for visualization. B,D. Enriched BTMs and IPA pathways for clusters from B. P and D. NP subjects. X axis represents–log10(FDR) generated from hypergeometric tests. Color indicates assignment of a BTM module or IPA pathway to a high-level annotation group. E. Circos plot showing overlap between P and NP clusters, numbered to match A. and C., and colored by cluster number (e.g: P_1, NP_1: green, P_2, NP_2:orange). F. Expression changes of type I interferon-associated genes in P and NP subjects over each interval.

https://doi.org/10.1371/journal.ppat.1010282.g002

Most gene clusters from both P and NP were strongly enriched for cell cycle-associated BTMs (Fig 2B and 2D). An important exception was P_1 and NP_1, both were associated primarily with various immune response modules including many type I IFN-associated modules. These two gene clusters have the largest total numbers of genes and associated BTMs. In total, P_1 and NP_1 have 317 genes in common (hypergeometric p = 4.2e-17) (Fig 2E), although numbers of enriched modules and gene expression dynamics are distinct (S3A, S3B and S3C Fig). Expression of most genes in these clusters was increased over the interval D0-1, with higher response magnitude in NP_1 (S3C Fig). Consistent with strong enrichment of type I IFN associated modules in P and NP-associated cluster 1, we observed that IFN-stimulated genes (ISGs) were strongly increased overall in NP compared with P and NP by day 1 after the first vaccination (Fig 2F, Mann Whitney U test, P < NP, p < 2.2e-16) [19]. Furthermore, we observed signficantly increased expression of genes related to sensing through Pattern Recognition Receptors (PRR) in NP participants by day 1, for both MyD88 dependent and independent pathways (S4A, S4B and S4C Fig). Altogether, these patterns indicate that overall innate sensing and inflammatory responses were highly elevated in NP compared to P subjects, strongly suggesting that they negatively affected the induction of adaptive immunity against SPZ challenge.

Gene set enrichment analysis identified early inflammation as a correlate of impaired immunity

To more broadly explore immunological processes and cell types associated with protection we performed gene set enrichment analysis (GSEA) separately for P and NP subjects at each time interval (Fig 3A). GSEA facilitated identifying important coherent response modules that our above analysis based on differentially expressed genes may not have revealed. GSEA revealed 80 BTMs significantly enriched at 1 or more timepoints (FDR < 0.05). Twenty of these BTMs overlapped with the 39 BTMs enriched in differentially expressed genes.

[Figure omitted. See PDF.]

Fig 3. Gene set enrichment analysis after the first immunization in P and NP subjects.

A. Heatmap of GSEA normalized enrichments scores (NES) derived from BTM expression changes over each time interval. Red represents activated BTMs and blue represents down-regulated BTMs. BTM clusters derived from hierarchical clustering are indicated by spaces between heatmap rows. B. Relationships between BTMs within the four identified clusters. Clusters are laid out left to right and top to bottom to match top-to-bottom order of clusters in A. BTMs with shared genes are connected by lines, with thicknesses and labels corresponding to numbers of shared genes, and node sizes corresponding to the number of BTM genes genes. Predominant cell-type module annotations for each cluster are shown.

https://doi.org/10.1371/journal.ppat.1010282.g003

Patterns of responses over time differed between P and NP. Two modules showed opposite responses: M4.0 cell cycle transcription over day 0–1, and M196 platelet activation + actin binding over day 3–7 were both increased in P and decreased in NP. Hierarchical clustering of the GSEA BTM normalized enrichment scores (NESs) identified four major clusters of BTMs (Fig 3A) and many of the BTMs within a cluster shared genes (Fig 3B). GSEA cluster 1 is associated with B and T cells, cluster 2 with cell cycle and potentially cell proliferation, cluster 3 with NK cells, and cluster 4 with monocytes, neutrophils and immune activation.

Cluster 1 contained BTMs associated with T and B cells. These BTMs decreased over day 0–1 specifically in NP, with most BTMs remaining unchanged at all other time intervals for both P and NP. BTMs in cluster 2 were generally associated with cell cycle and division. These BTMs were downregulated over day 0–1 and day 7–14 and upregulated from day 3–7 in both P and NP. Notably, most of these BTMs were activated in NP but not in P from day 3–7. Cluster 3 was the most heterogeneous cluster in terms of composition and time dynamics, and included NK cell, plasma cell and trancriptional regulatory BTMs among others. These modules tended to be increased at early time intervals, before decreasing at day 7–14 and day 14–28. Notably, NK BTMs were upregulated in P over day 0–1; however, their activation occurred late in NP over day 1–3 and day 3–7. Cluster 4 primarily represented innate inflammation, interferon, monocyte, neutrophil and dendritic cell BTMs. Consistent with our above gene-based analysis, cluster 4 responses in NP individuals increased sharply at day 0–1 and day 7–14 with an intermediate decrease from day 3–7, with relatively few changes in P individuals over the first four time intervals. Both P and NP individuals showed decreased responses in cluster 4 BTMs day 14–28. Cluster 4 also included platelet BTMs, and early activation of platelets in NP may reflect their influence on the elevated expression associated with immune activation and monocytes [20]. Activated neutrophils produce reactive oxygen species (ROS), and we observed an increase in ROS signalling in NP (S5A, S5B and S5C Fig) at day 0–1. The differences between P and NP subjects in this cluster highlights the greater magnitude of inflammatory responses by innate cells in NP versus P subjects early after the first immunization. Overall, this analysis supports our hypothesis that high levels of inflammatory responses and type I IFN are detrimental for protective immunity, and suggests these responses are associated with monocytes, DCs and NK cells.

Specific immune cell types associated with protection by trancriptomics and flow cytometry

GSEA and differential gene expression analyses indicated that specific immune associated transcriptional responses were induced at different times and to different extents between P and NP subjects. This suggests that specific immune cell populations responded differently in P vs NP, thus, we further investigated whether a cell type-specific signature could be identified after initial PfRAS immunization.

We generated heatmaps of expression changes by averaging expression of genes from significant cell-type associated BTMs identified by the GSEA analysis (Figs 4A, 4B and S6A–S6G). Monocyte, DC, neutrophil and platelet-associated responses were highly increased in NP over day 0–1, indicating strong innate immune responses in NP that were absent or significantly lower in P subjects (Figs 4A, 4B and S6A, S6B, S6C and S6D). In contrast, genes associated with T and B cells were strongly decreased in expression over day 0–1 in NP subjects (Figs 4A, S6E and S6F). Both P and NP subjects had a modest increase in expression of B cell associated genes in the last time interval (S6F Fig). Interestingly, P and NP subjects differed in the timing, magnitude and nature of active NK cell-related genes (Fig 4A). Once again, these results indicate that early innate immune activation strongly correlates with insufficient immunity to challenge after the completion of the vaccine regimen and transcriptome responses associated with adaptive cell types show the opposite expression pattern.

[Figure omitted. See PDF.]

Fig 4. Cell type-associated responses after the first vaccination.

A. Heatmap showing cell type specific BTM GSEA NESs in P and NP subjects. B. Heatmap showing expression changes of malaria responding neutrophil genes in P and NP subjects on the sampling day compared with the previous sampling day. Expression values were z-score transformed by row for visualization.

https://doi.org/10.1371/journal.ppat.1010282.g004

To validate this whole blood transcriptional analysis, we applied high parameter flow cytometry to characterize PBMCs isolated from immunized individuals after the first immunization. This was done at overlapping time points with the whole blood transcriptional analyses, although we lacked matching PBMC samples for the 1 day and 28 days post immunization timepoints (S4 Table and Fig 5A–5G). We assessed whether transcriptional responses associated with specific cell types were correlated with flow cytometry derived counts of the appropriate populations.

[Figure omitted. See PDF.]

Fig 5. Temporal changes of specific cell types across time intervals in P and NP subjects.

A-D. Flow-cytometry measured PBMC sub-populations percentages over time, stratified by protected (P: pink) and not-protected (NP: blue) status. A. non-classical monocytes and ILT3+ monocytes B. DCs C. T cells D. NK cells. Points represent cell proportion values in each subject. Line represents average values across P and NP subjects, with error bars indicating +/- SD. P-values derived from from ANOVA comparison of nested mixed models with and without protection term. E-G. Flow cytometry gating schemes E. DC subsets and monocytes F. Th1 and Th2 cells G. NK cells.

https://doi.org/10.1371/journal.ppat.1010282.g005

Consistent with the RNAseq data, NP subjects showed a trend towards increased proportions of circulating monocytes shortly after the first immunization, specifically non-classical monocytes (CD14+-CD16+) (Fig 5A, p = 0.06). These inflammatory cells have been implicated in several vaccination studies as being correlated with impaired immunity [21, 22]. We also observed increased proportions of ILT3+ monocytes, associated with IFN exposure [23, 24], after the first vaccination (Fig 5A, p = 0.0004). Furthermore, by applying a previously identified transcriptional signature of macrophage polarization [25], we observed that responses associated with classically activated macrophages were specifically increased in NP between day 0–1, but no change was observed for alternatively activated macrophage responses (S6H and S6I Fig). We compared the abundance and activation status of several DC subsets between P and NP subjects and found several trends that were consistent with RNA-seq results, albeit not statistically significant (Fig 5E). We observed a subtle increase in plasmacytoid DCs, and overall increased CD11c+ DC numbers 3 days after the first immunization in NP compared to P subjects (Fig 5B). Interestingly, the proportion of “mature”CD1c+ DCs expressing CD86 was higher in P than NP subjects whereas CD86 expressing CD11c+ DCs that lack CD1c or CD141 expression, and non-activated CD1c DCs were enriched in NP subjects. These findings suggest that CD11c+ DCs, non-activated CD1c+ DCs and non-classical monocytes contribute to the transcriptional signature in NP subjects where we observed an increase in DC and monocyte-associated genes (Figs 4A, 5B, S6B, S6H, and S6I). Because different types of DCs can activate and skew different T-helper cell responses, we also assessed the proportion of Th1 and Th2 CD4+ T cells [26] and observed greater proportions of Th2 CD4 T cells in NP, and a trend towards more Th1 cells in P subjects (Fig 5C). In addition, we found a trend of higher circulating numbers of CD8 T cells and T-bet expressing CD8 T cells in P subjects. In CD8 T cells, T-bet is an important transcription factor that is involved in memory cell formation [27]. These data suggest that differences in innate responses contributed to impaired protective immunity by skewing the CD4 T cells toward a type 2 phenotype, and protective responses are hallmarked by Th1 and CD8 T cell responses. The abundance of the NK cells in general, and CD38+ NK cells specifically, matched the transcriptome responses where NK-associated gene increases peaked at day 1–3 in P subjects and day 7 in NP subjects (Fig 5D). However CD8+ NK cells subsets did not show any early peak in P participants. Notably, CD8+ NK cells lacking FcRγ were more abundant in NP subjects across all measured timepoints [28–31].

To quantify the relationships between the transcriptional and cytometric data, we performed correlation analysis of temporal changes of BTMs and cell subsets from manually gated flow cytometry data from the NK and DC panels (Fig 6A). The NK and DC panels were used (without the the T and B cell panels) as these panels covered a comprehensive variety of innate immune cells as well as invariant T cells and basic T and B cell lineage markers. These analyses showed that transcriptional activation of the innate immune response (inflammatory responses, monocyte signatures, DC signatures, viral sensing & immunity, antigen processing & presentation, etc.) positively correlated with increases in the proportion of innate immune cells (DCs, monocytes) (Fig 6B). In contrast, decreased transcriptome T cell-related responses (T cell activation, T cell differentiation) paralleled cellular decreases of CD3+ T cells (Fig 6B), and transcriptional activation of several NK-associated genes positively correlated with increases in NK cells (Fig 6B). The activation of innate immune responses (increases in BTMs related to monocytes, DCs, inflammatory responses) in NP subjects correlated with decreases in total T cell counts (Fig 6C). The down-regulation of T cell modules (T cell activation, T cell differentiation) was associated with increases in the abundance of DCs and monocytes (Fig 6C). Additionally, the downregulation of B cell-associated genes was correlated with cell composition increases in NKT-like cells that co-express CD3 and CD56 (Fig 6C). These analyses indicate that the transcriptional data and the flow cytometry data align consistently and reveal cell abundance and activation changes that correlate with each other.

[Figure omitted. See PDF.]

Fig 6. Correlation of temporal changes between cell subsets and BTMs.

A. Heatmap showing correlations between flow-cytometry derived cell subset counts from NK and DC panels (columns) and BTM expression (rows) matched by participant and time-point. Red represents positive correlation, blue represents negative correlation, and white represents non-significant correlations. High level annotations of BTM and cell types are shown as colored bars above and to the right of the heatmap. Major regions of positive and negative correlation are labelled B-H, which corresponds to regions shown in other figure panels. B-H. Changes over time (relative to day 0) of BTMs and cell subsets that are B-E. positively correlated and F-G. negatively correlated with predominant cell types and high level BTMs named in the panel title. Each line represents average abundance changes (relative to day 0) of a cell subset or BTM in P and NP subjects.

https://doi.org/10.1371/journal.ppat.1010282.g006

RNAseq and flow cytometry findings were further validated with scRNAseq

To further validate the identity of cell subsets indicated by transcriptional analysis, we performed CITEseq with a panel of 14 antibodies using PBMC samples collected on day 0, day 3 and day 14 after the priming dose of PfSPZ from an IMRAS subject who dropped out of the trial (S5 Table). We obtained single cell RNA seq gene expression and surface protein marker profile data from a total of 12,442 cells. Mapping of these cells to a previously defined multimodal cell atlas based on reference clusters identified 29 cell types in the merged samples (Fig 7A) [32]. We then derived transcriptional cell-type associated signatures using genes highly expressed in CITEseq-identified cell types in our samples, and applied these signatures in GSEA analysis of the whole blood RNA-seq data (Fig 7B). Consistent with previous BTM-based GSEA analyses, GSEA with our CITEseq signatures indicated that CD14+ and CD16+ monocyte subsets were highly activated on day 1 in NP subjects. By contrast, genes associated with NK cells were highly activated on day 1 in P subjects, whereas the activation continued on day 3 in NP subjects (Fig 7B and 7C). Finally, comprehensive correlation analysis was performed to examine the relationship between cell subsets predicted by CITEseq and cell-specific BTMs identified in our previous GSEA (Figs 3 and 4). The high correlation between cell-specific BTMs and corresponding cell types identified in CITEseq data demonstrates that the cell subsets identified by both approaches are consistent (Fig 7D) and supports our interpretation of these responses. Thus, integration of whole blood RNA-seq, high parameter flow cytometry and scRNA-seq of individuals after initial PfRAS vaccination revealed a consistent picture where inflammatory transcripts and non-classical monocytes, and Th2 signalling are specifically increased in NP individuals in the first day after PfRAS vaccination while P individuals are marked by increased early NK-cell and Th1 signalling.

[Figure omitted. See PDF.]

Fig 7. Validation of transcriptional changes of specific cell types using CITEseq data.

A. Scatter plot illustrating cell types identified in CITEseq data, visualized using uniform manifold approximation and projection (UMAP). Each point represents a cell and is colored by cell type. B. Heatmap showing GSEA NES of CITEseq-derived cell signatures in whole-blood RNAseq in response to the first immunization in P and NP subjects. C. Lineplot showing kinetics of CITEseq-identified cell signatures in P and NP subjects. Each line represents the median gene expression levels of a CITEseq gene signature per-subject, and the black dashed line represents the median gene expression levels across all in P or NP subjects. D. Heatmap showing Spearman’s rank correlation between cell-type specific BTMs and CITEseq identified cell-specific signatures. The color row bar indicates cell-type annotation of BTMs.

https://doi.org/10.1371/journal.ppat.1010282.g007

Discussion

Whole blood RNAseq, high parameter flow cytometry and CITEseq systems analyses of the IMRAS trial identified numerous responses to the first vaccination, including early inflammatory responses correlated with a lack of protection. These early inflammatory responses were associated with myeloid cells, including neutrophils and increased levels of monocyte and DC subsets, and Th2 cells. In contrast, effective protection was correlated with early responses associated with NK cells, later responses associated with T cells and lower overall responses. These results underscore the influence of early innate and adaptive responses on the subsequent immunological trajectories as shown by different response profiles that do or do not ultimately lead to protection from infection (Figs 2 and 3). We hypothesize that differential priming of the adaptive compartment by PfRAS vaccinations impacts subsequent responses and ultimately protection from infection by protected vs not protected vaccinees. Thus, immune responses to the first vaccination may be decisive for the outcomes of vaccinations with attenuated sporozoites and may provide early biomarkers of protection.

Myeloid cell activation as a negative correlate of protection may seem counterintuitive since activation is required for antigen presenting cell (APC) priming of adaptive T cell responses in lymphoid tissues. The higher levels of inflammatory responses that we found to correlate with a lack of protection in this study imply that over-induction of inflammation negatively impacted the development of protective adaptive responses. Previously, innate pathways that suppress adaptive responses have been identified in the context of immune pathology [33, 34]. Similarly, the early increase in neutrophil-associated gene expression following vaccination correlated with a lack of protection. Although neutrophils are rapid innate responders to infection and function in pathogen clearance and immune modulation [35], a subset of neutrophils that are systemically induced upon acute inflammation suppress T cell responses through ROS production [36–38]. The increased expression of transcripts associated with ROS production observed in NP subjects implies that neutrophil suppression of T cell responses may also have negatively impacted development of protective adaptive immune responses after PfRAS vaccination.

Innate priming of naïve T cells into various T-helper cells can occur through the actions of neutrophils, basophils, and DC subsets via their intrinsic properties or their interplay with other innate cells, including the induction of Th2 cells in response to type 2 innate lymphoid cell (ILC2)-derived cytokines [39–43]. The reciprocal early relative increase in Th2 cells and decrease in Th1 cells in NP subjects implies that inhibition of polarization of Th1 cells by Th2 cytokines suppressed development of responses that protect against CHMI [39]. We also observed a trend towards higher numbers of Th1 and CD8 T cells in P vs NP subjects at early timepoints, suggesting that Th1 responses contribute to PfRAS vaccination-induced immunity. This is consistent with studies in animal models which concluded that sterilizing protection involves Th1 cells which secrete IFNγ and, in co-ordination with with macrophages and liver resident cytotoxic CD8 T cells, eliminate intracellular pathogens [13, 14, 44–47]. In addition, the correlation of CSP-specific, IFNγ-producing CD4 T cells with protection from infection following RTS,S/AS02 vaccination [48] also suggests that Th1 cells can contribute protective immunity in humans.

Increases in monocyte-associated transcriptome responses and of non-classical and ILT3-expressing monocyte counts in NP subjects suggest that these responses impair development of protection following Pf vaccination. Impairment of vaccine elicited immunity has been linked to non-classical or inflammatory monocytes in mice [21, 22, 49] perhaps involving ILT3 monocyte recruitment to lymph nodes and interference with T cell priming as implied by the enhanced T cell priming following monocyte depletion. GM-CSF production by monocytes and sequestration of cysteine have been suggested as mechanisms of T cell suppression during priming [22, 50]. The correlation between type I IFN expression and the lack of protection that observed in this work (Fig 2F) and in a vaccine trial of malaria SPZs given under chemoprophylaxis cover [51] indicates that these responses can hamper the development of adaptive responses. That ILT3 surface expression by monocytes can result from type I IFN stimulation [23, 24] and both type I IFN and ISGs increased expression in NP subjects early after vaccination suggests that the monocytes in NP subjects have encountered type I IFN stimulation. Type I IFNs are potent immune mediators that can directly activate DCs, NK cells and T and B cells and regulate immune responses to many pathogens. These IFNs signal through the interferon alpha receptor (IFNAR) that is present on almost all cells in the body and induce ISGs that function in the control of infection. However, type I IFN has also been linked to immune pathology in chronic viral diseases and in some bacterial infections [52]. The ultimate protective vs non protective effects of the type I IFN response may well depend on its timing, localization and magnitude of the specific IFN responses as has been suggested [53].

Protection elicited by vaccination with radiation attenuated SPZs requires an abortive liver infection [54]. Infection studies in animals have shown that malaria infected hepatocytes produce type I IFN, and IFNγ-secreting NK and NKT cells are recruited as liver stage parasites are eliminated [55]. However, type I IFN responses are also detrimental to long term immunity against infection [34, 55–57]. Caspase-mediated cell death of infected hepatocytes is required for the uptake and presentation of malaria antigens by innate phagocytic cells [58–60] but type I IFN can inhibit caspase activity and inflammasome activation thus potentially inhibiting the presentation of malaria antigens [61, 62]. IFNAR signaling during malaria liver infection may also impair the induction of Th1 and CD8 T cell responses and enhance exhaustion of liver resident CD8 T cells [34, 53, 63]. In addition, type I IFN can hamper protective immunity via inhibition of IFNγ responsiveness by monocytes and macrophages that in turn can limit the induction of Th1 responses [64, 65]. Overall IFNγ is an important mediator in anti-malarial immunity with a variety of downstream effects besides macrophage activation [66]. We did not find that the best known and potent producers of large quantities of type I IFN, namely plasmacytoid DCs (pDCs), had significantly higher levels of in NP subjects; however, other innate cells, e.g. neutrophils, monocytes and DCs can secrete type I IFN [67, 68].

NK cell-associated transcriptome responses (Figs 3A and 4A) and relative changes in NK cell subset numbers (Fig 6B) occur earlier in P than in NP subjects, which may indicate that they contribute to the development of immunity, perhaps via the balance between Th1 and Th2 responses. NK cells can be activated by neutrophils as well as inflammatory monocytes through type I IFN [69, 70] and they can have diverse functions. They can act as immune regulators that enhance or suppress adaptive responses and they can be innate effectors that rapidly respond to and eliminate infected or tumor cells [71, 72]. IFNy secretion by NK cells can support DC-mediated Th1 induction [71, 73], analogous to that we observed in P subjects, but NK cells can also limit adaptive responses by suppressing DCs, CD4 T cells and B cells and thus variably impact outcomes [73–76]. Furthermore, NK cells can reduce inflammation, e.g. as with COVID-19 related immune pathology [77–80]. Thus, the early NK responses in P subjects may support the priming of a Th1 polarized response and inhibit inflammation that in NP subjects primes a Th2 response. Interestingly, we identified a novel NK cell subset that expressed CD8 and lacked FcRγ expression and which is more abundant in NP subjects in the day 3–7 interval (Fig 5D). FcRγ-lacking NKs have previously been associated with an adaptive phenotype that is protective against seasonal malaria infection [31]. Further investigations into this phenotype could elucidate its functionality.

Overall, our analyses indicate that early innate responses to radiation-attenuated SPZ vaccination substantially impact the development of adaptive responses and ultimately protection from malaria infection. Understanding mechanisms of protection is complicated by the likelihood that protective effector processes are multi-functional, due to the large breadth of potential antigens, and protection may occur at multiple points between the introduction of SPZs and the establishment of a blood stage infection. Both antibody and cell-mediated mechanisms may contribute to protection: monoclonal antibodies derived from attenuated SPZ vaccination can protect humanized mice although antibody levels variably correlate with protection [14, 81] and liver resident CD8+ T-cells and IFNγ correlate with protection in non-human primates [82]. The complex balance of responses associated with protection are illustrated here by the differential early inflammatory responses between P and NP subjects and the potential effects on Th1 and Th2 responses. The events following the priming vaccination that correlated with protection are early NK associated responses in the context of limited inflammatory responses followed by CD4+ T cell responses and subsequently CD8+ T cell responses.

That this study was performed on blood samples, despite decisive immune events occuring in the liver, which is essentially experimentally inaccessible, limited us to indirect analysis of phenotypic differences between the P and NP subjects rather than functional assays. The sample size is relatively small, which reduced our sensitivity to discover more subtle changes associated with protection, especially given natural variation in the participants. However, the detection of robust correlates of protection despite the low numbers increases our confidence that we have identified meaningful responses. These findings can inform further studies to extend the understanding of protective immunity and its development. In addition, multiple factors may have influenced the differences that between P and NP subjects that we described here, and which influenced protection. These include intrinsic differences between subjects, such as HLA type and other genetic differences as well as baseline immune status at the time of immunization. Also vaccination via infected mosquito bites may have contributed to variability in the effective vaccine dose received by each participant, i.e. the number of liver cells infected by live attenuated SPZs. In addition, the point at which parasite development in the liver cells was arrested may have been variable since radiation damage is random which may have impacted the amount and type of parasite antigen available for presentation.

In conclusion, we show that a strong acute inflammatory response to a priming vaccination correlates with the ultimate lack of protection in this trial of malaria naive volunteers. We hypothesize that this results in skewing adaptive responses toward Th2-centered responses rather than protective Th1 responses and that this similarly impacts responses to subsequent immunizations. Thus, immune responses to the first vaccination can be decisive for the outcome of the trial.

Materials and methods

Ethics statement

The study was conducted at the Naval Medical Research Center (NMRC) Clinical Trials Center from 2014 to 2016; the CHMIs were conducted at the Walter Reed Army Institute of Research (WRAIR) secure insectary. The study protocol was reviewed and approved by the NMRC Institutional Review Board in compliance with all federal regulations governing the protection of human subjects. WRAIR holds a Federal-wide Assurance from the Office of Human Research Protections (OHRP) under the Department of Health and Human Services as does NMRC. NMRC also holds a Department of Defense/Department of the Navy Federal- wide Assurance for human subject protections. All key personnel were certified as having completed mandatory human subjects’ protection curricula and training under the direction of the WRAIR Institutional Review Board or the NMRC Office of Research Administration (ORA) and Human Subjects Protections Branch (HSPB). All potential study subjects provided written, informed consent before screening and enrollment and had to pass an assessment of understanding. This study was conducted according to the Declaration of Helsinki as well as principles of Good Clinical Practices under the United States Food and Drug Administration Investigational New Drug (IND) application BB-15767. This trial was performed under an IND allowance by the Food and Drug Administration (FDA).

Sample collection and RNA sequencing

Whole blood was collected from IMRAS trial participants directly into PAXgene blood RNA tubes (PreAnalytiX, Hombrechtikon, Switzerland) and stored at −20°C. RNA extraction and globin transcript depletion (GlobinClear, ThermoFisher Scientific, MA, USA) were performed prior to cDNA library preparation using the Illumina TruSeq Stranded mRNA sample preparation kit (Illumina, CA, USA). Globin transcript depletion, cDNA library preparation and RNA sequencing were performed by Beijing Genomics Institute (Shenzhen, China). A total of sixty-six RNA-seq samples were sequenced, with a target depth of 30 million reads per sample. Eleven of the samples were sequenced on Illumina (San Diego, CA) Hiseq2000 sequencers using 75 base-pair (bp) paired-end reads. The remaining one hundred and eighty-six samples were sequenced on BGI500 sequencers using 100 bp paired-end reads.

Quality control and processing of RNA-Seq data

RNAseq data were processed as previously described [83]. Read pairs were adjusted to set base calls with phred scores < 5 to ‘N’. Read pairs for which either end had fewer than 30 unambiguous base calls were removed, a method that indirectly removes pairs containing mostly adaptor sequences. Read pairs were aligned to the human genome (hg19) using STAR (v2.3.1d) [84]. Gene count tables were generated using htseq (v. 0.6.0) with the intersection-strict setting on and Ensembl gene annotations (GRCh37.74) used to link genomic locations to gene identifiers [85]. Log2-transformed TMM-normalized counts-per-million (CPM) expression matrices were computed using the cpm function of the edgeR package [86]. Batch correction for sequencer model (Hiseq2000, day 0 samples vs BGI5000, day 1–28 samples) was performed on log2-transformed counts using linear mixed-effects models with normally distributed errors and an unstructured covariance matrix. A total of 4 samples originally sequenced on Hiseq2000 were resequenced on the BGI5000 platform to facilitate batch correction unconfounded by time post vaccination. Mixed-effects models were fit using the R (https://www.r-project.org/) lme4 package [87]. The following formula was used:in which EXPRESSION represents the log2-transformed CPM value, SEQUENCER the sequencing platform, and including random intercepts for each PARTICIPANT. To create a final batch-corrected expression matrix, raw CPMs were adjusted by subtracting the fitted SEQUENCER coefficient.

Mixed-effects modeling to identify transcriptional signatures that were regulated by the primary vaccine and responded differentially in protected and non-protected immunized participants

Linear mixed-effects regression models (LMER) were used to model individual gene expression (EXPRESSION) as a function of sample collection time (TIME) and protection after CHMI (PROTECTION), with TIME and EXPRESSION as fixed effects, and PARTICIPANT as a random effect.

Mixed-models were fit as follows:

By contrasting the full model with reduced models lacking the TIME and PROTECTION terms, the significance of relationships between the TIME and PROTECTION variables and EXPRESSION were evaluated. ANOVA was used to compare the full model with reduced model 1, where P-values represent the significance of the improvement of fit associated with the TIME term in the analysis. FDR-adjusted P-values were computed using the Benjamini-Hochberg method. PROTECTION-associated genes were similarly identified within the TIME significant genes using ANOVA to compare the full model with reduced model 2. In a similar manner, the glmmSeq R package was used to calculate p-values using the same nested mixed model approach, however models were fit on raw RNAseq counts, as required by the glmmSeq statistical approach, rather than normalized transformed CPM values.

To identify genes with significant changes in EXPRESSION at specific time points relative to the pre-vaccination state, five full models were fit for each gene with two time points (each time point following vaccination 1 and its previous time point, i.e. time intervals) included in the TIME term. In addition, for each gene, models were fit that included all time points (Days 0,1,3,7,14,28) to identify transcriptional signatures that had temporal effects at any time point.

To determine the direction (UP/DOWN) of transcriptional responses relative to either pre-vaccination time point or the previous time point in all immunized subjects, 90% confidence intervals were estimated for the TIME coefficient of the reduced model 2 as above. Cases where the lower CI > 0 were considered UP genes, upper CI < 0 were considered DOWN genes.

Mixed-effects modeling to identify cell types that had significantly different cell proportion changes in P and NP immunized participants after the primary vaccination

Linear mixed-effects regression models (LMER) were used to model individual cell type proportion (PERCENTAGE) as a function of sample collection time (TIME) and protection after CHMI (PROTECTION), with TIME and PERCENTAGE as fixed effects, and PARTICIPANT as a random effect.

Mixed-models were fit as follows:

By contrasting the full model with reduced models lacking the PROTECTION terms, the significance of relationships between the PROTECTION variables and PERCENTAGE were evaluated. ANOVA was used to compare the full model with reduced model 1, where P-values represent the significance of the improvement of fit associated with the PROTECTION term in the analysis.

Gene set enrichment analysis (GSEA)

GSEA was performed for each vaccination time interval using the R fgsea package [88–90] with 500 permutations and whole blood transcriptional modules [16, 17, 51, 89]. Genes were ranked by average fold change across each time interval (day 1 to day 0, day 3 to day 1, day 7 to day 3, day 14 to day 7, day 28 to day 14) separately for samples from P and NP subjects. Normalized enrichment scores (NES) of non-significant modules (FDR-adjusted P-value > 0.05) were set to 0.

Statistical tests

The hypergeometric test was used to identify BTM modules enriched for subsets of genes. The resultant effect size (ES) was calculated as: (b/n)/(B/N), in which n: Number of genes of interest; N: Number of total mapped genes; b: Number of genes of interest from the given module; B: Number of genes from the module in the total mapped genes.

Ingenuity pathway analysis (IPA) was performed using the IPA software from Qiagen. P-values were calculated using Fisher’s Exact Test and FDR-adjusted P-values < 0.1 were considered significant.

Unsupervised clustering

Hierarchical clustering of summary measures representing gene expression/responses (average gene fold changes, GSEA NES, or ES scores) was computed by agglomerative complete linkage with 1 - (Pearson’s correlation coefficient) as the distance metric. The optimal number of clusters was determined by the “elbow” method [91].

Flow cytometry data

PBMCs were collected from IMRAS participants on day 0, and 3, 7 and 14 after the first immunization and frozen for later use. After thawing in RPMI supplemented with 10% FBS and benzonase nuclease (Millipore EMD 0.05 U/ml), the samples were incubated with LIVE/DEAD Fixable Blue Dead Cell Stain Kit and the Human BD Fc Block for 30 min at room temperature before being simultaneously stained with four phenotyping panels that have been previously described in OMIP-044 and OMIP-064 and further described in S4 Table [92, 93]. The cells were then acquired using a BD FACSymphony flow cytometer. The data were analyzed, and cellular populations gated and quantified using FlowJo Software (version 9.6.6). The percentage contribution of each manually gated cell subset was calculated using the counts of each defined cell subset divided by the total single live cells from that sample. Pearson correlation coefficients were calculated between cell type proportion changes and BTM mean expression level changes per-time interval for P an NP subjects separately.

CITE-seq single-cell RNA seq processing

Live frozen PBMCs were obtained from a single vaccinated individual in cohort 1 of the IMRAS trial at day 0, and three- and 14-days post first vaccination. Cells were thawed and washed with RPMI supplemented with 10% FBS and benzonase nuclease (Millipore EMD 0.05 U/ml). PBMCs were resuspended in 100 μl of PBS supplemented with 2% w/v Fetal Bovine Serum (FBS) and incubated with Fixable Viability Stain 510 and Human BD Fc Block for 30 minutes at room temperature. Cells were washed with 2% FBS PBS before incubating with a panel of previously titrated 14 barcoded oligo-conjugated antibodies (BioLegend TotalSeq-C), including FITC-anti-CD45. Stained PBMC samples were then sorted by fluorescence activated cell sorting (FACS) on a BD FACSMelody to enrich for live, hematopoietic cells. A standard viable CD45+ cell gating scheme was employed; FSC-A v SSCA (to exclude sub-cellular debris), two FSC-A doublet exclusion gates (FSC-W followed by FSC-H), dead cell exclusion gate (BV510 LIVE/DEAD negative) followed by CD45+ inclusion gate.

Sorted cells were resuspended in PBS supplemented with 1% BSA. Cells were loaded onto the 10X Chromium system, where we aimed for recovery of ~5000 cells per sample, and subjected to partitioning with barcoded 5’ V1.1 chemistry gel-beads (10X Genomics) to generate the Gel-Bead in Emulsions (GEMs). The RT reaction was conducted in the GEMs, barcoded cDNA extracted by post-GEM RT-cleanup, and cDNA and antibody barcodes amplified with 14 cycles. Amplified cDNA was subjected to SPRI bead cleanup at 0.6X. Amplified antibody barcodes were recovered from the supernatant and were processed to generate TotalSeq-C libraries as instructed by the manufacturers (10X Genomics and BioLegend, TotalSeq-C with 10x Feature Barcoding Protocol). The remaining amplified cDNA was subjected to enzymatic fragmentation, end-repair, A-tailing, adapter ligation and 10X specific sample indexing as per manufacturer’s protocol. Libraries were quantified using Bioanalyzer (Agilent) analysis. 10x Genomics scRNA-Seq and TotalSeq-C libraries were pooled and sequenced on an Illumina NovaSeq Sp100 flow cell using the recommended sequencing read lengths of 26 bp (Read 1), 8 bp (i7 Index Read), and 91 bp (Read 2), and depths of 50,000 and 5000 read pairs per cell for the 5’ Gene Expression and TotalSeq-C libraries respectively. Cell Ranger v3.1.0 (10x Genomics) was used to demultiplex raw sequencing data and quantitate transcript levels against the 10x Genomics GRCh38 reference.

Single-cell RNA seq processing and analysis

Raw count data were filtered to remove cells where 1) a mitochondrial RNA fraction greater than 7.5% of total RNA counts per cell, and 2) less than 200 or greater than 2500 genes were detected. The resultant count matrix was used to create a Seurat (v4.0.1) [32] object. Filtered read counts were normalized, scaled, and corrected for mitochondrial and rRNA read percentages with the SCTransform function. The ADT matrix was normalized per feature using center log normalization. Cell types in each sample were annotated by mapping to the annotated reference PBMC dataset provided in the Seurat v4 Azimuth workflow. Briefly, anchors between the query and reference datasets were identified using a precomputed supervised PCA on the reference dataset. Next, cell type labels from the reference dataset, as well as imputations of all measured protein markers, were transferred to each cell of the query datasets through the previously identified anchors. The query datasets were then merged and projected onto the UMAP structure of the reference. The genes expressed in each specific cell cluster were identified using the FindAllMarkers function from the Seurat4 package and filtered to include those with average log2 fold changes greater than 1 and FDR-adjusted P-values less than 0.05.

Supporting information

S1 Fig. Overview of the IMRAS trial.

A. Schematic indicating timing of vaccination and sampling. Black vertical arrows indicate immunizations, with the first (V1) and second (V2) immunization indicated. Blood samples taken between V1 and V2 indicated with vertical black lines. Time of CHMI shown as a red vertical arrow. B. Kaplan-Meier curve showing days to thick blood smear positivity for IMRAS subjects who were RAS-immunized (protected/non-protected) or mock immunized. C. Number of infectious mosquito bites received by each subject at each immunization. Circles and triangles indicate non-protected and protected subjects, respectively. Red lines indicate the median number of infectious mosquito bites.

https://doi.org/10.1371/journal.ppat.1010282.s001

(TIF)

S2 Fig. Vaccine induced genes relative to baseline.

Barplot showing numbers of vaccine-induced genes with increased (pink/red) or decreased (light green/dark green) expression relative to day 0 in all immunized subjects (TIME_UP+PROT_UP, TIME_DOWN+PROT_DOWN) (FDR < 0.2, p <0.05, 90% CI > 0 or < 0). Darker colors indicate genes that additionally differ significantly in expression between protected (P) and non-protected (NP) subjects (PROT_UP, PROT_DOWN).

https://doi.org/10.1371/journal.ppat.1010282.s002

(TIF)

S3 Fig. Comparison of protection associated clusters P_1 and NP_1.

A,B. Gene overlap of IPA pathways and BTMs enriched in A. cluster P_1 and B. NP_1, Node sizes indicate numbers of genes in a BTM or IPA pathway, and line thickness indicates the numbers of shared genes between two nodes. C. Heatmap showing expression of 317 genes common in cluster 1 of P subjects and cluster 1 of NP subjects. Expression values were z-score transformed in rows for visualization.

https://doi.org/10.1371/journal.ppat.1010282.s003

(TIF)

S4 Fig. Expression changes in pattern-recognition receptor pathways.

A,B. Heat maps of expression changes for A. MyD88-dependent toll-like receptor associated genes and B. MyD88-independent toll-like receptor associated genes in P and NP subjects. Genes shown were selected using Gene-ontology (GO) annotations GO:0002755 (MyD88 dependent TLR signalling pathway) and GO:0002756 (MyD88 independent TLR signalling pathway). Expression values were Z-score transformed in rows for visualization. C. Average gene expression of selected IPA pathways over time in P (green) and NP (red) subjects. Dots represent average gene expression values per-individuals and solid line represents the average gene expression of the IPA pathway across all participants.

https://doi.org/10.1371/journal.ppat.1010282.s004

(TIF)

S5 Fig. Expression changes in reactive-oxygen species pathway genes.

A. Heatmap showing expression of Hallmark Reactive Oxygen Species (ROS) Pathway genes. Expression values were z-score transformed in rows for visualization. B. STRING-DB derived protein-protein interaction networks of ROS genes, colored by expression changes on day 1 compared to day 0 seperately for P and NP. C. Average expression profiles of genes of Hallmark Reactive Oxygen Species Pathway in P and NP subjects. Dots represent average expression in individuals and solid line represents the average expression of the pathway over all P or NP subjects.

https://doi.org/10.1371/journal.ppat.1010282.s005

(TIF)

S6 Fig. Expression changes in pathways associated with specific immune cell types: A-I.

Heatmaps showing expression of genes of cell-type specific pathways: A. Neutrophils, B. Monocytes, C. DCs, D. Platelets, E. T-cells, F. B cells, G. NK cells, H. classically activated macrophages and I. alternatively activated macrophages. Expression values were z-score transformed in rows for visualization.

https://doi.org/10.1371/journal.ppat.1010282.s006

(TIF)

S1 Table. Correlation between lme4-determined p-values (based on fitting mixed models to TMM-normalized log2 transformed counts) and glmmSeq determined p-values.

https://doi.org/10.1371/journal.ppat.1010282.s007

(XLSX)

S2 Table. GSEA enrichment scores per-BTM per time interval, with high level BTM annotations.

https://doi.org/10.1371/journal.ppat.1010282.s008

(XLSX)

S3 Table. Genes and pathways associated with hierarchical clustering of protection-associated genes.

https://doi.org/10.1371/journal.ppat.1010282.s009

(XLSX)

S4 Table. List of antibodies used for high parameter flow cytometry.

Three separate staining panels for the detection of DC subsets, monocytes, T cells and NK cells are shown.

https://doi.org/10.1371/journal.ppat.1010282.s010

(XLSX)

S5 Table. List of CITEseq antibodies.

https://doi.org/10.1371/journal.ppat.1010282.s011

(XLSX)

S1 CONSORT checklist. Consort 2010 Checklist of information [9].

https://doi.org/10.1371/journal.ppat.1010282.s012

(DOC)

Acknowledgments

We would like to thank Michael Gale and Nana Minkah for their helpful comments on the manuscript. We would also like to acknowledge the contribution of all the IMRAS study volunteers.

Citation: Du Y, Hertoghs N, Duffy FJ, Carnes J, McDermott SM, Neal ML, et al. (2022) Systems analysis of immune responses to attenuated P. falciparum malaria sporozoite vaccination reveals excessive inflammatory signatures correlating with impaired immunity. PLoS Pathog 18(2): e1010282. https://doi.org/10.1371/journal.ppat.1010282

References

1. Vekemans J, Schellenberg D, Benns S, O’Brien K, Alonso P. Meeting report: WHO consultation on malaria vaccine development, Geneva, 15–16 July 2019. Vaccine. 2021; pmid:33931251

2. Weiss DJ, Bertozzi-Villa A, Rumisha SF, Amratia P, Arambepola R, Battle KE, et al. Indirect effects of the COVID-19 pandemic on malaria intervention coverage, morbidity, and mortality in Africa: a geospatial modelling analysis. Lancet Infect Dis. Lancet Publishing Group; 2021;21: 59–69. pmid:32971006

3. Datoo MS, Magloire Natama H, Somé A, Traoré O, Rouamba T, Bellamy D, et al. High Efficacy of a Low Dose Candidate Malaria Vaccine, R21 in 1 Adjuvant Matrix-M, with Seasonal Administration to Children in Burkina Faso. SSRN Electron J. 2021; pmid:33964223

4. RTSS Clinical Trials Partnership. Efficacy and safety of RTS,S/AS01 malaria vaccine with or without a booster dose in infants and children in Africa: Final results of a phase 3, individually randomised, controlled trial. Lancet. Lancet Publishing Group; 2015;386: 31–45. pmid:25913272

5. Coelho CH, Doritchamou JYA, Zaidi I, Duffy PE. Advances in malaria vaccine development: Report from the 2017 malaria vaccine symposium. npj Vaccines. Nature Publishing Group; 2017. pmid:29522056

6. Portugal S, Tipton CM, Sohn H, Kone Y, Wang J, Li S, et al. Malaria-associated atypical memory B cells exhibit markedly reduced B cell receptor signaling and effector function. Elife. eLife Sciences Publications Ltd; 2015;4. pmid:25955968

7. Obeng-Adjei N, Portugal S, Tran TM, Yazew TB, Skinner J, Li S, et al. Circulating Th1-Cell-type Tfh Cells that Exhibit Impaired B Cell Help Are Preferentially Activated during Acute Malaria in Children. Cell Rep. Elsevier B.V.; 2015;13: 425–439. pmid:26440897

8. Portugal S, Moebius J, Skinner J, Doumbo S, Doumtabe D, Kone Y, et al. Exposure-Dependent Control of Malaria-Induced Inflammation in Children. PLoS Pathog. Public Library of Science; 2014;10. pmid:24743880

9. Ishizuka AS, Lyke KE, DeZure A, Berry AA, Richie TL, Mendoza FH, et al. Protection against malaria at 1 year and immune correlates following PfSPZ vaccination. Nat Med. Nature Publishing Group; 2016;22: 614–623. pmid:27158907

10. Seder RA, Chang LJ, Enama ME, Zephir KL, Sarwar UN, Gordon IJ, et al. Protection against malaria by intravenous immunization with a nonreplicating sporozoite vaccine. Science (80-). 2013;341: 1359–1365. pmid:23929949

11. Trieu A, Kayala MA, Burk C, Molina DM, Freilich DA, Richie TL, et al. Sterile protective immunity to malaria is associated with a panel of novel P. falciparum antigens. Mol Cell Proteomics. Mol Cell Proteomics; 2011;10. pmid:21628511

12. Epstein JE, Tewari K, Lyke KE, Sim BKL, Billingsley PF, Laurens MB, et al. Live attenuated malaria vaccine designed to protect through hepatic CD8+ T cell immunity. Science. 2011;334: 475–80. pmid:21903775

13. Fernandez-Ruiz D, Ng WY, Holz LE, Ma JZ, Zaid A, Wong YC, et al. Liver-Resident Memory CD8+ T Cells Form a Front-Line Defense against Malaria Liver-Stage Infection. Immunity. Cell Press; 2016;45: 889–902. pmid:27692609

14. Hickey B, Teneza-Mora N, Lumsden J, Reyes S, Sedegah M, Garver L, et al. IMRAS—A clinical trial of mosquito-bite immunization with live, radiation-attenuated P. falciparum sporozoites: Impact of immunization parameters on protective efficacy and generation of a repository of immunologic reagents. Borrmann S, editor. PLoS One. 2020;15: e0233840. pmid:32555601

15. Lewis M, Goldmann K, Sciacca E. glmmSeq: General Linear Mixed Models for Gene-Level Differential Expression. 2021;

16. Chaussabel D, Quinn C, Shen J, Patel P, Glaser C, Baldwin N, et al. A modular analysis framework for blood genomics studies: application to systemic lupus erythematosus. Immunity. 2008;29: 150–64. pmid:18631455

17. Li S, Rouphael N, Duraisingham S, Romero-Steiner S, Presnell S, Davis C, et al. Molecular signatures of antibody responses derived from a systems biology study of five human vaccines. Nat Immunol. 2014;15: 195–204. pmid:24336226

18. Liberzon A, Birger C, Thorvaldsdóttir H, Ghandi M, Mesirov JP, Tamayo P. The Molecular Signatures Database (MSigDB) hallmark gene set collection. Cell Syst. 2015;1: 417–425. pmid:26771021

19. Kazmin D, Nakaya HI, Lee EK, Johnson MJ, Van Der Most R, Van Den Berg RA, et al. Systems analysis of protective immune responses to RTS,S malaria vaccination in humans. Proc Natl Acad Sci U S A. National Academy of Sciences; 2017;114: 2425–2430. pmid:28193898

20. Rolfes V, Ribeiro LS, Hawwari I, Böttcher L, Rosero N, Maasewerd S, et al. Platelets Fuel the Inflammasome Activation of Innate Immune Cells. SSRN Electron J. 2019;

21. George VK, Pallikkuth S, Pahwa R, De Armas LR, Rinaldi S, Pan L, et al. Circulating inflammatory monocytes contribute to impaired influenza vaccine responses in HIV-infected participants. AIDS. Lippincott Williams and Wilkins; 2018;32: 1219–1228. pmid:29683844

22. Mitchell LA, Henderson AJ, Dow SW. Suppression of Vaccine Immunity by Inflammatory Monocytes. J Immunol. The American Association of Immunologists; 2012;189: 5612–5621. pmid:23136203

23. Waschbisch A, Sanderson N, Krumbholz M, Vlad G, Theil D, Schwab S, et al. Interferon beta and vitamin D synergize to induce immunoregulatory receptors on peripheral blood monocytes of multiple sclerosis patients. PLoS One. Public Library of Science; 2014;9. pmid:25551576

24. Jensen MA, Yanowitch RN, Reder AT, White DM, Arnason BG. Immunoglobulin-like transcript 3, an inhibitor of T cell activation, is reduced on blood monocytes during multiple sclerosis relapses and is induced by interferon β2-1b. Mult Scler. Mult Scler; 2010;16: 30–38. pmid:20007427

25. Becker M, De Bastiani MA, Parisi MM, Guma FTCR, Markoski MM, Castro MAA, et al. Integrated Transcriptomics Establish Macrophage Polarization Signatures and have Potential Applications for Clinical Health and Disease. Sci Rep. 2015;5: 13351. pmid:26302899

26. Collin M, Bigley V. Human dendritic cell subsets: an update [Internet]. Immunology. Blackwell Publishing Ltd; 2018. pp. 3–20. pmid:29313948

27. Sullivan BM, Juedes A, Szabo SJ, Von Herrath M, Glimcher LH. Antigen-driven effector CD8 T cell function regulated by T-bet. Proc Natl Acad Sci U S A. 2003;100: 15818–15823. pmid:14673093

28. McKinney EF, Cuthbertson I, Harris KM, Smilek DE, Connor C, Manferrari G, et al. A CD8+ NK cell transcriptomic signature associated with clinical outcome in relapsing remitting multiple sclerosis. Nat Commun. Nature Research; 2021;12. pmid:33504809

29. Zambello R, Barilà G, Manni S, Piazza F, Semenzato G. NK cells and CD38: Implication for (Immuno)Therapy in Plasma Cell Dyscrasias [Internet]. Cells. NLM (Medline); 2020. pmid:32245149

30. Hwang I, Zhang T, Scott JM, Kim AR, Lee T, Kakarla T, et al. Identification of human NK cells that are deficient for signaling adaptor FcRγ and specialized for antibody-dependent immune functions. Int Immunol. Oxford University Press; 2012;24: 793–802. pmid:22962434

31. Hart GT, Tran TM, Theorell J, Schlums H, Arora G, Rajagopalan S, et al. Adaptive NK cells in people exposed to Plasmodium falciparum correlate with protection from malaria. J Exp Med. Rockefeller University Press; 2019;216: 1280–1290. pmid:30979790

32. Hao Y, Hao S, Andersen-Nissen E, Mauck WM, Zheng S, Butler A, et al. Integrated analysis of multimodal single-cell data. Cell. Cell Press; 2021; pmid:34062119

33. Wherry EJ, Ha S-J, Kaech SM, Haining WN, Sarkar S, Kalia V, et al. Molecular signature of CD8+ T cell exhaustion during chronic viral infection. Immunity. 2007;27: 670–84. pmid:17950003

34. Minkah NK, Wilder BK, Sheikh AA, Martinson T, Wegmair L, Vaughan AM, et al. Innate immunity limits protective adaptive immune responses against pre-erythrocytic malaria parasites. Nat Commun. Springer US; 2019;10: 1–14. pmid:30602773

35. Mantovani A, Cassatella MA, Costantini C, Jaillon S. Neutrophils in the activation and regulation of innate and adaptive immunity. Nat Publ Gr. 2011; pmid:21785456

36. Zemans RL. Neutrophil-mediated T-cell suppression in influenza: Novel finding raising additional questions [Internet]. American Journal of Respiratory Cell and Molecular Biology. American Thoracic Society; 2018. pp. 423–425. pmid:29717897

37. Pillay J, Kamp VM, Van Hoffen E, Visser T, Tak T, Lammers JW, et al. A subset of neutrophils in human systemic inflammation inhibits T cell responses through Mac-1. J Clin Invest. American Society for Clinical Investigation; 2012;122: 327–336. pmid:22156198

38. Kusmartsev S, Nefedova Y, Yoder D, Gabrilovich DI. Antigen-Specific Inhibition of CD8 + T Cell Response by Immature Myeloid Cells in Cancer Is Mediated by Reactive Oxygen Species. J Immunol. The American Association of Immunologists; 2004;172: 989–999. pmid:14707072

39. Kim B, Kim TH. Fundamental role of dendritic cells in inducing Th2 responses [Internet]. Korean Journal of Internal Medicine. Korean Association of Internal Medicine; 2018. pp. 483–489. pmid:29502361

40. Phythian-Adams AT, Cook PC, Lundie RJ, Jones LH, Smith KA, Barr TA, et al. CD11c depletion severely disrupts Th2 induction and development in vivo. J Exp Med. The Rockefeller University Press; 2010;207: 2089–2096. pmid:20819926

41. Tjota MY, Sperling AI. Distinct dendritic cell subsets actively induce Th2 polarization [Internet]. Current Opinion in Immunology. Elsevier Ltd; 2014. pp. 44–50. pmid:25290173

42. Otsuka A, Nakajima S, Kubo M, Egawa G, Honda T, Kitoh A, et al. Basophils are required for the induction of Th2 immunity to haptens and peptide antigens. Nat Commun. Nature Publishing Group; 2013;4: 1–9. pmid:23612279

43. Tang H, Cao W, Kasturi SP, Ravindran R, Nakaya HI, Kundu K, et al. The T helper type 2 response to cysteine proteases requires dendritic cell-basophil cooperation via ROS-mediated signaling. Nat Immunol. Nat Immunol; 2010;11: 608–617. pmid:20495560

44. McNamara HA, Cai Y, Wagle M V., Sontani Y, Roots CM, Miosge LA, et al. Up-regulation of LFA-1 allows liver-resident memory T cells to patrol and remain in the hepatic sinusoids. Sci Immunol. American Association for the Advancement of Science; 2017;2. pmid:28707003

45. Tse S-W, Radtke AJ, Espinosa DA, Cockburn IA, Zavala F. The chemokine receptor CXCR6 is required for the maintenance of liver memory CD8+ T cells specific for infectious pathogens. J Infect Dis. Oxford University Press; 2014;210: 1508–16. pmid:24823625

46. Cockburn IA, Amino R, Kelemen RK, Kuo SC, Tse S-W, Radtke A, et al. In vivo imaging of CD8+ T cell-mediated elimination of malaria liver stages. Proc Natl Acad Sci. National Academy of Sciences; 2013;110: 9090–9095. pmid:23674673

47. Perez-Mazliah D, Langhorne J. CD4 T-cell subsets in malaria: Th1/Th2 revisited [Internet]. Frontiers in Immunology. Frontiers Media S.A.; 2014. p. 671. pmid:25628621

48. Reece WHH, Pinder M, Gothard PK, Milligan P, Bojang K, Doherty T, et al. A CD4+ T-cell immune response to a conserved epitope in the circumsporozoite protein correlates with protection from natural Plasmodium falciparum infection and disease. Nat Med. Nature Publishing Group; 2004;10: 406–410. pmid:15034567

49. Li S, Nakaya HI, Kazmin DA, Oh JZ, Pulendran B. Systems biological approaches to measure and understand vaccine immunity in humans [Internet]. Seminars in Immunology. NIH Public Access; 2013. pp. 209–218. pmid:23796714

50. Serafini P, Carbley R, Noonan KA, Tan G, Bronte V, Borrello I. High-Dose Granulocyte-Macrophage Colony-Stimulating Factor-Producing Vaccines Impair the Immune Response through the Recruitment of Myeloid Suppressor Cells [Internet]. CANCER RESEARCH. 2004. Available: http://www.sigmaplot.com. pmid:15342423

51. Tran TM, Bijker EM, Haks MC, Ottenhoff THMM, Visser L, Schats Ret al. Whole-blood transcriptomic signatures induced during immunization by chloroquine prophylaxis and Plasmodium falciparum sporozoites. Sci Rep. Nature Publishing Group; 2019;9: 8386. pmid:31182757

52. McNab F, Mayer-Barber K, Sher A, Wack A, O’Garra A. Type I interferons in infectious disease. Nat Rev Immunol. Nature Publishing Group; 2015;15: 87–103. pmid:25614319

53. Nagai T, Devergne O, Mueller TF, Perkins DL, van Seventer JM, van Seventer GA. Timing of IFN-β Exposure during Human Dendritic Cell Maturation and Naive Th Cell Stimulation Has Contrasting Effects on Th1 Subset Generation: A Role for IFN-β-Mediated Regulation of IL-12 Family Cytokines and IL-18 in Naive Th Cell Differentiation. J Immunol. The American Association of Immunologists; 2003;171: 5233–5243. pmid:14607924

54. Doolan DL, Hoffman SL. The Complexity of Protective Immunity Against Liver-Stage Malaria. J Immunol. 2000;165: 1453–1462. pmid:10903750

55. Miller JL, Sack BK, Baldwin M, Vaughan AM, Kappe SHII. Interferon-Mediated Innate Immune Responses against Malaria Parasite Liver Stages. Cell Rep. Elsevier; 2014;7: 436–447. pmid:24703850

56. Liehl P, Zuzarte-Luís V, Chan J, Zillinger T, Baptista F, Carapau D, et al. Host-cell sensors for Plasmodium activate innate immunity against liver-stage infection. Nat Med. 2014;20: 47–53. pmid:24362933

57. Liehl P, Meireles P, Albuquerque IS, Pinkevych M, Baptista F, Mota MM, et al. Innate immunity induced by Plasmodium liver infection inhibits malaria reinfections. Infect Immun. American Society for Microbiology; 2015;83: 1172–1180. pmid:25583524

58. Kaushansky A, Metzger PG, Douglass AN, Mikolajczak SA, Lakshmanan V, Kain HS, et al. Malaria parasite liver stages render host hepatocytes susceptible to mitochondria-initiated apoptosis. Cell Death Dis. Cell Death Dis; 2013;4. pmid:23928701

59. Kurup SP, Anthony SM, Hancox LS, Vijay R, Pewe LL, Moioffer SJ, et al. Monocyte-Derived CD11c + Cells Acquire Plasmodium from Hepatocytes to Prime CD8 T Cell Immunity to Liver-Stage Malaria. Cell Host Microbe. Cell Press; 2019;25: 565–577.e6. pmid:30905437

60. Marques-Da-Silva C, Poudel B, Baptista RP, Peissig K, Hancox LS, Shiau JC, et al. Altered cleavage of Caspase-1 in hepatocytes limits control of malaria in the liver. Biorxiv. 2021;

61. Veeranki S, Duan X, Panchanathan R, Liu H, Choubey D. IFI16 protein mediates the anti-inflammatory actions of the type-I interferons through suppression of activation of caspase-1 by inflammasomes. PLoS One. PLoS One; 2011;6. pmid:22046441

62. Guarda G, Braun M, Staehli F, Tardivel A, Mattmann C, Förster I, et al. Type I Interferon Inhibits Interleukin-1 Production and Inflammasome Activation. Immunity. Cell Press; 2011;34: 213–223. pmid:21349431

63. Haque A, Best SE, De Oca MM, James KR, Ammerdorffer A, Edwards CL, et al. Type I IFN signaling in CD8-DCs impairs Th1-dependent malaria immunity. J Clin Invest. American Society for Clinical Investigation; 2014;124: 2483–2496. pmid:24789914

64. Rayamajhi M, Humann J, Penheiter K, Andreasen K, Lenz LL. Induction of IFN-αβ enables Listeria monocytogenes to suppress macrophage activation by IFN-γ. J Exp Med. J Exp Med; 2010;207: 327–337. pmid:20123961

65. de Paus RA, van Wengen A, Schmidt I, Visser M, Verdegaal EME, van Dissel JT, et al. Inhibition of the type I immune responses of human monocytes by IFN-α and IFN-β. Cytokine. Cytokine; 2013;61: 645–655. pmid:23299081

66. King T, Lamb T. Interferon-γ: The Jekyll and Hyde of Malaria. PLoS Pathogens. Public Library of Science; 2015. pmid:26426121

67. Rocha BC, Marques PE, Leoratti FM de S, Junqueira C, Pereira DB, Antonelli LR do V, et al. Type I Interferon Transcriptional Signature in Neutrophils and Low-Density Granulocytes Are Associated with Tissue Damage in Malaria. Cell Rep. Elsevier B.V.; 2015;13: 2829–2841. pmid:26711347

68. Ali S, Mann-Nüttel R, Schulze A, Richter L, Alferink J, Scheu S. Sources of type I interferons in infectious immunity: Plasmacytoid dendritic cells not always in the driver’s seat [Internet]. Frontiers in Immunology. Frontiers Media S.A.; 2019. p. 778. pmid:31031767

69. Lee AJ, Chen B, Chew M V., Barra NG, Shenouda MM, Nham T, et al. Inflammatory monocytes require type I interferon receptor signaling to activate NK cells via IL-18 during a mucosal viral infection. J Exp Med. Rockefeller University Press; 2017;214: 1153–1167. pmid:28264883

70. Costantini C, Cassatella MA. The defensive alliance between neutrophils and NK cells as a novel arm of innate immunity. J Leukoc Biol. Wiley; 2011;89: 221–233. pmid:20682626

71. Schuster IS, Coudert JD, Andoniou CE, Degli-Esposti MA. “Natural Regulators”: NK cells as modulators of T cell immunity [Internet]. Frontiers in Immunology. Frontiers Research Foundation; 2016. p. 235. pmid:27379097

72. Vivier E, Tomasello E, Baratin M, Walzer T, Ugolini S. Functions of natural killer cells [Internet]. Nature Immunology. Nature Publishing Group; 2008. pp. 503–510. pmid:18425107

73. Zhao L, Wang H, Thomas R, Gao X, Bai H, Shekhar S, et al. NK cells modulate T cell responses via interaction with dendritic cells in Chlamydophila pneumoniae infection. Cell Immunol. Elsevier; 2020;353: 104132. pmid:32446031

74. Zhang X, Ing S, Fraser A, Chen M, Khan O, Zakem J, et al. Follicular helper T cells: New insights into mechanisms of autoimmune diseases. Ochsner J. Ochsner Clinic, L.L.C. and Alton Ochsner Medical Foundation; 2013;13: 131–139. Available: /pmc/articles/PMC3603176/ pmid:23531878

75. Piccioli D, Sbrana S, Melandri E, Valiante NM. Contact-dependent stimulation and inhibition of dendritic cells by natural killer cells. J Exp Med. J Exp Med; 2002;195: 335–341. pmid:11828008

76. Hayakawa Y, Screpanti V, Yagita H, Grandien A, Ljunggren H-G, Smyth MJ, et al. NK Cell TRAIL Eliminates Immature Dendritic Cells In Vivo and Limits Dendritic Cell Vaccination Efficacy. J Immunol. The American Association of Immunologists; 2004;172: 123–129. pmid:14688317

77. van Eeden C, Khan L, Osman MS, Tervaert JWC. Natural killer cell dysfunction and its role in covid-19. Int J Mol Sci. 2020;21: 1–17. pmid:32883007

78. Mehta P, McAuley DF, Brown M, Sanchez E, Tattersall RS, Manson JJ. COVID-19: consider cytokine storm syndromes and immunosuppression [Internet]. The Lancet. Lancet Publishing Group; 2020. pp. 1033–1034. pmid:33098755

79. Li D, Chen Y, Liu H, Jia Y, Li F, Wang W, et al. Immune dysfunction leads to mortality and organ injury in patients with COVID-19 in China: insights from ERS-COVID-19 study. Signal Transduct Target Ther. Springer Nature; 2020;5: 62. pmid:32371949

80. Zheng M, Gao Y, Wang G, Song G, Liu S, Sun D, et al. Functional exhaustion of antiviral lymphocytes in COVID-19 patients. Cell Mol Immunol. Springer Nature; 2020;17: 533–535. pmid:32203188

81. Epstein JE, Paolino KM, Richie TL, Sedegah M, Singer A, Ruben AJ, et al. Protection against Plasmodium falciparum malaria by PfSPZ Vaccine. JCI insight. 2017;2: e89154. pmid:28097230

82. Pichyangkul S, Spring MD, Yongvanitchit K, Kum-Arb U, Limsalakpetch A, Im-Erbsin R, et al. Chemoprophylaxis with sporozoite immunization in P. knowlesi rhesus monkeys confers protection and elicits sporozoite-specific memory T cells in the liver. PLoS One. 2017;12: e0171826. pmid:28182750

83. Thompson EG, Du Y, Malherbe ST, Shankar S, Braun J, Valvo J, et al. Host blood RNA signatures predict the outcome of tuberculosis treatment. Tuberculosis. Churchill Livingstone; 2017;107: 48–58. pmid:29050771

84. Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, et al. STAR: Ultrafast universal RNA-seq aligner. Bioinformatics. Bioinformatics; 2013;29: 15–21. pmid:23104886

85. Anders S, Pyl PT, Huber W. HTSeq-A Python framework to work with high-throughput sequencing data. Bioinformatics. Oxford University Press; 2015;31: 166–169. pmid:25260700

86. McCarthy DJ, Chen Y, Smyth GK. Differential expression analysis of multifactor RNA-Seq experiments with respect to biological variation. Nucleic Acids Res. Nucleic Acids Res; 2012;40: 4288–4297. pmid:22287627

87. Bates D, Mächler M, Bolker BM, Walker SC. Fitting linear mixed-effects models using lme4. J Stat Softw. American Statistical Association; 2015;67: 1–48.

88. Mootha VK, Lindgren CM, Eriksson KF, Subramanian A, Sihag S, Lehar J, et al. PGC-1α-responsive genes involved in oxidative phosphorylation are coordinately downregulated in human diabetes. Nat Genet. Nat Genet; 2003;34: 267–273. pmid:12808457

89. Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, et al. Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A. Proc Natl Acad Sci U S A; 2005;102: 15545–15550. pmid:16199517

90. Korotkevich G, Sukhov V, Budin N, Shpak B, Artyomov M, Sergushichev A. Fast gene set enrichment analysis. bioRxiv. Cold Spring Harbor Laboratory; 2016; 060012.

91. Thorndike RL. Who belongs in the family? Psychometrika. Springer-Verlag; 1953;18: 267–276.

92. Mair F, Prlic M. OMIP-044: 28-color immunophenotyping of the human dendritic cell compartment. Cytom Part A. 2018;93: 402–405. pmid:29356334

93. Hertoghs N, Schwedhelm K V., Stuart KD, McElrath MJ, De Rosa SC. OMIP-064: A 27-Color Flow Cytometry Panel to Detect and Characterize Human NK Cells and Other Innate Lymphoid Cell Subsets, MAIT Cells, and γδ T Cells. Cytom Part A. 2020;97: 1019–1023. pmid:32415811