Gene signatures associated with adaptive humoral immunity following seasonal influenza A/H1N1 vaccination.

This study aimed to identify gene expression markers shared between both influenza hemagglutination inhibition (HAI) and virus-neutralization antibody (VNA) responses. We enrolled 158 older subjects who received the 2010-2011 trivalent inactivated influenza vaccine. Influenza-specific HAI and VNA titers and mRNA-sequencing were performed using blood samples obtained at Days 0, 3 and 28 post vaccination. For antibody response at Day 28 versus Day 0, several gene sets were identified as significant in predictive models for HAI (n=7) and VNA (n=35) responses. Five gene sets (comprising the genes MAZ, TTF, GSTM, RABGGTA, SMS, CA, IFNG and DOPEY) were in common for both HAI and VNA. For response at Day 28 versus Day 3, many gene sets were identified in predictive models for HAI (n=13) and VNA (n=41). Ten gene sets (comprising biologically related genes, such as MAN1B1, POLL, CEBPG, FOXP3, IL12A, TLR3, TLR7 and others) were shared between HAI and VNA. These identified gene sets demonstrated a high degree of network interactions and likelihood for functional relationships. Influenza-specific HAI and VNA responses demonstrated a remarkable degree of similarity. Although unique gene set signatures were identified for each humoral outcome, several gene sets were determined to be in common with both HAI and VNA response to influenza vaccine.

Introduction

Each year, seasonal influenza infection places a substantial burden on society, both medically and socioeconomically.[1] In the 2014–2015 season alone, nearly 18,000 hospitalizations and 5.0– 9.3% of all deaths were directly attributed to influenza infection and related complications, including pneumonia in the U.S.[2] Meanwhile, additional research suggests that medical expenses and loss of productivity associated with influenza infection place a financial burden of over $80 billion on the U.S. economy annually [1]. Yearly influenza vaccination is currently the most effective method available to reduce this morbidity and mortality.

Currently, the hemagglutination-inhibition (HAI) assay is the most widely accepted measure to assess the protection provided by influenza vaccine. An HAI titer of 1:40 is commonly accepted as a correlate of protection associated with ~50% protection against subsequent influenza virus challenge.[3-6] However, due to the inability of HAI titers to completely predict susceptibility to future infection, it is clear that HAI antibodies provide only a partial explanation of the immunity induced by influenza vaccination. For instance, a recent study noted the ability of patients with little or no HAI antibodies to resist influenza infection, suggesting that non-HAI neutralizing activity, such as broader virus neutralization antibodies (VNA), may contribute to protection.[7] Interestingly, HAI and VNA titers are often observed to be correlated, although the literature in this area is limited [8]. Another report suggests that VNA titers are a better correlate of protection against influenza infection than HAI titers after immunization with live attenuated influenza vaccine (LAIV).[9] Thus, it is important to consider both HAI and VNA titers when determining correlates of protection after influenza immunization.

In order to further explore this issue, it is critical to understand the genomic mechanisms underlying HAI and VNA responses to influenza vaccination. Currently, few biological or genetic markers[10-13] have been identified that characterize the adaptive immune response to influenza vaccine or predict vaccine failure, and additional research in this area is needed. In this study we utilized a systems biology approach to explore and identify such potential markers. Our aim for this study was to identify genomic markers (both individual genes and genesets) in common between both influenza A/H1N1-specific HAI and VNA responses following seasonal trivalent inactivated influenza vaccine (TIV) in older adults. Application of a systems biology approach will advance our knowledge by identifying novel mechanisms and generating hypotheses for variations in adaptive immune responses to vaccines.

RESULTS

Subject Demographics

In total, 159 subjects participated in this study,[14] and the median (IQR) age of the cohort was 59.5 (55.3, 66.3) years. Female participants represented 61.6% of the cohort while males represented 38.4%. Overall, a majority of the cohort self-identified as Caucasian (98.7%); the remaining 1.3% reported a race of “Other.” One subject’s sample failed mRNA-sequencing quality control; therefore, data from158 subjects were carried forward for analyses.

Antibody (Ab) Responses to Influenza A/H1N1 in Study Subjects

