Innate gene signature distinguishes humoral versus cytotoxic responses to influenza vaccination

Because innate immunity plays an important role in adaptive immune response to vaccination, researchers have been exploring novel strategies to shape the adaptive immune response by targeting innate immune cells through the use of adjuvants, viral vectors, or novel immunization routes (1, 2). Specifically, adjuvant formulation and alternative routes of immunization have sought to increase the immunogenicity of influenza vaccination and improve immune protection, or even induce cellular immunity in immunocompromised individuals (3–5). Besides the standard intramuscular (i.m.) vaccination route, which provides moderate and variable protection against influenza mainly through humoral response, intradermal (i.d.) and transcutaneous (t.c.) routes have been proposed to induce cellular immunity (1, 6, 7). Clinical trials assessing t.c. methods, such as that developed by our group targeting vaccine at hair follicles, have revealed that they induce a differential immune response and have the benefit of inducing CD8⁺ T cell responses against influenza and melanoma (8–12). Topical application facilitates reservoir formation of vaccine compounds in hair follicle openings, where they reach high densities of skin dendritic cells (DCs) in activated states and are capable of promoting CD8 immunity (13).

The i.d. vaccination route induces both humoral and CD8⁺ T cell responses due to Langerhans cells (LCs) and dermal dendritic cell (dDC) populations (7, 14, 15). The skin’s rich population of immune cells includes LCs in the epidermis and dDCs, both potent antigen-presenting cells that help to process and present pathogen-derived antigens to mount an adaptive immune response (16). Moreover, the i.d. administration of influenza vaccines using microneedles offers several advantages: it is painless, elicits stronger immune responses, and might improve responses in high-risk populations (4, 17, 18). Furthermore, although most vaccine strategies have focused solely on humoral response to induce protection against influenza, CD8⁺ T cells endowed with cytotoxic activity have also proved to be effectors for protective immunity against influenza (19–22). It is generally admitted that humoral response is one correlate of protective response to influenza virus. However, there is also evidence that cytotoxic T lymphocyte (CTL) responses are a surrogate of protection (21, 23). McMichael et al. brought insight to the role of cytotoxic CD8⁺ cells with a unique study, showing that in the absence of antibody response to influenza, the presence of CTL activity correlates with reduced viral shedding (23). However, coordinated response between adaptive immune compartments (both humoral and cellular immune responses) is needed to induce an effective immune response and long-term protection after vaccination (24, 25). Conventional i.m. immunization fails to induce CD8 cytotoxic responses against seasonal influenza (11).

Recent advances in our understanding of the innate immune system and the use of systems biology approaches are beginning to reveal the fundamental mechanisms by which the innate immune system orchestrates protective immune responses to vaccination (26). Studying the molecular signatures that are induced rapidly after vaccination may help pinpoint biomarkers of immune responses for early identification of impaired responses to vaccines (27). Until now, this mechanistic knowledge has come from studies of vaccine-induced immune regulation in the context of research using systems biology and data integration. Notable studies looked at yellow fever, influenza, pneumococcal, and Ebola vaccination (28–31). These systems vaccinology approaches have identified molecular signatures that correlate with the intensity of protective immune responses. Whether early molecular signatures enable researchers to determine the quality of immune responses and responder and nonresponder status remains to be studied.

Here, we report exploratory results of a randomized phase I/II clinical trial conducted in 60 healthy subjects 18–45 years old to compare the immunogenicity of the 2012–2013 seasonal trivalent inactivated influenza vaccine (TIV) administered by the t.c., i.d., and i.m. routes. We propose an integrative approach using systems vaccinology to improve our understanding of the early innate immunity induced by the route of administration and to better study the quality of adaptive immunity, i.e., the antibody and CD8 cellular responses after TIV vaccination.

