System vaccinology analysis of predictors and mechanisms of antibody response durability to multiple vaccines in humans

Abstract

We performed a systems vaccinology analysis to investigate immune responses in humans to an H5N1 influenza vaccine, with and without the AS03 adjuvant, to identify factors influencing antibody response magnitude and durability. Our findings revealed a platelet and adhesion-related blood transcriptional signature on day 7 that predicted the longevity of the antibody response, suggesting a potential role for platelets in modulating antibody response durability. As platelets originate from megakaryocytes, we explored the effect of thrombopoietin (TPO)-mediated megakaryocyte activation on antibody response longevity. We found that TPO administration enhanced the durability of vaccine-induced antibody responses. TPO-activated megakaryocytes also promoted survival of human bone-marrow plasma cells through integrin β1/β2-mediated cell-cell interactions, along with survival factors APRIL and the MIF-CD74 axis. Using machine learning, we developed a classifier based on this platelet-associated signature, which predicted antibody response longevity across six vaccines from seven independent trials, highlighting a conserved mechanism for vaccine durability.

Extended Data Fig. 1 |. Analysis of kinetics of transcriptional responses.

(a) Scatterplot of the mean log2 FC of all BTMs in adjuvanted participants on day 1 (x axis) and in nonadjuvanted participants on day 3 (y axis). The Pearson correlation coefficient and p value are reported. (b) Kinetics of differential BTMs between day 3 prime and boost. Lines represent average module fold change among adjuvanted participants. The 10 BTMs with the greatest fold change on day 24 are plotted (same as those labeled in Fig. 1g).

Extended Data Fig. 2 |. Analysis of circulating Tfh cells post-vaccination.

(a) Gating strategy for sorting of four different CD4+ CXCR5+ Tfh populations: quiescent Tfh1, quiescent Tfh2, activated Tfh1, and activated Tfh2. (b) Frequency of activated Tfh cells post-vaccination in adjuvanted (orange) and non-adjuvanted participants (green), defined as percentage of PD1+ICOS+ cells within the CXCR5+ CD4+ T cell population. n=34 (H5N1+AS03) and 16 (H5N1) participants. Orange/green p values represent post-vaccination changes within each group (two-sided Wilcoxon test), gray p values represent between group comparisons (two-sided Mann-Whitney test). (c) Correlation of the day 28/21 fold change in activated Tfh cell frequencies with the day 42/21 fold increase in MN titers. The Pearson correlation coefficient and p value are reported. (d) Average log2 fold change of genes in activated (PD-1+ICOS+) versus non-activated (PD-1- ICOS-) Tfh1 (x-axis) and Tfh2 (y-axis) cells. The top 20 genes with the greatest average absolute fold change are annotated in red. (e) Boxplot of estimated frequency of monocytes in non-activated and activated Tfh based on digital cytometry of transcriptional profiles using CIBERSORT. Lines represent means, shaded area represents 95% confidence interval, and lines represent standard deviation. n=27 (unactivated) and 28 (activated), unpaired two-sided t-test. (f) BTMs significantly enriched (FDR<0.05) in activated versus non-activated Tfh cells. CIBERSORTx was used to estimate CD4 T cell specific expression in sorted Tfh transcriptional profiles, and then GSEA was used to identify enriched BTMs using genes ranked by their fold change between activated and non-activated Tfh. See Methods section for further details. (g) Genes in BTM M219; each ‘edge” (gray line) represents a coexpression relationship; colors represent the fold change in activated versus non-activated Tfh. (h) Genes in BTM M4.2; each ‘edge” (gray line) represents a coexpression relationship; colors represent the fold change in activated versus non-activated Tfh. (i) Clustered heatmap of the top 40 genes by fold change between activated and non-activated Tfh. Colors represent row-wise z-scores. *, P < 0.05; **, P < 0.01; ***, P < 0.001; ****, P < 0.0001.

Extended Data Fig. 3 |. Relationship between peak and durable antibody responses.

