Seasonal Variability and Shared Molecular Signatures of Inactivated Influenza Vaccination in Young and Older Adults

Influenza is a major public health burden, particularly in high-risk populations such as older adults. The seasonal inactivated influenza vaccination is estimated to be 50–70% effective in randomized controlled trials of young adults (1–5), and efficacy is reduced to under 50% in adults over age 65 (6). Understanding the dynamics of vaccination-induced immune responses and the factors associated with immunological protection should provide insights important for improving vaccine design.

Systems vaccinology approaches using high-throughput immune profiling techniques have identified signatures of response to influenza vaccination (7–14). These include prevaccination transcriptional signatures of apoptosis-related gene modules (9) as well as B cell signaling and inflammatory modules (15). Postvaccination transcriptional signatures have also been identified, including an early IFN response one day postvaccination and a plasma cell response three and seven days postvaccination (13). IFN-stimulated genes were upregulated in both monocytes and neutrophils between 15 and 48 hours postvaccination and correlated with influenza-specific Ab responses (7, 12). In addition, the expression of genes enriched for proliferation and Ig production seven days postvaccination accurately predicted Ab response in an independent cohort (10). Studies of the influence of aging revealed that an early IFN response one to two days postvaccination as well as an oxidative phosphorylation and plasma cell response seven days postvaccination were correlated with the Ab response in young adults but were diminished or dysregulated in older adults (13, 14).

Notably, previous studies of influenza vaccine response studying the effects of aging used data from a single vaccine season (9) or from two consecutive seasons in which vaccine composition was identical (13, 14); consequently, the generalizability of these signatures is unknown. To date, no comprehensive characterization of vaccine response in both young and older adults has been reported to multiple influenza vaccines that vary in composition. To address this gap, we profiled young and older adults over five consecutive vaccination seasons (2010–2011, 2011–2012, 2012–2013, 2013–2014, and 2014–2015, hereafter referred to by the first year of each season). We developed a new automated metric to quantify Ab response while accounting for baseline titers and used this novel metric to identify predictive transcriptional signatures of vaccine response using postvaccination as well as baseline gene expression profiles.

A total of 317 subjects were recruited at Yale University over the five vaccination seasons between 2010 and 2014, and hemagglutination inhibition (HAI) titers pre- (day 0) and postvaccination (day 28) were available from the 294 subjects reported in Table I. Informed consent was obtained for all subjects under a protocol approved by the Human Subjects Research Protection Program of the Yale School of Medicine. Participants with an acute illness 2 wk prior to recruitment were excluded from the study, as were individuals with primary or acquired immune deficiency, use of immunomodulating medications, including steroids or chemotherapy, a history of malignancy other than localized skin or prostate cancer, or a history of cirrhosis or renal failure requiring hemodialysis. Blood samples were collected into Vacutainer sodium heparin tubes and serum tubes (Becton Dickinson) at four different time points: immediately prior to administration of vaccine (day 0) and on day 2 (2011, 2012, 2013, 2014) or day 4 (2010), day 7, and day 28 postvaccination.

To understand the transcriptional program underlying a successful vaccination response, we identified a subset of 134 subjects with strong or weak Ab responses to perform transcriptional profiling by microarrays. In the first three seasons, the selection criteria were a 4-fold increase to at least two strains (strong response) or no 4-fold increase to any strain (weak response), as described previously (14). In the fourth and fifth seasons, the adjusted maximum fold change (adjMFC) metric was used in addition to the fold change criteria to account for baseline titers (11). The maximum residual after baseline adjustment (maxRBA) response end point was developed after the study was completed; however, <10% (12/134) of subjects chosen for transcriptional profiling had indeterminate responses by maxRBA (neither high Ab responders [HR] nor low Ab responders [LR] using a 40% cutoff) (Table I). These 12 subjects were excluded from the predictive modeling of Ab response.

Blood samples were collected in sodium heparin tubes (Becton Dickinson) from volunteers with prior consent to an institutional review board–approved protocol. About 20 ml of fresh blood was layered on top of 20 ml of Histopaque-1077 solution (Sigma-Aldrich). Samples were carefully transferred to a centrifuge (Legend XT; Thermo Fisher Scientific) and centrifuged for 20 min at room temperature without brake at 2200 rpm. Buffy coats containing PBMCs at the interface were carefully collected to 50 ml of RPMI 1640 medium with 10% BSA. Samples were centrifuged to pellet cells, and the supernatant was discarded. Lymphocyte pellets were resuspended in 50 ml of RPMI 1640 containing 10% FBS. A 20-μl volume of suspended PBMC was mixed with equal volume trypan blue solution (Thermo Fisher Scientific) and incubated at room temperature for ∼2 min. A 20-μl cell suspension was immediately counted using a Haemocytometer under light microscope. Percent viability was determined by the formula (Number of total cells counted – Number of Blue cells counted) × 100. Samples with 95% or more viability were selected for all experiments including RNA extraction.

About 10 million cells (PBMC) in 1 ml of RPMI 1640 medium from samples with 95% or more viability were taken soon after Histopaque-1077 density gradient separation for RNA extraction. PBMC in RPMI 1640 medium were centrifuged for ∼10 min on a benchtop centrifuge at 10,000 rpm. Supernatant was removed carefully from each sample to ensure a clear pellet of PBMC without residual RPMI 1640 medium. To each pellet, ∼700 μl of QIAzol lysis reagent (QIAGEN) was added and mixed by pipetting at least 10 times to ensure proper cell lysis. Lysed cells were immediately frozen at −80°C until further extraction a using QIAcube instrument (QIAGEN).

