Comparisons of gene expression signatures in breast milk lymphocyte cells and PBMCs

Prior to differential gene analysis, we inspected global gene expression patterns to detect any strong outlying samples stratified by specimen type. The Day 0 sample for one subject (Subject B) in the TIV vaccine group showed a globally outlying gene expression pattern in BML samples (Supplemental Figures S21–S23). This sample was excluded from downstream analyses. Next, to assess differences in immune responses on the gene expression level, we first determined genes that were differentially expressed relative to pre-vaccination (FDR-adjusted p-value < 0.2 and a fold change of ≥1.5-fold) for each vaccine group, specimen type, and time point (Fig. 1, Supplemental Tables S10–S15). In BML, we identified a total of 382 DE genes in the LAIV group, most of which were up-regulated at Day 28 (Fig. 1A). For TIV, three DE genes were identified at two time points in BML, of which two up-regulated genes (IL1A and IL1B) overlapped with LAIV at Day 28. These genes, encoding for the Interleukin-1 alpha/beta cytokines, were strongly differentially up-regulated from pre-vaccination following administration of LAIV (10.9 and 18.9-fold, respectively) and TIV (11.8 and 9.2-fold, respectively). At Day 28, MTRNR2L8, SIGLEC14, CCL4, ST18, and RNVU1-19 showed the strongest differential responses between vaccine groups with 2.8-to-5-fold higher fold change responses in the LAIV group relative to the TIV group.

In panels (A, B), the Venn diagrams summarize differentially expressed (DE) genes for each time point post-vaccination versus pre-vaccination, the overall number of DE genes (left side), and both up-regulated (red font) and down-regulated (blue font) genes (right side) for each vaccine group. In (C, D), the heatmap summarizes subject-level log₂ fold changes for DE genes separately for each specimen type. DE differentially expressed, LAIV live attenuated influenza vaccine, TIV trivalent influenza vaccine, BML breast milk lymphocytes, PBMC peripheral blood mononuclear cells.

In PBMCs, at Day 28, 10 (2 up and 8 down-regulated) DE genes were identified for TIV (Fig. 1B) while one up-regulated gene was identified for LAIV but not TIV (MTRNR2L8) (Fig. 1A). This gene was also identified in BML LAIV at Day 28 (Fig. 1A). Up-regulation of MTRNR2L8 was 17.1-fold for the LAIV group and 1.04-fold for the TIV group compared to pre-vaccination in PBMCs as well as 5-fold higher in BML at Day 28 for LAIV (14.1-fold) vs. TIV (2.8-fold). The two up-regulated genes for TIV encoded for immunoglobulin proteins (IGHG1 and IGHG3). Both encoded for immunoglobulin (Ig gamma-1 and gamma-3 chain C region) proteins with 6.7-fold and 2.2-fold higher gene expression responses induced by TIV compared to LAIV vaccination, respectively.

Heatmap analysis of BML DE gene fold changes showed some initial up-regulation from baseline for DE genes at Day 2 and the strongest signal at Day 28 for the LAIV group with the majority of subjects showing increased responses for most of the DE genes (Fig. 1B). There was evidence of a weak Day 2 vaccine group effect, with LAIV subjects tending to have a higher fraction of subjects with up-regulated genes compared to subjects in the TIV group (Fig. 1B). At Day 8, a TIV vaccine effect was observed with 4 of 6 subjects showing a consistent upregulation of most DE genes (Fig. 1B). By Day 28, most subjects showed an up-regulation trend for most of the DE genes, but this pattern was stronger for LAIV compared to TIV with the highest consistent up-regulation identified in LAIV BML samples (Fig. 1B). In contrast, the Day 28 signal in PBMCs did not indicate strong overall effects for either vaccine (Fig. 1C).

Breast milk lymphocyte cells showed increased expression of interferon-inducible antiviral genes following LAIV vaccination

Next, to identify DE genes with similar responses following vaccination, we executed gene clustering analysis using log₂ fold changes of DE genes as input (Fig. 2, Supplementary Tables S17–S20 and Figure S12). Inspection of mean log₂ fold change time trends (Day 2, 8, 28) for 36 co-expressed DE gene clusters in BML showed that TIV generally induced responses ahead of LAIV (Day 8 vs. Day 28), albeit not as strongly, i.e., the peak was lower at Day 8 for TIV compared to that of LAIV at Day 28. A cluster of co-expressed interferon-inducible antiviral genes (IFIT3, OAS3, IFI44L, MX1, OAS2, IFIT1, IFI6) showed a statistically significant increase from pre-vaccination at Day 2 and an even stronger increase at Day 28 for the LAIV group. In both cases, the lower bound of the 95% bootstrap CI confidence interval (CI) was above 0 log₂ fold change. In contrast, the TIV group did not exceed pre-vaccination levels for this cluster in BLMC at any post-vaccination timepoint (Fig. 2). Interestingly, the Day 28 signal was reversed in PBMC samples with the TIV group showing statistically significantly higher than pre-vaccination levels, while the LAIV group remained at pre-vaccination levels.

The mean log₂ fold change from baseline and associated 95% bootstrap CIs for a co-expressed DE gene cluster of interferon-inducible genes (IFIT3, OAS3, IFI44L, MX1, OAS2, IFIT1, IFI6) are demonstrated by post-vaccination day, specimen type, and vaccine group. CI confidence interval, BML breast milk lymphocytes, LAIV live attenuated influenza vaccine, TIV trivalent influenza vaccine, PBMC peripheral blood mononuclear cells.

LAIV activated genes involved in innate immune signaling pathways in breast milk lymphocyte cells

