Distinct immune responses associated with vaccination status and protection outcomes after malaria challenge

JournalFeeds

Authors:
Damian A. Oyong , et al.

Abstract

Understanding immune mechanisms that mediate malaria protection is critical for improving vaccine development. Vaccination with radiation-attenuated Plasmodium falciparum sporozoites (PfRAS) induces high level of sterilizing immunity against malaria and serves as a valuable tool for the study of protective mechanisms. To identify vaccine-induced and protection-associated responses during malarial infection, we performed transcriptome profiling of whole blood and in-depth cellular profiling of PBMCs from volunteers who received either PfRAS or noninfectious mosquito bites, followed by controlled human malaria infection (CHMI) challenge. In-depth single-cell profiling of cell subsets that respond to CHMI in mock-vaccinated individuals showed a predominantly inflammatory transcriptome response. Whole blood transcriptome analysis revealed that gene sets associated with type I and II interferon and NK cell responses were increased in prior to CHMI while T and B cell signatures were decreased as early as one day following CHMI in protected vaccinees. In contrast, non-protected vaccinees and mock-vaccinated individuals exhibited shared transcriptome changes after CHMI characterized by decreased innate cell signatures and inflammatory responses. Additionally, immunophenotyping data showed different induction profiles of vδ2+ γδ T cells, CD56+ CD8+ T effector memory (Tem) cells, and non-classical monocytes between protected vaccinees and individuals developing blood-stage parasitemia, following treatment and resolution of infection. Our data provide key insights in understanding immune mechanistic pathways of PfRAS-induced protection and infective CHMI. We demonstrate that vaccine-induced immune response is heterogenous between protected and non-protected vaccinees and that inducted-malaria protection by PfRAS is associated with early and rapid changes in interferon, NK cell and adaptive immune responses.

Trial Registration: ClinicalTrials.gov NCT01994525.

Introduction

Malaria remains a significant threat to global health, causing an estimated 241 million cases and 627,000 deaths in 2020 [1]. Out of the five parasite species causing malaria, infection by Plasmodium falciparum accounts for the majority of global malaria cases and up to 95% of all cases in Africa [1]. The recent pandemic caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV2) infection has also exacerbated malaria health burdens by disrupting treatment services, increasing the number of malaria-related mortality [1]. The development of an effective anti-malaria vaccine is critically needed to prevent malaria and reduce its health burdens that disproportionally affect young children and pregnant women. To date, there is only one malaria vaccine recommended by the WHO for widespread use among sub-Saharan African children, which is the RTS,S/AS01 subunit vaccine. Nonetheless, vaccine efficacy for RTS,S/AS01 is very modest, at about 26–36% when assessed in children during a phase 3 field trial [2]. Although another similarly designed subunit vaccine showed better efficacy in a phase 2b trial, a wider age range, larger number of participants, and assessments across different regions and transmission settings are required to fully assess its effectiveness [3,4]. A promising alternative to subunit vaccines are attenuated whole P. falciparum sporozoite vaccines which include radiation-attenuated sporozoite (RAS); wild-type sporozoite administration under drug cover, also known as CPS (chemical prophylaxis with sporozoite) or ITV (infection treatment vaccination) [5]; and genetically attenuated sporozoite (GAP) approaches [6,7]. RAS arrests early in liver-stage development and, historically, have been highly effective in inducing sterile protective immunity in animal models and humans, with a >90% efficacy rate [5,811]. However, the efficacy of RAS vaccines in malaria endemic regions may be less strong, as reduced efficacy has been observed in pre-exposed individuals [1214]. The human immune system is critical in determining the outcome of sterile immunity against malaria, and our incomplete understanding of the dynamics and mechanisms of natural or vaccine-induced immune responses to malaria is hampering the development of effective vaccines.

Identifying the mechanistic correlates of malaria protection have remained elusive, in part due to the complexity of the parasite’s multi-stage life cycle and host-parasite interactions [15]. Studies of RAS vaccines have associated sterilizing protection with multiple adaptive and innate immune cell types, including CD4+ T cells [16], CD8+ T cells [9,1719], γδ T cells [20,21], NK cells [22], and dendritic cells (DCs) [23]. For example, abrogation of CD4+ T cells [16] and γδ T cells [21] in RAS-immunised mice resulted in a loss of malaria protection. Also, a higher frequency of CD8+ T cells secreting IFN-γ has been previously observed in individuals protected against controlled human malaria infection (CHMI) challenge compared to those who were non-protected [18]. Indeed, higher interferon-related transcription signatures have been associated with malaria protection [24]. Although the practicality of using RAS for primary vaccinations is uncertain due to logistical and large-scale manufacturing challenges [25,26], it remains a valuable tool for the analysis of immunity to malaria. For example, suboptimal dosing with RAS can produce different outcomes of protection or no protection when assessed with CHMI. This can be harnessed as an intentional strategy to investigate correlates of malaria protection.

In this study, we examined whole blood immune responses after malaria challenge in volunteers receiving P. falciparum RAS vaccine. Trial participants either received 5 RAS immunizations, dosed to achieve ~50% protection outcome, or uninfected mosquito bites (mock vaccination), followed by CHMI to assess malaria protection. To characterize vaccine-induced protection, we analyzed blood samples at sampling time points prior to and after CHMI and then utilized systems biology analyses of whole blood transcriptomic and cellular immune profiles. We provide evidence that malaria protected vaccinees induce distinct immune responses, including interferon and adaptive responses, and with different kinetics following CHMI compared to non-protected vaccinees and mock-vaccinated individuals.

Results

Study participants

Blood samples from volunteers enrolled in Cohort 1 of the Immunization by Mosquito bite with RAS (IMRAS) trial [27] were the focus of this study. Briefly, individuals in the IMRAS trial received 5 doses, at approximately monthly intervals, P. falciparum radiation attenuated sporozoites (PfRAS) by bites of infected mosquitoes (PfRAS vaccinees, n = 11) or non-infected mosquitoes (mock, n = 5). The number of infectious bites delivered was selected to achieve incomplete protection in the study cohort, with a target of ~50% protection outcome, to facilitate assessment of correlates of protection. An average of 1,027 infectious bites were administered to Cohort 1 over five vaccination doses. Protection from malaria was assessed with CHMI by bites from 5 mosquitoes infected with homologous P. falciparum strain (Fig 1A). Of the PfRAS vaccinated group, 5 individuals developed blood-stage parasitemia following CHMI (non-protected, n = 5) while 6 others remained negative (protected, n = 6). All mock-vaccinated individuals developed blood-stage infection after CHMI. Blood parasitemia was initially assessed with blood smear microscopy and retrospectively confirmed using qPCR (Fig 1B). Blood-stage infection was detected between 7–11 days after CHMI in mock and non-protected individuals (Fig 1B).

thumbnail

Fig 1. Immunization by Mosquito bite with RAS (IMRAS) trial.

A) IMRAS trial timeline and sampling time points. Volunteers received 5 vaccine doses (V1-V5), followed by CHMI. The two vaccination groups are shown: Immunized (n = 11) = PfRAS vaccination, mock (n = 5) = non-infected mosquito bites. Time points indicate sampling times used for this study. B) Kaplan-Meier analysis showing parasitemia-free survival curves based on qPCR detection of parasites in peripheral blood.


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

Whole blood transcriptome profiles were obtained for samples taken at pre-vaccination baseline (B0), day 0 pre-CHMI (C0), day 1 post-CHMI (C1), and day 7 post-CHMI (C7). High dimensional flow cytometry and single-cell RNA sequencing (scRNA-seq) for immune cell profiling was carried out on PBMC samples obtained via leukapheresis at pre-vaccination baseline (B0), day 6 post-CHMI (C6), and day 112 post-CHMI (C112) (Fig 1A and S1 Table).

Single-cell RNA-sequencing elucidates the dynamics of immune cell profiles responding to CHMI

We first investigated responses to Plasmodium infection by applying CITE-seq profiling of single cell transcriptomics and selected protein markers to PBMC samples from 4 mock-vaccinated individuals at pre-vaccination baseline (B0), and 6 and 112 days (C6, C112) post-CHMI. Unsupervised clustering of the samples revealed a total of 41 distinct cell clusters (Fig 2A). Surface protein and gene expression profiles were used to assign cell types to clusters (S1 Fig and S1 Data). Longitudinal changes in cluster frequency were then calculated using mixed effect linear regression analysis, controlled for within-individual variation. Five clusters significantly changed in frequency over time (FDR<0.01), responding to CHMI (Fig 2B). Cluster 20 corresponds to CD4+ Treg with elevated expressions of CD4 and CD25 surface proteins and FoxP3 gene (Fig 2C). Cluster 21 was identified as a recently activated CD4+ effector T cell (Teff) with high CD4 and CD5 and low or negative expression of CD25, CD62L, and CD127 surface markers and elevated cell proliferation and growth genes such as RPL4 and EEF1A1 (Fig 2C). Both cluster 27 and 34 expressed CD56, and therefore were annotated as NK cells. Proliferating state on cluster 34 and NKT cell annotation were determined with high expression of CD38 and other proliferating genes such as PCNA, PLCLAF, MKI67, and TOP2A, while also expressing the CD8 (Fig 2C). Lastly, cluster 28 was identified as plasmacytoid dendritic cell (pDC) with high surface expression of HLA-DR, CD45RA and CD123 accompanied by HLA-DRB1, JCHAIN, and FCER1A gene expression (Fig 2C).

