Longitudinal immune profiling after radiation-attenuated sporozoite vaccination reveals coordinated immune processes correlated with malaria protection

Ethics statement

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

IMRAS study design

The IMRAS trial (ClinicalTrials.gov NCT01994525) design has been described previously [13–15]. In this work, we analyzed IMRAS Cohort 1, which consisted of 11 individuals who received 5 RAS immunizations via approximately 200 mosquito bites each. Protective efficacy was determined by CHMI with bites from 5 NF54 P. falciparum infected mosquitos. Protected individuals were defined as those that remained free of detectable blood stage malaria infection via PCR and thick-blood smear post-CHMI. These vaccinees were accompanied by 3 mock immunized individuals (mosquito bites only) who served as infectivity controls, all of whom were not protected after CHMI.

Serum antibody profiling via P. falciparum antigen protein array

Serum samples were taken from the 11 IMRAS RAS vaccinees and 3 mock-immunized study participants at day 0 (pre-vaccination) and day 140 (3 weeks after the fifth immunization and immediately pre-CHMI). Serum was profiled on P. falciparum whole-proteome microarrays, (PfWPM: Antigen Discovery Inc, Irvine CA), consisting of 7,455 P. falciparum full or partial proteins, representing 4,805 unique genes (91% of the P. falciparum proteome), to quantify levels of P. falciparum antigen-specific serum antibodies. P. falciparum array antigens were produced using either a cell-free ORF-expression clone library followed by in vitro transcription/translation (IVTT). After profiling, per-spot microarray signal intensity (MFI) was normalized by taking median foreground intensity minus median background fluorescence intensity based on IVTT control spots subtracted, and log2 transformed. To identify vaccine-induced serum antibody responses, the normalized MFI response was calculated as the increase in normalized MFI per-participant from day 0 to day 140 for RAS-immunized individuals. Significantly induced antibody responses were identified as those increased at least 2-fold from day 0 to day 140 with an FDR-adjusted single tailed t-test significance threshold of 0.05.

Whole blood RNAseq

Whole blood RNAseq was performed as previously described [13] on 11 IMRAS Cohort 1 RAS vaccinees. Whole blood was collected directly into PAXgene blood RNA tubes (PreAnalytiX, Hombrechtikon, Switzerland) and stored at −20°C. RNA extraction and globin transcript depletion (GlobinClear, ThermoFisher Scientific, MA, USA) were performed prior to cDNA library preparation using the Illumina TruSeq Stranded mRNA sample preparation kit (Illumina, CA, USA). Globin transcript depletion, cDNA library preparation and RNA sequencing were performed by the Beijing Genomics Institute (Shenzhen, China). A total of sixty-six RNAseq samples were sequenced, with a target depth of 30 million reads per sample. Samples from four timepoints (day 0, 35, 119, 140) were sequenced on Illumina (San Diego, CA) Hiseq2000 sequencers using 75 base-pair (bp) paired-end reads. The remaining samples were sequenced on BGI500 sequencers using 100 bp paired-end reads.

RNAseq differential expression and gene-set enrichment analysis

Raw RNAseq fastq files were filtered for quality control by adjusting base calls with phred scores < 5 to ‘N’, and read pairs for which either end had fewer than 30 unambiguous base calls were removed. Read pairs were aligned to the human genome (hg19) using STAR (v2.3.1d) [22]. Gene count tables were generated using htseq (v. 0.6.0) [23] with the intersection-strict setting on and Ensembl gene annotations (GRCh37.74) used to link genomic locations to gene identifiers. Raw aligned RNAseq counts were normalized to account for different sequencing platforms (Hiseq2000 vs BGI500) using ComBat-seq [24], taking advantage of 4 samples resequenced on both sequencers to ensure an unbiased data transformation. Principal component analysis was performed on log-transformed counts-per-million gene expression data.

The R limma and voom packages [25,26] were used to identify differentially expressed genes augmented by the dream package [27] to allow the use of a mixed modelling approach. Timepoint and protection were set as fixed effects, including interaction terms, with a random intercept for each participant. For each gene, models were fit using the formula:

To identify genes significantly changed in expression relative to the day of most recent immunization, i.e. vaccine-induced genes, contrasts were constructed of the form ‘timepoint – most recent immunization timepoint’, e.g. ‘D003 – D000’ represents changes at day 3 relative to day 0 (day of first immunization). An FDR cutoff of 0.05 was used to identify significant genes at all timepoints. Protection-associated genes were identified by F-test, as implemented in the limma topTable function, to identify genes that where any combination of protection:timepoint interaction coefficients were significant.

Gene-set enrichment analysis (GSEA) was performed using the R fgsea package [28], using the BloodGen3 coherently expressed blood transcriptional modules (BTMs) [29]. Gene rankings used as input for GSEA were based on estimated log2 fold changes obtained from the limma-dream model. The limma topTable function was used to extract fold changes for each ‘timepoint – most recent immunization induced timepoint’ contrast, which were input to fgsea to calculate vaccine-induced BTM enrichment. Similarly each protection:timepoint interaction coefficient was used to input fold changes into fgsea to to identify enriched protection-associated BTMs.

BTM expression was calculated as the 25% trimmed mean of log-transformed counts-per-million gene expression of each BloodGen3 BTM gene, per sample. BTMs with reliable expression in less than half of their genes were excluded. Non-coherent BTMs, defined as BTMs where less than half of their genes were strongly correlated (Pearson’s r > 0.5) with the trimmed mean of BTM expression were also excluded.

PBMC flow cytometry profiling

Flow cytometry profiling was carried out using previously published antibody panels and gating strategies [13,30,31]. PBMCs were collected from IMRAS participants and frozen for later use. After thawing in RPMI supplemented with 10% FBS and benzonase nuclease (Millipore EMD 0.05 U/ml), the samples were incubated with LIVE/DEAD Fixable Blue Dead Cell Stain Kit and the Human BD Fc Block for 30 min at room temperature before being simultaneously stained with four phenotyping panels as shown in OMIP-044 [30], OMIP-064 [31] and Table S 1. The cells were then profiled using a BD FACSymphony flow cytometer. The data were analyzed, and cellular populations gated and quantified as previously using FlowJo Software (version 9.6.6). The abundance of each manually gated cell subset was calculated as a percentage, using the counts of each defined cell subset divided by the total single live cells from that sample.

Statistical analysis of flow-cytometry cell counts

Mixed models were used to identify vaccine-induced (i.e. significantly changed in abundance over time) and protection-associated cell type subsets using the R lme4 package [32]. A nested approach was taken where a ‘full’ model including protection and timepoint as fixed effects, as well as interaction terms, was compared to protection-only and timepoint-only fixed effect models. In all cases, per-participant intercepts were fit as random effects.

ANOVA was used to compare the FULL model to the TIMEONLY model (to determine the effect of protection on abundance) and the PROTONLY model (to determine the effect of sampling timepoint on abundance) and calculate associated p-values, and an FDR threshold of 0.05 was applied to determine significance.

The IMRAS study

The IMRAS study design has previously been described [13–15]. Briefly, study participants were immunized five times with radiation attenuated P. falciparum NF54 sporozoites (RAS) delivered by mosquito bite, with roughly 200 bites per-immunization. Immunizations were spaced 4-5 weeks apart with protection determined by controlled human malaria infection (CHMI) three weeks after the final immunization. In this study, we analyzed serum antibody levels, whole blood transcriptomics and PBMC cell phenotypic profiles from the 11 RAS-immunized (6 protected, 5 non-protected) participants in IMRAS Cohort 1. The cohort also included three mock-immunized individuals who received sporozoite-free mosquito bites followed by CHMI on the same schedule as RAS vaccinees. Serum, PBMCs and whole-blood PAXgene samples were obtained for up to 5 timepoints per-immunization per-participant for analysis (Table 1)

RAS immunization induced detectable serum antibody responses

Serum samples were taken from all 11 vaccinees and 3 mock-immunized individuals immediately prior to immunization (day 0) and challenge (day 140). Each individual serum was screened on P. falciparum antigen protein arrays [Antigen Discovery Inc.] to detect the presence of P. falciparum-specific antibodies. Strongly significant antibody responses to three P. falciparum proteins: circumsporozoite protein (CSP), Apical membrane antigen 1 (AMA1) and thrombospondin related adhesion protein (TRAP) were seen in RAS immunized individuals (Figure 1 A). All three of these proteins are known to be expressed on the surface of on the apical membrane of P.falciparum sporozoites, and are involved in hepatocyte invasion and binding [33–35]. Antibody responses did not significantly differ between protected and non-protected individuals (Figure 1 B).