Influenza A/H1N1-specific Ab responses were measured in pre-vaccination serum samples (Day 0) and at Days 3 and 28 post-vaccination. Day 0 influenza-specific median HAI and VNA titers (median 1:80; IQR 1:40–1:320) demonstrated the presence of pre-existing antibodies in all subjects. Both HAI and VNA titers increased by Day 28 (1:320; IQR 1:160–1:640, p<0.001 for both outcomes). As expected, we found a strong positive linear relationship between HAI titer and VNA titer for Days 0 (r=0.924, p= 7.07x10−67), 3 (r=0.934, p= 1.65x10−71), and 28 (r=0.940, p=8.62x10−75), which suggests that Ab responses measured by the HAI and VNA assays are similar (Figure 1). The intra-class correlation coefficient (ICC) was 0.91 for the HAI assay (Day 0) and 0.83 for the VNA assay (Day 0). There were no significant differences in HAI (p=0.22, p=0.17 and p=0.79; for Day 0, 3, and 28, respectively) or VNA (p=0.37, p=0.31, and p= 0.41; for Day 0, 3, and 28, respectively) titers between male and female subjects at any timepoint (Figure 2A). Figure 2B illustrates that there were no significant differences at any timepoint in HAI (p=0.56, p=0.44 and p=0.15; for Day 0, 3, and 28, respectively) or VNA (p=0.86, p=0.51 and p=0.32; for Day 0, 3, and 28, respectively) responses.

Figure 1

Comparison of HAI and VNA responses (Day 28 vs Day 0) (n=158 subjects)

Heatmap of the overlap between the log2 fold-change of HAI (rows) and VNA (columns) for Day 28 relative to Day 0. The color scheme is determined by the percent of the total subjects in each cell with white indicating no overlap, pink indicating a small percentage of overlap, and blue indicating the largest overlap. The majority of the results either fall on the diagonal or off-diagonal, indicating strong concordance in these assays.

Figure 2

Distribution of HAI and VNA responses (Day 0, Day 3, and Day 28) by sex and age group

A) Scatterplot of HAI and VNA titers by male (M) and female (F) at Day 0, Day 3 and Day 28. B) Scatterplot of HAI and VNA titers by age: 50–64 years old, and 65 and older.

Common Genesets (Log2 Day 28 vs Day 0) Associated with HAI and VNA Responses

Externally defined genesets with significant mRNA expression changes from Day 0 to Day 28 (p < 0.005) were used in cross-validated penalized regression models to predict HAI or VNA change from Day 0 to Day 28 (adaptive immune response). We identified 7 and 35 genesets containing genes with changes in expression that demonstrated as association with variability in the response and had at least one gene entering the model for the HAI and VNA responses, respectively. In contrast, no genesets entered the models to predict early HAI or VNA change from Day 0 to Day 3. The genes associated with HAI are RNA transcription factors (TTF2, MAZ); chemokine/cytokine/receptor (CCR9, IL10RA, IFNG); transferase activity markers (GSTM1,2, SMS, RABGGTA ); carbonic-anhydrases (CA2,6,8,11,14); oxidoreductases (NDUFS3, CYB5R2,3, CYB561, NQO1,2); and the endosome transport dopey family marker (DOPEY2) (Supplemental Table 1). Among genesets associated with VNA response are TNF ligand TNFSF11; cytokines/receptors (IFNG, IL7, IL27); interferon-inducible transcription factors (IRF7, IRF9); and other genes with unknown role in the regulation of humoral immunity (Supplemental Table 2). Five genesets were in common with both HAI and VNA (Table 1). These genes/genesets were involved in cellular protein processes (MAZ, TTF2); transferase (GSTM1, GSTM2, RABGGTA, SMS) and carbonate dehydratase (CA2, CA6, CA8, CA11, CA14) activities; KEGG pathway (IFNG), and endosome transport (DOPEY2).

Table 1

Common Genesets with genes entering regression models for HAI and VNA Responses with the log2 Day 28 vs Day 0 fold-change in gene expression as the explanatory variables.

	HAI		VNA