thumbnail

Fig 2. Immune cell profiling to CHMI in mock-vaccinated individuals using single-cell RNA sequencing.

A) Uniform manifold approximation and projection (UMAP) plot indicating clusters of distinct cell subsets. Numbers indicate cell clusters. High annotation of cell subsets is indicated by the dashed circled lines. B) Cell clusters that significantly changed in frequency overtime. Bar plots indicate significantly different cell clusters. Line plots indicate frequency of the significant clusters overtime at baseline pre-vaccination (B0), 6 days post-CHMI (C6), and 112 days post-CHMI (C112). Statistical tests were performed using linear mixed-effect regression with FDR-adjusted P value <0.01. C) Gene and protein expression profile for each of the significant clusters. Dot colors indicate average expression level while dot sizes indicate percent of marker expression in the corresponding cluster. D) Heatmap plot showing GSEA on comparisons between time points, individually tested in each of the significant clusters. Each square indicates a gene set. Color represents statistically significant gene sets and normalized enrichment score (NES) obtained from GSEA. All gene expression changes were used as inputs for GSEA. Threshold for significant gene sets was set to an FDR-adjusted P value of <0.01. Four (n = 4) mock-vaccinated individuals were selected.


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

A decrease in cluster frequency following CHMI at the C6 time point was observed for all significant clusters with the exception of the pDC cluster, which had an increase in frequency (Fig 2B). Frequency of all five significant clusters was higher at 112 days post-CHMI (C112) compared to B0 and C6 (Fig 2B), possibly indicating a new baseline immune state post-CHMI or a long-term memory response. We also observed an increased frequency in the γδ T cell cluster (Cluster 15) at C6 and C112 compared to B0 time point, albeit not statistically significant (FDR<0.01). We next performed gene set enrichment analysis (GSEA) [28] on gene expression changes between time points in each of the significant clusters. We determined enrichment of MSigDB Hallmark gene sets comprising sets of coherently expressed genes in blood that represent well-defined biological process [29]. Overall, there was a significant enrichment in gene sets associated with TNFα, IFNα, IFNβ, and inflammatory responses in the cell clusters, most notably in the activated CD4+ Teff at C6 vs B0 time comparison (Fig 2D). Both NK cell clusters, CD56hi NK cell (at C6 vs B0) and proliferating NK/NKT cell (at C112 vs B0) were significantly enriched for IFNa and IFNγ responses (Fig 2D). Additionally, TNFa and IFNγ gene sets were also significantly enriched in the pDC cluster at both C6 vs B0 and C112 vs B0 time comparisons. Overall, scRNA-seq showed that both innate and adaptive immune cells were responding to CHMI in mock-vaccinated individuals and that these immune cells were characterized with interferon and inflammatory responses.

Differential gene expressions in response to PfRAS vaccinations and CHMI

We next investigated transcriptomic response to PfRAS vaccinations and CHMI by selecting both PfRAS-vaccinated and mock-vaccinated individuals. Whole blood samples were profiled with RNA sequencing and differentially expressed genes (DEGs) were quantified by examining gene expression changes over time in each vaccination group; protected, non-protected, and mock. DEGs were identified using mixed effect linear regression to account for within-individual variation with an FDR threshold set of 0.2. The pre-CHMI (C0) time point (3 weeks after the fifth RAS vaccination) was compared to the pre-vaccination baseline (B0) time point to determine changes induced by PfRAS vaccination. A total of 52 genes (17 increased, 35 decreased) and 106 genes (15 increased, 91 decreased) were differentially expressed in protected and non-protected vaccinees, respectively (Fig 3A). The greatest differential expression at C0 compared to B0 was observed in non-protected vaccinees, possibly indicating immune status change in response to vaccination. 21 genes (11 increased, 10 decreased) were differentially expressed significantly at C0 compared to B0 in the mock group. The transcriptomic changes induced by CHMI compared to C0 were determined at the 1 day (C1) and 7 day (C7) post-challenge time points. Within protected vaccinees, 52 genes were differentially expressed between C1 and C0 while 43 DEGs were observed at C7 relative to C0 (Fig 2C). Similarly, a total of 85 and 45 DEGs were observed in non-protected vaccinees at C1 and C7 relative to C0 (Fig 3A). In contrast to PfRAS vaccinees, a small number of genes were differentially expressed in mock individuals following CHMI (Fig 3A); 12 DEGs and 18 DEGs over C1 and C7 relative to C0.

thumbnail

Fig 3. Gene expression changes to PfRAS vaccination and CHMI.

A) Number of differentially expressed genes (DEGs) calculated using mixed effect linear regression between two time points and for each sample group. The number of DEGs for each group and time point comparisons is indicated above plotted bars. B) Volcano plot illustrating DEGs between sample groups at individual time points. Comparisons were determined between P and NP vaccinees and between mock individuals and PfRAS vaccinees, P and NP. Colored dots indicate DEGs after multiple testing correction while black dots indicate non-statistically significant (NS) comparisons. -Log10 P values shown on the plots are prior to multiple testing correction. Horizontal line indicates P = 0.05. Threshold for DEGs was set at an FDR-adjusted P value of <0.2. P (n = 6) indicates PfRAS-vaccinated protected individuals, NP (n = 5) indicates PfRAS-vaccinated non-protected individuals, Mock (n = 3) indicates non-infected mosquito bites vaccinated individuals.


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

Further analysis was done to explicitly compare 1) protected versus non-protected vaccinees and 2) mock versus vaccinated individuals at each time point. The greatest number of DEGs (40 genes) at C7 was observed when comparing protected and non-protected vaccinees (Fig 3B). This differential response coincides with the first detection of parasites in peripheral blood in non-protected vaccinees (Fig 1B), which likely indicates transcriptomic changes associated with the onset of blood-stage infection. At C0 and C1, only 1 and 2 DEGs were observed between the two PfRAS vaccinee groups, respectively. For comparisons between mock and the PfRAS vaccinee groups, the largest differences were observed at C1 with 129 and 87 DEGs observed on comparisons against non-protected and protected vaccinees, respectively (Fig 3B). Notably, more DEGs were seen comparing mock and PfRAS vaccinee groups than between protected and non-protected vaccinees at any time point, suggesting that vaccination was influencing immune response to an infection in non-protected vaccinees. Of note, the genes B4GALNT3 and TUBB2A were consistently reduced in expression in non-protected vaccinees while TREML4 expression was higher in mock individuals compared to the other groups, across all time points, pre- and post-CHMI (Figs 3B and S2), suggesting inherent differences at the genetic level that are not induced by vaccination nor CHMI. These genes are involved in protein N-glycosylation (B4GALNT3), intracellular transport (TUBB2A), and modulation of inflammation and immune responses (TREML4).

Coherently expressed gene sets associated with PfRAS vaccination and protection from malaria

To identify functional responses and biological processes represented by the gene changes observed following PfRAS vaccination and CHMI, we performed an independent GSEA on all gene expression changes using published coherently expressed blood transcriptional modules (BTMs) [30,31] (S2 Data). Transcriptional responses to PfRAS vaccination, determined by comparing gene expression at C0 relative to B0, were associated with decreased expression of cell cycle and T cell annotated BTMs, and increased expression of erythrocyte and monocyte BTMs in both protected and non-protected vaccinees (Figs 4A and S3). However, increases in interferon and NK-cell-associated responses were observed specifically in protected vaccinees (Fig 4B). Interferon signaling has been previously associated with sterile protection after PfRAS vaccination [9]. In non-protected vaccinees, inflammation and monocyte responses were amongst the top upregulated BTMs in the C0 to B0 time interval (S3 Fig). In contrast to the PfRAS vaccinees, a lesser number of significantly enriched BTMs were observed in the mock group, with changes in immune responses at C0 compared to B0 associated with increased platelet and erythrocyte-related BTM expression. The observed higher number of significantly enriched BTMs in protected (60 BTMs) and non-protected (86 BTMs) vaccinees compared to mock (33 BTMs) individuals prior to CHMI may indicate transcriptomic changes in response to PfRAS vaccinations. Overlap in the directionality of the BTM response was moderately distributed between sample groups (Fig 4C), indicating unique immune response patterns in each group prior to CHMI.

thumbnail

Fig 4. Blood transcriptional module profile on responses to PfRAS vaccination and CHMI.

A) Gene set enrichment analysis (GSEA) on comparisons between time points, individually tested in each sample group. Each square indicates enrichment of a single of blood transcriptional module (BTM). Colored row annotations represent high-level annotations of the BTMs. Square color represents statistically significant gene sets and normalized enrichment score (NES) obtained from GSEA. All gene expression changes were used as inputs for GSEA. Threshold for significant gene sets was set at an FDR-adjusted P value of <0.01. B) Average gene expression of representative high-level BTM annotation. Lines indicate each individual. Dashed lines indicate locally estimated scatterplot smoothing (LOESS) regression with 95% confidence interval shown in the highlighted color. Transcript abundance is in log-transformed counts per million (logCPM). C) Venn diagrams showing proportions of shared significant BTMs, with identical NES direction, between sample groups at each time interval. P (n = 6) indicates PfRAS-vaccinated protected individuals, NP (n = 5) indicates PfRAS-vaccinated non-protected individuals, Mock (n = 3) indicates non-infected mosquito bites vaccinated individuals.


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

