Systems Analysis of Immunity to Influenza Vaccination across Multiple Years and in Diverse Populations Reveals Shared Molecular Signatures

We vaccinated with TIV a total of 212 subjects across 5 influenza vaccine seasons from 2007 to 2011, among whom 54 were elderly (>65 years old in 2010 and 2011 cohorts) (Figures S1A and S1B). The 2011 cohort also included 17 subjects diagnosed with type 2 diabetes (T2D). Overall, there was a 65:35 ratio of females to males, and the majority of subjects were of European descent (Figure S1B). Blood samples were collected at baseline (day of vaccination) and at the time points indicated (Figure 1A). We performed microarray analyses on the peripheral blood mononuclear cells (PBMCs) as indicated (Figure 1A). For the 2008 and 2009 seasons, we included in our analyses published data from a previous study by Franco et al. (2013)) containing 218 additional subjects. For a subset of the 2010 season cohort, fluorescence-activated cell sorting (FACS) measurements as well as miRNA profiling were also performed (Figure 1A). The entire dataset is publicly available online at Immport (Bhattacharya et al., 2014) (http://immport.niaid.nih.gov/).

We evaluated the plasma antibody responses of all vaccinated individuals via hemagglutination-inhibition (HAI) antibody titer assays. The individual antibody response for each of the three influenza strains included in the vaccine was calculated as the fold change between the HAI titer at day 28 relative to the baseline titer. We then defined the magnitude of the HAI response as the maximum fold change among the three influenza strains (Figure 1B; Nakaya et al., 2011). Subjects were classified as “high responders” if (1) their HAI response had at least a 4-fold increase at 28 days (Sullivan et al., 2010) and (2) the day 28 antibody titer was 1:40 or more for at least one strain; subjects were classified as “low responders” otherwise. In the 2010 and 2011 seasons, in which elderly subjects were included, antibody responses showed a significant decrease with age (Figures 1B, 1C, and S1C), consistent with previous studies (Seidman et al., 2012). Vaccines from different seasons induced varying levels of HAI responses, and there was a significant decrease in the response with age (Figure 1C). Subjects vaccinated during the 2010 season, the year after the H1N1 2009 pandemic, generated the strongest HAI responses. There was no significant difference in antibody response based on gender (Figure S1D), diabetic status (Figure S1E), or race (data not shown). There was no significant difference in the response between males and females in any of the seasons examined (data not shown).

In order to identify transcriptional pathways associated with the antibody response to vaccination, we performed gene set enrichment analysis (GSEA) (Subramanian et al., 2005) on genes correlated with the HAI response in each influenza season (Figure 2A). For this, we used a set of blood transcriptional modules (BTMs) previously identified by our group through large-scale network integration of publicly available human blood transcriptomes (Li et al., 2014). BTMs related to the induction of interferons as well as the activation of dendritic cells (Figures 2A and 2B) were enriched on days 1 and 3 after TIV vaccination, whereas modules related to T cells at these time points were negatively associated with the antibody response. These results were validated by analysis of both seasons of the Franco et al. (2013)) dataset. On day 7, there was a robust enrichment of ASC and cell cycle-related modules, consistent with our original study (Nakaya et al., 2011).

Additionally, we compared transcriptional responses based on gender, age, race, or diabetic status. We first performed differential expression analysis between these groups (male versus female, type-2-diabetes-positive versus -negative, and European versus African and/or Asian) using fold change expression data from day 3 and day 7. There were very few differentially expressed genes based on gender (Figure S2A) or race (data not shown), indicating similar responses after vaccination. Comparison of gene expression between type-2-diabetes-positive and -negative subjects revealed a modest number of differentially expressed genes (Figure S2B), but GSEA on genes ranked by correlation with the day 28 HAI response in these two groups showed a high degree of overlap in the modules associated with a high antibody response (Figure S2C). Based on the similarity of these transcriptional responses, we decided not to segregate subjects based on gender, race, or diabetic status in subsequent analyses.

Given that we could identify pathways that consistently correlated with the magnitude of HAI response (Figure 2A), we sought to identify transcriptional signatures of immunogenicity across multiple influenza seasons. In our previous study, we have used DAMIP (discriminant analysis via mixed integer programming) (Brooks and Lee, 2010) to find sets of 3–5 discriminatory genes that could predict the antibody response of subjects from the 2007–2009 seasons (Nakaya et al., 2011). Here we applied DAMIP to a larger number of subjects across a greater number of influenza vaccine seasons (Figure S3), as described in the Experimental Procedures and Supplemental Experimental Procedures. Here DAMIP was able to correctly predict the HAI response in the blind testing group with 67.6%–70.0% accuracy (Tables S1 and S2). Consistent with our previous work, the most frequently selected genes included a large number of ASC genes (Table S1).