Geneset [64]	MSE*	Genes	MSE*	Genes	Coefficient **	Median log2 Fold-Change
CPCD: (CELLULAR_PROTEIN_COMP LEX_DISASSEMBLY)	2.245	MAZ	2.290	MAZ	−	−0.005
				SUPT16H	−	0.026
		TTF2		TTF2	−	−0.009
				UPF1	+	0.037
ALKYL: (TRANSFERASE_ACTIVITY_T RANSFERRING_ALKYL_OR_ ARYLOTHER_THAN_METHY LGROUPS)	2.279	GSTM1	2.334	GSTM1	−	−0.029
		GSTM2		GSTM2	−	−0.039
				MAT2A	+	0.070
		RABGGTA		RABGGTA	+	0.004
		SMS		SMS	−	−0.005
CDA: (CARBONATE_DEHYDRATAS E_ACTIVITY)	2.286	CA11	2.346	CA11	+	−0.080
		CA14		CA14	−	−0.101
		CA2		CA2	−	−0.038
		CA6		CA6	−	−0.100
		CA8		CA8	−	0.197
RA: (KEGG_REGULATION_OF_AU TOPHAGY)	2.290	IFNG	2.403	IFNG	+	0.113
ET: (ENDOSOME_TRANSPORT)	2.330	DOPEY2	2.357	DOPEY2	+	0.023

Cross validated mean squared error (MSE). Genesets presented had genes remain in the penalized regression models for both HAI and VNA. The geneset name provides the abbreviation that is used for simplicity in the text, a brief description and geneset name from the MSigDB [64] and the actual gene.

A positive (+) coefficient from the regression models indicates that as the log2 fold change for the gene increases from Day 0 to Day 28 then the estimated response increases (upregulated with respect to the change). If the coefficient is negative (−) as the log2 fold change for the gene increases the estimate response decreases.

Common Genesets (Log2 Day 28 vs Day 3) Associated with HAI and VNA Responses

Utilizing the approach detailed above for genesets having statistically significant changes between Days 3 (innate immune response) and 28 post-influenza immunization, we identified 13 genesets associated with HAI titers. Additionally, 41 genesets were associated with variation in VNA response to vaccine. Genes in the HAI regression model include endoplasmic reticulum-associated enzyme (MAN1B1), DNA polymerase (POLL), and many other genes with unknown immune function (Supplemental Table 3). Genes in the VNA regression model include signal transducers (STAT1, STAT3); tyrosine kinase (TYK2); Golgi vesicular transport protein (GOSR1); and others. For additional genes and genesets, see Supplemental Table 4. Overall, 10 genesets and related genes were associated with both HAI and VNA responses. Included are immune-related genes MAN1B1, POLL, CEBPG, FOXP3, IL12A, TLR3, and TLR7 (Table 2).

Table 2

Common Genesets with genes entering regression models for HAI and VNA Responses, withlog2 Day 28 vs Day 3 fold-change in gene expression as the explanatory variables.

	HAI		VNA
Geneset [64]	MSE*	Genes	MSE*	Genes	Coefficient **	Median log2Fold- Change
ERCC: (REACTOME_N_GLYCAN_TR IMMING_IN_THE_ER_AND_C ALNEXIN_CALRETICULIN_C YCLE)	2.222	MAN1B1	2.333	MAN1B1	−	−0.022
NHEJ: (KEGG_NON_HOMOLOGOUS _END_JOINING)	2.288	POLL	2.374	POLL	−	−0.026
IFNG: (INTERFERON_GAMMA_PRO DUCTION)	2.302	CEBPG	2.419	CEBPG	−	−0.022
		EBI3			−	0.004
		FOXP3		FOXP3	−	−0.041
		IL12A		IL12A	+	−0.026
		INHBA		INHBA	+	0.005
		TLR3		TLR3	+	0.151
		TLR7		TLR7	+	0.058
		TLR9			+	0.147
BER: (KEGG_BASE_EXCISION_REP AIR)	2.302		2.313	LIG3	+	0.047
				NEIL1	+	−0.016
		PARP2		PARP2	−	−0.001
		POLL		POLL	−	−0.026
				SMUG1	+	0.015