Study population. Randomized subjects received A/California/7/2009 (H1N1)pdm09–like strain, A/Victoria/361/2011 (H3N2)–like strain, and B/Wisconsin/1/2010–like strain for the 2012–2013 season by the allocated route of administration on d0 from October 15, 2012, to November 10, 2012, with follow-up until April 25, 2013 (Figure 1). One subject in the t.c. group was lost to follow-up at d21. Subjects in the 3 groups had similar demographics, and the sex ratio was around 60% to 40% women to men in all groups (Supplemental Table 1; supplemental material available online with this article; https://doi.org/10.1172/JCI125372DS1). The percentage of subjects who had previously received 1 or more (A)H1N1pdm09 pandemic influenza vaccine (P = 0.0177) or seasonal TIV (P = 0.0340) injection was significantly higher in the i.d. group (70% and 75%, respectively) than in either the t.c. (30% and 35%) or the i.m. group (30% and 40%). Blood and serum samples were collected from all subjects at visits on d0, d1, d21, and month 5 (M5), for the measurement of gene expression at early time points and humoral and cellular responses at d0, d21, and M5. Subjects used diary cards for 5 days after vaccination to report daily local reactions, systemic signs or symptoms, or influenza-like illnesses (ILIs). Neither any ILI nor any serious adverse event was reported during the study. The significantly higher local skin reactogenicity in individuals vaccinated i.d. compared with t.c. and i.m. confirmed previous reports of the safety of these modes of administration (Supplemental Table 2) (4, 10, 11, 32).

Figure 1

Flow diagram of 60 participants through the phases of the FLUWAY randomized phase I/II clinical trial. Flow of participants through the FLUWAY clinical trial, according to Consolidated Standard of Reporting Trials (CONSORT). Sixty participants were enrolled, randomized in 3 arms (1:1:1), to receive TIV vaccine by t.c., i.d., or i.m. routes on d0.

Different administration routes induce different qualities of influenza-specific humoral and T cell responses. In an exploratory analysis of our randomized clinical study, a multiparametric analysis was used to define the quality of immune responses after each route of TIV administration. Anti–influenza virus–specific humoral responses were measured by hemagglutinin inhibition (HI) and microneutralization (MN) antibody titers (Supplemental Tables 3 and 4, respectively). TIV-specific CD4⁺ and CD8⁺ T cells producing cytokines and/or Cd154 and cytotoxic granzyme B (GRZ) were also measured according to the gating strategy given in Supplemental Figure 1. Figure 2 summarizes the ratio at d21 relative to baseline (d0) for each immune parameter. Of note, baseline HI titers were similar in all subjects, regardless of whether they had or had not received (A)H1N1pdm09 vaccine or seasonal TIV before 2012 (Supplemental Table 1 and Supplemental Figure 2). Moreover, we didn’t find any differences in the TIV-specific immune response ratios (d21/d0) between individuals who previously received or did not receive at least 1 pandemic H1N1 or seasonal vaccination (Supplemental Figure 3). After TIV administration by the i.d. and i.m. routes, geometric mean titers of HI antibodies to all 3 influenza strains increased significantly (d21 seroprotection rates of 100% for A/H3N2, 80% for A/H1N1, and 35%–50% for influenza B; P ≤ 0.0001 for all 3 viral strains) (Supplemental Table 3). As expected, t.c. vaccination with TIV did not induce higher HI or MN (Figure 2, A and B) antibody titers against A/H1N1 and A/H3N2 strains than did i.d. and i.m. vaccination (11). Intensity of immune responses at d0 and d21 are shown in Supplemental Figure 4.

Figure 2

Influenza vaccine immunogenicity following t.c., i.d., and i.m. TIV vaccination. Adaptive immune responses before (d0) and after (d21) vaccination for each administration route: t.c. (blue, n = 19), i.d. (green, n = 20), and i.m. (black, n = 20). Data analysis presents d21/d0 ratios in box-and-whisker plots with the minimum to maximum showing all points (from the bottom up: the minimum 25th percentile Q1, median, 75th percentile Q3, and maximum values). Antibody titers specific for H1N1 (A) and H3N2 (B) virus by hemagglutination inhibition assay (left panels) and microneutralization (right panels). (C) TIV-specific CD8⁺ (left) and CD4⁺ (right) T cells expressing intracellular IL-2, IFN-γ, TNF-α, and/or cell-surface Cd154 (CD40 ligand), analyzed by flow cytometry. (D) Vaccine-specific CD8⁺GRZ⁺ T cells. Vaccine responders are defined with d21/d0 ratio of MN or HI antibody titers ≥ 4 and with d21/d0 ratio of TIV-specific CD4⁺, CD8⁺, or CD8⁺GRZ⁺ T lymphocytes ≥ 2. Therefore, dotted lines represent the threshold between responders (above) and non-responders (below) for each parameter. The Kruskal-Wallis test compared the 3 routes (*P < 0.05; ****P < 0.0001) and the X² test was used for the qualitative variables (^##P < 0.01).

We monitored TIV-specific CD4⁺ and CD8⁺ T cells expressing intracellular IL-2, IFN-γ, TNF-α, and/or cell-surface Cd154 (CD40 ligand) as well as TIV-specific CD8⁺ GRZ⁺ T cells (Figure 2, C and D). The d21/d0 ratios of the TIV-specific cytokine-producing CD8⁺ T cells were similar in all 3 arms and highly heterogeneous between individuals (Figure 2C) (Kruskal-Wallis test; P not significant). We found significant differences (P = 0.0369) among the 3 arms for TIV-specific cytokine-producing CD4⁺ T cells. This was confirmed by a χ² test applied to responders and nonresponders as defined by a d21/d0 ratio of at least 2 for vaccine-specific CD4⁺ T cells. We found significant differences (P = 0.0204) among the 3 routes; response rates were higher for i.d. (8/20) than i.m. (3/20) or t.c. (1/19).

Interestingly, the d21/d0 ratios of the TIV-specific CD8⁺GRZ⁺ T cells were similar in all 3 arms according to a Kruskal-Wallis test, although after application of the χ² test to responders and nonresponders we found a significantly higher response rate observed in the t.c. (7/18) and i.d. (6/20) groups compared with the i.m. group (0/20) (P = 0.0098; Figure 2D). Overall, the intensity of vaccine-specific humoral, cytokine-producing T cell and CD8⁺GRZ⁺ T cell responses differed by administration route, with higher humoral responses (HI and MN antibody titers) after i.m. and i.d vaccination and higher vaccine-specific CD8⁺GRZ⁺ T cells after t.c. and i.d. vaccination.

Early innate signatures discriminate immunological outcome. We used transcriptomic approaches to compare the quality of early innate immunity regardless of adaptive immunity outcomes. Comparison of the expression profiles of blood genes at 24 hours after vaccination (d1) and at baseline (d0) for each route of immunization showed 389 genes differentially expressed in the i.m. arm, 127 in the i.d. arm, and 4 in the t.c. arm (significance for the Bayesian t test was set at a corrected P < 0.1), as represented in a Venn diagram (Figure 3A). We noticed that 24 genes are common between the routes. Of note, comparison of gene expression at d0 (baseline) of all samples showed no significant difference among study arms (Supplemental Figure 5).

Figure 3

Early innate signature discriminates immunological outcome, rather than route of immunization. (A) Venn diagram of the genes differentially expressed between paired d1 and d0 for t.c. (blue, n = 17), i.d. (green, n = 18), and i.m. (black, n = 17). Significance for the Bayesian t test was set at a corrected P less than 0.1. The number of genes differentially expressed between d1 and d0 is indicated for each arm. (B) PCA of these 496 genes differentially expressed between d1 and d0. Vaccine administration routes are indicated. (C) Gene and sample hierarchical clustering is based on the expression profiles of these 496 genes. Two sample clusters named C1 (gray) and C2 (red) are observed. For C1, n = 28: t.c. n = 2, i.d. n = 16, i.m. n = 10. For C2, n = 24: t.c. n = 15, i.d. n = 2, i.m. n = 7. The color-gradient from green (–1, low) to red (1, high) indicates fold change d1/d0 gene expression (log²). (D) A major cluster of gene expression, with expression significantly different between C1 and C2, corresponds to the IFN-signaling pathway. Of note, IPA analysis indicated master genes (IRF7 and STAT1) at the center of the network. The red gradient indicates the differential expression between C1 and C2 (1.08- to 4.5-fold change gene expression C1/C2). And the box-and-whisker plot (from the bottom up: the minimum 25th percentile Q1, median, 75th percentile Q3, and maximum values) with minimum to maximum showing all sample points (y axis of graph in log¹⁰ scale). It represents the d1/d0 ratio of the CXCL10 serum level from the C1 (gray) and C2 (red) samples. Box-and-whisker plots represent (E) H1N1-specific MN antibody titers (d21/d0 ratios) and (F) TIV-specific CD8⁺GRZ⁺ T cells (d21/d0 ratios) according to C1 and C2 clusters (left panels). Respective histograms in E and F show immune response intensities for each individual in C1 (gray) and C2 (red) clusters (right panels). The Mann-Whitney t test compares C1 and C2 for these responses; *P < 0.05; **P < 0.01; ****P < 0.0001.

To study whether these 496 genes could differentiate immunization routes, we performed principal component analysis (PCA) of the data for the 496 genes from each individual (Figure 3B) and a hierarchical clustering (HCL) of the individuals (Figure 3C) (t.c., blue; i.d., green; and i.m., black). As the HCL analysis shows (Figure 3C), the samples were distributed into 2 major clusters: C1 (gray) and C2 (red). Of note, samples from individuals are distributed in both clusters. C1 includes 2/17, 16/18, and 10/17 t.c.-, i.d.-, and i.m.-vaccinated individuals, respectively. C2 includes the other 15/17, 2/18, and 7/17 t.c-, i.d.-, and i.m.-vaccinated individuals, respectively. Using Ingenuity pathway analysis (IPA) on 496 genes, we found that this C1-C2 dichotomy is due mainly to the genes involved in multiple pathways, such as those for antigen-presentation (IPA: P = 1.32 × 10^–5), DC maturation (IPA: P = 7.14 × 10^–4), and IFN signaling (IPA: P = 9.18 × 10^–15). Figure 3D displays the top pathway IFN signaling, which includes STAT1, IRF9, IFI35, IRF7, PSMB8, and FCGR1A from a cluster of significantly upregulated genes in C1 compared with C2, among the 496 genes. It has been previously shown that CXCL10 is also one of the key mediators of immune response (33). We thus measured CXCL10 production in subjects’ serum samples at d1 compared with d0. We found that d1/d0 ratios of CXCL10 levels were significantly higher in the C1 than in the C2 cluster (P < 0.0001) (Figure 3D), CXCL10 serum levels being correlated with influenza-specific MN antibody titers (strains H1N1: P = 0.0002, r = 0.49; and H3N2: P < 0.0001, r = 0.61).

Based on the 496 genes differentially expressed at d1 and the distribution of samples according to the HCL, we analyzed adaptive immunity (d21) in C1 and C2 clusters. We observed significant differences in the immune responses (d21/d0 ratio) between C1 and C2 (Figure 3, E and F). C1 individuals had significantly higher influenza-specific MN antibody titers (H1N1: P = 0.0013, Figure 3E; and H3N2: P = 0.0004, Supplemental Figure 6A) than C2 subjects. Similar data were observed for HI Ab titers in C1 and C2 (data not shown). In contrast, TIV-specific CD8⁺GRZ⁺ T cell responders were more frequent in the C2 than the C1 group (P = 0.0406) (Figure 3F). However, we did not observe any significant difference in distribution of TIV-specific CD4⁺ and CD8⁺ T cells expressing intracellular IL-2, IFN-γ, TNF-α, and/or cell-surface Cd154 (CD40 ligand) in C1 and C2 (Supplemental Figure 6, B and C). This was confirmed by the immunological assays performed at d0 and d21, according to C1 and C2 sample clusters (Supplemental Figure 7). In conclusion, the expression profile of the whole-blood genes (d1/d0) allowed us to distinguish H1N1- and H3N2-specific antibody responses from TIV-specific CD8⁺GRZ⁺ T cell responses.

Early innate gene signature and serum levels of cytokines allows identification of humoral and cytotoxic T cell responders. Among the 496 differentially expressed genes (d1/d0), 80 were correlated with MN antibody titers (P < 0.05, 0.3, less than r < 0.6) mainly involved in interferon signaling (IPA: P = 9.18 × 10^–5) and antigen presentation pathways (IPA: P = 4.11 × 10^–4) and 31 genes with TIV-specific CD8⁺GRZ T cell responses (P < 0.05, 0.3, less than r < 0.4) involved in metabolic pathways (IPA: P = 5.88 × 10^–3). We then looked at whether these genes were differentially expressed in responders and nonresponders, as defined in Figure 2. The top 5 genes with significantly different gene expression levels between humoral responders and nonresponders (MN antibody titers) are presented in Figure 4A and Supplemental Table 5, as are the top 4 genes with significantly different gene expression levels among subjects with and without TIV-specific CD8⁺GRZ⁺ T cell responses.

Figure 4

Minimal innate signature correlates with immunological outcome independently of route of immunization. (A) Box-and-whisker plots (minimum to maximum) are depicted for 9 gene expression ratios between d1 and d0 in nonresponders (NR; gray) and responders (R; brown) (d21/d0 ratio ≥ 4 for A/H3N2 or A/H1N1 MN titers [left panel] and ratio ≥ 2 for TIV-specific CD8⁺GRZ⁺ T cell responses [right panel]). The Mann-Whitney t test compares responders and nonresponders (*P < 0.05; **P < 0.01; ***P < 0.001). (B) Gene and sample hierarchical clustering of 9 gene expression profiles revealed 2 clusters of subjects: S1 (purple, n = 33), with 3/17 t.c., 13/18 i.d., and all i.m. subjects (n = 17/17), and S2 (yellow, n = 19), with 14/17 t.c. and 5/18 i.d. subjects. Samples are color-coded according to vaccination route: t.c. (blue), i.d. (green), and i.m. (black). The color gradient from green (–0.4, low) to red (0.4, high) indicates the intensity of gene expression. Box-and-whisker plots with the minimum to maximum points of d21/d0 immune responses (from the bottom up: the minimum 25th percentile Q1, median, 75th percentile Q3, and maximum values) for (C) A/H1N1 MN titers, (D) A/H3N2 MN titers, and (E) TIV-specific CD8⁺GRZ⁺ T cell responses in the S1 (purple) and S2 (yellow) clusters (left panels). The Mann-Whitney t test compared the S1 and S2 clusters (***P < 0.001). Respective histograms show immune response intensities for each individual in the S1 (purple) and S2 (yellow) clusters (right panels). t.c., blue; i.d., green; and i.m., black.

Hierarchical clustering based on the expression profile of these 9 genes revealed 2 clusters of subjects: S1 (purple) including 3/17 (9%) vaccinated by the t.c. route, 13/18 (39%) by i.d., and 17/17 (52%) by i.m., and S2 (yellow), which included 14/17 (74%) vaccinated by t.c. and 5/18 (26%) by i.d. (Figure 4B). We observed significantly higher influenza-specific MN antibody responses in S1 than in S2 (A/H3N2: P = 0.0002; A/H1N1: P = 0.0005) and the same distribution is observed for the d1/d0 ratios of CXCL10 concentrations (Figure 4, C and D and Supplemental Figure 8). In contrast, the S2 cluster contains significantly higher numbers of individuals with higher TIV-specific CD8⁺GRZ⁺ T lymphocyte counts (P = 0.0007) (Figure 4E). This is confirmed by the immunological assays at d0 and d21 independently according to S1 and S2 samples clusters (Supplemental Figure 9). Moreover, it has been shown that i.d. vaccination induces both TIV-specific MN titers and CD8⁺GRZ⁺ T cell responses, whereas t.c. vaccination preferentially induces CD8⁺GRZ⁺ T cell responses. However, we did not observe any significant difference in distribution of TIV-specific CD4⁺ and CD8⁺ T cells expressing intracellular IL-2, IFN-γ, TNF-α, and/or cell-surface Cd154 (CD40 ligand) in S1 and S2 (Supplemental Figure 8, A and B).

PCA performed on 9 genes (genes correlated with immune responses and were differentially expressed between responders and nonresponders) showed sample distribution according to the administration route (t.c., blue; i.d., green; i.m., black; Figure 5A).

Figure 5

Minimal innate signature and serum level of CXCL10 allows identification of humoral and cytotoxic T cell responders after TIV vaccination by different routes of administration. (A) Principal component analysis of the 9 gene expression profiles (d1/d0). Sample vaccine administration routes are color coded: i.d. in green (n = 18), i.m. in black (n = 17), and t.c. in blue (n = 17). (B) ROC curves showing the specificity and sensitivity of the logistic regression models, i.e., the proportion of correctly anticipated responders and nonresponders, respectively. The curves correspond to the d21/d0 ratio of immune responses: H1N1-specific MN antibody titers (blue; AUC: 0.9676), H3N2-specific MN antibody titers (green; AUC: 0.9585), and TIV-specific CD8⁺GRZ⁺ T cells (red; AUC: 0.9354). The 9 gene expression profiles and the CXCL10 serum levels significantly explain each of these variables. (C) The histogram represents the number of responders for each immune response: H1N1/H3N2 MN antibodies (d21/d0 ratio ≥ 4) and TIV-specific CD8⁺GRZ⁺ T lymphocytes (d21/d0 ratio ≥ 2), in regard to provided status (black) and correctly predicted status (gray) based on the early expression of the 9 genes and on CXCL10 serum levels (d1/d0).

We then used logistic regression analysis to consider whether this minimal gene signature, together with the CXCL10 serum levels, might supposed response by either TIV-specific CD8⁺GRZ⁺ T lymphocytes or MN antibody titers. The receiver operating characteristics (ROC) curve showed that the 9 genes and the d1/d0 ratios for CXCL10 serum levels were the most accurate prognosticators of responses by TIV-specific CD8⁺GRZ⁺ T cells (AUC: 0.9354) and by H1N1- and H3N2-specific MN antibody titers (AUC: 0.9676 and 0.9585, respectively) (Figure 5B). Figure 5C shows the number and percentage of individuals correctly foreseen to be responders by their status for each immune parameter. The percentage of individuals predicted as responders for H1N1- and H3N2-specific MN antibody and TIV-specific CD8⁺GRZ⁺ T cell responses is 100%, 84.21%, and 83.33%, respectively.

We believe this work is the first phase I/II clinical trial that aims to assess the safety and immunogenicity of TIV administered by t.c., i.d., and i.m. routes, and provides important insights into the dichotomy shown by early innate markers for the induction of either TIV-specific humoral or CD8 responses.

We found that TIV immunogenicity varied according to delivery routes. Stronger humoral responses (HI and MN antibody titers) were observed with i.m. and i.d. vaccination compared with t.c. vaccination, and vaccine-specific cytotoxic CD8⁺GRZ⁺ T lymphocyte responses were higher after t.c. and i.d. vaccination, compared with i.m. vaccination. These results confirm our previous demonstration that a t.c. route, preferentially targeting LCs, induces CD8⁺ T cell responses (6, 7, 11, 13). It also confirms that the i.d. vaccination route induces both humoral and CD8⁺ T cell responses due to the presence of LC and dDC populations. These cells include LCs in the epidermis and dDCs in the dermal layer, both involved in processing and presenting pathogen-derived antigens to T and B cells from draining lymph nodes to initiate the antigen-specific cellular and humoral responses. Influenza vaccine by the i.d. route (a) elicits stronger immune responses in preclinical studies than those induced by the i.m. route (34–36), (b) may improve response in high-risk populations such as the elderly (37–40), and (c) requires a lower vaccine dose than the i.m. route (17, 41). We proposed that immunization routes dictate the immunological outcomes of either humoral or CD8⁺ cellular responses in mice and humans (6, 11). In our study, we also took into account the impact of previous vaccination status on innate and adaptive immunity. We did not find any significant difference at baseline in both transcriptomic profile and adaptive immunity in previously vaccinated individuals compared with individuals who were not vaccinated in previous years.

We analyzed whole-blood gene expression at d1 relative to d0. Although previous work by others has studied later time points (D3–D7) (29, 33), we chose our design to study early innate immunity only with samples collected at baseline and d1 after vaccination. Because the route of immunization dictates innate immune response, we expected approximately 3 gene expression clusters representing the 3 modes of administration. Surprisingly, however, we observed that subjects were distributed into 2 major clusters (C1 and C2). This dichotomy was mainly due to the IFN-signaling pathway and CXCL10 production in serum 1 day after vaccination. C1 and C2 each displayed 1 of the 2 major immunological outcomes: humoral (influenza-specific MN antibodies) and CD8⁺GRZ⁺ T cells, respectively, at d21 relative to d0. Analysis of HI antibody titers produced similar results (data not shown). Overexpression of the IFN-signaling pathway (at D3 after vaccination) has previously been correlated with a higher TIV-specific HI antibody response (28, 33). Similarly, we found that the IFN-signaling pathway (at d1) and the CXCL10 serum level characterized individuals with high humoral responses at d21. Other studies have identified this chemokine, known to mediate antiviral immunity, in serum 3 days after influenza immunization by the i.m. route (where it was negatively correlated with HI antibody response 28 days afterward) (33) and 3 days after recombinant Vesicular Stomatitis Virus-Zaire Ebola virus (rVSV-ZEBOV) vaccination administration, where it was finally found to be an independent correlate of antibody response (31). On the other hand, early upregulation of IFN-signaling and antigen-presentation pathways after i.m. TIV immunization is known to be associated with higher influenza-specific HI antibody response (28, 33) and to result in overexpression of STAT1, IRF9, IFI35, IRF7, PSMB8, and FCGR1A, as found in our analysis for both the i.d. and i.m. routes.

Systems vaccinology is helping to increase knowledge about the molecular networks that orchestrate immunity in response to vaccination in humans (42, 43) and to identify signatures used to predict vaccine-induced immunity, such as humoral and cellular responses (44, 45). Studies of the innate blood transcriptome have helped to predict the magnitude of later antibody and CD8⁺ T cell responses against the live-attenuated yellow fever vaccine (30, 43, 46). Predicting CD8⁺ T cell responses against influenza vaccines is more challenging because of the poor induction of these immune responses by i.m. immunization. In circumstances when humoral immunity is not protective, or is inadequately protective, against the influenza strain causing an outbreak, cytotoxic CD8⁺ T cells might provide some degree of protection against infection and attenuate disease severity and mortality (21). This is especially the case for the majority of CD8⁺ T lymphocytes directed to epitopes that are cross-reactive between viral subtypes (47, 48). The t.c. route of immunization gave us the opportunity to investigate the induction of CD8⁺ responses without humoral responses. Among the 496 genes differentially expressed (on d1 compared with d0), we found 5 genes correlated with A/H1N1 and A/H3N2 MN antibody titers and 4 with TIV-specific CD8⁺GRZ⁺ T cell responses. They constitute so-called minimal gene signature. We found an innate gene signature induced 1 day after TIV vaccination that distinguished patients with late MN antibody responses from those with CD8⁺GRZ⁺ T lymphocyte responses, regardless of delivery route. In our study, CXCR2P1, C2, and CKS1B, induced on d1 early after TIV immunization, were positively correlated with A/H1N1 MN response at d21, whereas PRKAA1 and TMEM8B were negatively correlated with it. Among the genes correlated to A/H1N1 MN antibodies, CKS1B is involved in cell cycle and proliferation, JAK/STAT3 signaling (through its role in cytokine receptor signaling and inflammation), and the MEK/ERK pathway (where CKS1B regulates monocyte autophagy) (49). C2 is involved in the complement system; it triggers inflammation and plays a central role in promoting humoral immune response through antigen uptake and presentation to B and T cells (50). CXCR2P1 is the pseudogene of the IL-8 receptor and is mainly involved in neutrophil migration (51). Protein kinase PRKAA1, the catalytic subunit of AMPK, is expressed principally in macrophages and lymphocytes and plays a pivotal role in inducing proinflammatory signals through TAK1 and NFkB activation (52). AMPK promotes inflammation, innate immunity, and antiviral defenses through modulation of the host protein, a stimulator of IFN gene signaling (STING) (53). TMEM8B is involved in cell signal transduction and controls cell-cycle progression (54).

Our study also identified the following gene signature related to CD8⁺GRZ⁺ T cell response: MAP2K5, PVRL1, SARM1, and CXCR4. Interestingly, MAP2K5 is a cofactor required for innate immune responses dependent on TLR7 and TLR9, which are toll-like receptors of DCs, macrophages, and natural killer cells (55). CXCR4 is a chemokine receptor involved in migration and proliferation of white blood cells and plays a chemotactic role for neutrophils, monocytes, and lymphocytes (56, 57). Of note, PVRL1 is involved in cell movement and proliferation and is expressed on DCs. SARM1 is a TLR adaptor carrying a TIR (toll interleukin receptor) domain to ensure efficient immune responses against bacterial infections by negatively regulating TLR signaling (58). Another study has also proposed biomarkers related to apoptotic pathways that reflect impaired TIV responses in young compared with older individuals; however, these biomarkers are related to the individuals’ chronological ages (59).

Franco et al. identified 20 genes exhibiting a transcriptional response to vaccination, and significant genotype effects on gene expression, which are correlated to the intensity of antibody responses (60). A number of loci that have the strongest evidence of genetic variation might also influence the intensity of humoral response to vaccine. These results suggest that genetic variation in identified genes could also influence immunological outcomes.

Several inflammatory cytokines were measured in the sera at d1. We found slightly but nonetheless significantly more serum level of IL-6 production following vaccination (defined by a d1/d0 ratio ≥ 2) in S2 (7/19) compared with S1 samples (4/33) (χ², P = 0.0356) (Supplemental Figure 10). Of note, the adjusted ROC curve showed that the 9 genes and the d1/d0 ratios for IL-6 and CXCL10 serum levels were the most accurate predictors of responses by TIV-specific CD8⁺GRZ⁺ T cells (AUC: 0.9417, predicted responders 91.66%) (Supplemental Figure 10). IL-6 production by DCs enhances CD8⁺ T cell proliferation, triggers TCR activation, and is dispensable for memory CD8⁺ T cell generation in mice (61). In addition, the interaction of keratinocytes, LCs, and dDCs after cutaneous vaccination tends to promote inflammatory cytokine production, such as IL1-α, IL-6, IL-8, and CXCL10, that might contribute to the overall inflammatory reaction and further promote the attraction of inflammatory cells at the immunization site (62, 63). Moreover, vaccine adjuvant monophosphoryl lipid A (MPLA) enables effective CD8⁺ T cell differentiation into cytotoxic cells in the presence of DC-induced IL-6. MPLA and alum with an epitope from a nucleoprotein highly conserved between influenza strains generate cytotoxic T lymphocyte differentiation and long-lived memory T cells that protect mice against influenza A (64). Additional studies are necessary to determine whether the genes identified have a causal relationship to either the antibody response or the cell-mediated immune response to the vaccine.

In conclusion, we identified potential blood biomarkers at the level of gene expression and CXCL10 in serum that are good indicators of either humoral or cytotoxic T cell response after seasonal influenza vaccination. It is clear that further studies and validation in clinical trials are necessary to identify both a larger group of biomarkers induced after vaccination and predictors of the quality of immune responses. These studies must consider the variability, sensitivity, and specificity of these biomarkers in humans. Unfortunately, we did not have access to such studies comparing routes of immunization since our proposed work is the only clinical trial comparing 3 routes of administration. Although the transformation of data into predictive markers of immunity remains a challenge (65, 66), this study provides what we believe is new insight into the impact of the route of immunization and innate signature in the quality of immune responses.

Clinical study protocol. The FLUWAY (Routes of Immunization and Flu Immune Responses) randomized phase I/II clinical study enrolled 60 volunteers aged 18–45 years from October 15, 2012, to November 10, 2012, with follow-up until April 25, 2013; 20 were allocated to each vaccine administration group in a 1:1:1 ratio after randomization at the Cochin-Pasteur Center for Clinical Investigation in Vaccinology at Cochin Hospital, Paris, France. The study has been registered at ClinicalTrials.gov (NCT01707602). Exclusion criteria were as follows: pregnancy (positive urine test); HIV, hepatitis B, or hepatitis C virus infections; Guillain Barré syndrome; immunosuppressive treatment or other immunodeficiency; allergy to a vaccine component; medical history of skin cancer; acute infection or vaccination within 4 weeks of enrollment; body mass index lower than 21 or greater than 26; skin phototype V–VI; excessive terminal hair growth; and planned sun exposure 6 weeks before or during the study.

The number of subjects to be included in the research is not based on a statistical hypothesis but on the feasibility of carrying out this pilot study for systems vaccinology. Subjects were assigned to groups centrally through a balanced-block randomization process (random block size of 3, 6, or 9). The allocation schedule was computer generated in SAS 9.3 by an independent statistician who was not otherwise involved in planning or analysis. All subjects received their allocated vaccine on d0 (Figure 1) by the allocated route of administration to assess its immunogenicity and safety against A/California/7/2009 (H1N1)pdm09-like strain, A/Victoria/361/2011 (H3N2)-like strain, and B/Wisconsin/1/2010-like strain for the 2012–2013 season. The i.m. group received an injection of Vaxigrip (Sanofi Pasteur). The i.d. group received Intanza15 (Sanofi Pasteur) delivered via the BD Soluvia microinjector device. The t.c. group received Intanza15 by a t.c. needle-free method targeting hair follicles, according to the previously described standard operating procedure (10). All vaccines were administered in the deltoid region of the nondominant arm. Briefly, the t.c. vaccination area (4×4 cm) was gently shaved, covered with a thin layer of cyanoacrylate (Superglue, UHU GmbH & Co. KG), and then stripped with adhesive tape to open hair follicles. Vaccine was applied to the stripped skin inside a silicone barrier to limit the spread of the liquid, allowed to dry for 20 minutes, and then covered with a Comfeel adhesive bandage (Coloplast) for 24 hours. To measure humoral and cellular responses, blood and serum samples were collected from all subjects at visits on d0, d1, and d21.

The primary outcome was the CD8⁺ T cell response against the specific vaccine strain, measured at baseline and d21 after vaccination by 3 routes of immunization. Secondary outcomes were safety at each visit, measurement of microneutralizing antibodies and cellular responses at d0, d21, and M5, and transcriptomic analysis at d0 and d1. Immunomonitoring was performed blinded. Subjects used diary cards for 5 days after vaccination to report daily local reactions, systemic signs or symptoms, or ILIs. Data were collected in the electronic case report form. Neither any ILI nor any serious adverse event was reported during the study.

Hemagglutination inhibition and microneutralization assays. Serum antibodies against the 3 influenza viral strains contained in the 2012–2013 influenza vaccine were measured at d0, d21, and M5 with microtiter HI and MN assays. The method is described in Supplemental Materials and Methods.

Transcriptomic analysis. For transcriptomic analysis, 2.5 ml whole blood was collected in PAXgene tubes (Qiagen) from each volunteer before and 1 day after vaccination. RNA was amplified with the Ovation Pico WTA system version 2 (NuGEN) and the resulting cDNA hybridized on Illumina whole human genome oligo microarrays (HT12 version 4.0). Genome Studio (Illumina) was used to analyze the data. Nonspecific background signal was filtered; quantile normalization of samples and log² transformation were applied to the gene raw data. In all, 47,323 probes were analyzed.

Transcriptomic data were analyzed with R software 3.4.1 (http://www.r-project.org/) and the limma package. Raw data were analyzed according to a linear model fit with a moderated t test by empirical Bayes statistics. Adjusted P values for simple multiple testing procedures were generated including Benjamini & Yekutieli (2001) and Benjamini & Hochberg (2006) procedures for the control of the false discovery rate (FDR). Corrected P values less than 0.1 were considered to indicate statistical significance from transcriptomic data analyses for the 3 administration routes, and comparisons between d0 and d1 (d1–d0). Statistical analyses and graphic representations were performed with Prism 6.0 (GraphPad Software Inc.). Comparisons between routes or between identified clusters used an unpaired Mann-Whitney t test. All tests were 2-sided. P values less than 0.05 were used to identify genes with significant differential expression among vaccine responders (for the d21/d0 ratio ≥ 4 for HI and MN antibody titers and the d21/d0 ratio ≥ 2 for CD8⁺GRZ⁺ T lymphocytes) and nonresponders. Correlations were identified with the Spearman coefficient, with significance defined by P values less than 0.05. Heatmaps and hierarchical clusters were generated with Multi Experiment Viewer MeV 4.9.0 (which is part of TM4 microarray software suite developed at The Institute for Genomic Research) (67). The Pearson coefficient of correlation with the complete linkage method generated hierarchical clusters, represented with dendograms, and R software was used to perform principal component analyses (PCA). Ingenuity pathway analysis (IPA) was used to perform functional enrichment analyses and identify new targets or candidate biomarkers within the context of biological systems. It provided the canonical pathways, molecular/cellular functions, and networks that were statistically overrepresented in the gene signatures.

Statistics. Baseline characteristics of the study population were expressed as percentages for qualitative variables and as medians (with interquartile intervals) for continuous variables. Geometric means and 95% CIs were calculated for HI and MN antibody titers. Intention-to-treat analyses were performed and included all randomized subjects. Characteristics were compared between the 3 groups with χ² or Fisher’s exact tests for qualitative variables and with Kruskal-Wallis or Mann-Whitney tests for quantitative variables. Immunological data analyses used samples from 60 subjects. Appropriate comparative statistics are stated in the text. Data were analyzed with SAS 9.3 and Graphpad Prism software version 6.0, with a 2-sided P value of less than 0.05 considered statistically significant. Analysis of the difference between week 0 and week 14 within each individual group used a paired nonparameterized t test (Wilcoxon matched-pairs signed rank test), with statistical significance set at P less than 0.05. The logistic regression analyses and ROC curves were performed and generated with R.

Data availability statement. All data analyzed in this study are included in the published article. The normalized microarray data that support the finding of this study have been deposited in ArrayExpress with the accession code E-MTAB-7741.

Study approval. The FLUWAY original studies were conducted in accordance with the Declaration of Helsinki and the International Conference on Harmonization Good Clinical Practice guidelines and approved by relevant regulatory and independent ethics committees (Comité de Protection des Personnes, CPP Ile de France III Paris). Each participant provided written informed consent before study entry.