(a) Kinetics of HAI titers in response to H5N1+AS03 and IIV vaccination. Lines represent the geometric mean, and shaded areas represent the geometric standard deviation. IIV titers are from young adults (<65 years old) vaccinated with the 2010 and 2011 seasonal influenza vaccine (Nakaya et al.13). n=34 (H5N1+AS03) and 42 (IIV). (b) Scatterplot of the day 100 and day 42 HAI titers. The day 100/42 HAI residual is defined here as the vertical distance between a given point and the regression line, with the regression line representing the average day 100 titer expected given a particular day 42 response. Participants above/below the regression line are considered ‘persistent’ and ‘temporary’ responders, respectively. The Pearson correlation coefficient and p value are reported. (c) Heatmap of the Pearson correlation coefficients between expression of plasma cell and cell cycle BTMs and peak antibody titers (prime - day 21, boost - day 42) in adjuvanted participants. (d) Scatterplot of the day 28/day 21 mean log2 FC of M156.0 (x axis) versus the day100/day42 HAI residual in adjuvanted participants. The Pearson correlation coefficient and p value are reported.

Extended Data Fig. 4 |. Platelet-associated signatures of antibody response durability in a NHP vaccination model.

(a) Study design for AS03-adjuvanted COVID-19 vaccination in NHPs. (b) Kinetics of pseudotyped lentivirus neutralization antibody titers following vaccination in NHPs. (c) Bar plot of BTMs associated with antibody persistence on day 7 post-boost in the NHP COVID-19 vaccine study. GSEA was performed on genes ranked by correlation of their day 7 post-boost expression with the day 180/42 antibody residual. Modules shown are those with FDR < 0.05, NES ≥ 2, platelet-associated modules are highlighted in red.

Extended Data Fig. 5 |. CITEseq quality control data.

(a) Per-cluster cell proportions from day 21 and 28 samples before QC filtering. (b) Per-cluster cell proportions from each subject before QC filtering. (c) Scatterplots of day 28/21 FCs among day 28 DEGs via microarray (x axis) and pseudobulk estimates via CITEseq (y axis) for each subject. The Pearson correlation coefficient and p value are reported. (d) CITE-seq antibody abundance in each cell before QC filtering. (e) DEGs in each cluster compared to all other clusters before QC filtering. (f) Percell total reads and number of detected transcripts by cluster before QC filtering.

Extended Data Fig. 6 |. Effect of TPO administration on the bone marrow compartment.

(a) Study design for TPO administration in mice. Recombinant mouse TPO was injected intraperitoneally at 12.5 μg/kg daily for 5 days. Bone marrow cells were analyzed by flow cytometry on days 4, 7 and 11 following the initial TPO injection. (b-d) Flow cytometry gating strategy (b), megakaryocytes frequency in bone marrow live cells (c), and megakaryocytes ploidy stages (d) after TPO injection were analyzed. n = 6–15, Dunnett’s multiple comparisons test for panel C, unpaired t test for panel D. (e-f) Flow cytometry gating strategy and frequencies of bone marrow immune cell populations. n = 6–15, Dunnett’s multiple comparisons test. (g-h) Flow cytometry gating strategy and frequencies of bone marrow stem and progenitor cells. n = 6–15, Dunnett’s multiple comparisons test. *, P < 0.05; **, P < 0.01; ***, P < 0.001; ****, P < 0.0001.

Extended Data Fig. 7 |. Analysis of platelet-associated durability signature across vaccines and cell types.

(a) Scatterplots of Day 7 vs Day 0 log2 fold changes (relative to last vaccination) of BTM M85 versus antibody residual in each vaccine dataset. Pearson correlation coefficient and p values are reported. (b) Heatmap of expression of platelet-associated BTMs across cell clusters in the CITE-seq data (Fig. 4). Colors represent row-wise z-scores of average expression of all genes in each module. (c) Scatterplots of Day 7 vs Day 0 log2 fold changes of platelet-expressed BTMs from PBMC (x-axis) and Paxgene (y-axis) samples in a cohort of healthy adults vaccinated with seasonal influenza vaccine. Pearson correlation coefficients and p values are reported.