FAB: (REACTOME_FATTY_ACYL_ COA_BIOSYNTHESIS)	2.306	ACSL1	2.381	ACSL1	+	0.063
		ACSL5		ACSL5	+	0.020
				ELOVL7	+	−0.012
		FASN		FASN	+	0.013
		HSD17B12			−	0.034
		SLC25A1		SLC25A1	+	0.009
SRA: (SECRETIN_LIKE_RECEPTOR _ACTIVITY)	2.309	BAI1	2.398	BAI1	+	−0.029
		TAPT1		TAPT1	−	−0.057
KREB: (BIOCARTA_KREB_PATHWA Y)	2.310		2.360	ACO2	+	0.017
		CS		CS	−	−0.002
				FH	−	0.016
				IDH2	−	0.002
		OGDH		OGDH	−	0.055
				SDHA	+	0.036
		SUCLA2		SUCLA2	+	0.043
PNS: (PERIPHERAL_NERVOUS_SY STEM_DEVELOPMENT)	2.318	BAI1	2.422	BAI1	+	−0.029
		GSTM3		GSTM3	−	−0.067
MCC: (MEIOTIC_CELL_CYCLE)	2.325	LIG3	2.363	LIG3	+	0.047
				RAD51B	+	−0.032
		RAD54B		RAD54B	−	−0.060
		REC8		REC8	+	−0.018
		STAG3		STAG3	−	−0.013
		TUBG1		TUBG1	−	0.030
LM: (LEARNING_AND_OR_MEMO RY)	2.326	ACSL4	2.354	ACSL4	+	0.015
		CRHBP		CRHBP	+	−0.025
		DLG4		DLG4	−	−0.011
				FYN	−	0.024
		GALR2		GALR2	−	0.059
				NF1	+	0.069
		S100B		S100B	−	0.030
				VLDLR	−	0.071

Cross validated mean squared error (MSE). Genesets presented had genes remain in the penalized regression models for both HAI and VNA. The geneset name provides the abbreviation that is used for simplicity in the text, a brief description and geneset name from the MSigDB [64] and the actual gene.

A positive (+) coefficient from the regression models indicates that as the log2 fold change for the gene increases from Day 3 to Day 28 then the estimated response increases (upregulated with respect to the change). If the coefficient is negative (−) as the log2 fold change for the gene increases the estimate response decreases.

Biologic Functions Implicated by Genes Associating with HAI and VNA Responses

Network interactions among all genes within significant genesets were extracted, revealing a number of highly interconnected modules. Our statistically prioritized genes were typically representatives from these modules (Supplementary Figure S2). Within this network, we calculated the extent of interconnectivity between genes prioritized by our statistical approach (Figure 3) and identified 51% as directly connected to at least one other prioritized gene. Randomly selected genesets (see Methods) exhibited a mean interconnectivity of 29 ± 4%, indicating a significant level of interrelationships among our prioritized genes; none of the randomly generated sets of genes showed this extent of interconnectivity. Focusing on statistically prioritized genes, many of those associated with both HAI and VNA (including INFG and TLR genes) were found to have known interactions with genes that were uniquely associated with HAI or VNA. Further, a high degree of network interactions were identified between the genesets prioritized at different timepoints, indicating a high likelihood for functional relationships between the genes. Thus, the distinct genes/genesets prioritized for each outcome and timepoint tend to directly interact with one another, indicating probable participation in common biologic mechanisms.

Figure 3

Statistically prioritized genes exhibit a high degree of network interactions

A) Comparing to randomly selected genesets, our prioritized genes have a significant level of direct interactions. B) Visualizing the full network of all genes within prioritized genesets reveals the presence of network modules (Supplementary Figure S2). C) The network interactions between our statistically prioritized genes demonstrate significant interactions across time points and outcomes. Edges are bundled to increase legibility.

To add further annotation-based evidence for functional relationships between these genes, we performed GO term enrichment. The most significantly enriched VNA-specific terms were purine salvage and metabolism. The most significantly enriched terms shared by both outcomes are the biosynthetic processes “long chain fatty-acyl-CoA,” “positive regulation of IFNG,” and “triglyceride.” While many of the statistically enriched terms matched prior expectation (e.g., IFNG regulation),[10] some present more novel hypotheses for future study. See Supplementary Figure S1 for the full list of significantly enriched GO terms, by outcome and timepoint.