Following CHMI, immune profiles represented by the BTMs were strikingly different between protected and non-protected vaccinees. Between C0 and C1, the protected group showed increases in several interferon-related BTMs and decreases in adaptive immune BTM expression including B and T cells (Figs 4A and 4B and S3). In contrast, responses to CHMI in the non-protected group during this time were marked by increases in NK cell-associated BTMs and decreases in BTMs associated with inflammation, neutrophils, and monocytes (Figs 4A and 4B and S3). Notably, BTM response profiles in non-protected vaccinees following CHMI seemed to overlap more closely with mock individuals versus protected vaccinees in the C1 vs C0 time interval, with 13 BTMs sharing identical response directionality between non-protected vaccinees and mock individuals, while only 1 BTM response was shared between non-protected and protected vaccinees (Fig 4C). This similarity in BTM responses following CHMI between non-protected vaccinees and mock individuals continued at time point C7, when blood-stage infection was detected and with shared BTMs representing increased expression in T cells, NK cells, cell cycle, and other biological processes and decreased expression in inflammation, interferon, neutrophils, and monocytes (Fig 4A). These immune responses were consistent with our scRNA-seq data, where cell clusters that were reduced in frequency at C6 also primarily mediate interferon and inflammatory responses (Fig 2B). Only a few BTMs were significantly enriched at C7 compared to C0 in protected vaccinees, and none of these BTMs shared responses of similar directionality to non-protected vaccinees or mock individuals (Fig 4C). These included erythrocyte and cell cycle associated BTMs, perhaps indicating a return to the pre-CHMI state following parasite clearance in protected vaccinees. Additional GSEA using the Hallmark gene sets revealed that both type I (IFNα) and type II (IFNγ) were increased in protected vaccinees between C1 and C0 time point while gene sets associated with inflammation were decreased in non-protected vaccinees and mock individuals at C7 relative to C0 (S4 Fig).

To identify protection-associated genes, we selected leading-edge genes that were present in more than half of the significant gene sets from the interferon and NK cell-associated BTMs in protected vaccinees. A total of 5 and 3 genes were present in more than half of the significant interferon BTMs at C0 vs B0 and C1 vs C0 time comparison, respectively (S5A Fig). These genes include HERC5, IFIH1, IFIT1, IRF7, and RSAD2 and were all associated with interferon response. Consistently, expression level of the interferon genes was elevated following PfRAS vaccination and peaked at 1 day following CHMI in the protected vaccinees (S5B Fig). NK cell-associated BTMs were significant only at C0 vs B0 time comparison in protected vaccinees, of which IL2RB and TGFBR2 gene were present in 4 of the 5 significant BTMs (S5A Fig). Eight other genes including FASLG, GPR56, KLRB1, NKG7, S1PR5, and the KIR gene family (KIR2DL1, KIR2DL3, and KIR3DL1), were present in 3 significant BTMs. IL2RB and KIR gene expression levels peaked at C0 in protected vaccinees (S5B Fig). Interestingly, expression level of the regulatory gene, TGFBR2, was reduced at C0 and peaked at C1 in protected vaccinees.

Changes in biological pathways and deconvoluted cell responses following PfRAS vaccination and CHMI

To further understand the pathways that lead to immune processes observed in the GSEA, we selected leading-edge genes from the significantly enriched BTMs and functionally profiled them using Ingenuity Pathway Analysis (IPA) [32]. Comparing changes over the full course of vaccinations, i.e. between B0 and C0, IPA predicted activation in pathways related to interferon signaling and NK cells in protected vaccinees (Fig 5A and S3 Data). By comparison, pathways associated with innate cells, including phagocytes and monocytes were activated in non-protected vaccinees, and EIF2 signaling, a pathway related to protein synthesis, was activated in mock individuals in the same time interval. One day following CHMI, at C1 relative to C0, pathways associated with protein synthesis (EIF2) and mitochondrial activity (mTOR) were significantly activated in protected vaccinees. In non-protected and mock groups, innate and phagocyte related pathways, including phagosome formation, NK cell signalling, Fcy receptor-mediated phagocytosis in macrophages and monocytes, and production of nitric oxide and reactive oxygen species in macrophages were inactivated between C1 and C0 as well as C7 and C0 time intervals. Curiously, NK cell-related modules which were upregulated in non-protected vaccinees in GSEA (Figs 4A and S3), were inactivated in IPA at C1 vs C0 (NK cell signaling) (Fig 5A and S3 Data). At C7 relative to C0, pathways related to erythropoiesis and cell cycle were activated in protected vaccinees. Overall, both IPA and GSEA showed concordant results that indicate immune responses associated with interferon in protected vaccinees and innate responses associated with non-protected vaccinees and mock individuals.

thumbnail

Fig 5. Biological pathway and transcriptome cell subset response.

A) Ingenuity pathway analysis (IPA) on leading edge genes from significantly enriched gene sets. Leading edge genes were obtained from significantly enriched gene set modules in Fig 3A. Top 10 most significantly predicted pathways are shown in each group and time point. Bar width indicates the strength of IPA prediction following FDR P value adjustment. Colors indicate zScore prediction state; activated (purple), inactivated (green), and undetermined (grey) or z-score = 0 (white). B) ImmunoStates cell deconvolution analysis on comparisons between time points, individually tested in each sample group. Each square indicates a cell type. Color represents mean difference in cell proportion between time points. Squares outlined in purple represent significant statistical comparisons of FDR-adjusted P value <0.2. P (n = 6) indicates PfRAS-vaccinated protected individuals, NP (n = 5) indicates PfRAS-vaccinated non-protected individuals, Mock (n = 3) indicates non-infected mosquito bites vaccinated individuals.


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

To determine the cell types that contribute to the bulk transcriptomic responses, we performed cell-mixture deconvolution using a previously described reference expression matrix, immunoStates [33]. Consistent with GSEA, significant reductions in cell proportion were observed for the aggregate T and B cell response in protected vaccinees while significant reductions in proportion of granulocytes and neutrophils were also indicated in mock individuals after CHMI at C1 vs C0 (Fig 5B).

Distinct dynamics of immune cell subsets between immune status group following CHMI

To quantify changes in immune cell frequency after CHMI, we performed a comprehensive immunophenotyping by high-dimensional flow cytometry with 4 different antibody panels on PBMC samples collected at pre-vaccination baseline (B0) and 6 days (C6) and 112 days (C112) after CHMI. Unsupervised clustering using Seurat [34] identified a total of 50 clusters with distinct marker expression profiles (S5 Fig). Information on antibody marker expressions was used to annotate the clusters (S5 Fig and S4 Data). Mixed effect linear regression analysis was used to calculate changes in cluster frequency over time and to identify statistically significant differences between groups. We identified 3 clusters that significantly changed in frequency overtime and were significantly different between protected and non-protected vaccinees (Fig 6A). Cluster 6 of the invariant T cell panel expressed TCRγδ and vδ2, indicating a vδ2+ γδ T cell subset (Figs S5 and S4 Data). Cluster 3 of the DC panel had high expression of CD16, HLA-DR, CD11c, and CX3CR1 and low expression of CD14, suggesting that this is a non-classical monocyte subset. CD8 and CD56 surface markers were highly expressed on cluster 4 of the T cell panel while CD45RA and CCR7 marker were absent, and therefore, this cluster was designated as CD56+ CD8+ T effector memory cells (Tem). Consistent with other studies that reported an increase in γδ T cell frequency following CHMI [35,36], vδ2+ γδ T cell cluster frequency was higher at C6 compared to B0 in all groups (Fig 6B). Interestingly, the frequency of the vδ2+ γδ T cell cluster was increased at C112 relative to C6 in non-protected vaccinees and mock individuals, long after blood-stage parasitemia had resolved. Of note, we did not observe significant changes in the frequency of non- vδ2+ γδ T cell (Cluster 9 of the invariant T cell panel) following PfRAS vaccination and CHMI, as reported in other studies [37]. The decrease in non-classical monocyte cluster frequency in non-protected vaccinees and mock individuals at C6 was consistent with our GSEA results, which showed a decreased expression of gene sets associated with innate cells at C7 compared to C0 (Figs 6B and 4A and S3), albeit the sampling time points differed by a day between the two datasets. Similar to the vδ2+ γδ T cell cluster, the CD8+ Tem cluster expressing CD56 surface protein was also expanded over time in non-protected vaccinees and mock individuals after CHMI (Fig 6B), indicating the effect of exposures to blood-stage parasites.

thumbnail

Fig 6. Distinct cell subset responses following PfRAS and CHMI.

