Systems analysis of protective immune responses to RTS,S malaria vaccination in humans

Challenge Model for the RTS,S/AS01 and AdVac Malaria Vaccines.

The clinical trial (NCT01366534) was conducted at Walter Reed Army Institute of Research, as described (18). Forty-six healthy malaria-naïve volunteers, randomized to two study arms, participated in this study testing the efficacy of RTS,S and AdVac malaria vaccine candidates (Fig. 1), as described (18). Study arm 1 (hereafter referred to as ARR), comprised of 25 volunteers who received the AdVac vaccine composed of Ad35 vector expressing full-length CSP, as a primary immunization, was followed by two doses of RTS,S/AS01 vaccine. The subjects in the second arm, consisting of 21 volunteers, received three doses of RTS,S/AS01 (RRR regimen). Participants in both study arms were vaccinated at 28-d intervals, and subjected to CHMI 21 d following the final immunization. Parasitemia was monitored for 28 d, and immunomonitoring continued for 159 d following challenge. The study also included 12 nonvaccinated subjects as infectivity controls. Vaccine efficacy was 44% and 52% in ARR and RR arms, respectively, and was not statistically different between the two arms (18). All subjects in the control group developed parasitemia (18).

Adaptive Immune Responses.

The RRR regimen induced significantly greater antibody titers against CSP than ARR regimen at all time points before or on the day of challenge (Fig. S1A and ref. 18). Similar results were also seen for antibody titers against HBsAg, the protein fused to CSP, although the differences at later time points were modest (Fig. S1B). Two doses of RTS,S/AS01 following the ARR were not able to induce as high a magnitude of antibody titers as two doses of RTS,S/AS01 vaccine in the RRR arm (Fig. S1A). We also assessed the number of antibody secreting cells (ASCs) induced after immunization, using ELISPOT. Both vaccines induced similar frequencies of CSP-specific (Fig. S1C), or HBsAg-specific ASCs, 6 d after the second and third immunizations. This was surprising because RRR vaccination induced a greater magnitude of CSP-antibody titers compared with ARR vaccination (Fig. S1A). This discordance may reflect differences in kinetics of the ASC response induced by ARR versus RRR. Alternatively, a different population of ASC (which was not sampled in this study) may contribute to enhanced antibody response in the RRR group.

CSP-specific CD4⁺ and CD8⁺ T-cell responses to vaccination were also assessed (18). There was negligible induction of CD8⁺ T-cell responses by RRR and a modest induction by ARR. In contrast, there was a significant induction of CSP-specific CD4⁺ T-cell response by ARR, and to a much weaker degree by RRR (18). The functionality of T cells was monitored by FACS analysis using a panel including four markers: CD40L, IL-2, TNFα, and IFNγ. In the ARR vaccine group, there was a markedly enhanced frequency of polyfunctional (expressing three or four functions) CSP-specific CD4⁺ T cells at D14, D42, and D77, and postchallenge (Fig. S1D).

Immunologic Correlates of Protection.

As reported (18), in the RRR vaccine arm, individuals who did not develop parasitemia within 28 d after challenge (referred hereafter as “protected”) had higher concentration of anti-CSP antibodies at the time of challenge than nonprotected individuals (Fig. 2A). In the ARR arm, the concentration of anti-CSP antibodies was substantially lower than that in the RRR arm (Fig. S1A), and there was no statistically significant difference in the titers on the day of challenge, between the protected versus nonprotected subjects. Anti-HBsAg antibody concentrations were not significantly different between protected and nonprotected individuals (Fig. 2B). In the ARR arm, it was in fact the frequency of CSP-specific polyfunctional CD4⁺ T cells at day 14 that significantly correlated with protection (Fig. 2C). The frequencies of CSP-specific polyfunctional T cells were similar between protected and nonprotected groups at all later time points (Fig. 2C). The frequency of polyfunctional CD4⁺ T cells did not correlate with protection in the RRR arm (Fig. 2C).

Transcriptional Signatures Induced by Vaccination.

Vaccination with ARR or RRR induced potent transcriptional responses in PBMCs, with several thousands of genes being induced or repressed (Fig. 3A). The transcriptional responses at day 1 and day 6 after vaccinations, which include signatures of inflammation/TLR signaling and cell cycle genes in proliferating ASCs, correspond to the early innate and the later ASC responses (Fig. S2 A and B). ARR and RRR induced a small number of differentially expressed genes at D1 after primary vaccination (Fig. S2C). The genes more strongly induced by ARR included genes associated with the type I IFN antiviral response and innate immunity, such as IFI27, IFI44L, IFI6, and HESX1 (27), consistent with previous studies (28). We then investigated the regulation of known IFN type I response-associated genes (27), as well as genes up-regulated in response to the live attenuated virus YF-17D vaccine (24). Both Ad35.CS and RTS,S/AS01 primary vaccinations induce potent expression of IFN type I and YF-17D signatures (Fig. S3), suggesting that the virus-like particles and the AS01 adjuvant contained in RTS,S/AS01 induces a potent antiviral type I IFN response, similar to that observed with viruses such as Ad35 or YF-17D. The identity of such genes induced by Ad35.CS and RTS,S/AS01 was largely overlapping (Fig. S3).