All RNA samples were extracted using a miRNeasy Kit (catalog no. 217004; QIAGEN) following the instructions provided for using a QIAcube (QIAGEN). Briefly, samples lysed in QIAzol reagent were incubated for 5 min at room temperature (15–25°C). To each sample, ∼140 μl of chloroform was added and shaken vigorously and left at room temperature for ∼2–3 min. Subsequently, samples were centrifuged at 4°C at 12,000 × g for 15 min. The upper aqueous phase containing the RNA species were carefully transferred to a 2-ml collection tube (catalog no. 990381; QIAGEN) without touching the interphase and placed in the QIAcube for extraction. For every sample, a rotor adapter was prepared with an RNeasy Mini Spin Column at position L1 and a 1.5-ml collection tube at position L3 and placed in the QIAcube rotor. All reagents were prepared by adding the proper amount of 100% ethanol (44 ml to buffer RPE and 30 ml to buffer RWT) prior to extraction and placed in respective positions in reagent rack in QIAcube (100% ethanol in position 3, Buffer RWT in position 4, buffer RPE in position 5 and RNase-free water in position 6). RNA extraction was carried out by executing recommended protocol (FIW-50-001-J_FW_MB and PLC program version FIW-50-002-G_PLC_MB) available from the QIAcube web portal. RNA samples with RIN value above 7.0 were used for gene expression analysis.

PBMCs isolation and RNA preparation on freshly isolated PBMCs was performed as previously described (1). For cell sorting, frozen PBMCs were thawed, washed, and stained using Abs directed against: CD19(HIB19), CD4(SK3), CD20(2H7), CD8(SK1), and CD3(UCHT1) (from BD Biosciences) and using a BDFortessa instrument. T and B cells were sorted from a subgroup of samples in seasons 2010–2012. We evaluated the purity of 268/270 sorts for CD4⁺ T cells with median purity of 97%, 254/254 sorts for CD8⁺ T cells with median purity of 97%, and 252/256 sorts for B cells with median purity of 99%. 90% of the samples had a postsort purity of at least 90% for all three cell types. Very few samples had a purity of <80% (0.4, 1.1, and 5.6% for CD4 T cells, CD8 T cells, and B cells, respectively).

Serum samples were collected prevaccination (day 0) and 28 d postvaccination. HAI assays were performed as previously described (2). Virus neutralization assays (VNA) were performed as described elsewhere (3, 4). Briefly, 2-fold dilutions (50 μl) of the receptor destroying enzyme–treated sera in sterile Opti-MEM (Invitrogen, Carlsbad, CA) were mixed with 200 PFU of influenza virus (5 μl). The serum–virus samples were then incubated at room temperature for 60 min to allow any hemagglutinin (HA)-specific Abs present in the serum to neutralize the influenza virus. The serum–virus samples (55 μl) were then transferred to Madin–Darby canine kidney cell cultures cultured in 96-well flat-bottom plates. Following virus absorption for 60 min, the serum–virus inocula were removed, and the Madin–Darby canine kidney cells were cultured for 4 d in Opti-MEM supplemented with 1 μg/ml tosylsulfonyl phenylalanyl chloromethyl ketone–trypsin (Sigma-Aldrich). Virus production was determined by HA assay. The neutralization titer was defined as the reciprocal of the highest dilution of serum that neutralizes 200 PFU of influenza virus.

To adjust for inverse correlations between HAI titer fold changes and baseline titers, we developed an automated metric: maxRBA. First, young and older cohorts were separated, and endpoints were calculated in each season and each age group separately. Any fold changes <1 were set to 1 because we did not expect the number of Abs in the blood to decrease by 2-fold in the 28 d after vaccination and this was likely due to measurement error. Next, the baseline and fold changes were log₂ transformed, and an exponential curve was fit to the fold change versus baseline titers for each strain. Next, the residuals were calculated, and for each subject, the maximum residual across all strains was selected as the maxRBA. Finally, HR and LR were defined as the top and bottom 40th percentile of maxRBA, respectively. The maxRBA values presented in Fig. 1 and Supplemental Fig. 3 were calculated by combining young and older adults within each season to allow for comparison across age groups. The code to calculate maxRBA is available in the Calculate_maxRBA() function from the titer R package (https://bitbucket.org/kleinstein/titer).

Each RNA sample was quantified, and integrity was assessed by the Agilent 2100 Bioanalyser (Agilent). Samples were processed for cRNA generation using the Illumina TotalPrep cRNA Amplification Kit and subsequently hybridized to the Human HT12-V4.0 BeadChip (Illumina). For gene expression analyses, samples were processed and hybridized to HumanHT-12v4 Expression BeadChip (Illumina). Arrays were processed at Yale’s Keck Biotechnology Resource Laboratory, and raw expression data were output using Illumina GenomeStudio software. Samples from each season were processed in batches, and all samples from each subject were run on the same chip to mitigate batch effects. Data from each season are available via ImmPort (https://www.immport.org) under accession numbers SDY63, SDY404, SDY400, SDY520, and SDY640. Microarray data are available through the Gene Expression Omnibus Database (https://www.ncbi.nlm.nih.gov/geo/) with accession numbers GSE59635, GSE59654, GSE59743, GSE101709, and GSE101710 for PBMC data and GSE65440, GSE65442, and GSE95584 for B and T cell data. Data from seasons 2010 and 2011 (GSE59635 and GSE59654) were previously published (1) as well as data from day 0 and day 7 time points from the 2012 season (GSE59743) (5). The day 0 expression in PBMC (GSE59635, GSE59654, GSE59743, GSE101709, and GSE101710) of a single gene, MINCLE, was also recently published (6). The remainder of the data in this work has not, to our knowledge, been previously published. The data were quantile normalized and log₂ transformed within each season. Multiple probes were collapsed to unique Entrez Gene IDs by selecting the probe with the highest average expression. The Bioconductor package illuminaHumanv4.db version 1.26 was used to map the probes (7).

Differentially expressed modules (DEMs) were defined at each postvaccination time point using Quantitative Set Analysis for Gene Expression (QuSAGE) with a false discovery rate (FDR) <0.05 (8). Differentially expressed genes (DEGs) were defined at each postvaccination time point using two criteria: 1): an absolute fold change of at least 1.25 relative to the prevaccination time point and 2) a significant change in expression by limma (version 3.30.7) after correction for multiple hypothesis testing (FDR <0.05) (9). Clusters of DEGs were determined by hierarchical clustering using Ward minimum-variance method, with the distance between any two genes defined by one less the Pearson correlation coefficient. Enrichment of DEG clusters was performed using Enrichr v1.0 with the following databases: GO_Cellular_Component_2018, GO_Biological_Process_2018, GO_Molecular_Function_2018, KEGG_2016, and Reactome_2016 (10). Enrichment of principal components (PC) was performed using the geneSetTest function in the limma R package v3.24.15 (9). Control of the FDR was performed according to the procedure in Ref. 11 unless otherwise noted.