Cell subsets were identified using unsupervised clustering cells based on normalized marker fluorescence intensities, and annotated as described in S6 Fig and S4 Data. A) Bar plot showing significantly different cell subsets comparing protected vs. non-protected groups and mock vs. protected and non-protected groups. Significant clusters were identified using FDR adjusted P-values derived from linear mixed-effect model where FDR<0.05. B) Line plots showing significant cluster frequencies, calculated relative to pre-vaccination baseline (B0), for protected, non-protected, and mock individuals. Each solid black line represents an individual while dashed lines represent locally estimated scatterplot smoothing (LOESS) regression for the group, with 95% confidence interval indicated by the shaded areas. P (n = 6) indicates PfRAS-vaccinated protected individuals, NP (n = 5) indicates PfRAS-vaccinated non-protected individuals, Mock (n = 3) indicates non-infected mosquito bites vaccinated individuals.


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

To better understand the effect of PfRAS vaccinations on immune cell development and immune response to CHMI, comparisons were made between mock individuals and protected and non-protected vaccinees. A total of 4 and 7 clusters were significantly different between mock individuals and non-protected vaccinees and between mock individuals and protected vaccinees, respectively. These clusters represent CD56+ CD8+ Tems (cluster 4 –T cell panel), CD8+ T cells expressing CD57+ (cluster 4 –invariant t cell panel), atypical memory B cells (MBCs) (cluster 6 –B cell panel), naïve B cells expressing FCRL5+ (cluster 5 –B cell panel), non-classical monocytes (cluster 3 –DC panel), double negative (DN) CD3+ T cells (cluster 12 –invariant T cell panel), and IgG MBCs (cluster 3 –B cell panel) (Fig 6A). The increase in the frequency of the CD56+ CD8+ Tem cluster over the course of the trial was significantly higher in the mock group compared to PfRAS vaccinated groups (Fig 6B), possibly indicating a response to higher parasitemia in the mock group consistent with parasite qPCR data (Fig 1B). A similar increase in frequency in mock-vaccinated individuals was also observed for a CD8+ T cell cluster expressing CD57, suggesting a senescent and highly cytotoxic CD8 T cell. Curiously, the FCRL5 marker that was previously associated with impaired B cell responses in malaria [38], was highly expressed in a naïve B cell and atypical MBC cluster, and in a higher frequency in mock individuals compared to PfRAS vaccinated groups (Fig 6B). Clusters that were exclusively different between mock individuals and protected vaccinees, but not non-protected vaccinees included non-classical monocytes, DN CD3+ T cells, and IgG MBCs. There was a similar reduction in frequency in the non-classical monocyte cluster in mock individuals and non-protected vaccinees at C6 relative to B0 whereas DN CD3+ T cell and IgG MBCs had a larger increase in frequency following CHMI in mock individuals compared to protected vaccinees (Fig 6B). Of note, our observations in cell phenotyping data from flow cytometry were consistent with the scRNA-seq data, where we observed an overall increase in frequency across many cell types at convalescent C112 time point. Both flow cytometry and scRNA-seq also identified a cell cluster that was associated with CD56+ CD8+ T cell, i.e. cluster Cluster 4 (T panel–flow cytometry) and Cluster 34 (scRNA-seq) (Fig 6B). Overall, cell phenotyping data largely support transcriptomic observations of distinct immune responses among study individuals.

Discussion

In this study, we observed that systemic immune responses to P. falciparum infection are distinct between protected and non-protected PfRAS vaccinees as well as mock-vaccinated individuals. Blood transcriptome profiles in the PfRAS vaccinees diverged following CHMI in which non-protected vaccinees and mock individuals shared a more similar profile than protected vaccinees, before and after they developed blood-stage parasitemia. Transcriptomic responses in protected vaccinees were associated with early increase in NK cell and interferon (type I and II) responses and decreased expression of adaptive immunity, including T and B cells. In contrast, non-protected vaccinees and mock individuals had transcriptomic responses associated with decreased inflammation and innate cell-related signatures such as neutrophils and monocytes. Data from immunophenotyping revealed that PfRAS vaccinations and CHMI were associated with an increased frequency of vδ2+ γδ T cells in peripheral blood which remained elevated at convalescence, regardless of protection status. Non-classical monocytes were reduced in mock individuals and non-protected vaccinees but not protected vaccinees following CHMI, consistent with the whole blood transcriptomic response. Immune signatures related to malarial infection from CHMI in mock individuals were characterized by higher expression of the senescence and cytotoxic CD57 protein on CD8+ T cells and intriguingly, by higher expression of the dysfunctional marker FCRL5 on atypical MBCs and naïve B cells. Our study has revealed important differences in immunological response patterns among protected and non-protected PfRAS vaccinees and mock-vaccinated individuals following malarial infection providing additional insights into immune mechanisms required to achieve sterilizing immunity.

Interferon responses, including type I IFN and IFNγ, have been associated with the suppression of liver-stage infection in mice [39]. Conversely, we found that both type I and type II interferon responses were largely unique to protected PfRAS vaccinees, prior to and one day post CHMI. Such responses were largely absent in both non-protected vaccinees and mock individuals. Other clinical studies analyzing transcriptomic profiles of malaria-specific PBMCs have also shown increases of interferon responses in protected compared to non-protected individuals vaccinated with RTS,S or CPS vaccine [24]. Nonetheless, type I IFN responses have also been associated with reduced protective efficacy via impairment of CD8+ T cell responses in GAP-immunized mice [40]. It is likely that optimal interferon responses are required to induce sterilizing immunity, and that insufficient or excessive response could lead to a lack of protection. Our previous studies have demonstrated that the dynamics of transcriptome responses were crucial in determining protection outcome, in which strong and early increased interferon and inflammatory-associated responses following the first dose of PfRAS vaccination were associated with lack of protection [4143]. Although it is unclear which cell types are mediating interferon responses in protected vaccinees, our scRNA-seq data suggest that in mock-vaccinated individuals, T and NK cells were the primary cell subsets that responded to CHMI and induced interferons.

In addition to IFN responses, we observed early decreased expression in T and B cell-related BTMs in protected vaccinees one day after CHMI. Given that such observations are lacking in the non-protected and mock-vaccinated individuals, we speculate that these adaptive immune modules represent immune cells that responded to parasitic infection and trafficked into the liver to target sporozoites and activate effector functions. Different subsets of CD8+ T cells can contribute to sterilizing immunity in malaria, including the memory and liver resident subsets [17,44,45]. These CD8+ T cells are thought to mediate anti-parasitic effector mechanisms through IFNγ secretion [9,18,4648]. Although it is widely perceived that non-circulating liver resident CD8+ T cells are the main subset targeting liver stage malaria, recent work has demonstrated that circulating antigen-specific CD8+ Tem could rapidly infiltrate liver cells during infection and mediate parasite clearance [49]. It is possible that our observation of decreased expression of T cell-related BTMs in the blood in protected vaccinees were indicative of parasite-specific CD8+ Tems that may also secrete IFNγ. Along with T cells, B cells could also mediate sporozoite clearance by generating antibodies that function through multiple mechanisms, including complement-fixation [50] and interactions with Fcγ receptors on various immune cells [51]. Nonetheless, other immune cells are also likely to contribute to the development of sterilizing immunity to malaria, including the γδ T cells [20,21], CD4+ T cells [16], DCs [23,52], and NK cells [22]. In contrast to protected vaccinees, an increase in BTMs associated with the adaptive immune response was only observed 7 days post CHMI in non-protected vaccinees and mock individuals and correspond with blood-stage parasitemia, suggesting the expansion of adaptive immune cells in response to malarial infection and parasite replication.

Neutrophils and monocytes are innate immune cells that provide the first line of defense against infection. Early activation of these cells during malarial infection have been described in CHMI and cross-sectional malaria studies [5355]. Unexpectedly, while most previous studies in P. falciparum and P. vivax reported an increase in innate cell frequency and activated inflammatory pathways after malarial infection [5358], we observed decreased expression in innate-cell associated BTMs and reduced frequency of non-classical monocytes specifically in non-protected and mock groups after CHMI. BTMs associated with inflammation were also decreased in these CHMI-susceptible groups. Different sampling times post-CHMI or after natural infection and past malaria exposures could explain the discrepancies in the direction of innate cell frequency or enrichment score between studies. It is possible that the decreased innate cell response observed in our cohort indicates trafficking activity of these cells into the tissue or liver where sporozoite infection occurred and that these innate cells were also responsible for mediating inflammatory responses. Data from mice showed that neutrophils and monocytes can rapidly expand in the skin within hours after intradermal injection of irradiated sporozoites [59], demonstrating rapid recruitment of these cells to the sites of infection in the tissue. Other clinical studies have also reported transcriptomic activation in inflammatory responses post-CHMI in malaria naïve individuals, albeit with an opposite enrichment score direction to our study [58]. Thus, we hypothesize that the non-protected and mock groups in our cohort lacked sufficient parasite-specific adaptive immune cells that target and suppress parasite replication after CHMI, which subsequently resulted in the activation of innate cells in order to further capture and present parasite antigens to the adaptive immune cells. In line with this hypothesis, BTMs associated with adaptive immune cells were significantly enriched while innate cells were largely absent in the protected group after CHMI.