• Download figure
• Open in new tab

Figure 1:

Protein array determined serum antigen responses. A. Volcano plot. X-axis shows increase in protein array normalized mean fluorescent intensity (MFI) from pre-immunization (day 0) to immediately pre-challenge (day 140) timepoints for RAS-immunized participants. Y-axis shows matching log-transformed false-discovery rate (FDR) adjusted single-tailed t-test p-values. Significant (FDR < 0.05) and highly increased (MFI increase > 1) antibody responses are highlighted with red circles and labelled with the corresponding protein. B. Box and dot-plots of pre-vaccine and pre-challenge MFIs for each study participant for labelled antibody responses from A. shown separately for protected (P) and non-protected (NP) RAS vaccine recipients. Mock immunized individuals are included for reference. T-test p-values comparing pre-challenge antibody binding MFIs between P and NP individuals are indicated.

Differential expression analysis revealed vaccine-induced and protection-associated genes

To identify vaccine-induced responses that correlate with protection, we carried out longitudinal whole blood transcriptomic profiling of IMRAS RAS vaccinees over the entire course of repeated immunizations. A mixed modelling approach was applied (see Methods) to identify genes that were responsive to vaccination and/or associated with protection that also accounted for gene expression levels primarily attributable to per-participant effects. Principal component analysis of transcriptome profiles showed that whole blood expression profiles clustered closely together by participant, rather than by longitudinal timepoint (Figure S1), underscoring the need to incorporate transcriptional per-participant effects as a factor in our differential expression analysis.

Vaccine-induced differentially expressed genes (DEGs) were identified as those significantly changed in expression in blood relative to the most recent previous immunization. Hundreds of vaccine-induced DEGs were identified (FDR < 0.05) at the majority of timepoints across the entire course of immunizations (Figure 2 A, Datafile S 1). Most strikingly, we observed large numbers of genes decreased in expression after the second immunization. In contrast, there were relatively few vaccine induced DEGs seen after the third immunization. We quantified the overlap of these DEGs with the recently published BloodGen3 [29] coherent blood transcriptional modules (BTM) to interpret biological processes induced post-immunization. BTMs significantly enriched in vaccine-induced genes that were uniformly increased or decreased per-timepoint are shown in Figure 2 B (FDR-adjusted hypergeometric p-value < 0.05). Genes increased in expression in the two weeks after the first immunization were enriched for interferon, oxidative phosphorylation and erythroid cell BTMs. Genes reduced in expression after the second immunization were enriched for protein synthesis, protein modification, gene transcription and oxidative phosphorylation BTMs at multiple timepoints. Fewer significant BTMs were identified in DEGs after the 3^rd – 5^th immunizations, although BTMs representing cell death, and inflammation were increased 3 days after the 3^rd immunization. Numbers and expression patterns of DEGs and accompanying enriched BTMs were highly distinct after each immunization.

• Download figure
• Open in new tab

Figure 2:

Differentially-expressed genes associated with vaccine-induced and/or protection-associated responses to immunization. A. Bar plot showing counts of genes significantly increased (positive) or decreased (negative bar) relative to the day of the most recent vaccination. B. Bar plot showing counts of BloodGen3 modules enriched in genes shown in A., per-timepoint and direction (positive bars represent modules enriched in genes with increased expression relative to the day of the most recent immunization, negative bars represent enrichment for genes with decreased expression). Module functional annotation is shown by the color key. C. Line plots showing per-individual normalized expression of each protection-associated gene over time colored by protection status. Group-averaged expression is shown as a ribbon. For all panels, vertical black lines indicate times of immunizations.