Given the robustness of the module-level responses across different TIV seasons, we wanted to investigate whether module-level features were capable of accurately predicting vaccine-induced immunity in a multi-year dataset and in the elderly. Therefore, we first generated BTM normalized enrichment scores for each subject using single-sample GSEA (ssGSEA) (Barbie et al., 2009) and then used these scores as inputs to an artificial neural network classifier (Dreiseitl and Ohno-Machado, 2002) (see Supplemental Experimental Procedures for details). In brief, the young subjects from all seasons were randomly divided into training and testing groups in an 85%/15% ratio. For each bootstrapping cycle, we split the young subjects (from all years) as 85% (young training set) and 15% (young testing set) and put all elderly aside (from both 2010 and 2011 years, the elderly testing set). Then, we used the young training set to select the features (in this case the BTMs) and to check their performance in the same young training set. Next, we checked the performance of these “predictive BTMs” in the young testing set and elderly testing set. This process was repeated for 100 trials in order to ensure that the performance was robust across many random divisions of the data. For each of the 100 times that this 85%/15% split in young subjects was done, we tested the prediction on the same elderly testing set. With this approach, we identified BTMs that predicted HAI responses in both the young and elderly groups with accuracies ranging between 79.0% and 80.2% and between 64.7% and 72.3%, respectively, using both day 3 and day 7 signatures (Figure 2C). Examination of the modules most frequently selected by the algorithm in the 100 randomized trials revealed pathways related to innate immune cell responses as well as leukocyte differentiation and antigen presentation on day 3 and B cell/immunoglobulin production and cell cycle pathways on day 7 (Table S3).

Interestingly, we observed similar classification accuracy in both the young and elderly test sets when using the predictive signatures trained on the young dataset (Figure 2C). This suggests that although the elderly had lower antibody responses to vaccination, the molecular signatures that were predictive of high HAI response remained the same. One potential caveat to this was that in the analysis in Figure 2C we had defined young as <65 years and elderly as >65 years, a narrow chronological distinction that could hamper quantitative analysis of age-related effects. To address this issue, we performed additional analysis by removing the intermediate-age subjects (between 40 and 65 years) from our analysis (Figure 2D). Thus here we defined young as subjects younger than 40, and elderly as subjects older than 65. The results are similar to our original work (Figure 2C versus Figure 2D). In addition, as the 2010 cohort has a sharper distinction between the young and the elderly (Figure S1A), we performed GSEA on genes ranked by correlation with HAI response (similar to the analysis ofFigure 2A) for young and elderly subjects in the 2010 cohort (Figure S4). There were a number of consistently enriched modules in both the young and elderly, in particular those relating to the interferon response and activation of dendritic cells on day 1 and the plasmablast response on day 7. These results suggest that there were consistent transcriptional signatures associated with the day 28 antibody response in both groups. Together these results demonstrate that module-level features are able to successfully predict vaccine immunogenicity in TIV seasons in both young and elderly vaccinees.

We also examined whether there were pre-vaccination transcriptional signatures that were associated with antibody responses. Although there has been previous work aimed at identifying baseline cellular and molecular predictors of influenza vaccine response (Frasca et al., 2010, Furman et al., 2013, Tsang et al., 2014), these studies examined subjects within only a single influenza season, and the robustness of these signatures in predicting immunity to influenza vaccination across multiple seasons has not been determined. First, we performed GSEA (Subramanian et al., 2005) in each season by using genes ranked by correlation between baseline expression and the HAI response. However, this approach resulted in little overlap between seasons, with no module consistently enriched in more than three out of five seasons (Figure S5). This might have been due to a reduced signal-to-noise ratio and increased batch effects owing to the inability to use fold-change expression at baseline. Although we attempted to adjust for batch effects in the day 0 expression data across seasons by using the ComBat software (Johnson et al., 2007), these effects might have not been completely removed.