DISCUSSION

The goal of this study was to identify genomic markers in common with both HAI and VNA responses following influenza A/H1N1 vaccination. Although the magnitude of HAI and VNA titers to influenza vaccine are often correlated, as they are in this study (Figure 1), the two measures of humoral immunity are different; therefore, it was expected that gene expression models would vary somewhat between the two immune outcomes. Thus, we chose to identify common genesets/genes based on similarities in HAI and VNA responses.

Through our modeling, we identified five genesets containing genes with expression changes from Day 0 to Day 28 that demonstrated association with both HAI and VNA titers (Table 1). Similarly, we identified 10 genesets containing genes with expression pattern changes from Day 3 to Day 28 (Table 2) that were associated with variation in both HAI and VNA responses. One of the genesets associated with both HAI and VNA for Day 28 vs Day 0 expression was CPCD. Genes of particular interest within this geneset include MAZ, a transcription factor, which was prioritized in both models and TTF2, which is a transcription termination factor for RNA polymerase II [15]. Often referred to as Pur1, MAZ encodes for a protein that regulates inflammation-induced expression of serum amyloid A (SAA) proteins.[16, 17] Data suggest that SAA is involved in at least two aspects of immune regulation. SAA proteins are released from hepatocytes during acute inflammation where they typically collaborate with high-density lipoprotein (HDL) [18]. This newly formed complex prompts the synthesis of several cytokines and acts as a chemotactic agent for both neutrophils and mast cells. Additionally, SAA is capable of binding to, and activating, TLR4 on B lymphocytes to initiate maturation.[18, 19]

We also identified the RA geneset, with significant change between Day 28 vs Day 0. This geneset includes the immune gene IFNG. Studies have shown that IFNγ is capable of inhibiting proliferation of pre-activated B lymphocytes or stimulating activated B lymphocyte proliferation and isotype switching.[20] Thus, MAZ or IFNG gene expression signature may contribute to differential influenza-specific Ab expression post-immunization.

Upon analysis of changes in gene expression between Day 28 vs Day 3 post- influenza vaccination, we identified 10 genesets associated with both HAI and VNA titers (Table 2). There was one gene in the ERCC geneset (MAN1B1) that entered the regression model. MAN1B1encodes the enzyme endoplasmic reticulum mannosyl-oligosaccharide 1,2-alpha-mannosidase (ER mannosidase I, Mnl1). MAN1B1 is commonly referenced in disorders of lysosomal storage [21] due to the enzyme’s role in cleaving mannose monomers from newly synthesized peptides within the ER. Because Mnl1 is an essential component of quality control for proper protein folding, secretion, and function,[22, 23] polymorphisms in MAN1B1 could potentially result in an alteration of the protein profile of immune response to TIV antigens. Recent studies utilizing kifunensine (Mnl1 inhibitor) suggest that Mnl1 affects immune regulation by preventing naïve CD4+ T cell activation,[24] leading to a lack of B cell co-stimulation and Ab production. Although this phenomenon has been more widely studied for MAN1A1-encoded Mnl1, the observed actions could potentially be generalized to MAN1B1-encoded Mnl1.

When evaluating the change in early-activation gene expression from Day 3 and Day 28 as predictors in the model, the NHEJ geneset had genes enter the regression model for both HAI and VNA. Within this geneset is the POLL gene, which codes for DNA polymerase λ (Pol λ), an enzyme with a critical role in both DNA replication and repair [25]. Previous studies suggest that Pol λ is involved in DNA base excision repair in the wake of oxidative damage.[26, 27] One study shows that Pol λ-deficient mice experience a decreased rate of germinal center B cell receptor somatic hypermutation (SHM).[27] This observed decline in SHM is thought to be caused by a reduction in the prevalence of B cells due to their inability to successfully repair oxidative DNA damage. It can be proposed that mutations in the POLL gene may affect the ability of germinal center B cells to diversify Ab populations, thereby altering the humoral immune response to influenza.