We also found 21 genes significantly associated with protection when considering all timepoints together (Figure 2 C). The 21 protection-associated genes significantly overlapped (FDR-adjusted hypergeometric p = 3 × 10⁻⁸) with a single BTM: six protection-associated genes (EPSTI1, IFI44L, IFI6, ISG15, LY6E, RSAD2) are members of BTM M8.3 representing Type I Interferon responsive genes. These genes were also increased in expression immediately post-immunization, but whereas non-protected individuals had high and sustained increased expression of these genes after the first immunization, protected individuals had stronger increases after the remaining 4 immunizations (Figure 2 C). This suggests that while interferon responses are broadly associated with protection during RAS immunization, the effect of increased interferon gene expression on protection may change over the course of vaccination. 13 of these 21 total protection-associated genes and 5 of 6 type I interferon protection-associated genes were also significantly vaccine-induced to for at least one timepoint.

Gene-set enrichment analysis revealed complex timing of vaccine responses

Modules identified as having enriched overlap with DEGs were limited by the binary significance cutoff used to identify DEGs. Therefore, we applied gene-set enrichment analysis (GSEA) to identify BTMs that were significantly enriched in the full rank-ordered gene list. Genes were ranked based on the above determined linear model coefficients and normalized enrichment scores (NES) calculated to classify BTMs as being increased or decreased relative to the day of the most recent immunization (Figure 3 A) or associated with protection at a specific timepoint (Figure 3 B). Overall, vaccine-induced responses identified as enriched in DEGs (Figure 2 B) were also captured in the GSEA analysis; for example, increased interferon and erythroid cell BTMs after the first immunization, and decreased gene transcription and protein synthesis BTMs after the second immunization. However, the vaccine-induced GSEA captured many events after the third, fourth and fifth immunizations not seen in the DEG-based analysis. This includes increases in erythroid cell and gene transcription and decreased inflammation after the third immunization as well as increased inflammation, interferon and protein phosphorylation after the fourth and fifth immunization. Many significantly enriched BTMs fall into the ‘TBD’ category, where the overall functional annotation was not determined. These show enrichment patterns similar to inflammation and cytokine/chemokine BTMs, suggesting co-ordination with defined immune responses. Concordant with the individual DEG analysis, we saw highly distinct response patterns induced after each immunization.

• Download figure
• Open in new tab

Figure 3:

Gene set enrichment analyses showing patterns of vaccine-induced and protection-associated responses profiled using blood transcription modules (BTMs). Heatmaps show significant (FDR < 0.05) normalized enrichment scores (NES). A. BTMs significantly enriched for genes changing in expression relative to day of the most recent immunization. B. BTMs significantly enriched for genes with expression differences between protected and non-protected individuals at each timepoint. Vertical black lines indicate times of immunizations in both panels

Focusing on protection-associated BTMs, we saw more consistent patterns after each immunization. The majority of protection-associated processes comprised BTMs annotated as associated with specific immune cell types (adaptive and innate immune cells), cytokines, erythroid cell and inflammation and proliferation (protein synthesis/phosphorylation/modification) that were consistently reduced in protected individuals relative to non-protected individuals 1-3 days after each immunization (Figure 3 B). As with the vaccine-responsive BTMs above, ‘TBD’ modules mimicked the protection-associated enrichment patterns of inflammation associated BTMs. Interferon and prostanoid BTMs which were increased in protected participants immediately after immunizations 2-5, but not the first immunization. These interferon responses are consistent with what was observed in the DEG analysis above.

Hierarchical clustering revealed distinct patterns of protection-associated responses

In order to decipher the different patterns of protection-associated responses, we calculated average BTM expression and hierarchically clustered the modules based on their correlation matrix. BTMs that shared a functional annotation tended to cluster closely (Figure 4 A) with Interferon BTMs and erythroid cell BTMs forming highly distinct clusters. We also observed a distinct cluster that represented inflammation/neutrophil/cytokine BTMs. This indicated that BTMs representing similar functional annotations shared coherent expression patterns. Representative BTMs were then selected from these sub-clusters and their averaged expression was plotted (Figure 4 B). Alongside the previously noted type I Interferon BTM, we observed a variety of protection-associated expression patterns, including both transient responses to vaccination and sustained protection-associated differences over the study. Notably, M12.2 monocyte BTM expression showed strong increases in protected individuals shortly after immunizations 2-5. In contrast, the M11.3 erythroid cells and M8.2 prostanoid BTMs showed sustained increased expression in non-protected individuals throughout the study. Thus, GSEA revealed numerous functional gene sets and expression dynamics that were correlated with protection.