To identify functional pathways perturbed by the two regimens, we used Gene Set Enrichment Analysis (GSEA) using blood transcription modules (BTMs) (29) as gene sets. Transcripts were ranked according to a fold change differences relative to the D0 baseline. The functional responses elicited by the two vaccines were broadly similar (Fig. S4 and Dataset S1). Both vaccines induced strong innate responses (Fig. S4), including inflammatory/TLR/chemokines BTMs, following each vaccination (Fig. S4). Enrichment of cell cycle and plasma and B-cell–related BTMs was observed 6 d after each vaccination (Fig. S4). A noticeable difference was the contraction of B-cell and plasma cell BTMs at D2 following prime immunization, observed in the ARR regimen, but absent in the RRR arm (Fig. S4). Interestingly BTMs related to cell cycle were enhanced even at D14 after primary vaccination, suggesting the persistence of cycling cells. Furthermore, in both arms, we observed a repression of BTMs related to NK cells at D1 following each vaccination (Fig. S4).

Molecular Signatures of Immunogenicity.

We next analyzed the transcriptional signatures that correlated with immunogenicity of vaccination. In the case of RRR vaccination, we assessed transcriptional correlates of CSP-specific antibody titers on D77 (the day of challenge). Following RRR vaccination, the expression of BTMs related to plasmablasts at D1 after each vaccination was positively associated with the antibody titers on the day of challenge (Figs. S5A and S6A). This observation was surprising, given that the plasmablast response in humans to vaccination with other vaccines such as influenza (30) has been shown to peak at day 7 after vaccination, and with ARR and RRR vaccination robust, plasmablast responses were observed 6 d after each boost (Fig. S1C). The observed correlation between BTMs related to B cells and plasmablasts, at day 1 after each boost, and immunogenicity might reflect a transient burst of genes related to B-cell activation within a day of vaccination, but this hypothesis needs further exploration. Additionally, cell division BTMs showed positive correlation to the antibody titers even at later time points (D14, 28, and 56), suggesting the persistence of cycling cells (Figs. S5A and S6A). Furthermore, the expression of several innate immunity modules (antigen presentation M95, dendritic cell activation M165), including many antiviral and type I IFN-related modules at day 6 post primary and secondary vaccinations, were positively correlated with the antibody titers on the day of challenge (Figs. S5A and S6A and Dataset S2). Most strikingly, on the day of first and second boosts, gene modules relevant to NK cells showed strong negative correlation to the antibody titers at the day of challenge (Fig. S5A and Dataset S2). Indeed, we observed that the majority of genes included in these NK cell-related BTMs showed negative association with antibody titers.

In the case of ARR, the frequency of polyfunctional CD4⁺ T cells at day 14 was associated with protection (Fig. 2C). We thus assessed whether early transcriptional signatures correlated with the polyfunctional CD4⁺ T-cell response at day 14. At day 1, several modules representative of innate immune activation (antigen presentation M71, M95.1; activated dendritic cells and monocytes M168, M11; TLR and inflammatory responses M16, M25, M146) were strongly associated with polyfunctional CD4⁺ T-cell response at day 14 (Figs. S5B and S6B). Interestingly, modules representative of respiratory electron transport were strongly associated with the response. In contrast, modules representative of NK cells and T cells were negatively associated with the response. Similar, but weaker, associations were observed at D2 after prime Ad35.CS immunization. By day 6, the landscape of correlates changed, with many modules representing DC markers becoming negatively enriched, whereas NK cell modules continued to be negatively associated (Figs. S5B and S6B and Dataset S2).

Association of Molecular Signatures with Protection.

We then analyzed transcriptional signatures associated with protection. In the RRR group, BTMs related to plasma cells, B cells, and cell cycle at D1, D29, and D57 (i.e., 1 d after each vaccination) were positively associated with protection (Fig. 4 and Fig. S7A), consistent with the correlations between the expression of such BTMs at 1 d after each booster vaccination, and antibody titers (Figs. S5A and S6A). We also observed positive associations of multiple innate immunity modules at 6 d after primary and secondary vaccination, similar to the observed correlations with CSP-specific antibody titers (Fig. 4 and Fig. S7A). Additionally, several NK cell modules at D56 (day of the second boost) negatively associate with protection (Fig. 4 and Fig. S7A), consistent with their correlation with antibody titers (Figs. S5A and S6A). Strikingly, there were negative correlations of the expression of almost all of the genes contained within the NK cell-related BTMs, at D56, and protection (Fig. 5).