The immune-associated geneset IFNG also had genes that entered the model for both HAI and VNA when the change in gene expression for Day 28 vs Day 3 was evaluated. Genes of interest in this geneset are FOXP3, INHBA, TLR3, and TLR7; others include CEBPG, EBI3, IL12A and TLR9. However, TLR9 entered models only for the HAI, and not for the VNA response.

Studies have shown that FOXP3 serves as the master regulator of T regulatory cell (Treg) production. Tregs, which are critical for dampening the immune response, have also been shown to secrete cytokines, such as IFNγ, that have important immunoregulatory functions.[28] IFNγ is essential for the induction of FOXP3 in CD4+ T cells.[29, 30] Less research has been conducted on the influence of Tregs on humoral immunity; however, recent studies suggest that Tregs also inhibit B cell responses to antigens through IFNγ-related mechanisms. It has been suggested that alterations in the ratio of Tregs to T follicular helper cells within germinal centers exert control over Ab responses.[31]

Of further importance in the IFNG geneset is INHBA. This gene encodes the peptide inhibin beta A, a subunit of both activin and inhibin; these proteins operate antagonistically within the immune system.[32] Activin A is produced by activated B cells before it directly stimulates naïve B cells to increase production of IgG and indirectly stimulates activated B cells to increase IgG and IgE secretion.[33]

The Toll-like receptor (TLR) genes, TLR3 and TLR7, also entered the models for both the HAI and VNA response to influenza antigen models when using Day 28 vs Day 3 post-vaccination expression. TLRs are well known for their involvement in the innate immune response; however, it has been shown that B cells also express TLRs that, when bound to ligand, enhance cell survival and provoke increased Ab production.[34] Other genes in this geneset include CEBPG, a B cell transcription factor; [35] EBI3, a gene encoding a subunit of both IL-27 and IL-35;[36] and IL12A, a gene encoding a subunit of the Treg inducer IL-35.[37]

We have demonstrated that the genes prioritized by our statistical methods are involved in diverse cellular functions that are important biologically. Our statistical approach filters out genes with highly correlated expression (i.e., members of protein complexes or proteins mediating similar biologic functions). Thus, the somewhat diverse set of genes prioritized by our methods are likely to be representatives of broader classes of genes or of specific protein complexes (see Figure S2).

To interpret the coordinated activity of prioritized genes, we utilized two types of prior information: network resources and GO term enrichment. These two types of resources provide complementary information: evidence of physical or molecular interactions and participation in the same molecular or biochemical processes, respectively. Note that statistically significant GO term enrichment in this context is a descriptive measure about the gene’s functions: the genes were prioritized by our statistical approach and we use term enrichment as a means to describe what biologic activities to which they contribute. Through mapping all of the genes identified by our statistical approach onto biologic networks and inspecting the degree of known relationships among them, we find that many share direct interactions. The set of common genes, those prioritized by models of both HAI and VNA, have many direct connections with genes that were prioritized by one outcome, indicating participation in similar functions. This is also true across timepoints. In order to further identify the functions shared by these genes, we performed GO term enrichment for the entire set and also for each subset (by timepoint and outcome). Some functions are identified across the individual outcome comparisons, but some are only identified when considering genes across our outcome-associated results. Specifically, multiple lipid metabolism terms, including phosphatidylinositol (PI) signaling (including PIP3), lipid (cholesterol, triglyceride, PI) metabolism, and lipid transport, are implicated by genes across our comparisons. This broad network-based summary highlights the multifaceted nature of immune response.

Our study design utilized rigorous randomization approaches in order to minimize the influence of experimental artifacts. Our data reduction steps were agnostic to the outcome being modeled and utilized externally available immune system knowledge. All modeling involving the outcome was cross validated to ensure reproducibility. The novelty of this study is evident in that it is the first of its kind to identify genesets and genes whose changes in expression influence variation in both HAI and VNA responses to influenza A/H1N1 immunization. However, there are limitations to this study. First, the utilization of a larger cohort would enhance statistical power. The examination of gene expression underlying the innate immune response would be beneficial to this study; in order to do this, sample collection on Day 1 post-immunization would have been ideal. Lastly, the scope of this study is limited to A/H1N1 influenza and may or may not be generalizable to other components of the TIV.