To increase the power of this analysis so as to detect the more subtle transcriptional signatures at baseline, we repeated the procedure with all five seasons combined. In addition, publicly available data from previously published influenza vaccine studies by Franco et al. (2013)) and Furman et al. (2013)) were included as independent validation sets (see Supplemental Experimental Procedures for details). This approach revealed several B-cell- and T-cell-related modules whose expression pre-vaccination was positively correlated with an increased antibody response to vaccination in all studies (Figure 3A). In contrast, modules related to monocytes were negatively correlated with antibody responses (Figures 3B–3G). Interestingly, contained within the monocyte-related modules were genes such as TLR4, TLR8, NOD2, and ASGR2, which encode proteins involved in innate sensing, and genes encoding IFN-γR, IL-13R, TIMP2, LYN, SYK, and other molecules involved in inflammation. This supports the concept that inflammatory responses at baseline might be detrimental to the induction of vaccine-induced antibody responses (Haq and McElhaney, 2014,Pawelec et al., 2014). Investigation of the correlation with antibody response at a gene level (Figures 3B–3G) highlights the strength of pathway-level analyses: although the module enrichment was consistent across datasets, the individual genes contributing to this enrichment varied from study to study.

We next investigated the influence of age on the vaccine-induced transcriptional signatures. To this end, we first identified differentially expressed genes after vaccination in the young versus elderly subjects from the 2010 season (Figure 4A) (in which there is a sharp chronological separation between the young and the elderly [Figure S1A]). By far the largest difference in overall expression between the two groups occurred on day 1, where the young exhibited a much greater number of both up- and downregulated genes. Weighted correlation network analysis (WGCNA) (Langfelder and Horvath, 2008) was used to find clusters of highly correlated genes among the TIV-regulated genes in young (clusters Y1–Y6) and elderly (clusters E1–E5) (Figure 4B). This method compares the correlation in expression patterns among genes across all young or elderly subjects and utilizes hierarchical clustering with dynamic tree cut to define modules of genes that are temporally co-expressed within each group. Genes in common to clusters Y4 and E3 were associated with antibody-secreting cells (ASCs), whereas the overlap between clusters Y1 and E1 contained several interferon-related genes (Figure 4C). Although the genes shared by clusters from young and elderly vaccinated individuals had similar temporal expression profiles, the magnitude of the expression of interferon-related genes differed, being higher in young individuals (Figure 4C).

We then performed GSEA using genes ranked by their correlation with age to identify modules with enriched expression in either young or elderly vaccinees (Figure 4D). It is important to note that in this analysis we used age as a quantitative variable, rather than arbitrarily splitting the cohorts into young versus elderly. In order to ensure the robustness of our results, we performed this analysis separately on the two trials containing elderly subjects (2010 and 2011) and identified modules that were consistently enriched in both seasons. Natural killer (NK) cell-related modules showed increased expression with age on both day 3 and day 7 after vaccination, and several monocyte modules showed increased expression with age on day 7. In contrast, many B cell modules showed decreased expression in older subjects on both day 3 and day 7. These results indicate that elderly subjects might be mounting a significant innate response, but that their adaptive B cell response is diminished. One of the modules showing increased expression in the elderly was BTM M61.0 (Figure 4E), which contains many killer cell immunoglobulin-like receptor (KIR) and lectin-like receptor (KLR) genes. These are inhibitory receptors that suppress the cytotoxic activity of NK cells when bound to MHC I molecules (KIRs) and cadherins and other ligands (KLRs) (Long et al., 2013).

Next, we examined whether or not any of the age-related transcriptional differences were also associated with impaired antibody responses. By plotting the age-based enrichment score of each module against its antibody response-based enrichment score on day 7, we were able to identify five modules whose expression varied with age and were associated with either a high or low antibody response (Figure 4F). Two B cell modules, M47.0 and M69, showed increased expression in younger subjects and in subjects with high antibody responses, whereas three monocyte- and myeloid-cell-related modules, S4, M4.3, and M11.0, were enriched in the elderly and were also associated with a decreased antibody response in the 2010 and 2011 seasons. Interestingly, two of these modules, M11.0 and M4.3, showed association with an increased antibody response in the 2007 and 2008 seasons (Figure 2A), indicating that this association is not consistent across all five seasons on day 7. These results demonstrate that TIV induces distinct but overlapping transcriptional responses in the young versus elderly. The early innate response at day 1, comprised of antiviral and type IFN-related genes, seems to be impaired in the elderly. In contrast, several transcriptional modules related to NK cells and monocytes appear to have enhanced expression at baseline and after vaccination, in the elderly relative to the young.