Figure 1.. AS03 induces potent early transcriptional signatures which are enhanced after a booster vaccination

(A) Study overview. A total of 50 healthy participants aged 21–45 years old were randomized 2:1 to receive two doses 21 days apart of a monovalent, split-virion, inactivated H5N1 clade 2.1 A/Indonesia/05/2005 influenza vaccine, administered with (n=34) or without (n=16) the AS03 adjuvant. Biologic samples were collected and analyses performed at regular intervals (gray squares) as illustrated in the diagram. (B) Number of DEGs (p <0.01 and log₂ FC >0.2) post-vaccination in adjuvanted (orange) and non-adjuvanted (green) participants. Two-sided paired t-test was used to compute p values in each group. (C-D) Average blood transcriptional module (BTM) enrichment scores by cell type/pathway on day 1 (dark) and day 7 (light) after prime (C) and boost (D) in each group. (E) Enrichment scores of interferon-related BTMs on days 1–7 after prime (top) and boost (bottom) in each group. (F) Scatter plot of day 1 fold changes (x – prime, y – boost), for BTMs differentially expressed (FDR<0.03) in adjuvanted participants on day 1 between prime and boost. BTMs are color-coded as indicated in the legend. (G) Genes in BTM M111.1; each “edge” (gray line) represents a coexpression relationship; colors represent the day 1 fold change after prime (left) and boost (right) in adjuvanted participants. (H) Scatter plot of day 3 fold changes (x – prime, y – boost), for BTMs differentially expressed (FDR<0.03) in adjuvanted participants on day 3 between prime and boost. BTMs are color-coded as indicated in the legend. (I) Genes in BTM 89.0; each “edge” (gray line) represents a coexpression relationship; colors represent the day 3 fold change after prime (left) and boost (right) in adjuvanted participants. n=34 (H5N1+AS03 group) and 16 (H5N1 group). Panels F-I represent results from the H5N1+AS03 group only.

Figure 2.. Adjuvanting with AS03 enhances the magnitude, affinity, and breadth of antibody responses to H5N1 vaccination.

(A-B) Microneutralization (MN, A) and hemagglutination inhibition (HAI, B) titers against the H5N1 A/Indonesia vaccine strain in adjuvanted (orange) and non-adjuvanted (green) subjects. Geometric means are presented in thick lines, while shades are for geometric standard deviations (SD). (C) MN titers against heterologous clade 2 H5N1 strains. Geometric means are presented in thick lines, while error bars represent geometric standard deviations (SD). (D) H5 head and stem IgG binding capacity in resonance units (RU) as measured by surface plasmon resonance (SPR). Median values and interquartile ranges are shown in boxplots, violin plots show sample distributions. (E) Day 42 post-vaccination fold change in IgG affinity against the H5 head (left) and stem (right) as measured by SPR. (F-G) Scatterplot of the day 42/21 fold change in IgG antibody binding (F) and affinity (G) against the H5 head and MN titers. The Pearson correlation coefficient and p value are reported. n=34 (H5N1+AS03) and 16 (H5N1) participants. All between group comparisons were performed using two-sided Mann-Whitney tests. *p<0.05; **p<0.01; ***p<0.001; ****p<0.0001.

Figure 3.. Durability of antibody responses to influenza vaccination is associated with a transcriptional signature of cellular migration