The unique data and results generated from this study could potentially be used to predict humoral response to influenza vaccine in the future and inform the development of an individualized vaccine schedule paradigm.[38-40] The identification of gene signatures associated with humoral immunity may provide a better understanding of the genetic markers of immune response, and may assist with the design of better vaccines and adjuvants.

MATERIALS AND METHODS

Study Subjects

The methodology used for the selection and recruitment of study subjects has been previously reported elsewhere.[14, 41–43] Briefly, the study included 159 healthy adults, ranging in age from 50 to 74 years old, who were immunized with an intramascular single dose (0.5 mL) of the 2010–2011 seasonal TIV Fluarix (GlaxoSmithKline), containing A/California/7/2009 (H1N1), A/Perth/16/2009 (H3N2), and B/Brisbane/60/2008 viral strains.[14, 41–43] TIV is prepared from split-virion influenza viruses propagated in embryonated chicken eggs and contains a ratio of 15 ug of the hemagglutinin antigen for each of the three influenza strains.[44]

Prior to the onset of the study, each participant provided written, informed consent. All subjects reported stable health and provided detailed vaccination histories. Subjects were excluded from the study if they already received the 2010–2011 TIV. Exclusion criteria included the display of flu-like signs or symptoms at any point throughout the duration of the study. Subjects were also excluded from enrollment if they were diagnosed with influenza, or exhibited symptoms consistent with influenza, at any time from the beginning of the influenza season in Minnesota (as defined by the first reported cases over the period of the individual subject’s participation).[45] Additionally, one participant was excluded due to cDNA library preparation malfunction. Blood samples (90 ml) from each subject were obtained at three separate timepoints: pre-vaccination (Day 0), and Days 3 (innate immune response) and 28 (adaptive immune response) post-vaccination.[14] This study was approved by the Mayo Clinic Institutional Review Board.

Hemagglutination-Inhibition and Virus Neutralization Antibody Assays

We have previously described the HAI and VNA assays.[14, 41, 43] As published elsewhere,[46-48] the standard WHO protocol [49] was utilized to determine influenza H1N1-specific (virus strain A/California/07/2009) Ab titers from each subject’s serum at all three timepoints. The HAI titer was defined as the highest dilution of serum that inhibits red blood cell (RBC, 0.5%) hemagglutination. The VNA titer was defined as the reciprocal of the highest dilution of serum that neutralizes 200 plaque-forming units (PFU) of influenza A/H1N1 virus.

Seroconversion to influenza vaccine antigens was defined as a four-fold increase in serum Ab titers between Day 0 (prior to vaccination) and Day 28, or as an increase in serum Ab titers from <10 to ≥40 from Day 0 to Day 28.[50] The average coefficients of variation for the assays performed in this study were 2.9% and 4.7% for HAI and VNA, respectively.

Separation of peripheral blood mononuclear cells (PBMCs)

Methodology used for PBMC preparation is identical to what we have previously published.[42] PBMCs were isolated from blood samples (100 mL) from each subject at each timepoint pre- and post-vaccination using cell preparation tubes with sodium citrate (CPT™), as previously published.[42] Purified PBMCs were resuspended at a concentration of 1×107/mL in RPMI-1640 medium containing L-glutamine (Invitrogen), supplemented with 10% dimethyl sulfoxide (DMSO) and 20% fetal calf serum (Hyclone). Cells were frozen overnight at 80°C in freezing containers (Thermo Fisher Scientific) to reach an optimal rate of cooling and then transferred for storage to liquid nitrogen.

mRNA-Sequencing

To carry out mRNA next generation sequencing, we used protocols identical to those utilized and published for our transcriptomic studies.[51, 52] To summarize, the RNeasy Plus mini Kit and the RNAprotect reagent from Qiagen (Valencia, CA) were used to extract total RNA from 1x106 PBMCs. Manufacturer protocols were used to create full-length cDNA libraries using the mRNA-Seq 8 Sample Prep Kit by Illumina (San Diego, CA). The DNA 1000 Nano Chip kits were run on an Agilent 2100 Bioanalyzer (Agilent; Palo Alto, CA) for library validation and quantification, and cDNA libraries (5–7pM) were loaded onto individual flow cell lanes. The Illumin HiSeq 2000 (Illumina; San Diego, CA), in addition to Illumina’s Single Read Cluster Generation kit (v2) and 50 Cycle Illumina Sequencing Kit (v3), was used to perform single-end read sequencing, and the sequencing reads were aligned to the human genome build 37.1 using TopHat (1.3.3) and Bowtie (0.12.7). Gene counting was performed using HTSeq (0.5.3p3), and the reads mapping to individual exons were counted using BEDTools (2.7.1).[53-55]