• Download figure
• Open in new tab

Figure 4:

Expression dynamics of protection-associated BTMs. Expression of protection-associated BTMs identified by GSEA was calculated as the average of BTM gene expression and used to create a BTM correlation matrix. A. The BTM correlation matrix, with BTM-BTM correlation coefficients (Pearson’s r) indicated, and BTMs ordered identically on the x- and y-axes by hierarchical clustering of the correlation matrix. B. Averaged expression of representative BTMs from subclusters in A. plotted with lines representing individuals BTM expression, colored by protection status. Protection group average expression shown as thick ribbons

Monocyte, Dendritic cell and invariant T cell populations were responsive to RAS vaccination

To better understand the cellular RAS vaccination response profile, we performed high-parameter flow cytometry profiling of PBMCs taken from the IMRAS vaccinees at timepoints overlapping with whole blood gene expression (Table 1). Four protein marker panels were used, consisting of 21-28 markers each, specific for B cells, T cells, antigen-presenting cells, and NK/ILC/ γδT T cells (Table S 1). An extensive gating strategy was applied (see Methods) to determine counts for hundreds of different cell types and corresponding activation states (Datafile S 2). As with whole-blood transcriptomics, cell abundance profile variability was most strongly predicted by person-to-person variability (Figure S 2) rather than by vaccine-induced changes, which necessitated a mixed-modelling approach to account for this variability while identifying vaccine-induced cell subsets.

Mixed-modelling revealed vaccine-induced cell types (Figure 5 A, Datafile S 3) including Vδ2 γδ T (Figure 5 B) cell subsets, classical CD14+ monocyte subsets (Figure 5 C), non-classical (NC) monocyte subsets (CD11c+ HLA-DR+ CD16+) (Figure 5 D) MAIT cell subsets, (Figure 5 E), and CD4+ and CD8+ T cell subsets. We observed few or no changes in B cell, NK cell or NKT cell subsets. Several cell type groupings showed coherent changes in abundance over the course of vaccination: NC monocyte subsets were rapidly and transiently increased after immunizations 1 and 3 and 5 in both protected and non-protected individuals. Larger increases were apparent in non-protected individuals after the first immunization with smaller increases after the third immunization, while the opposite pattern was observed in protected individuals (Figure 5 B). In contrast, Vδ2 γδ T cell subsets continuously increased over the course of vaccination and this increase was more apparent in non-protected individuals (Figure 5 C). Similar to NC monocytes, classical CD14+ monocytes subsets showed transient increases after each immunization (Figure 5 D). CD272- (a.k.a. BTLA: B-and T-lymphocyte attenuator) and CD85k+ (a.k.a. LILRB4, leukocyte immunoglobulin-like receptor B4) CD14+ monocytes represented the most abundant induced monocyte populations. For both CD272- and CD85k+ monocytes, increased abundance was apparent at the earliest post immunization 1 timepoint (day 3) in non-protected individuals, but not until day 7 post-immunization for protected individuals. A distinct pattern was observed for MAIT cell abundance. Total MAIT cell abundance was generally increased after each immunization for both protected and non-protected individuals. However, increases in total and CD8+ MAIT cells were stronger and more rapid in protected individuals after the first immunization. (Figure 5 E). Vaccine-induced T cell subsets, both CD4+ and CD8+, were not confined to a single memory or effector phenotype. However, since malaria sporozoite-specific T cells represent a small (or zero) proportion of total T cells it is unclear what is driving changes in circulating T cell proportion. These changes may simply reflect changes in the total number of circulating myeloid cells in response to vaccination.

• Download figure
• Open in new tab

Figure 5:

Cell type responses to RAS vaccination. A. Counts of significantly (FDR < 0.05) vaccine-induced cell types grouped by lineage and/or function. Red columns indicate total numbers of gated subsets in the corresponding group, adjacent teal columns show proportion of gated subsets that were vaccine-induced. B-E Averaged counts of vaccine-induced cell types showing vaccine induced cell types from specific lineage/function groupings in A., over time. Vertical black lines indicate vaccination times. Legend shows color for each gated cell type subset, ordered from most abundant (top) to least abundant (bottom). B. Vδ2 γδT cells, C. Classical CD14+ monocytes, D. Non-classical (NC) monocytes, E. MAIT cells. F. Averaged counts of the protection-associated (FDR < 0.05) cell type: CD1c+ CD40+ DCs.