The transcriptional signatures of protection for the ARR arm were different from those for RRR. Here, multiple innate immunity modules positively associate with protection at D1 and 2 after the prime, and D28 (day of the first boost) (Fig. S7 B and C), similar to the transcriptional correlates of polyfunctional CD4⁺ T cells described above (Figs. S5B and S6B). Again, NK modules display strong negative association with protection at multiple time points (D2, D28, D29, D56) (Fig. S7 B and C and Dataset S3).

We then determined the overlap between the molecular signatures of protection and immunogenicity. For the ARR arm, at D1 and 2 BTMs related to antigen presentation, TLR signaling and dendritic cells were associated with both immunogenicity (i.e., polyfunctional CD4⁺ T cells at D14) and protection (Fig. S7C). For RRR, there was considerable overlap between signatures of protection and immunogenicity. Several BTMs that correlate with protection were also correlated with immunogenicity (Fig. 4). BTMs related to plasma and B cells, and the cell cycle were correlated with both protection and immunogenicity at 1 d after the primary and secondary vaccinations (Fig. 4). In contrast, several innate immunity modules, at day 6 after prime and day 6 after boost (D34), were correlated with both anti-CSP–specific antibody response and protection (Fig. 4). Strikingly, at D56 (the day of the final boost), we observed negative correlations of several NK cell-related BTMs with protection and immunogenicity (i.e., CSP-specific antibody titers at D77) (Fig. 4). A full description of all common associations with immunogenicity and protection is provided in Dataset S4.

Predictive Modeling of Protection.

We then developed predictive signatures of protection based on the transcriptional response to RRR vaccination. To achieve this goal, a discriminant analysis via mixed integer programming (DAMIP) (31) was used to generate the candidate predictive signatures. The two groups to classify are group 0 (protected) versus group 1 (not protected). For signature validation, we used a transcriptional dataset from an independent malaria challenge study with RTS,S/AS01 (32), NCT00075049, hereby referred to as “Vahey data set.” Responses in the two studies were broadly similar (Fig. S8 and Datasets S5 and S6). The baseline-normalized expression values for the RRR cohort in the present study was used as a training set, and candidate signatures that passed a minimum accuracy threshold in the training set were then applied to the independent validation set for blind prediction. The outline of the predictive modeling experiments is provided in Fig. 6A. A full list of signatures and their performance metrics is provided in Datasets S7 and S8. Analysis of the transcripts included in the successful predictive signatures revealed a high prevalence of transcripts that were commonly found in a large number of signatures at D56 (Table S1). One such transcript, KIR2DS1 (an NK cell marker), was found in 57 of 99 successful predictive signatures. This observation is consistent with the fact that NK cell-related BTMs are negatively associated with protection in both datasets at D56 (Fig. 6B). To validate these signatures, we used the set of transcripts identified in RRR to generate and train signatures in the Vahey data set (Dataset S9). Many of the mRNAs that were highly represented in signatures trained in RRR arm of this study are also highly represented in signatures trained in the Vahey data set (Datasets S10 and S11). Up-regulated mRNA common to both RRR- and Vahey-generated signatures include several NK markers KIR2DS1, KIR2DL2, and KIR3DL1. Notably, many of the mRNAs that were included in the predictive signatures in 10-fold cross validation (10× CV) in RRR and the Vahey data set individually were also included in predictive signatures that were trained by using RRR, and were shown to blind predict outcome in Vahey data set (Dataset S11). Therefore, we conclude that a small number of mRNAs with high prevalence in predictive signatures are likely to be determining factors that distinguish protected versus nonprotected individuals.

Finally, we illustrated the ability of the generated signatures of protection to segregate the samples in their respective protection groups. For this analysis, we used signatures consisting of features that include the high-prevalence genes noted above and plotted the distribution of protected and nonprotected subjects as a function of baseline normalized expression values of genes contained in these signatures. The results for three representative signatures are shown in Fig. 6C. Notably, although the overall accuracy of prediction was 80% or higher in both protected and unprotected groups, there were specific individuals in the Vahey validation set that were consistently misclassified. One of these individuals was misclassified by more than 90% of signatures, whereas four others were classified into one or the other group with nearly random probability. These subjects are indicated by black dots in Fig. 6C. In summary, we confirmed that expression of genes included in representative predictive signatures is sufficient to segregate protected and nonprotected subjects.