Statistical Methods

Demographic results are presented as a percent of the total sample for discrete variables or as the median and interquartile range (IQR) for continuous variables. Spearman’s correlation was used to calculate associations between Ab titers. The Wilcoxon rank-sum test was used to test for difference in HAI and VNA titers between male and female subjects, and between subjects who were 50–64 years of age and those who were 65 and older.

Conditional quantile normalization, which adjusts for gene length and GC content, was used to calculate a normalization offset for use in subsequent analyses of the counts from the mRNA-sequencing data.[56] Moderated dispersion estimates were estimated using edgeR, assuming variance was a non-linear function of the mean, with the tagwise dispersion calculated from the trended dispersion.[57, 58] To identify genes that changed significantly over time, we utilized generalized linear models, assuming the negative binomial distribution with generalized estimating equations in order to account for correlation between multiple observations within a subject to obtain per gene p-values.[59-62] Self-contained geneset tests were performed using the gamma method with soft truncation threshold of 0.15, an extension of Fisher’s method of combining per-gene p-values.[63] The genesets tested were externally defined genesets from the Molecular Signatures Database [64, 65] that were in the Biocarta, GO, KEGG, Reactome, Sigma-Aldrich, Signaling Gateway Signal Transduction KE pathways or had “virus,” “infect,” “pathogen,” or “innate” as a keyword in the geneset description.

Penalized regression methods were used to build regression models with the goal of understanding the biological processes that may explain variation in vaccine response.[66] Two response variables were evaluated: 1) Log2 fold-change of the Day 28 HAI titer relative to the Day 0 HAI titer; and 2) log2 fold-change for Day 28 VNA titer relative to Day 0 VNA titer. Predictor variables were the individual genes from genesets having statistically significant changes over time (p< 0.005, n=339). Specifically, genesets with significant changes in Day 28 vs Day 0, and Day 28 vs Day 3 were used in the models (log2 difference in normalized expression between timepoints). Models were fit to genes from one geneset at a time. The predictor variables were first filtered using redundancy analysis, [67] with an R2 < 0.75 as the cutoff. The remaining genes were included as independent variables in the ‘glmnet’ function in R, which was used for model selection.[68] Gender was also included as a covariate in the models. Ten-fold cross validation was used to select the lambda parameter (based on the minimum cross validated mean squared error, or MSE), governing the selection of genes to be incorporated in the final model. The tuning parameter alpha was set to 0.9, reflective of the elastic net penalty (a combination of the L1 LASSO and L2 ridge penalties). The R statistical software version 3.0.2 was used for all analyses (www.r-project.org).

Genes prioritized by the above statistical models were used in network analysis. In order to be comprehensive and also focus on high-confidence interaction data, multiple network resources were combined including HPRD,[69] CCSB,[70] the Pathway Interaction Database,[71] and the subset (7.8%) of STRING,[72] where all interactions had a confidence score of at least 70%. We evaluated the significance level of network connections using random sampling. To do this, we generated a network interconnectivity metric defined as the fraction of genes that are first neighbors of each other, using all statistically significant genes that map to the network (n=166). We then randomly selected 10,000 sets of the same size from the network, limited to genes that were detected by RNA-Seq in our study (n=15,708). Network operations were performed using the igraph R package, version 0.7.1. Networks were visualized using Cytoscape [73] version 3.2.1 and layouts refined using AllegroLayout v.2.2.1.[74] The biologic functions of genesets were evaluated using gene ontology (GO) term enrichment.[75, 76] GO terms were extracted from the human Gene Ontology Annotation (GOA) database, [77] and hypergeometric tests were used to determine enrichment. We reported FDR-corrected p-values (q-values) for terms significant at the q < 0.01 level.