First, genes were filtered to those that were detected in at least 20% of samples with a p value <0.05. For the gene module meta-analysis, QuSAGE activity distributions were taken from the single-season analysis and convoluted to create a probability density function for the meta-analysis (12). Modules with an FDR <0.10 in the meta-analysis were chosen as significant. For the single gene meta-analysis, a random effects model was fit for every gene at each time point in each age group separately. This model allows for variation in effect sizes across seasons. Effect sizes (mean differences) were calculated for every gene in each season separately and then combined using the rma() function from the metafor R package (13). The restricted maximum likelihood estimator for the amount of heterogeneity was used (14). Genes with an FDR <0.05 in the meta-analysis were chosen as significant.

Because females tended to respond better than males, genes located on the X and Y chromosomes were removed to avoid selection of sex-linked genes that may be confounded with vaccine response. For baseline predictors, the 1000 genes with the largest coefficient of variation were selected as the initial feature set. For postvaccination predictors, the log fold change from day 0 was calculated for each gene, and the 1000 genes with the largest fold change magnitudes were selected as the initial feature set. The baseline or fold changes were standardized by subtracting the mean and dividing by the SD. Finally, this preprocessed data from each individual season were combined to form the discovery data. The young adult models were tested on GSE47353 at baseline, day 1, day 7, and day 70, whereas the older adult models were tested on GSE41080 at baseline and GSE74813 at baseline, day 1, day 7, and day 14 postvaccination. The logistic multiple network–constrained regression (LogMiNeR) framework was performed as previously described (5). Briefly, 5-fold cross-validation was used to select the optimal tuning parameters, and 50 iterations of cross-validation were performed on different splits of the discovery data set. The prior knowledge networks were defined for Reactome (15, 16), gene ontology (17), blood transcriptional modules (BTM) (18) and cell-specific signatures (19) by connecting all pairs of genes within each gene set. The Kyoto Encyclopedia of Genes and Genomes (KEGG) network incorporated pathway topology and was built using the KEGGgraph R package v1.26.0 (20). ImmuneGlobal (ImmuNet) and ExpOnly_ImmuneGlobal (ImmuNet_Exp) networks were obtained from ImmuNet, and edges were restricted to those with confidence of at least 0.1 (21). The Search Tool for Retrieval of Interacting Genes/Proteins network incorporated all experimental evidence from the Search Tool for Retrieval of Interacting Genes/Proteins database v10.0 (22).

Influenza vaccination signatures were manually curated from several publications (23–28). DEGs were defined for each contrast using a significant change in expression by limma (version 3.30.7) (p < 0.001). DEMs were defined as described above using QuSAGE but with a p value cutoff of p < 0.001. Single gene meta-analysis and gene module meta-analysis were performed as described above for each published signature without filtering lowly expressed genes. HAI-associated comparisons were discretized and tested using maxRBA-defined HR versus LR, as described above.

We evaluated 294 healthy young (21–30 y old, n = 147) and older (≥65 y old, n = 147) adults over five consecutive influenza vaccination seasons from 2010 to 2014. All subjects received the standard dose trivalent (2010–2012) or quadrivalent (2013, 2014) seasonal inactivated influenza vaccine. We measured influenza-specific HAI titers prevaccination (day 0) and 28 d postvaccination. Over the course of our study, the vaccine composition changed, relative to the previous season in three of five seasons (Table I).

In all seasons, prevaccination titers were negatively correlated with the increase in titers postvaccination (Supplemental Fig. 1). Previous work defined an adjMFC end point that removes the nonlinear correlation between fold change and baseline titers (11). However, adjMFC separates subjects into manually defined bins, making it difficult to perform high-throughput analysis. Furthermore, adjMFC does not allow for information sharing between bins, as each bin is adjusted independently. To address these limitations, we developed maxRBA, which corrects for the dependence on baseline titers for each strain by modeling titer fold changes as an exponential function of prevaccination titers and selecting the maximum residual across strains (Fig. 1A). All vaccine strains were approximately equally responsible for the maximum residual in any given season. HR and LR were defined as the top and bottom 40th percentiles of the residuals, respectively. maxRBA can be interpreted as the maximum change from expected fold change given the initial titer; it is fully automated, is strain agnostic, and is correlated with plasmablast frequencies 7 d postvaccination (Supplemental Fig. 2A, 2B). Thus, maxRBA allows a completely automated assessment of the relative strength of each subject’s Ab response independent of pre-existing Ab titers.

Older adults had significantly lower prevaccination titers than young adults for three of five seasons (Fig. 1B). The maximum fold change to any vaccine strain showed an increasing trend in young adults compared with older adults (Supplemental Fig. 3C). Because of the inverse relationship between baseline titers and fold change (Supplemental Fig. 1), we adjusted for baseline titers using maxRBA and found that the difference in vaccine response between young and older adults was statistically significant in more seasons (Fig. 1C). Males and females had similar prevaccine geometric mean titers (Supplemental Fig. 3A). However, the Ab response calculated by maxRBA showed a trend toward stronger Ab responses in females compared with males with similar baseline titers in both age groups (Fisher combined p = 0.02 (young), p = 0.12 (older); Supplemental Fig. 3B). We did not detect any significant difference in baseline titers or titer responses across seasons when stratifying subjects by body mass index, smoking history, aspirin use, or diabetes medication use (p > 0.05, two-sided Wilcoxon rank-sum test (discrete) or simple linear regression [continuous]).