Next, to assess the functional context of DE genes, we carried out pathway enrichment analysis (Supplementary Tables S20–S34). Two days following LAIV vaccination, several innate immune signaling pathways were significantly enriched in DE genes in BMLs at Day 2 including IL1 signaling and Signaling by interleukins as well as NF-kappa B signaling, Cytokine-cytokine receptor interaction, Toll-like receptor signaling, and TNF signaling pathways. Two of the top three enriched MSigDB immunologic gene sets were from a TIV/LAIV vaccine study by Nakaya et al. (PMCID:PMC3140559) (GSE29618 MONOCYTE VS MDC UP, MONOCYTE VS MDC DAY7 FLU VACCINE UP)¹⁹. Similar innate immune signaling pathways were enriched at Day 28 following LAIV vaccination. The top significantly enriched pathways included Osteoclast differentiation, Cytokine-cytokine receptor interaction, TNF signaling, NF-kappa B signaling, and Toll-like receptor signaling pathways. In contrast to Day 2, Influenza A, Cytokine signaling in immune system, Signaling by interleukins, Interferon signaling, and more specifically Interferon alpha-beta signaling were significantly enriched at Day 28. The association of DE genes and cell signaling was further confirmed by a significant enrichment of cell membrane-associated proteins (74 of 382 (19%) DE genes, source: GO Cellular Components). At Day 28, the best overlap with the LAIV/TIV study by Nakaya et al. was observed for GSE29618 MONOCYTE VS MDC UP, GSE29618 MONOCYTE VS PDC UP, and GSE29618 MONOCYTE VS MDC DAY7 FLU VACCINE UP with 41–47 overlapping up-regulated DE genes¹⁹. For the TIV group at Day 28, IL1 SIGNALING (Reactome) was significantly enriched in up-regulated DE genes (2 DE genes: IL1A and IL1B). Time trends of average pathway fold change responses and associated 95% confidence intervals confirmed significant activation of key innate immune signaling pathways at Day 2 and then later at Day 28 in BML following LAIV vaccination and activation of IL1 SIGNALING at Day 28 following TIV vaccination (Fig. 3). Pathway maps for enriched KEGG pathways color-coded by vaccine effect are included in the Supplemental Text (Figures S28–S80).

The median average log₂ fold change from baseline and associated 95% bootstrap Cis for a subset of enriched KEGG and Reactome pathways is demonstrated by post-vaccination (post-treatment) day, specimen type, and vaccine group. Each panel demonstrates the median average Log₂ fold change in a pathway over time: complement and coagulation cascades (A), NF-kappa B signaling (B), TNF signaling pathway (C), and IL-1 signaling pathway (D). LAIV live attenuated influenza vaccine, TIV trivalent influenza vaccine, BML breast milk lymphocytes, PBMC peripheral blood mononuclear cells, TNF tumor necrosis factor, IL1 interleukin 1.

Correlation between ELISA IgG/IgA and IgG/IgA gene expression responses

Next, we assessed whether IgG immunoglobulin heavy chain gene expression (IGHG1 and IGHG3) in PBMCs could serve as a marker for serologic response. In our cohort, as for the main study¹⁷, the directionality of the geometric mean titers (GMTs) was the same for all 5 antigens tested with higher GMTs for the TIV group (Table 1). Antibody responses to TIV at Day 28 were approximately three times higher than the antibody responses to LAIV, with the exception of the vaccine antigen for B/Brisbane/60/2008 (2011/2012 strain), in which there was a minimal difference between vaccine groups. In our study, a significant increase in gene expression from pre-vaccination for IGHG1 and IGHG3 genes in PBMCs was observed in the TIV group with the lower bounds of the CIs exceeding 0 log₂ fold change (Fig. 4A). The effect was stronger for IGHG1 which encodes the constant region of IgG1. The Day 28 increase in IGHG1 gene expression in PBMCs showed weak positive correlations with fold changes in ELISA IgG titers against 4 vaccine antigens: A/H3N2 (Perth and Victoria), A/H1N1 (California), and influenza B (Brisbane) (Supplemental Figures S2–S5). However, we observed a moderate positive association between IGHG1 gene expression and IgG responses for the influenza B antigen – Wisconsin, which is presented in Fig. 4B (r = 0.483, p = 0.0584). Sensitivity analysis following the removal of one strongly outlying observation (highlighted in red) resulted in a statistically significant correlation of r = 0.805 (p < 0.0003).

In (A, C), the mean log₂ fold change from baseline and associated 95% bootstrap CIs for co-expressed IgG genes (IGHG1 and IGHG3) in PBMC at Day 28 is demonstrated by vaccine group. For BML, gene expression levels for these two genes did not pass the low expression cut off for both vaccine groups. In (B, D), the scatterplot displays per-subject log₂ fold change in PBMC IGHG1 and IGHG3 gene expression versus log₂ fold change in ELISA IgG titers for the influenza B antigen (Wisconsin) 28 days following TIV or LAIV vaccination (n = 16, r: Pearson correlation, in green: linear regression fit, in blue: locally weighted regression fit). Sensitivity analysis following the removal of the outlying observation highlighted in red resulted in a correlation of r = 0.805 and p = 0.0003 for (B). LAIV live attenuated influenza vaccine, TIV trivalent influenza vaccine, PBMC peripheral blood mononuclear cells, ELISA enzyme-linked immunosorbent assay.

Following these results, we assessed whether immunoglobulin heavy constant alpha 1 and 2 gene expression (IGHA1 and IGHA2) in breast milk lymphocyte cells could serve as a marker for IgA response in breast milk. We observed a moderate positive association for the influenza B Wisconsin (r = 0.406, p = 0.1191) and the Swine H1N1 (0.486, p = 0.0565) antigens, which is presented in Supplemental Figure S9, S10, respectively; however, neither correlation was statistically significant.