(A) Heat map of BTMs whose post-vaccination activity is associated with antibody persistence in both H5N1+AS03 (orange) and IIV (pink) vaccine responses. Gene set enrichment analysis (GSEA, FDR < 0.05; 1,000 permutations) (Subramaniam et al. 2004) was used to identify positive (red), negative (blue), or no (gray) enrichment of BTMs within pre-ranked gene lists, where genes were ranked according to their correlation with the residual of the day 100 versus day 42 (H5N1+AS03) or day 180 versus day 28 (IIV) antibody response. (B) Genes in BTM M196; each “edge” (gray line) represents a coexpression; colors represent the correlation of the day 7 gene expression with the day 100 (H5N1+AS03, B) or 180 (IIV, C) antibody response residual (positive – red, negative – blue) vaccination. (C) Workflow for prediction of antibody response durability. Single sample GSEA was used to generate participant-level enrichment scores based on gene day 7/0 fold changes for all BTMs. BTMs whose scores were most correlated with the antibody residual (p<0.25) were selected as inputs to a linear regression model using forward feature selection. See Methods for additional details. (D) Scatterplot of actual versus predicted day 100 antibody response residuals in the training (H5N1+AS03), validation (IIV) and blind test (CHI, H5N1+AS03) datasets. The Pearson correlation coefficient and p value are reported. (E) Bar plot of genes from BTMs selected by the predictive model (D) whose expression on day 7 significantly correlates with antibody longevity (day 100 or day 180 response residual) in the H5N1+AS03 and IIV datasets. Bars represent the combined Pearson correlation coefficient from both datasets.

Figure 4.. CITEseq analysis reveals a platelet origin for transcriptional signature of antibody persistence

(A) UMAP representation of PBMCs from all analyzed samples (n = 12, 3 ‘persistent’ and 3 ‘waning’ responders, day 21 and day 28 samples from each subject) colored by manually annotated cell type. (B) UMAP representation of PBMCs from all analyzed samples showing the per-cell sum of expression for all genes in the predictive signature of antibody persistence (Figure 3C). (C) Left panel: boxplot of pseudobulk fold change of the antibody persistence signature among persistent (red) and waning (blue) responders. Right panel: line graph of changes in pseudobulk expression of the antibody persistence signature when a given cell cluster is removed from the pseudobulk calculation. n=3 (persistent group) and 2 (waning group). One-sided, unpaired t-test was used to compare groups. Boxes show median and 25th–75th percentiles, and whiskers show the range. (D) Heatmap of average expression of predictive model genes correlated with antibody residuals (Figure 3E) across cell clusters. Colors represent row-wise z-scores. (E) Bar chart of significantly enriched BTMs (FDR<0.05) via overrepresentation testing of DEGs (two-sided Wald p ≤ 0.05 & FC ≥ 0.25) between persistent and waning responders within plasmablasts at day 28 after vaccination. (F) Heatmaps of genes within M238, M219, and M216 (left), M4.1 (center), and M250 (right) modules. Colors represent row-wise z scores of average expression within plasmablasts.

Figure 5.. Antibody response durability is modulated by thrombopoietin receptor signaling.

(A) Study design for TPO administration to mice during spike+AS03 immunization. (B) Kinetics of anti-Spike IgG titers in Mock and TPO-treated wild-type mice post-boost vaccination. n=10–14, two-way ANOVA with post hoc Bonferroni test. Data pooled from three independent experiments. (C) Representative flow cytometry plots (left), percentage of CD41hiCD61+ megakaryocytes (CD3-CD11b-CD19-) among CD45+ bone marrow cells (middle), and megakaryocyte counts in femurs in (right) from WT and cMplKO mice. n = 9–10, two-sided unpaired t test. (D) Total IgA (left) and IgA (right) levels in the serum of naïve WT and cMplKO mice. n = 27–28, two-sided unpaired t test. (E) Representative flow cytometry plots (left) and percentage of CD138+TACI+ plasma cells (CD3-CD19-) (right) among live CD45+ bone marrow cells in naïve WT and cMplKO mice. n = 8–10, Two-sided Mann-Whitney test. (F) Fold change in anti-spike IgG titers (relative to WT mock-treated) at day 21 post-boost in WT and cMplKO mice following prime/boost spike+AS03 vaccination with or without TPO injections at boost. n=4–5, Two-sided one-way ANOVA with post hoc Bonferroni test. (G) Schematic illustration of mouse splenic B cell (CD45.1) and bone marrow megakaryocyte (CD45.2) co-culture experiment. (H) Left, representative flow cytometry plots showing the percentage of live Annexin V- cells among CD45.1+CD138+ cells co-cultured with (MK) or without megakaryocyte (No MK). Right, the percentage of live Annexin V-CD138+ cells among total CD45.1+ cells. (I) IgM concentration in the supernatant at day 2 after co-culture measured by ELISA. (J-K) Percentage of live Annexin V-CD138+ cells among total CD45.1+ cells (J) and the IgM concentration of supernatant (K) from in vitro differentiated plasmablasts co-cultured with megakaryocytes for two days in the presence of anti-IL-6 (aIL-6, 5 μg/ml) or anti-APRIL (aAPRIL, 25 μg/ml). For panels H and I, data pooled from 6 independent experiments. n = 15, paired t test. For panels J and K, data pooled from 3 independent experiments. n = 7, one-way ANOVA with post hoc Dunnett test.