We also examined the dynamics of viral titers over the course of the five seasons (Supplemental Fig. 3D). The A/California 7/2009 H1N1 strain was introduced into the seasonal vaccine in 2010 and remained through the 2014 season; however, prevaccine titers to this strain were consistently lower in older versus young adults for 2011–2014. Although we did not follow the same subjects across multiple seasons, 50–80% of young and 80–98% of older adults self-reported receiving influenza vaccine in the previous year. Taken together, these results support existing evidence that the capability for Ab persistence is reduced with age (16).

To identify correlates and predictors of vaccine response, we selected a subset of individuals (20–40 subjects per season) from young and older adult cohorts who had strong or weak Ab responses according to HAI titers and performed longitudinal transcriptional profiling prevaccination (baseline) and 4 (2010 cohort) or 2 (all other cohorts), 7, and 28 d postvaccination (Table I; 2Materials and Methods). We first performed differential expression analysis independently in each season without differentiating subjects by Ab response. We compared each postvaccination time point to baseline and found a vaccine-induced signature that comprised a total of 2462 significantly DEGs over all five seasons (FDR <0.05, absolute fold change >1.25; Supplemental Table I).

Most of the DEGs were from the first two seasons, whereas vaccination in the latter three seasons induced relatively weak changes (Fig. 2A, Supplemental Fig. 4E). In fact, a substantial fraction of DEGs were unique to a single season and not differentially expressed at any time point in another season (young: 38%, older: 75%). In young adults, there were 1330 DEGs shared across two or more seasons, whereas in older adults there were 265 shared DEGs. In both young and older adults, a substantial fraction of these shared genes was differentially expressed 28 d postvaccination (Supplemental Fig. 4F). To assess whether vaccine-induced changes were consistent between seasons, we divided the 2462 DEGs into seven clusters by hierarchical clustering (Fig. 2A, Supplemental Table II) and tested for their activity in every season using QuSAGE (17) (Supplemental Fig. 5). In young adults, three of the clusters (B, F, G) had significant but opposite activity during the 2010 and 2011 seasons, whereas these clusters were relatively consistent across seasons in older adults. Genes in cluster A were induced strongly in the 2011 season in both age groups and notably enriched for multiple pathways related to mitochondria, including mitochondrial inner membrane, oxidative phosphorylation, respiratory electron transport, citric acid cycle and respiratory electron transport, and mitochondrial respiratory chain complex assembly (FDR <0.05; Supplemental Table II). These findings reflect our previous identification of a mitochondrial biogenesis signature associated with influenza vaccine Ab response (14). Cluster D was only significantly induced in the 2013 season at 7 and 28 d postvaccination and was not significantly enriched for any gene sets tested (FDR >0.05; Supplemental Table II). The cluster with the most consistent expression pattern across the five seasons was cluster C, which was enriched for pathways related to TLR signaling, B and T cell signaling, NF-κB signaling, MAPK signaling, cell senescence or proliferation, and apoptosis (Supplemental Table II). Interestingly, cluster C contains three genes (DUSP1, DUSP2, CCL3L3) that were significantly downregulated 28 d postvaccination in four of five seasons. CCL3L3 is a ligand for CCR1, CCR3, and CCR5, known to be chemotactic for monocytes and lymphocytes (18). DUSP1 and DUSP2 are dual-specificity phosphatases; DUSP2 dephosphorylates STAT3, leading to inhibition of survival and proliferation signals (19–21), and an age-associated decrease in DUSP1 function contributed to inappropriate IL-10 production in monocytes before and after influenza vaccination (22). To determine whether downregulation of these three genes was a result of changes in cell subset composition or observed in subpopulations of cells, we performed transcriptional profiling on sorted B and T cells in a subset of individuals from three seasons. DUSP1 and DUSP2 but not CCL3L3 were significantly downregulated 28 d postvaccination over multiple seasons in CD4 and CD8 T cells of young adults (p < 0.01, one-sided t test; Fig. 2C, Supplemental Fig. 4C, 4D). Furthermore, whereas DUSP2 was only significantly decreased in PBMCs of older individuals in the 2011 season, expression of DUSP2 was significantly decreased 28 d postvaccination in sorted CD4 and CD8 T cells from older individuals in multiple seasons (Fig. 2C). Thus, the downregulation of DUSP2 28 d postvaccination is observed in the T cell compartment of both young and older adults.

To further assess shared patterns in vaccine-induced changes across five seasons, we performed a PC analysis (PCA) on gene expression fold changes postvaccination for all DEGs. The first two components together explained 38% of the variation in young adults’ and 46% of the variation in older adults’ transcriptional changes postvaccination (Fig. 2B, Supplemental Fig. 4B). Notably, in young adults, the 2011 and 2014 seasons (both with vaccine composition identical to the previous year) had similar trajectories, increasing along PC1 by day 28 postvaccine. Examining the genes contributing to PC1 reveals that four of the top 10 genes (SLMAP, MATR3, MBNL3, RANBP3) increase in expression postvaccination more in the 2011 and 2014 seasons than in any other season. The shared trajectories along PC1 are not significantly enriched for any BTMs (23), KEGG pathways (24), or cell subset signatures (25) (FDR >0.05; Supplemental Table III). The trajectory of the 2010 season was quite distinct from the other seasons in young adults. This season is consistently elevated on PC2, which is significantly enriched for monocytes, TLRs, and inflammatory signaling (FDR <0.05; Supplemental Table III). The 2012 and 2013 seasons also appear to have similar trajectories, both decreasing in PC2 over time. The vaccines used in these two seasons each introduced multiple new strains and also retained the A/California/7/2009 strain. Five of the top 10 genes (ZNF493, ZNF652, OCIAD1, C21orf58, IL-11RA) contributing to PC2 increased in expression 28 d postvaccination in the 2012 and 2013 seasons, whereas they decreased in expression in the other seasons. This differential expression analysis shows that there are large variations in vaccine-induced transcriptional signatures between seasons that in young adults might be explained, in part, by vaccine composition.