The aforementioned transcriptional signatures of NK cells and monocytes in the elderly subjects could have been caused by de novo transcriptional induction of genes expressed in these cell types or by changing representations of these cell types within the PBMC compartment. In order to distinguish between these two possibilities, we performed FACS analysis on PBMCs in a subset of subjects from the 2010 cohort. Our analysis showed that proportions of total NK cells in elderly subjects were higher than those of young subjects at baseline as well as at all the time points studied (days 0–14). Of note, vaccination induced an increase in the percentage of NK cells in the elderly subjects; in contrast, in the young there was a reduction in the percentage of NK cells on day 1 after vaccination (Figure 5A). Similar trends were observed in the absolute numbers of NK cells (data not shown). To examine this trend in the 2011 cohort, where FACS data were unavailable, we used a deconvolution method (Abbas et al., 2009, Gaujoux and Seoighe, 2013, Shen-Orr and Gaujoux, 2013) to estimate cell frequencies based on gene expression of a mixed cell population. This approach revealed increased estimated NK cell frequencies in the elderly on day 3 after vaccination when compared to young, agreeing with the 2010 FACS data (data not shown). These results confirm our findings at a molecular level (Figure 4D) and show that the increase in NK-cell-related expression in the elderly is due in part to enhanced representation of this cell population after vaccination.

We also explored the distribution of the NK compartment into three subpopulations: CD56++ NK (CD56bright), CD56++CD16+ NK, and CD56dimCD16++NK (Figure 5B). These NK cell subsets have distinct functions: CD56dimCD16++cells, the most abundant population of NK cells in the blood, have significantly higher cytotoxic activity, whereas CD56++ cells are characterized by increased cytokine production and can exhibit immunoregulatory properties under certain conditions (Poli et al., 2009). We observed a higher frequency of CD56dimCD16++cells coincident with a lower frequency of CD56++ cells in the elderly at baseline, consistent with previous studies (Figure 5C; Solana et al., 2014). The distribution of these subpopulations was relatively stable during the response, with CD56dimCD16++ cells showing a small but significant increase in the elderly on day 1 after vaccination and CD56++ cells increasing in the young on day 7. Additionally, we examined the activation of these subsets through immunostaining for CD69, a marker associated with increased cytotoxicity and IFN-γ production (Figure 5D; Clausen et al., 2003, Gorski et al., 2006). Elderly subjects showed increasing activation after vaccination in all three subsets, with highest expression on day 3. In contrast, the young exhibited activation of the CD56dimCD16++ subset as measured by CD69 expression only on day 1 and at a lower level than in the elderly.

In addition to NK cells, we also observed higher proportions of monocytes among the elderly at baseline and during the entire duration of the study compared to young subjects (Figure S6A). In both groups, the percentage of total monocytes in PBMCs increased substantially on day 1 after vaccination, then returned to or below baseline levels. We further examined the distribution of monocytes into three subsets: “classical” (CD14+CD16−), “intermediate” (CD14+CD16+), and “non-classical” (CD14dimCD16++) (Figure S6B). The distribution of monocyte subsets was similar in both groups, with the elderly having a slight increase in intermediate monocytes compared to young (Figure S6C). The activation of these subsets was also analyzed by immunostaining for CCR5, a receptor for a number of inflammatory cytokines, and CD86, a costimulatory molecule involved in T cell activation (Figure S6D). In both the young and elderly, CCR5 expression peaked in classical and intermediate monocytes on day 1 after vaccination, corresponding with the inflammation associated with the innate immune response. Elderly subjects showed increased CCR5 expression in classical and intermediate monocytes compared to the young, whereas young subjects had higher CD86 expression in these subsets. These results are consistent with the transcriptional changes observed (Figure 4) and support the concept that persistent and excessive inflammatory responses might be detrimental to the induction of vaccine-induced antibody responses (Frasca et al., 2014, Haq and McElhaney, 2014, Pawelec et al., 2014).