GeneArray Probe Level Processing.

Total mRNA was isolated from frozen PBMCs provided by Walter Reed Army Institute for Research by using Quagen RNAeasy kit (Qiagen) and stored at −80 °C. RNA was then quantified and checked for integrity by using Agilent BioAnalyzer (Agilent Technologies). For label preparation, 50 ng of total mRNA were amplified and labeled by using NuGEN Ovation WB target labeling kit (NuGEN) and hybridized to HU-133 plus 2.0 GeneChip (Affymetrix). Samples were processed in batches of 96; special care was taken to ensure that all time points from a particular subject were processed in the same batch. Subjects were assigned to batches in a way that balanced the distribution according to arm and protection status. Two pooled references were included in each batch to warrant against any batch effects.

Posthybridization QC analysis was performed by using Bioconductor packages arrayQualityMetrics (www.bioconductor.org/packages/release/bioc/html/arrayQualityMetrics.html) and AffyQCReport. Spike controls were evaluated following RMA normalization. Based on probe level QC analysis and RNA quality control, 56 of 639 samples were removed from analysis. Background correction on the remaining samples was performed by GC-RMA, with probes summary by median polish, followed by log₂ transform.

Procedures were implemented to detect the contribution of batch to the overall variance (batch effect). No batch effect was observed.

Microarray data have been submitted to National Center for Biotechnology Information GEO (accession no. GSE89292).

Significance Testing.

Significance testing comparing the gene expression at each time point to the baseline at D0 was performed by paired two-sided t test, followed by Benjamini–Hochberg FDR correction as described (24). Genes with q values <0.05 and absolute fold change greater than 1.5 were considered to be significantly regulated.

GSEA Analysis.

For GSEA analysis (37), either fold change with respect to baseline (D0) or correlation coefficients were used as a ranking metric. For genes represented by multiple probe sets, the probe set with the greatest extreme value of the ranking metric was used. GSEA was run with 10,000 permutations by using Java interface (Broad Institute, software.broadinstitute.org/gsea/index.jsp) or standalone R code (37). When BTM modules were used, modules with fewer than 10 gene identifiers were ignored.

SPICE Analysis.

Analysis and presentation of distributions was performed by using SPICE version 5.3, downloaded from https://exon.niaid.nih.gov (37). Comparison of distributions was performed by using a Student's t test and a partial permutation test as described (37). Briefly, a χ² like test is used to test the null hypothesis that the distribution of measurements of any of all possible combinations of cytokines (for four cytokines that would equal to 15, not including all-negative subset) is drawn from the same distribution for two groups (for example, ARR vs. RRR at time = D14), and is not significantly different between the two comparison groups. Once the measurement of χ² metric is obtained, we need to test whether this value represents a significant difference or can be obtained stochastically. To test this possibility, we use a partial permutation test. In this approach, all samples from the two groups being compared are pooled into one set, and groups to be compared are reformed by randomly assigning each sample to one or the other group to be compared. The χ² statistic is calculated for each of these random permutations. After a large number (1,000–10,000) of permutations, it becomes possible to estimate whether the values of χ² statistic obtained on the original (nonpermuted) combination of samples could be obtained with certain frequency by randomly permuting the samples, that is, stochastically. This frequency forms the P value for the test (for instance, P = 0.05 indicates that there is no more than 5% probability that a value of χ² test obtained in this comparison could be achieved by a stochastic (random) assignment of samples to the comparison groups.

Enrichment Analysis.

Enrichment analysis was performed by using Gather (changlab.uth.tmc.edu/gather/gather.py) or Reactome (www.reactome.org).

Predictive Modeling.

Predictive models were built using DAMIP. Application of DAMIP to vaccine clinical studies has been described earlier (24, 25). Briefly, subjects in each arm were separated into training and blind test sets. The criteria for subject selection were (i) no fewer than 10 subjects must have nonmissing data at each time point (because of some samples being removed after QC) and (ii) the ratio of protected and nonprotected subjects in the training set is closely matched the ratio in the overall set of subjects. Feature selection and model training were performed through 10-fold cross-validation (10× CV) loops in the training set. The predictive models/rules/signatures were selected by using a predefined accuracy cutoff for predicting the protection status (typically in the ≥80% accuracy range). Models passing filtering criteria in the training set were then evaluated in the blind test set. Models passing filtering criteria in the blind test set were reported as candidate signatures of protection.

Analysis of ASC Responses.

Analysis of ASC responses was performed as described (30).

Analysis of CSP-Specific T-Cell Responses.

FACS analysis of CSP-specific CD4⁺ and CD8⁺ T-cell responses were performed as described (18).