While we also observed increases in NK cell-related BTMs in non-protected vaccinees and mock individuals after CHMI, the same BTMs were also increased in expression earlier prior to CHMI in protected vaccinees (S3 Fig), suggesting that early NK cell response may be pivotal to protection outcomes against malaria. It is unknown if these NK cells were specific to Plasmodium parasites or if they indicate the evidence of adaptive NK cells via trained immunity. Studies in mice have demonstrated that hepatic NK cells can acquire memory-like and antigen-specific capacity to viral antigens [60]. In malaria, the protective capacity of the CD8 T cell response is dependent on the NK cell response [22]. Our results showed that several of the KIR genes including KIR2DL1, KIR2DL3, and KIR3DL1, were directing the overall NK cell BTM response in the protected vaccinees prior to CHMI. These KIR genes are classified as inhibitory KIRs, due to long intracytoplasmic tails [61]. Coincidentally, CD56dim NK cells express diverse inhibitory KIRs for HLA class I and can inhibit P. falciparum growth through antibody-dependent cellular cytotoxicity (ADCC) [62]. It is possible that CD56dim NK subset was mediating malaria immunity in protected vaccinees in our study, albeit it is unclear if this was induced by PfRAS vaccination or due to genetic variation of the highly polymorphic KIR genes between individuals [63]. An in-vitro study exposing cells from different malaria-naïve donors with parasitized RBCs showed that NK cell activation profile was heterogenous between individuals [64]. Thus, the early and elevated NK cell response in protected vaccinees in this study could also be explained by genetic variation between individuals.

Immunophenotyping by flow cytometry revealed an increase in vδ2+ γδ T cells in all sample groups following CHMI. This increase was sustained out to 112 days post-CHMI in non-protected vaccinees and mock individuals and was higher compared to protected vaccinees at convalescence. Similar longitudinal expansion in vδ2+ γδ T cells in mock individuals was also observed from the scRNA-seq data, although this was not statistically significant. Expansion of γδ T cells in subsets expressing vγ9+ vδ2+ chain, has been consistently reported in different malaria settings, including in vaccine trials [18,20], CHMI [37,65,66], and natural infection [37,67,68]. Late expansion of vδ2+ γδ T cells following acute parasitemia has also been previously reported [37,69,70], up to 35 days post-CHMI. While vδ2+ γδ T cells is often considered as innate-like, responding to phosphoantigens (Pags) [71,72], we observed long-lived and late expansion of vδ2+ γδ T cells 112 days after CHMI and after parasite clearance. A similar observation was also reported in macaques, with expansion of vδ2+ γδ T cells that persisted up to 7 months after a second BCG vaccination [73]. While data from mice and humans suggest that γδ T cells play a role in mediating malaria protection [20,21], our results suggest that vδ2+ γδ T cells expand along with parasite exposures, from either or both sporozoites and infected RBCs. This expansion is also likely mediated from a TCR-mediated clonal expansion, as previously demonstrated in mice [69]. Contrary to a previous study that saw an increase in vδ1+ γδ T cell frequency and its differentiation into effector subsets after CHMI [37], we did not observe any significant changes in the non-vδ2+ γδ T cell fraction. This could possibly be due to lack of vδ1 marker in our study or due to differences in study design, that is, PfRAS vaccination followed by CHMI versus repeated CHMI exposures. Specific to our study and PfRAS vaccination strategy, our results suggest that the vδ2+ γδ T cells did not play a primary role in mediating protection and that sterile protection was mainly mediated by the adaptive immune cells.

When comparing immunophenotyping profile between PfRAS vaccinees (both protected and non-protected) and mock individuals, we observed a higher frequency of cell subsets expressing the senescent marker CD57 on CD8 T cells and the dysfunction marker FCRL5 on atypical MBCs and naïve B cells in the mock group. CD57 expression on T cells has been implicated in reduced malaria immunity, in which the frequency of CD57+ T cells was higher in symptomatic malaria individuals compared to asymptomatic and healthy individuals [74]. Additionally, FCRL5+ atypical MBCs with reduced function have been associated with increased exposure to malaria [38]. Increased frequency in these cell subsets may indicate dysfunctional cell responses due to over-replication of parasitemia.

Our study has several limitations. First, we were unable to directly test for correlations between gene expression and cell type profiling data due to a lack of matching CHMI time points between blood and PBMC samples. Further, the relatively low sample size and missing sampling time points in the PBMC data are limitations that reduce the statistical power and sensitivity of the study. We utilized linear mixed-effect regression on the overall time points to min this limitation and to investigate the overall change in transcriptome cell type responses in the cohort. Using this approach, we were nonetheless able to identify genes and cell subsets that responded to PfRAS and CHMI and were associated with malaria protection. Inherent differences between study volunteers could also contribute to differences in immune response to PfRAS vaccination and CHMI as well as to protection outcomes in our study. Indeed, elevated inflammatory responses states prior to any immunization have been associated with malaria protection [75]. Given that volunteers in our study were given a suboptimal dose, it may be that higher vaccination dose is required to overcome these variations and induce a higher protection rate. Although our data have indicated that early transcriptome interferon responses post-CHMI are associated with protection, we could not demonstrate if the mRNA transcripts are translated into protein synthesis and secretion in the blood. Lastly, it is important to note that the immune responses described in our cohort do not reflect antigen-specific nor tissue-specific immune responses, rather, they reveal total responses in the peripheral blood. Further studies may employ techniques that can characterize antigen-specific responses such as in-vitro stimulation or antigen-coupled tetramers as well as comparing tissue-specific to circulating immune responses.

Materials and methods

RNA sequencing and data generation

Total RNA was extracted from whole blood using PAXgene Blood RNA Kit (PreAnalytiX), followed by treatment with GLOBINclear™ kit (ThermoFisher Scientific) to remove unwanted globin mRNA. The remaining RNA product was used to prepare cDNA library using Illumina TruSeq Stranded mRNA library preparation kit (Illumina). RNAseq was performed by Beijing Genomics Institute using the Illumina Hiseq2000 with 75 base-pair (bp) paired-end reads to a depth of at least 30 million reads per sample.

RNAseq output was processed as previously described [76]. Read pairs were pre-processed to adjust base calls with phred <5 to ‘N’ and to remove read pairs where either end had fewer than 30 unambiguous base calls. This method indirectly removes read pairs containing mostly adaptor sequences. Read pairs were aligned to the human genome (hg19) using STAR (v2.3.1d) [77]. Gene counts were tabulated using htseq (v0.6.0) with the intersection-strict setting turned on and Ensembl gene annotations (GRCh37.74) used to map genomic locations to gene identifiers [78]. The edgeR (v3.36.0) package function, cpm, was used to calculate TMM-normalized counts-per-million (CPM) expression matrices [79].

Single cell RNA sequencing and CITE-seq

PBMC samples from mock individuals (n = 4) in cohort 1 (subjects S049 and S054) and 2 (subjects S104 and S137) of the IMRAS trial were obtained at day 0 (B0), 6 days post-CHMI (C6), and 112 days post-CHMI (C112). Cells were thawed and counted and processed over three batches in combination with other samples entering the same processing pipeline (batch 1 = S049 all timepoints; batch 2 = S054 and S104 all timepoints; batch 3 = S137 all timepoints). 5 x 105 cells from each sample were stained for combinatorial hashing using oligonucleotide-conjugated Hashtag antibodies targeting the ubiquitous cell surface antigens Beta-2-Microglobulin and CD298 (Biolegend Totalseq A Hashtag antibodies). Hashed samples were pooled and stained together with a cocktail of 73 oligonucleotide-conjugated antibodies targeting cell surface antigens, including principal lineage antigens, that also includes isotype control antibodies (Biolegend Totalseq A). The stained pool of cells was washed, and non-viable cells removed using a magnetic Dead Cell Removal kit (StemCell). Each batch of pooled samples was then encapsulated across 5 lanes of a 10X Genomics v3.0 chip with a target cell capture of 3 X 104 cells per lane. An average of ~4600 cells (minimum 3844, maximum 6404) were encapsulated per IMRAS trial sample analyzed in this study.

Separate sequencing libraries were generated for gene expression (GEX), oligonucleotide-conjugated antibodies (ADT), and hashtags (HTO). Hashtag libraries were sequenced using the Illumina NextSeq 500/550 High Output (75 cycle) and NovaSeq S1 (100 cycle) kits at a target depth of 1,000 reads per cell, and gene expression and antibody libraries were sequenced using the Illumina NovaSeq S2 (100 cycle) kit at a target depth of 25,000 and 10,000 reads per cell, respectively. The binary base call (BCL) sequence files were base-called and demultiplexed using cellranger mkfastq (v3.1.0). Alignment to the human reference GRCh38 and feature counts were performed using cellranger count (v3.1.0). GEM wells for each sample were aggregated and depth normalized using cellranger aggr (v3.1.0).

Clustering of the CITE-seq data using weighted nearest neighbor (WNN) analysis

A Seurat (v4.0.1) [34,83] object was created for each sample following debris removal and barcode mapping to samples. Raw count data was pre-processed to remove low quality cells with high mitochondrial RNA fraction (25%), very high UMI counts (nCount_RNA >25,000), and very low or very high gene counts (nFeature_RNA <200 and >4000). The resulting per sample RNA and ADT data were independently normalized using standard or center log normalization respectively, then each integrated across all samples using the canonical correlation analysis (CCA) method to correct for any batch effects and control for between-sample technical variation. To allow for the integration of multimodal datasets, RNA and ADT weighted similarities were used to construct the weighted nearest neighbor (WNN) graph and to perform WNN-based uniform manifold approximation and projection (UMAP) plot as described [84]. Manual annotation of cellular identity was performed by finding differentially expressed genes and ADTs for each cluster using Seurat’s implementation of the Wilcoxon rank-sum test (“FindMarkers” function) and comparing those markers with known cell type–specific genes and surface marker proteins (S4 Fig and S4 Data).