Figure 6.. Thrombopoietin receptor signaling in human megakaryocytes promotes plasma cell survival and antibody production in through APRIL, MIF/CD74 axis and integrin-dependent contact.

(A) Schematic illustration shows isolation of plasma cells (PC) and megakaryocytes from human bone marrow aspirates. Human bone marrow plasma cells then were co-cultured with autologous bone marrow megakaryocytes at a ratio of 1 to 10 in the presence of recombinant human TPO. (B) Live plasma cell counts at 3.5 days after co-culture. n=7, paired two-tailed t-test. (C) IgM and total IgG concentration in the supernatant measured by ELISA. Donors are represented by dots, and dots from the same donor are connected with dashed lines. Means of the data points were connected with solid lines. n=8, paired two-tailed t-test. (D-E) Plasma cells were co-cultured with autologous megakaryocytes for 3.5 days in the presence of isotype (mouse IgG) or various neutralizing antibodies. Live plasma cell counts (D), IgM and total IgG concentration in the supernatant (E) are shown. n=6, Two-sided Dunnett’s multiple comparisons test. *p ≤ 0.05; **p ≤ 0.01; ***p<0.001; ****p<0.001.

Figure 7.. Summary of proposed mechanism for involvement of megakaryocytes in promotion of durability of antibody responses to vaccination.

When activated by certain inflammatory mediators following vaccination, megakaryocytes promote survival of plasma cells in the bone marrow through integrin β1 and β2 mediated cell-cell contact and production of the cytokines APRIL and MIF. As megakaryocyte activation also causes increased thrombopoiesis, platelet-associated transcriptional signatures in blood are a marker of megakaryocyte activity and predictive of antibody response durability. Diagram created in BioRender.

Figure 8.. Platelet signature association with antibody response durability is conserved across multiple vaccines.

(A) Kinetics of antibody titers following vaccination across vaccines. Included vaccines were those for which day 7 transcriptomic data and peak (∼day 21) + long-term (> day 90) antibody titers were available. (B-C) Heatmap of associations of (B) platelet and (C) cell adhesion/migration BTM expression with antibody durability residuals across vaccines. Color represents the Pearson correlation t test statistic, which accounts for differences in sample sizes between studies. (D) ROC curves of prediction of high/low antibody response durability in training (Emory H5N1+AS03), validation (IIV) and test (CHI H5N1+AS03, Pfizer COVID-19 mRNA, MPSV4, MCV4, malaria RTS,S) datasets. High versus low antibody response durability predictions were generated using an elastic net model using day 7 expression of platelet and cell adhesion/migration BTMs as input features. One-sided p values were estimated via permutation test with n=10,000 permutations. See Methods section for further details. (E) Heatmap of average expression of genes within BTMs selected in the predictive model (D) across cell clusters in the CITE-seq data (Figure 4). Colors represent row-wise z-scores.