Given the substantial seasonal variation in the number of DEGs, we next performed an analysis of differential expression of gene modules using QuSAGE to quantify the gene module activity of 346 previously defined BTMs (23). There were 262 DEMs (FDR <0.05; Supplemental Fig. 4A, Supplemental Table IV). Similar to the gene level analysis, no significant changes were identified in the 2014 season, but six modules (cell cycle and growth arrest [M31], chemokines and inflammatory molecules in myeloid cells [M86.0], enriched for TF motif TTCNRGNNNNTTC, leukocyte differentiation [M160], putative targets of PAX3 [M89.1], and signaling in T cells [I] [M35.0]) were significantly downregulated in young adults at day 28 in four of five seasons (Fig. 2D). These changes were largely driven by decreases in DUSP1/2, EGR1/2, JUN/JUNB, FOS/FOSB, TNF, CD83, and IL-1B. Thus, whereas there was substantial variability in the signatures induced by vaccination across multiple seasons, there is a shared signature consisting of three genes and six modules that was downregulated at day 28 in four of five seasons.

The differential expression approach is limited by fixed fold change and significance cutoffs that may vary between seasons. To increase our power to identify shared signatures across seasons and in older adults, we performed a meta-analysis at the individual gene and gene module level. We identified 338 genes with significantly altered expression postvaccination (FDR <0.05; Supplemental Table V). In young adults, we identified significant genes at day 2, day 7, and day 28, with little overlap among genes on each day. Genes induced on day 2 were moderately enriched for innate immune genes from InnateDB (http://innatedb.com/), including MYH9, TYK2, GLRX, and IP6K1 (p = 0.12, hypergeometric test). Some of the genes consistently induced at day 7 included IGLL1, CD38, ITM2C, TNFRSF17, MZB1, and TXNDC5. We previously identified TNFRSF17, B cell maturation Ag, as induced 7 d following influenza vaccination (26), and it was also identified as a predictive marker gene of Ab response to multiple vaccines, including influenza, meningococcal conjugate (MCV4), and yellow fever (YF17D) vaccines (11, 23, 27–29). Consistent with the individual season analysis, the majority of genes identified by the meta-analysis were altered at day 28; these day 28 DEGs included DUSP1, DUSP2, and CCL3L3, identified in the single-season analysis, and many other downregulated genes, including IL-1B, CCL3, and JAK1. Thus, there are consistent changes identified across all seasons in young adults at every time point measured.

In older adults, we identified 125 genes with significantly altered expression at day 28, but no genes with significantly altered expression at day 2 or day 7 (Supplemental Table V). The most significantly increased gene at day 28 is XRN1, the primary 5′ to 3′ cytoplasmic exonuclease involved in mRNA degradation (30). XRN1 plays a critical role in the control of RNA stability in general but, in addition, appears to regulate the response to viral infection at several levels [for example, by targeting viral RNAs for degradation (31)], or regulating levels of potential activating ligands such as dsRNA (32). Notably, XRN1 has also been reported to facilitate replication of influenza and other viruses by inhibiting host gene expression (33, 34), suggesting that dysregulated expression of XRN1 in older adults could influence host response to vaccination. We identified three genes shared between both age groups: ARRDC3 and USP30 were downregulated, whereas TNPO1 was upregulated, all at day 28. ARRDC3 encodes a member of the arrestin protein family that regulates G protein–mediated signaling and is implicated in regulating metabolism (35). USP30 is a ubiquitin-specific protease that acts as a mitochondrial deubiquitinating enzyme (36). TNPO1 encodes Transportin-1 that serves to import proteins into the nucleus (37). The effect sizes of all genes at day 28 were positively correlated between young and older adults with weak positive associations at day 2 and day 7 (Fig. 3). These results provide additional evidence that transcriptional changes are broadly similar in young and older adults at day 28 postvaccine.

We carried out a gene set level meta-analysis using QuSAGE to combine probability density estimates of gene module activity for each season (38). We identified 186 BTMs significantly altered postvaccination across five seasons (FDR <0.05; Supplemental Table IV). The module with the largest increase in activity was plasma cells, Igs (M156.1), which peaked on day 7 with a combined fold change of 1.17 in young adults and 1.08 in older adults at day 7. Most BTMs showing significant changes were identified in young adults and, unlike the individual gene level, there was a large overlap between sets at each time point, suggesting the same module changes were sustained over the 28 d following vaccination (Supplemental Fig. 4A). Indeed, a heatmap of module activity shows that in young adults transcriptional changes continued to intensify at day 28 for many modules rather than returning to the baseline state (Supplemental Fig. 6). Older adults showed a qualitatively similar pattern to young adults on day 2 and day 28 but not day 7. The majority (40/59) of the modules significantly altered in older adults on day 28 were also significantly altered in young adults at day 28 (Supplemental Fig. 4A). The modules downregulated on day 28 in both young and older adults were annotated with Ag processing and presentation (M95.0, M95.1, M28, M71, M200, M5.0) and T cell activation (M36, M44, M52). The modules upregulated on day 28 included Golgi membrane (II) (M237), enriched in DNA-interacting proteins (M182), and chaperonin-mediated protein folding (I, II) (M204.0, M204.1). Taken together, the high correlation between individual gene changes and overlap of many BTMs suggests a convergence toward a common transcriptional program in young and older adults at day 28.

A meta-analysis across all five seasons revealed markedly different baseline transcriptional profiles in young versus older adults, with 1072 genes significantly altered (FDR <0.05, Supplemental Table VI). Of these age-associated genes, 204 genes were also significantly induced by the vaccine in young adults and 125 genes were induced in older adults. We tested whether age-associated genes were enriched for vaccine-induced genes at each time point and found that the overlap was significantly more than expected by chance for the six age-associated genes induced on day 7 in young adults (p = 0.017, hypergeometric test). Of these six overlapping genes, five genes (ITM2C, MZB1, IGLL1, TNFRSF17, and TXNDC5) exhibited decreased basal expression in older adults, whereas one (SELENOS) exhibited increased basal expression compared with young adults. Although these genes were induced in young adults, they were not significantly induced in older adults on day 7. Notably, MZB1 and TNFRSF17 are B cell–associated genes, suggesting that older adults have decreased B cell activity prevaccination and fail to induce the same B cell response as young adults at day 7. SELENOS encodes selenoprotein S, which is involved in degrading misfolded endoplasmic reticulum (ER) proteins and influences inflammation via the ER stress response (39, 40). Our results show that age-associated genes are significantly overrepresented in the set of genes altered in young adults 7 d postvaccination.

We next performed a meta-analysis of BTMs between age groups at baseline and identified 120 modules significantly altered with age (FDR <0.05, Supplemental Table VII). Most of the modules that were decreased with age were associated with adaptive immunity, whereas those that had increased expression with age were mostly innate and inflammatory modules (reflecting age-associated inflammatory dysregulation; Supplemental Fig. 7B). Of these 120 modules, 52 were also significantly altered postvaccination; however, the overlap at each time point was not significantly more than expected by chance (hypergeometric p > 0.05, Supplemental Fig. 7A). Thus, age-related genes are enriched among the genes induced at day 7 in young adults, whereas no gene modules were significantly overrepresented.

We next asked whether any transcriptional changes postvaccination could discriminate HR from LR. Regularized logistic regression models with an L1 (lasso) or with L1 and L2 (elastic net) penalties were fit to identify genes predictive of Ab response. In addition, to identify biologically interpretable predictors we used the LogMiNeR framework (26) that facilitates the generation of predictive models with improved biological interpretability over standard methods. We combined the fold changes in gene expression data postvaccination from five seasons and trained LogMiNeR to predict HR versus LR in young and older cohorts separately. At each time point, models were trained on all five seasons of data (except for day 2, which was not available in the 2010 season; see 2Materials and Methods). Publicly available data sets from independent groups were used to validate the models. For the models built from expression changes at day 2 or day 28, no studies at identical time points were available, so we attempted to validate these models on studies with similar time points (day 1 or 3 in Ref. 11 and day 14 in Ref. 13). Although we could build predictive models on our data (median area under the curve [AUC] ≥ 0.75), they did not validate on other data sets at the (different) time points available (median AUC ≤ 0.55).

For day 7 postvaccine, direct validation data were available in independent datasets. In young adults, day 7 models were predictive for HR in the discovery and validation (11) datasets (Fig. 4A). Another MAP kinase phosphatase acting on ERK1/2, DUSP5, was 1 of 37 genes selected by the lasso model whose expression was increased in HR (Fig. 4C). DUSP5 is expressed in multiple immune cell types such as B cells (including plasma cells), T cells, dendritic cells, macrophages and eosinophils (41). In murine T cells, DUSP5 appears to promote the development of short-lived effector CD8⁺ T cells and inhibit memory precursor effector cell generation in a lymphocytic choriomeningitis virus infection model (42); whereas optimizing memory precursor cell generation would be the goal of vaccination, the upregulation of DUSP5 in HR could reflect regulation of the balance between short-lived versus memory precursor effector CD8⁺ T cells. A sensitivity analysis of the maxRBA cutoff shows that the average expression of predictive genes is consistent across a range of definitions for HR and LR (20th–40th percentile; Supplemental Fig. 2C, 2D). Using LogMiNeR, the models were consistently enriched for the B cell signature as well as the KEGG chemokine signaling pathway (Supplemental Table VIII).

In older adults, models predicting Ab responses built from day 7 gene expression were highly predictive in the discovery dataset but did not validate on an independent dataset (13) (Fig. 4B, 4D). Expression of the solute carrier family 25 gene SLC25A20 of mitochondrial transporters contributes to predicting HR versus LR in older adults. SLC25A20 is the carrier for carnitine and acylcarnitine (43) and so would be expected to be crucial for the transport of fatty acids into mitochondria. The models of response in older adults were significantly enriched for several BTMs of monocyte signatures as well as TLR and inflammatory signaling (M16), which positively predicted vaccine response; together with previous studies linking age-associated impairments in TLR function to influenza vaccine Ab response (44, 45), these findings provide additional support for the crucial role of innate immune function in vaccination (Supplemental Table VIII).

Notably, none of the models built in young adults at any time point are predictive in older adults (AUC ≤ 0.5). In fact, models built on transcriptional changes at day 28 in young adults had a median AUC near 0.8 in young adults but no more than 0.3 in older adults, suggesting that the same genes predictive of HR in young adults predicted LR in older adults (Supplemental Fig. 8E). The lasso models making these predictions often chose a single gene, killer cell lectin-like receptor B1 (KLRB1, also known as CD161), which was driving this inverse pattern (Fig. 4E). KLRB1 is an inhibitory receptor on NK cells (46, 47) and is also a biomarker of Th17 cells (48–50). Notably, changes in KLRB1 expression in sorted CD4 and CD8 T cells at day 28 closely mirrored the changes in PBMCs for young but not older adults (Supplemental Fig. 8A, 8B). We confirmed this inverse correlation between age groups on a genome-wide scale by performing a meta-analysis comparing HR versus LR (Supplemental Table IX). We observed a weak negative correlation in effect sizes between young and older adults at day 28 (r = −0.27; Fig. 4F). We confirmed this negative correlation in effect sizes between young and older adults using a VNA in a test sample of blood from seasons 2011 and 2012 (r = −0.32; Supplemental Fig. 8D). Thus, expression changes of many genes at day 28 have opposing signs between age groups for the effect size comparing HR versus LR, and a single gene, KLRB1, predicts response with AUC >0.7 in opposing directions in young versus older adults.

We next sought to identify baseline transcriptional predictors of Ab response. In young adults, LogMiNeR models were predictive above random on discovery and validation (11) (Fig. 5A) datasets. Lasso models included the gene VASH1, known as an angiogenesis inhibitor and mediator of stress resistance in endothelial cells, which was expressed at lower levels in HR (Fig. 5C); notably, the KEGG gene set leukocyte transendothelial migration was significantly enriched in over 50% of the models when LogMiNeR was used with ImmuNet as prior knowledge (51). Another predictive gene, EIF4E, a translation initiation factor important in type I IFN production, was decreased in HR. A sensitivity analysis of the maxRBA cutoff shows that the average expression of predictive genes is consistent across a range of definitions for HR and LR (20th – 40th percentile; Supplemental Fig. 2E, 2F). Finally, the BTMs cell adhesion (M51) and B cell surface signature (S2) were consistently enriched in the models (Supplemental Table VIII). In older adults, LogMiNeR models were also predictive on the discovery and one validation dataset (9) (Fig. 5B) but not another (13) (Supplemental Fig. 8C). Two of the individual genes that predict response, ALDH1A1 and ALDH3B1, are aldehyde dehydrogenases that metabolize vitamin A to retinoic acid (Fig. 5D). Recently, aldehyde dehydrogenases were implicated in antiviral innate immunity as mediators of the IFN response through their role in the biogenesis of retinoic acid (52). Multiple monocyte gene sets are enriched in the predictive genes, including the BTM enriched in monocytes (II) (M11.0), which negatively predicts vaccine response (Supplemental Table VIII). Thus, these baseline predictive models built from five seasons of transcriptional profiling data provide further evidence for functional distinctions present in subjects prior to vaccination that influence the immunologic response to influenza vaccine in young and older adults.

To link our findings to previously identified influenza vaccine signatures, we performed a comprehensive assessment of the behavior of 1603 previously published individual gene and gene module signatures in our data set. We manually curated published signatures from studies that carried out transcriptional profiling on adult cohorts after influenza vaccination (9, 11, 13, 15, 27, 53). We further limited the signatures to shared time points postvaccination. This set of findings describe 935 response-associated and 653 temporal signatures in B cells and PBMCs as well as 15 age-associated signatures (Supplemental Table X).

Most of the previously published signatures we validated in our data were single genes induced 7 d postvaccination in PBMCs or B cells (Supplemental Fig. 9). Of the 135 signatures that showed significant differential expression (p < 0.001), 103 changed in the same direction as the published signature. In PBMCs we validated 26 day 7 vaccine-induced genes, including four genes independently discovered in our meta-analysis: CD38, ITM2C, TNFRSF17, and SPATS2 (Supplemental Fig. 9B) (11). CD38 is upregulated on the surface of Ab-secreting cells, and TNFRSF17, or B cell maturation Ag (BCMA) is a receptor for BAFF expressed on memory B cells and plasma cells (54). Notably, validated vaccine-induced genes in B cells include several associated with mitochondrial function whose expression was upregulated at day 7, including ubiquinol cytochrome C reductase, complex III, subunit VII (UQCR), NAD-dependent malic enzyme (ME2), transaldolase 1 (TAL), and glycine decarboxylase (GLDC) (Supplemental Fig. 9A). We validated several modules significantly associated with Ab response at baseline in young and older adults (Supplemental Fig. 9D) (13). Of these modules, one positively associated with Ab response (enriched in B cells [I] [M47.0]) is enriched in our baseline predictive model of young adults, and three negatively associated with Ab response are enriched in our baseline predictive model of older adults (monocyte surface signature [S4], myeloid cell-enriched receptors and transporters [M4.3], enriched in monocytes [II] [M11.0]). Interestingly, these latter three modules are also enriched in predictive models of HR versus LR from day 7 fold changes. Finally, there are seven validated single genes whose fold change at day 7 is positively associated with Ab response in young adults (Supplemental Fig. 9C) (11). One of these genes, HSP90B1, or gp96 (an ER-based chaperone protein implicated in innate and adaptive immune function) is also selected as a predictive gene of Ab response (55, 56).

This study is the first, to our knowledge, to evaluate the transcriptomic response to influenza vaccination in young and older adults in five consecutive vaccine seasons with three different vaccine compositions. We sought to address whether common signatures of vaccine response or transcriptional predictors of Ab response could be elucidated despite differences in seasonal vaccine composition.

To adjust for the inverse relationship between baseline Ab titers and vaccine-induced Ab production, we developed a novel vaccine response end point, maxRBA, to automatically correct for variation in baseline titers; this allowed us to demonstrate an age-associated decrease in Ab response in gender-matched participants. Comparing the transcriptional profiles across five seasons revealed substantial seasonal variability in both the magnitude as well as direction of response. For example, the vaccines administered in the 2010 and 2011 seasons elicited large changes in gene expression, but no statistically significant DEGs were found in the 2014 season, despite a comparable sample size. Potentially, the large transcriptional changes observed in 2010 and 2011 could reflect the introduction of the A/California/7/2009 viral pandemic strain to the seasonal vaccine (as well as a change in the H3N2 vaccine strain beginning in 2010, the only year of the five studied when both influenza A strains changed). Notably, a PCA revealed similar vaccine-induced signatures in the 2011 and 2014 seasons and in the 2012 and 2013 seasons. The similarities between the 2011 and 2014 seasons are intriguing because in both seasons the composition of the vaccine was identical to that in the preceding year, perhaps suggesting that these gene signatures reflect a relatively recent recall response. In contrast, the 2012 and 2013 vaccines each contained two strains that had not been present in the previous year’s vaccine. However, we did not observe the same trends in older adults, and the design of this study precluded the derivation of strain-specific transcriptional signatures. Nonetheless, our results indicate that changes in vaccine composition, influencing factors, such as vaccine strain immunogenicity, and the effects of previous vaccination or infection can alter the transcriptional response to influenza immunization.

Despite substantial interseason variability, we identified shared vaccine-induced signatures in both young and older adults at day 28. We expected day 28 expression profiles to be similar to baseline; however, there were numerous transcriptional changes at day 28 that were consistent across seasons with different vaccine compositions. Some of the most significant changes identified from single-season differential expression analysis in four out of five seasons were in DUSP1, DUSP2, and CCL3L3; moreover, DUSP2 expression was also decreased in sorted CD4⁺ and CD8⁺ T cells from both young and older adults at day 28. It is notable that a basal age-related alteration in phosphorylation of DUSP1, a negative regulator of IL-10 production, was associated with increased expression of IL-10 in monocytes from older adults (seen before and after influenza vaccination) (22) and that increased DUSP6 expression was associated with impaired TCR signaling in CD4⁺ T cells from older adults (57). These results emphasize the importance of modulation of MAP kinase function, such as through phosphatases of the DUSP family, in the regulation of influenza vaccine response. Surprisingly, early response signatures at day 2 and day 7 postvaccination were not as consistent across seasons as day 28 signatures in a meta-analysis of genes and gene modules. One potential hypothesis that explains this observation is that temporal variations in early responses across seasons were not captured at the time points used and that responses at day 28 are less variable and thus were captured in every season. It is possible that this common transcriptional program at day 28 reflects a convergence toward resolution of the vaccine response in both young and older adults. However, a substantial number of BTMs showed upregulated activity at day 28 without evidence of resolution to baseline, particularly in young adults; notably, we previously found evidence of enhanced TNF-α and IL-6 production in monocytes 28 d post–influenza immunization (22) that was blunted in monocytes from older adults. Thus, it remains possible that the transcriptional signature we observed also reflects elements of an ongoing immune-activated state several weeks after vaccination.

We built predictive models of Ab response from postvaccination transcriptional responses that were successfully validated in an independent cohort of young adults. Although transcriptional changes were correlated between age groups at day 28, models of Ab response built in young adults did not validate in older adults. Strikingly, we identified a genome-wide inverse correlation between the effect size of genes discriminating HR and LR at day 28 and confirmed this finding with both HAI and VNA titers. A similar inverse correlation related to age was recently reported using baseline (day 0) gene expression signatures (15). We identified a single gene, KLRB1, whose expression alone predicted response in both age groups but in opposite directions. In young adults, changes in KLRB1 expression were also observed in sorted CD4 and CD8 T cells, perhaps reflecting the finding that KLRB1 expression is increased in populations of memory T cells (58). Furthermore, KLRB1^hi CD8⁺ T cells are self-renewing memory cells that are able to reconstitute the memory T cell pool after chemotherapy (59). Thus, KLRB1 induction in young adults may reflect an increase in memory T cell populations. In older adults, these expression patterns were not observed in sorted T cells, implying that KLRB1 expression in another cell type, perhaps NK cells or Th17 cells, was the basis for the predictive performance.

We also built and validated predictive models of Ab response in young and older adults from day 0 gene expression data. One of the predictive genes in young adults, VASH1, showed evidence of genetic regulation of gene expression in a previous study of influenza vaccination, suggesting that genotype may have predictive power to explain the Ab response (8). Leukocyte migration and a B cell surface signature were enriched in the predictive models. This is consistent with a recently reported meta-analysis that included baseline transcriptional profiles from the 2010, 2011, and 2012 seasons of the current study and validated a temporally stable BCR signaling gene module that positively predicted response at baseline (15). Whereas the B cell surface signature (S2) module we identified was not the same one identified in the previous study, our findings further support the implication of B cell transcriptional signatures as prevaccine biomarkers of Ab response in young adults. In older adults, we incorporated prior knowledge on gene coexpression using LogMiNeR to identify monocyte signatures that were enriched in the predictive models and were negatively associated with Ab response. Our model validated on one older adult cohort (9) but not another (13); this may reflect substantial variability in cohorts of older adults, which would be expected to be more heterogeneous in terms of comorbid medical conditions or medication use compared with young adults. Finally, we linked our findings to previously identified influenza vaccination signatures by performing a comprehensive assessment of 1603 previously published individual gene and gene module signatures. We present the signatures that validate in any season or a meta-analysis of all seasons of our data to highlight the most consistent set of genes and gene modules associated with vaccination or Ab response in PBMC and B cells.

One limitation of this study is that it was not placebo controlled. Thus, we cannot rule out that some of the changes in gene expression postvaccination are unrelated to Ag-specific vaccine responses. Participants were specifically recruited from among those presenting to the clinic to receive an influenza vaccination, and placebo vaccination was not a practical consideration, especially in older adults who are at higher risk of morbidity and mortality from influenza infection. Nevertheless, the fact that transcriptional changes are occurring postvaccination means that they have the potential to influence the ultimate quality of the response.

In summary, we recruited nearly 300 young and older adults across five vaccination seasons and, despite substantial seasonal variability in vaccine-induced transcriptional signatures, identified a core transcriptional signature shared between seasons and across age groups 28 d postvaccination. Our results suggest that vaccine composition, in concert with differences in pre-existing immunity and other individual factors, dramatically influences immune response to inactivated influenza vaccination. A deeper understanding of the cause of the shared transcriptional signatures we consistently identified in both young and older adults 28 d postvaccination may help to better understand the mechanisms by which the current vaccines function. In addition, we defined a new end point (maxRBA) to capture Ab response relative to baseline titer and were able to predict response in young and older adults separately using baseline transcriptional profiles. These signatures, such as the monocyte signature negatively predicting vaccine response in older adults, could be further developed and used to predict the subgroup of individuals who would most likely benefit from vaccination or, in contrast, the subgroup of predicted nonresponders to current vaccines who should be monitored more closely during flu season and enrolled in studies of investigational influenza vaccines. Baseline predictive signatures also open the possibility of modulating the immune system of individuals to mimic the immune state of predicted HR before vaccination or as part of a combination therapy with vaccination. This work represents a step toward better understanding of how the immune system responds to influenza vaccination in young and older adults and may be beneficial for rationally designing more effective vaccines.

We thank Dr. Randy Albrecht and Dr. Adolfo Garcia-Sastre at the Icahn School of Medicine at Mount Sinai, who led the Human Immunology Project Consortium core for influenza viral neutralization assays.

S.A. was employed by the Janssen Pharmaceutical Companies of Johnson & Johnson while writing this manuscript. The other authors have no financial conflicts of interest.