Statistical and bioinformatics analysis

To calculate differentially expressed genes (DEGs) over time points in the longitudinal data, linear mixed-effect regression model (LMER) was fit through the R package glmmSeq (v.0.1.0) [85]. Gene expression (CPM) was fit as a dependent variable while sample time point was fit as a fixed effect, and individual identifier as a random effect. Dispersion was calculated using the edgeR (v3.36.0) package function, estimateDisp. For DEG calculation on comparisons between sample groups; protected, non-protected, and mock, at each time point, the R package DESeq2 (v1.34.0) was used [86].

Gene set enrichment analysis (GSEA) was performed using the R package fgsea (v1.20.0) [87]. Genes were ranked by the average fold change across different time interval comparisons and separately for each sample groups. Normalized enrichment score (NES) for non-statistically significant gene sets were adjusted to 0. The blood transcriptomic module (BTM) gene sets [30,31] were obtained from tmod R package (v0.46.2) [88] while the MsigDB Hallmark gene set [29] was obtained from download site http://www.gsea-msigdb.org/gsea/msigdb/genesets.jsp?collection=H.

Cell deconvolution analysis from bulk RNA transcriptomic data was performed using the immunoStates matrix [33] as part of MetaIntegrator R package (v2.1.3) [89]. For Ingenuity Pathway Analysis (IPA) [32], content version 68752261 was used to generate Canonical Pathway enrichment analysis results.

For flow cytometry data, LMER was used to calculate cell subsets that significantly change over time points and were significantly different between group comparisons. The R lme4 package (v1.1–27.1) [90] was used to fit the model using the lmer function. Cell cluster frequency generated by the unsupervised clustering was fit as a dependent variable while sample time points and group were fit as fixed effects, and individual identifier as a random effect. The mixed-effect models were fit as follow:

Full model : Cell frequency ~ 1 + Time point * Group + (1|Subject ID)

Reduced model 1 : Cell frequency ~ 1 + Time point + (1|Subject ID)

Reduced model 2 : Cell frequency ~ 1 + Group + (1|Subject ID)

Comparisons between the full model and reduced models were evaluated to determine significant relationships between cell frequency and time point and group variables. Specifically, ANOVA was used to compare the full model and reduced model 2, in which the P values represent statistical significance of improved fit associated with time variable, i.e. response-associated cell subsets. Within the statistically significant response-associated cell subsets, full model and reduced model 1 were compared using ANOVA to calculate significantly different cell subsets between groups.

For single-cell RNA sequencing data, significant changes in cell subsets overtime and GSEA were analyzed as above. Only “reduced model 1” was used to fit the mixed model, followed by ANOVA to test for statistical significance of the fixed effect, i.e. time point variable.

The FDR-corrected P values were calculated using the Benjamini-Hochberg method. For DEG, deconvoluted cell proportion, and flow cytometry cell subset frequency data, the FDR threshold was set at <0.2. For GSEA (both bulk and single-cell transcriptomics), IPA, and scRNA-seq cell subset frequency data, the FDR threshold was set at <0.01.

Supporting information