It is important to note that none of the protection-associated trends listed in the vaccine-induced cell types described above reached statistical significance. We identified a single significantly protection-associated cell type: CD40+ CD1c+ DCs (Figure 5 F). These showed a similar overall expression pattern as the NC monocytes, but with a more pronounced increase in abundance in non-protected individuals after the first immunization and a more pronounced spike in protected individuals after the third immunization. In summary, we observed that myeloid cell types including classical CD14+ monocytes, non-classical monocytes and DCs tended to show similar vaccine response and protection dynamics as type I interferon-associated gene expression, while other vaccine responsive lymphoid cell types such as MAITs and Vδ2 γδT cells had distinct response patterns over time.

Combined analysis of cell populations and gene expression revealed cell types correlated with vaccine-induced and protection-associated transcriptional responses

To better understand the relationships between cell type abundances and transcriptional responses after RAS immunization, we created a comprehensive correlation map relating RAS induced BTM expression to cell type abundances. Correlations between BTM expression and vaccine-induced cell type abundances were calculated by matching transcriptomic profiles with flow-cytometry derived cell type abundanced from each participant at shared timepoints (Table 1). Multiple significant cell type-BTM correlated pairs were observed (Figure 6 A), and these correlations often reflected known functions of the cell type in question. For example, BTMs with the ‘inflammation’ functional annotation were correlated with CD14+ monocyte subsets and BTMs with the ‘cytotoxic lymphocyte’ annotation were correlated with γδ T cells. Many BTM functionally annotated as ‘antigen presentation’, ‘lymphocytes’; ‘T cells’, and ‘cell cycle’ were specifically correlated with naïve CD4+ and CD8+ T cell subset abundance (Figure 6 B), possibly reflecting initiation of adaptive responses post immunization.

• Download figure
• Open in new tab

Figure 6:

Correlations between vaccine-induced cell type abundances and vaccine-induced BTM expression. A. Heatmap showing Spearman correlations between vaccine-induced cell type abundance and average BTM expression. Only BTMs with at least one significant cell type correlation are shown. Strong positive and negative (|π| > 0.5 & FDR < 0.05) correlations are highlighted with red and blue boxes respectively. B. Subset of heatmap A. with T cell subsets labelled. C. Subset of heatmap A. with DC and Monocyte subsets labelled.

Several monocyte subsets were significantly (FDR < 0.05) and strongly positively correlated (Spearman’s π > 0.5) with Interferon BTM expression, specifically the CD16hi, CD272- and CD38hiCD163+ subsets of CD14+ monocytes, and CD38+CD32+ NC monocytes (Figure 6 C). Vaccine-induced NC monocyte subsets were uniformly positively correlated (all NC monocyte – interferon π > 0.2) with all 4 interferon BTMs. While most CD14+ monocyte populations were positively correlated with interferon BTMs, this was not always the case for CD38lo CD32+ or CD16hi populations of CD14+ monocytes. Therefore, it appeared that NC monocyte abundance was linked to interferon expression in the blood across the entire study period regardless of NC monocyte activation/maturity status, while this was not the case for CD14+ monocyte populations.

To further explore the temporal linkage between interferon transcriptional responses and DC/monocyte abundance, we investigated how the strength of their correlation varied over time. After immunization 1 (day 0-28) specifically, 2 CD14+ monocyte subsets were significantly positively correlated with interferon BTM expression. In contrast only 1 NC monocyte subset and no CD1c+ DCs (Table 2) were correlated with interferon BTMs. This changed after immunization 3, when 9 NC monocyte subsets, including total NC monocyte counts, were significantly correlated with interferon along with 3 CD14+ monocyte subsets and one CD1c+ DC subset (the protection associated CD1c+ CD40+ DC subset). This suggests that non-protection associated interferon responses after immunization 1 are primarily associated with subsets of CD14+ monocytes, while protection-associated interferon responses after immunization 3 were correlated with increased abundance of a wider spectrum of both CD14+ and NC monocyte populations.