A desirable feature of a good vaccine is the ability to induce long-lasting protection from infection. To investigate how effectively the influenza vaccine generates a long-lived response and identify the mechanisms responsible for this response, we measured antibody titers in a subset of subjects at 180 days (D180) after vaccination. The antibody responses peaked at day 28 (D28), with a significant decline in most subjects by D180 (Figure 6A). By applying the same criteria for seroconversion used for the D28 high- versus low-responder classification to the D180 antibody responses, we saw that just over half of D28 high responders maintained their responder status on day 180 (“persistent” responders), whereas the remainder no longer met the criteria for seroconversion on day 180 (“temporary” responders). Because the D180 antibody titer shows significant dependence on the original antibody response generated on day 28 (Figure 6B), examining the molecular profile associated with the D180 response would identify many of the mechanisms responsible for generating the D28 response. In order to study the relative persistence of the response, we calculated residuals from a linear fit between the D28 and D180 responses, thereby removing the effect of the D28 response on the D180 response. We observed that subjects with a higher D28 titer had a greater decrease in HAI response between D28 and D180, indicating increased difficulty to maintain higher levels of antibodies. The linear fit therefore represents the expected D180 response given a particular initial D28 response, and the residual can be considered an “adjusted” D180 response taking into account the original D28 response. Subjects with positive residuals have a more persistent response, whereas a negative residual represents a waning response. As expected, correlation analysis of plasmablast module BTM S3 shows a positive association with the D180 antibody response due to the effect of D28 antibody responses, generated through plasmablast expansion (Figure 6C). However, when we performed correlation between the plasmablast module activity and the D180/D28 residual, we no longer saw an association with the plasmablast expression (Figure 6D). This result suggests that factors other than the expansion of plasmablasts are responsible for maintaining a persistent antibody response.

To identify the pathways associated with a persistent or waning response, we performed GSEA on genes correlated with the D180 versus D28 residual. Among the BTMs most enriched in subjects with a persistent response on both days 3 and 7 were several modules associated with cell movement and adhesion, such as BTM 51 (Figure 6E). In fact, one of the genes in module M51 that was highly correlated with the response persistence was P-selectin (SELP), which plays an important role in the adhesion and extravasation of leukocytes out of the blood circulation and into organs or tissues (Wang et al., 2007). In addition, the day 7 expression of a number of T-cell-related modules was negatively associated with the D180/D28 residual, indicating that increased T cell responses might potentially result in more transient antibody responses (Figure 6F). The plasmablast signature at day 7, which was predictive of the HAI titers at day 28, did not correlate with the D180/D28 residual (Figure 6D).

MicroRNAs (miRNAs) have been identified as key regulators of gene expression at a post-transcriptional level (Filipowicz et al., 2008). Although there is some knowledge about miRNA regulation of immune pathways (Lodish et al., 2008), their role in responses to vaccination remains unaddressed. To assess whether or not miRNAs contribute to the differences in response we saw between young and elderly vaccinees, we measured the miRNA expression profiles of 672 human miRNAs on days 1, 3, and 7 after vaccination from a subset of subjects in the 2010 cohort (Figure 7A). We saw a significant difference in the miRNA profiles between these two groups. Whereas a majority of the differentially expressed miRNA in the elderly were upregulated after vaccination, particularly on days 1 and 7, the young exhibited predominantly downregulation of miRNA. We also identified miRNA whose expression correlated with the day 28 antibody response (Figure 7B). We again saw differences between the young and elderly subjects, with the elderly having a larger number of miRNA negatively correlating with the antibody response when compared with the young. These results suggest that miRNA might indeed be important regulators of the immune response to influenza vaccination. miRNA expression reduces the translation of its target genes through mRNA degradation or silencing via the RNA-induced silencing complex, so the increase in miRNA expression in the elderly after vaccination suggests a possible mechanism for the impaired responses in this population that merits further exploration.

To identify the miRNA that regulate particular immune pathways during vaccine response, we implemented a strategy to integrate the mRNA and miRNA omics-level measurements (Figure 7C). Because miRNA often have multiple related target genes, we first generated normalized enrichment scores for 252 BTMs on a per subject basis by performing ssGSEA on the transcriptomic data. We then performed correlation analysis between all miRNA-BTM pairs, revealing groups of modules that appear to be regulated together by several miRNA (Figure 7D) on day 1 after vaccination.

Finally, we used the TargetScan miRNA database (Garcia et al., 2011) to identify the genes in a given module predicted to be targets of negatively correlating miRNAs. Of particular interest was the regulation of module M75-antiviral interferon signature (Figure 7E) on day 1 after vaccination, because the interferon pathway plays a key role during the innate immune response. This miRNA-mRNA network suggests potential novel regulators of the interferon response after vaccination, such as miR-424. This miRNA was predicted to target OAS3, a 2′-5′-oligoadenylate synthase involved in viral RNA degradation (Samuel, 2001), as well as CXCL10, an important chemokine induced by IFN-γ that serves as a chemoattractant for lymphocytes (Dufour et al., 2002).