References

  1. 1.
    WHO. World malaria report 2021. Geneva; 2021.
  2. 2.
    RTS SCTP. 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. 2015/04/29 ed. 2015;386: 31–45. pmid:25913272
  3. 3.
    Datoo MS, Natama MH, Somé A, Traoré O, Rouamba T, Bellamy D, et al. Efficacy of a low-dose candidate malaria vaccine, R21 in adjuvant Matrix-M, with seasonal administration to children in Burkina Faso: a randomised controlled trial. Lancet. 2021/05/09 ed. 2021;397: 1809–1818. pmid:33964223
  4. 4.
    Birkett A, Miller RS, Soisson LA. The Importance of Exercising Caution When Comparing Results from Malaria Vaccines Administered on the EPI Schedule and on a Seasonal Schedule. Am J Trop Med Hyg. 2022; tpmd220544a. pmid:36410329
  5. 5.
    Roestenberg M, McCall M, Hopman J, Wiersma J, Luty AJ, van Gemert GJ, et al. Protection against a malaria challenge by sporozoite inoculation. N Engl J Med. 2009/07/31 ed. 2009;361: 468–77. pmid:19641203
  6. 6.
    Murphy SC, Vaughan AM, Kublin JG, Fishbauger M, Seilie AM, Cruz KP, et al. A genetically engineered Plasmodium falciparum parasite vaccine provides protection from controlled human malaria infection. Sci Transl Med. 2022;14: eabn9709. pmid:36001680
  7. 7.
    Mueller A-K, Labaied M, Kappe SHI, Matuschewski K. Genetically modified Plasmodium parasites as a protective experimental malaria vaccine. Nature. 2005;433: 164–167. pmid:15580261
  8. 8.
    Chattopadhyay R, Conteh S, Li M, James ER, Epstein JE, Hoffman SL. The Effects of radiation on the safety and protective efficacy of an attenuated Plasmodium yoelii sporozoite malaria vaccine. Vaccine. 2008/12/17 ed. 2009;27: 3675–80. pmid:19071177
  9. 9.
    Epstein JE, Tewari K, Lyke KE, Sim BK, Billingsley PF, Laurens MB, et al. Live attenuated malaria vaccine designed to protect through hepatic CD8+ T cell immunity. Science. 2011/09/10 ed. 2011;334: 475–80. pmid:21903775
  10. 10.
    Hoffman SL, Goh LM, Luke TC, Schneider I, Le TP, Doolan DL, et al. Protection of humans against malaria by immunization with radiation-attenuated Plasmodium falciparum sporozoites. J Infect Dis. 2002/04/04 ed. 2002;185: 1155–64. pmid:11930326
  11. 11.
    Roestenberg M, Teirlinck AC, McCall MB, Teelen K, Makamdop KN, Wiersma J, et al. Long-term protection against malaria after experimental sporozoite inoculation: an open-label follow-up study. Lancet. 2011/04/26 ed. 2011;377: 1770–6. pmid:21514658
  12. 12.
    Jongo SA, Shekalaghe SA, Church LWP, Ruben AJ, Schindler T, Zenklusen I, et al. Safety, Immunogenicity, and Protective Efficacy against Controlled Human Malaria Infection of Plasmodium falciparum Sporozoite Vaccine in Tanzanian Adults. Am J Trop Med Hyg. 2018/06/27 ed. 2018;99: 338–349. pmid:29943719
  13. 13.
    Oneko M, Steinhardt LC, Yego R, Wiegand RE, Swanson PA, Kc N, et al. Safety, immunogenicity and efficacy of PfSPZ Vaccine against malaria in infants in western Kenya: a double-blind, randomized, placebo-controlled phase 2 trial. Nat Med. 2021/09/15 ed. 2021;27: 1636–1645. pmid:34518679
  14. 14.
    Sissoko MS, Healy SA, Katile A, Omaswa F, Zaidi I, Gabriel EE, et al. Safety and efficacy of PfSPZ Vaccine against Plasmodium falciparum via direct venous inoculation in healthy malaria-exposed adults in Mali: a randomised, double-blind phase 1 trial. Lancet Infect Dis. 2017/02/22 ed. 2017;17: 498–509. pmid:28216244
  15. 15.
    Rénia L, Goh YS. Malaria Parasites: The Great Escape. Front Immunol. 2016/11/23 ed. 2016;7: 463. pmid:27872623
  16. 16.
    Carvalho LH, Sano G, Hafalla JC, Morrot A, Curotto de Lafaille MA, Zavala F. IL-4-secreting CD4+ T cells are crucial to the development of CD8+ T-cell responses against malaria liver stages. Nat Med. 2002/02/01 ed. 2002;8: 166–70. pmid:11821901
  17. 17.
    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. 2016/10/21 ed. 2016;45: 889–902. pmid:27692609
  18. 18.
    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. 2013/08/10 ed. 2013;341: 1359–65. pmid:23929949
  19. 19.
    Weiss WR, Jiang CG. Protective CD8+ T lymphocytes in primates immunized with malaria sporozoites. PLoS One. 2012/02/23 ed. 2012;7: e31247. pmid:22355349
  20. 20.
    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. 2016/05/10 ed. 2016;22: 614–23. pmid:27158907
  21. 21.
    Zaidi I, Diallo H, Conteh S, Robbins Y, Kolasny J, Orr-Gonzalez S, et al. γδ T Cells Are Required for the Induction of Sterile Immunity during Irradiated Sporozoite Vaccinations. J Immunol. 2017/10/29 ed. 2017;199: 3781–3788. pmid:29079696
  22. 22.
    Doolan DL, Hoffman SL. IL-12 and NK cells are required for antigen-specific adaptive immunity against malaria initiated by CD8+ T cells in the Plasmodium yoelii model. J Immunol. 1999/07/08 ed. 1999;163: 884–92. pmid:10395683
  23. 23.
    Montagna GN, Biswas A, Hildner K, Matuschewski K, Dunay IR. Batf3 deficiency proves the pivotal role of CD8α(+) dendritic cells in protection induced by vaccination with attenuated Plasmodium sporozoites. Parasite Immunol. 2015/08/19 ed. 2015;37: 533–543. pmid:26284735
  24. 24.
    Moncunill G, Scholzen A, Mpina M, Nhabomba A, Hounkpatin AB, Osaba L, et al. Antigen-stimulated PBMC transcriptional protective signatures for malaria immunization. Sci Transl Med. 2020/05/15 ed. 2020;12. pmid:32404508
  25. 25.
    Duffy PE, Patrick Gorres J. Malaria vaccines since 2000: progress, priorities, products. NPJ Vaccines. 2020/06/23 ed. 2020;5: 48. pmid:32566259
  26. 26.
    Itsara LS, Zhou Y, Do J, Grieser AM, Vaughan AM, Ghosh AK. The Development of Whole Sporozoite Vaccines for Plasmodium falciparum Malaria. Front Immunol. 2019/01/09 ed. 2018;9: 2748. pmid:30619241
  27. 27.
    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. PLoS One. 2020/06/20 ed. 2020;15: e0233840. pmid:32555601
  28. 28.
    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 A. 2005/10/04 ed. 2005;102: 15545–50. pmid:16199517
  29. 29.
    Liberzon A, Birger C, Thorvaldsdóttir H, Ghandi M, Mesirov JP, Tamayo P. The Molecular Signatures Database (MSigDB) hallmark gene set collection. Cell Syst. 2016/01/16 ed. 2015;1: 417–425. pmid:26771021
  30. 30.
    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. 2013/12/18 ed. 2014;15: 195–204. pmid:24336226
  31. 31.
    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/07/18 ed. 2008;29: 150–64. pmid:18631455
  32. 32.
    Krämer A, Green J, Pollard J, Tugendreich S. Causal analysis approaches in Ingenuity Pathway Analysis. Bioinforma Oxf Engl. 2014;30: 523–530. pmid:24336805
  33. 33.
    Vallania F, Tam A, Lofgren S, Schaffert S, Azad TD, Bongen E, et al. Leveraging heterogeneity across multiple datasets increases cell-mixture deconvolution accuracy and reduces biological and technical biases. Nat Commun. 2018;9: 4735. pmid:30413720
  34. 34.
    Stuart T, Butler A, Hoffman P, Hafemeister C, Papalexi E, Mauck WM, et al. Comprehensive Integration of Single-Cell Data. Cell. 2019;177: 1888–1902.e21. pmid:31178118
  35. 35.
    Lyke KE, Ishizuka AS, Berry AA, Chakravarty S, DeZure A, Enama ME, et al. Attenuated PfSPZ Vaccine induces strain-transcending T cells and durable protection against heterologous controlled human malaria infection. Proc Natl Acad Sci U S A. 2017;114: 2711–2716. pmid:28223498
  36. 36.
    Rutishauser T, Lepore M, Di Blasi D, Dangy J-P, Abdulla S, Jongo S, et al. Activation of TCR Vδ1+ and Vδ1-Vδ2- γδ T Cells upon Controlled Infection with Plasmodium falciparum in Tanzanian Volunteers. J Immunol Baltim Md 1950. 2020;204: 180–191. pmid:31801816
  37. 37.
    von Borstel A, Chevour P, Arsovski D, Krol JMM, Howson LJ, Berry AA, et al. Repeated Plasmodium falciparum infection in humans drives the clonal expansion of an adaptive γδ T cell repertoire. Sci Transl Med. 2021;13: eabe7430. pmid:34851691
  38. 38.
    Sullivan RT, Kim CC, Fontana MF, Feeney ME, Jagannathan P, Boyle MJ, et al. FCRL5 Delineates Functionally Impaired Memory B Cells Associated with Plasmodium falciparum Exposure. PLoS Pathog. 2015/05/21 ed. 2015;11: e1004894. pmid:25993340
  39. 39.
    Miller JL, Sack BK, Baldwin M, Vaughan AM, Kappe SHI. Interferon-mediated innate immune responses against malaria parasite liver stages. Cell Rep. 2014/04/08 ed. 2014;7: 436–447. pmid:24703850
  40. 40.
    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. 2019/09/04 ed. 2019;10: 3950. pmid:31477704
  41. 41.
    Du Y, Hertoghs N, Duffy FJ, Carnes J, McDermott SM, Neal ML, et al. Systems analysis of immune responses to attenuated P. falciparum malaria sporozoite vaccination reveals excessive inflammatory signatures correlating with impaired immunity. PLoS Pathog. 2022;18: e1010282. pmid:35108339
  42. 42.
    Duffy FJ, Du Y, Carnes J, Epstein JE, Hoffman SL, Abdulla S, et al. Early whole blood transcriptional responses to radiation-attenuated Plasmodium falciparum sporozoite vaccination in malaria naïve and malaria pre-exposed adult volunteers. Malar J. 2021;20: 308. pmid:34243763
  43. 43.
    Duffy FJ, Hertoghs N, Du Y, Neal ML, Oyong D, McDermott S, et al. Longitudinal immune profiling after radiation-attenuated sporozoite vaccination reveals coordinated immune processes correlated with malaria protection. Infectious Diseases (except HIV/AIDS); 2022 Sep. pmid:36591224
  44. 44.
    Schmidt NW, Podyminogin RL, Butler NS, Badovinac VP, Tucker BJ, Bahjat KS, et al. Memory CD8 T cell responses exceeding a large but definable threshold provide long-term immunity to malaria. Proc Natl Acad Sci U A. 2008/09/11 ed. 2008;105: 14017–22. pmid:18780790
  45. 45.
    Schmidt NW, Butler NS, Badovinac VP, Harty JT. Extreme CD8 T cell requirements for anti-malarial liver-stage immunity following immunization with radiation attenuated sporozoites. PLoS Pathog. 2010/07/27 ed. 2010;6: e1000998. pmid:20657824
  46. 46.
    Jobe O, Lumsden J, Mueller AK, Williams J, Silva-Rivera H, Kappe SH, et al. Genetically attenuated Plasmodium berghei liver stages induce sterile protracted protection that is mediated by major histocompatibility complex Class I-dependent interferon-gamma-producing CD8+ T cells. J Infect Dis. 2007/07/13 ed. 2007;196: 599–607. pmid:17624847
  47. 47.
    Nganou-Makamdop K, van Gemert GJ, Arens T, Hermsen CC, Sauerwein RW. Long term protection after immunization with P. berghei sporozoites correlates with sustained IFNγ responses of hepatic CD8+ memory T cells. PLoS One. 2012/05/09 ed. 2012;7: e36508. pmid:22563506
  48. 48.
    Seguin MC, Klotz FW, Schneider I, Weir JP, Goodbary M, Slayter M, et al. Induction of nitric oxide synthase protects against malaria in mice exposed to irradiated Plasmodium berghei infected mosquitoes: involvement of interferon gamma and CD8+ T cells. J Exp Med. 1994/07/01 ed. 1994;180: 353–8. pmid:7516412
  49. 49.
    Lefebvre MN, Surette FA, Anthony SM, Vijay R, Jensen IJ, Pewe LL, et al. Expeditious recruitment of circulating memory CD8 T cells to the liver facilitates control of malaria. Cell Rep. 2021;37: 109956. pmid:34731605
  50. 50.
    Kurtovic L, Atre T, Feng G, Wines BD, Chan J-A, Boyle MJ, et al. Multi-functional antibodies are induced by the RTS,S malaria vaccine and associated with protection in a phase I/IIa trial. J Infect Dis. 2020 [cited 16 Apr 2020]. pmid:32236404
  51. 51.
    Feng G, Wines BD, Kurtovic L, Chan JA, Boeuf P, Mollard V, et al. Mechanisms and targets of Fcγ-receptor mediated immunity to malaria sporozoites. Nat Commun. 2021/03/21 ed. 2021;12: 1742. pmid:33741975
  52. 52.
    Jobe O, Donofrio G, Sun G, Liepinsh D, Schwenk R, Krzych U. Immunization with radiation-attenuated Plasmodium berghei sporozoites induces liver cCD8alpha+DC that activate CD8+T cells against liver-stage malaria. PloS One. 2009;4: e5075. pmid:19347042
  53. 53.
    Berens-Riha N, Kroidl I, Schunk M, Alberer M, Beissner M, Pritsch M, et al. Evidence for significant influence of host immunity on changes in differential blood count during malaria. Malar J. 2014/04/25 ed. 2014;13: 155. pmid:24758172
  54. 54.
    Dobbs KR, Embury P, Vulule J, Odada PS, Rosa BA, Mitreva M, et al. Monocyte dysregulation and systemic inflammation during pediatric falciparum malaria. JCI Insight. 2017/09/22 ed. 2017;2. pmid:28931756
  55. 55.
    van Wolfswinkel ME, Langenberg MCC, Wammes LJ, Sauerwein RW, Koelewijn R, Hermsen CC, et al. Changes in total and differential leukocyte counts during the clinically silent liver phase in a controlled human malaria infection in malaria-naïve Dutch volunteers. Malar J. 2017/11/12 ed. 2017;16: 457. pmid:29126422
  56. 56.
    Antonelli LR, Leoratti FM, Costa PA, Rocha BC, Diniz SQ, Tada MS, et al. The CD14+CD16+ inflammatory monocyte subset displays increased mitochondrial activity and effector function during acute Plasmodium vivax malaria. PLoS Pathog. 2014/09/19 ed. 2014;10: e1004393. pmid:25233271
  57. 57.
    Rojas-Peña ML, Vallejo A, Herrera S, Gibson G, Arévalo-Herrera M. Transcription Profiling of Malaria-Naïve and Semi-immune Colombian Volunteers in a Plasmodium vivax Sporozoite Challenge. PLoS Negl Trop Dis. 2015/08/06 ed. 2015;9: e0003978. pmid:26244760
  58. 58.
    Tran TM, Jones MB, Ongoiba A, Bijker EM, Schats R, Venepally P, et al. Transcriptomic evidence for modulation of host inflammatory responses during febrile Plasmodium falciparum malaria. Sci Rep. 2016;6: 31291. pmid:27506615
  59. 59.
    Mac-Daniel L, Buckwalter MR, Berthet M, Virk Y, Yui K, Albert ML, et al. Local immune response to injection of Plasmodium sporozoites into the skin. J Immunol. 2014/07/02 ed. 2014;193: 1246–57. pmid:24981449
  60. 60.
    Paust S, Gill HS, Wang BZ, Flynn MP, Moseman EA, Senman B, et al. Critical role for the chemokine receptor CXCR6 in NK cell-mediated antigen-specific memory of haptens and viruses. Nat Immunol. 2010/10/26 ed. 2010;11: 1127–35. pmid:20972432
  61. 61.
    Tukwasibwe S, Nakimuli A, Traherne J, Chazara O, Jayaraman J, Trowsdale J, et al. Variations in killer-cell immunoglobulin-like receptor and human leukocyte antigen genes and immunity to malaria. Cell Mol Immunol. 2020;17: 799–806. pmid:32541835
  62. 62.
    Arora G, Hart GT, Manzella-Lapeira J, Doritchamou JY, Narum DL, Thomas LM, et al. NK cells inhibit Plasmodium falciparum growth in red blood cells via antibody-dependent cellular cytotoxicity. eLife. 2018;7: e36806. pmid:29943728
  63. 63.
    Middleton D, Gonzelez F. The extensive polymorphism of KIR genes. Immunology. 2010;129: 8–19. pmid:20028428
  64. 64.
    Korbel DS, Newman KC, Almeida CR, Davis DM, Riley EM. Heterogeneous human NK cell responses to Plasmodium falciparum-infected erythrocytes. J Immunol Baltim Md 1950. 2005;175: 7466–7473. pmid:16301654
  65. 65.
    de Jong SE, van Unen V, Manurung MD, Stam KA, Goeman JJ, Jochems SP, et al. Systems analysis and controlled malaria infection in Europeans and Africans elucidate naturally acquired immunity. Nat Immunol. 2021/04/24 ed. 2021;22: 654–665. pmid:33888898
  66. 66.
    Mordmüller B, Surat G, Lagler H, Chakravarty S, Ishizuka AS, Lalremruata A, et al. Sterile protection against human malaria by chemoattenuated PfSPZ vaccine. Nature. 2017/02/16 ed. 2017;542: 445–449. pmid:28199305
  67. 67.
    Ho M, Tongtawe P, Kriangkum J, Wimonwattrawatee T, Pattanapanyasat K, Bryant L, et al. Polyclonal expansion of peripheral gamma delta T cells in human Plasmodium falciparum malaria. Infect Immun. 1994/03/01 ed. 1994;62: 855–62. pmid:8112855
  68. 68.
    Taniguchi T, Md Mannoor K, Nonaka D, Toma H, Li C, Narita M, et al. A Unique Subset of γδ T Cells Expands and Produces IL-10 in Patients with Naturally Acquired Immunity against Falciparum Malaria. Front Microbiol. 2017/08/05 ed. 2017;8: 1288. pmid:28769886
  69. 69.
    Mamedov MR, Scholzen A, Nair RV, Cumnock K, Kenkel JA, Oliveira JHM, et al. A Macrophage Colony-Stimulating-Factor-Producing γδ T Cell Subset Prevents Malarial Parasitemic Recurrence. Immunity. 2018/02/11 ed. 2018;48: 350–363.e7. pmid:29426701
  70. 70.
    Teirlinck AC, McCall MBB, Roestenberg M, Scholzen A, Woestenenk R, de Mast Q, et al. Longevity and composition of cellular immune responses following experimental Plasmodium falciparum malaria infection in humans. PLoS Pathog. 2011;7: e1002389. pmid:22144890
  71. 71.
    Tanaka Y, Morita CT, Tanaka Y, Nieves E, Brenner MB, Bloom BR. Natural and synthetic non-peptide antigens recognized by human gamma delta T cells. Nature. 1995;375: 155–158. pmid:7753173
  72. 72.
    Hintz M, Reichenberg A, Altincicek B, Bahr U, Gschwind RM, Kollas AK, et al. Identification of (E)-4-hydroxy-3-methyl-but-2-enyl pyrophosphate as a major activator for human gammadelta T cells in Escherichia coli. FEBS Lett. 2001;509: 317–322. pmid:11741609
  73. 73.
    Shen Y, Zhou D, Qiu L, Lai X, Simon M, Shen L, et al. Adaptive immune response of Vgamma2Vdelta2+ T cells during mycobacterial infections. Science. 2002;295: 2255–2258. pmid:11910108
  74. 74.
    Frimpong A, Kusi KA, Adu-Gyasi D, Amponsah J, Ofori MF, Ndifon W. Phenotypic Evidence of T Cell Exhaustion and Senescence During Symptomatic Plasmodium falciparum Malaria. Front Immunol. 2019/07/19 ed. 2019;10: 1345. pmid:31316497
  75. 75.
    Neal ML, Duffy FJ, Du Y, Aitchison JD, Stuart KD. Preimmunization correlates of protection shared across malaria vaccine trials in adults. NPJ Vaccines. 2022;7: 5. pmid:35031601
  76. 76.
    Thompson EG, Du Y, Malherbe ST, Shankar S, Braun J, Valvo J, et al. Host blood RNA signatures predict the outcome of tuberculosis treatment. Tuberc Edinb Scotl. 2017;107: 48–58. pmid:29050771
  77. 77.
    Dobin A, Gingeras TR. Mapping RNA-seq Reads with STAR. Curr Protoc Bioinforma. 2015;51: 11.14.1–11.14.19. pmid:26334920
  78. 78.
    Anders S, Pyl PT, Huber W. HTSeq—a Python framework to work with high-throughput sequencing data. Bioinforma Oxf Engl. 2015;31: 166–169. pmid:25260700
  79. 79.
    McCarthy DJ, Chen Y, Smyth GK. Differential expression analysis of multifactor RNA-Seq experiments with respect to biological variation. Nucleic Acids Res. 2012;40: 4288–4297. pmid:22287627
  80. 80.
    Hertoghs N, Schwedhelm KV, 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 J Int Soc Anal Cytol. 2020;97: 1019–1023. pmid:32415811
  81. 81.
    Mair F, Prlic M. OMIP-044: 28-color immunophenotyping of the human dendritic cell compartment. Cytom Part J Int Soc Anal Cytol. 2018;93: 402–405. pmid:29356334
  82. 82.
    Hahne F, LeMeur N, Brinkman RR, Ellis B, Haaland P, Sarkar D, et al. flowCore: a Bioconductor package for high throughput flow cytometry. BMC Bioinformatics. 2009;10: 106. pmid:19358741
  83. 83.
    Butler A, Hoffman P, Smibert P, Papalexi E, Satija R. Integrating single-cell transcriptomic data across different conditions, technologies, and species. Nat Biotechnol. 2018;36: 411–420. pmid:29608179
  84. 84.
    Hao Y, Hao S, Andersen-Nissen E, Mauck WM, Zheng S, Butler A, et al. Integrated analysis of multimodal single-cell data. Cell. 2021;184: 3573–3587.e29. pmid:34062119
  85. 85.
    Lewis M, Goldmann K, Sciacca E, Cubut C, Surace A. glmmSeq: General Linear Mixed Models for Gene-level Differential Expression. R Package Version 001. 2021; https://github.com/KatrionaGoldmann/glmmSeq.
  86. 86.
    Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15: 550. pmid:25516281
  87. 87.
    Sergushichev AA. An algorithm for fast preranked gene set enrichment analysis using cumulative statistic calculation. bioRxiv. 2016; 060012.
  88. 88.
    Weiner J 3rd, Domaszewska T tmod: an R package for general and multivariate enrichment analysis. 2016.
  89. 89.
    Haynes WA, Vallania F, Liu C, Bongen E, Tomczak A, Andres-Terrè M, et al. EMPOWERING MULTI-COHORT GENE EXPRESSION ANALYSIS TO INCREASE REPRODUCIBILITY. Pac Symp Biocomput Pac Symp Biocomput. 2017;22: 144–153. pmid:27896970
  90. 90.
    Bates D, Mächler M, Bolker B, Walker S. Fitting Linear Mixed-Effects Models Using lme4. J Stat Softw. 2015;67: 1–48.

Source link