Systems biological assessment of human immunity to BNT162b2 mRNA vaccination

The emergency use authorization of two COVID-19 mRNA vaccines in less than a year since the emergence of SARS-CoV-2 represents a landmark in vaccinology1,2. Yet, how mRNA vaccines stimulate the immune system to elicit protective immune responses is unknown. Here we used a systems biological approach to comprehensively profile the innate and adaptive immune responses of 56 healthy volunteers vaccinated with the Pfizer-BioNTech mRNA vaccine. Vaccination resulted in robust production of neutralizing antibodies (nAbs) against the parent strain and a variant of concern, B.1.351, and significant increases in antigen-specific polyfunctional CD4 and CD8 T cells after the second dose. There was also a robust innate response induced within the first 2 days of the booster vaccination, compared to the first dose. Specifically, there were strongly enhanced: (i) frequency of CD14+CD16+ inflammatory monocytes; (ii) concentration of IFN-γ in the plasma, which correlated with enhanced pSTAT3 and pSTAT1 levels in monocytes and T cells; and (iii) transcriptional signatures of innate responses characteristic of antiviral responses, within 2 days following booster vaccination, compared to the primary response. Consistent with these observations, single-cell transcriptomics analysis of 242,479 leukocytes demonstrated a ~100-fold increase in the frequency of a myeloid cell cluster containing monocytes and dendritic cells, enriched in interferon-response transcription factors (TFs) and reduced in AP-1 TFs, only after the second immunization. Finally, we identified distinct molecular pathways of innate activation that correlate with CD8 T cell and nAb responses, and identify an early monocyte-related signature that was associated with the breadth of the nAb response against the B1.351 variant strain. Collectively, these data provide insights into the cellular and molecular responses induced by mRNA vaccines and demonstrate their capacity to prime the immune system to mount a more potent innate immune response following booster immunization.

responses were primarily Th1-type expressing IL-2, IFN-g or TNF-a although there were low levels of IL-4 induction (Fig. 1e). On the other hand, IFN-g and TNF-a were the dominant responses in CD8 T cells with 25 of the 31 participants responding with IFN-g response atleast 3 times higher than the baseline (Fig.   1f, g). Of note, three individuals with no known exposure to SARS-CoV-2 or clinical symptoms associated with COVID-19 demonstrated ~0.2% of Spike-specific CD8 T cell responses at the baseline consistent with previous studies indicating that 10% of healthy individuals have cross-reactive CD8 T cell responses 9 .
There was no significant correlation between T cell response and age or nAbs against the parent or B.1.351 strains (Extended Data Fig. 1b -d).

BNT162b2 vaccination does not cause new onset autoantibodies or anti-cytokine antibodies (ACA).
Multiple studies have demonstrated the presence of serum autoantibodies 10 11,12 13 14 and ACA 5 15 in patients infected with SARS-CoV-2, as well as development of new-onset antibodies in a subset of hospitalized COVID-19 patients 16 . We screened serum samples from 31 participants and compared mean fluorescence intensity (MFI) of IgG autoantibodies and ACA on days 0, 21, and 42 using a 55-plex antigen array and a 58-plex cytokine array. We included prototype serum samples from 17 patients with autoimmune and immunodeficiency disorders as positive controls. Five vaccinated subjects had preexisting autoantibodies (three suggestive of autoimmune thyroiditis, one low level PDC-E2+ associated with primary biliary cirrhosis, and one subject positive for connective tissue disease antigens RPP25 (Th/To), PM/Scl75, and SSB (La), (Extended Data Fig. 2). Anti-cytokine antibodies were largely absent or were observed at low MFI (Extended Data Fig. 3). Two subjects, including one who was also TPO+, had anti-IL-21 autoantibodies, and two additional subjects had anti-IL-1 antibodies (Extended Data Fig. 4, 5).
Importantly, in subjects with pre-existing autoantibodies or ACA, none had adverse events, nor did levels of pre-existing autoantibodies or ACA change in response to vaccination. New-onset autoantibodies or ACA were not observed in any of the vaccinated subjects (Extended Data Fig. 2, 3).

Innate immune responses to mRNA vaccination
While previous studies have analyzed adaptive immune responses to BNT162b2 vaccination 3,17 , little is known about the innate immune responses induced by BNT162b2. To analyze innate immune responses, we first assessed whole blood samples of 27 individuals collected across multiple time points after vaccination using a 38-parameter mass cytometry (Cytometry by Time of Flight, CyTOF) panel containing an assortment of cytokines and phospho signaling markers (Extended Data Table 3). Unsupervised clustering identified 14 major cell types ( Fig. 2a and Extended Data Fig. 6a) which were further subtyped manually (Extended Data Fig. 6b). The frequency of intermediate monocytes (CD14 + CD16 + monocytes), key orchestrators of innate immunity, significantly increased 2 days after the first immunization. Strikingly, the frequencies were substantially higher on day 23, 2 days post-secondary vaccination, compared to the frequencies at day 2 post primary vaccination (Fig. 2b, and Extended Data Fig. 6c). There was no correlation with age (Extended Data Fig. 6d). In addition, the CyTOF analysis revealed greatly enhanced phosphoSTAT3 (pSTAT3) and phosphoSTAT1 (pSTAT1) in multiple cell types on day 1 after secondary immunization, relative to their modestly increased expression on day 1 post primary immunization (Fig.   2c, d). These data suggested that the BNT162b2 mRNA vaccination induced a heightened innate immune response following secondary immunization, relative to the response after the primary immunization.
To further investigate this phenomenon, we measured 92 different cytokines in plasma collected at various time points from 31 vaccinees using the Olink platform. Of the 67 cytokines that were detected within the dynamic range of the assay, the concentration of two cytokines, IFN-g and CXCL-10, were significantly increased on days 1 and 2 after primary immunization (Fig. 2e, left panel). Similar to the observations on intermediate monocytes and pSTAT3 signaling, the concentrations of these cytokines were increased even further after the secondary immunization (Fig. 2e, right panel). IFN-g in particular rose 11.3fold between day 1 and 22 (Fig. 2f). CXCL-10, on the other hand, peaked on day 2 suggesting a response driven by IFN-g (Extended Data Fig. 6d). Interestingly, the anti-inflammatory cytokine IL-10 also showed a similar pattern of response to that of IFN-g although this trend did not reach statistical significance (Extended Data Fig. 6e). We also measured type I IFN (IFN-a and IFN-b) using ELISA, which were below detection limit (<2 pg/ml) at any time point after vaccination. Furthermore, there was a strong correlation between plasma IFN-g levels and pSTAT1/3 expression levels across several cell types (Fig.   2g, h). Of note, the concentration of cytokines returned to baseline levels by day 28 (i.e., 7 days postsecondary vaccination), suggesting origin from an innate immune cell type. Collectively, these data demonstrate that vaccination with BNT162b2 stimulates low levels of innate immune responses after primary immunization, which strikingly increase after the secondary immunization.

Transcriptional signatures induced by vaccination
We next investigated the transcriptomic changes induced by BNT162b2 vaccination. We performed bulk mRNA sequencing of 185 samples obtained from 31 participants across 7 time points. Six of 185 samples did not pass quality control and were removed from the analysis (Extended Data Fig. 7a, b). Strikingly, secondary immunization generated a much greater transcriptional response compared to primary immunization, with nearly a four-fold increase of DEGs found at day 22 compared to day 1 (Figure 3a).
This was consistent with the increased markers of innate immunity demonstrated after secondary immunization by both CyTOF and Olink (Figure 2b-h).
In order to explore the specific transcriptional pathways altered in response to mRNA-based vaccination, we performed gene set enrichment analysis (GSEA) 18 using a set of previously defined blood transcriptional modules (BTMs) 19 at each post-vaccination timepoint. Both doses of BNT162b2 induced upregulation of antiviral and interferon response modules, including M75, M127 and M111.0 (Figure 3b).
However, booster immunization led to a significantly broader innate response. In addition to induction of antiviral pathways, the boost dose led to increases in dendritic cell activation and upregulation of TLR signaling, monocyte, and neutrophil modules on days 22-23, which were previously decreased post-prime ( Figure 3c). We compared fold changes of the genes within these modules and found that many antiviral genes showed a greater increase following the boost dose (Figure 3d), and other inflammatory genes switched from downregulation to upregulation between prime and boost ( Figure 3e). These results were consistent regardless of the baseline timepoint used (Extended Data Fig. 7c, d). A complete list of enriched BTMs is presented in Supplementary Table 1.
As mortality to COVID-19 is highest among the elderly and older populations are known to mount inferior responses to many vaccines 20,21 , an important question we sought to address was whether or not there were age-associated differences in response to mRNA vaccination. We correlated the per-person fold changes of each gene with participant age and used GSEA to identify age-dependent response pathways.
On day 22, we observed that younger subjects tended to have greater changes in monocyte, inflammatory response, and platelet-related expression, whereas older subjects had increased response in B and T cell modules (Figure 3f).
Finally, given that our serum cytokine analysis revealed that IFN-g responses were also significantly higher following booster immunization, we asked whether there was any association between the level of IFN-g and the increased innate responses following the boost. Indeed, both interferon response and other inflammatory modules were significantly enriched by GSEA when using genes ranked by correlation with IFN-g on day 22 (Figure 3g). Furthermore, the average fold changes of these modules also correlated with IFN-g (Figure 3h), suggesting that IFN-g may play an important role in driving enhanced innate and antiviral expression post-boost.

Single-cell transcriptional response to BNT162b2 vaccination
We used cellular indexing of transcriptomes and epitopes by sequencing (CITE-seq) to determine the cellular origin of the enhanced antiviral and inflammatory gene signatures, and to more broadly characterize transcriptional signatures induced by vaccination at the single-cell level. To this end, we analyzed 45 PBMC samples from 6 individuals across seven time points (days 0, 1, 2, 7, 21, 22, 28 and 42). We used an "enrichmix" strategy in which we enriched dendritic cells (DCs) and mixed with total PBMCs at a ratio of 1:2 to capture minor populations such as plasmacytoid dendritic cells (pDCs) sufficiently in the CITE-seq 4 . After preprocessing and quality control, we obtained 242,479 high quality transcriptomes, which were segregated into 18 cell clusters (Fig. 4a, Extended Data Fig. 8a, b). Strikingly, one cluster C8 (annotated C8_CD14 + BDCA1 + PD-L1 + ), expressing CD14, VNN, CD1C, FCGRIA, and CD274 and other myeloid markers at gene and/or protein level, emerged on day 22, 1-day after secondary vaccination (Fig. 4b). These cells also expressed several ISGs including IFI30, IFITM3, WARS and GBP1, and constituted only ~0.01% of the Lin -HLA-DR + population on day 1 after primary immunization but increased almost 100X to ~1% one day after secondary immunization ( Fig. 4c). Notably, the emergence of C8 correlated with plasma IFN-g levels measured by Olink or an independent ELISA assay (performed because data of 2 participants were unavailable in Olink due to technical reasons), with the participant 2053 demonstrating a delayed increase in IFN-g as well as C8 on day 28 (Extended Data Fig. 8c, d). Furthermore, iterative removal of each cluster from a pseudobulk score showed that C8 contributed to IFN and monocyte BTMs observed in the bulk transcriptomics data (Extended Data Fig. 8e). To further delineate the cellular composition of C8, we reembedded C8 with UMAP, using harmony to correct for participant-specific biases. Using Louvain clustering, we resolved seven distinct clusters within the original C8 cluster (Fig. 4d). The cluster C8 proved to a heterogeneous mix of classical monocytes (C8_0, C8_1 and C8_3), cDC2 (C8_2) and intermediate monocytes (C8_4) as evidenced by the proximity to the original clusters measured by Euclidean distance (Fig. 4e). More interestingly, two sub clusters, C8_1 and C8_2, expressed significantly higher ISGs compared to their parent clusters (Fig. 4f) 22 . We asked if C8 represents an analogous cell type at the transcriptional level.
We found that the C8 has a relatively higher expression of TFs IRF1, STAT1, STAT2, STAT3, IRF8 and reduced levels of AP-1 TFs FOS, JUNB, JUND and ATF3, the same TFs that defined the monocyte population in the previous study (Fig. 4g). We also confirmed this using an extended set of genes for which the chromatin accessibility profile was higher 21 days after H5N1/AS03 vaccination (Fig. 4h).
Next, we set out to identify transcriptional changes induced by vaccination more broadly within each cell type. Given that we observed a higher level of pSTAT1/3 in multiple cell types using CyTOF and a higher magnitude of IFN response after secondary immunization by bulk RNAseq, we asked if there is an enhanced IFN response in multiple cell types or it is primarily driven by the emergence of cluster C8.
Interestingly, the IFN signature was induced in all cell types present on day 1 and day 22, and the higher magnitude of response on day 22 was more evident (Fig. 5a). The decrease in NK cell signatures on day 22 in bulk RNAseq was another intriguing feature observed in the bulk RNAseq. We have previously shown that the NK cells decrease in frequency following TIV vaccination, especially in young adults 23 . In line with this, we observed a significant reduction in the frequency of NK cells on day 22 in the CyTOF dataset ( Fig. 5b). However, the genes in the NK cell modules were also significantly downregulated within NK cells present on day 22 (Fig. 5c). Conversely, the NK cells on day 22 had a higher activation status and higher levels of AP-1 transcription factors known to be driven by IL-2-mediated activation of NK cells 24 (Fig. 5d).

Comparison of transcriptional responses with other vaccines
As mRNA vaccines have only recently received approval for use in humans, the degree to which these vaccines induce similar or distinct immune responses compared to other vaccine types, such as inactivated or live attenuated vaccines, is unknown. To address this, we utilized a set of previously published vaccine trials from our group as well as several publicly available datasets to perform a comparative analysis with Pfizer-BioNTech BNT162b2 (see Extended Data Table. 4 for a summary of included vaccine datasets).
In order to compare the relative similarity in transcriptional responses between vaccines, we generated similarity matrices through pairwise correlations of mean gene fold changes between vaccines at days 1 and 7 post-vaccination. While the day 1 response to the prime dose of BNT162b2 showed little overlap with other vaccines, the day 1 boost response was broadly similar to a group of vaccines containing either potent adjuvants (H5N1+AS03), live viral vectors (Ebola and HIV), or inducing a recall response (inactivated influenza) (Fig. 6a). Meanwhile, day 7 responses to both the prime and boost dose exhibited very weak correlation both between themselves and with other vaccines, suggesting little commonality in induced transcriptional signatures (Fig. 6b).
To identify the specific pathways induced by BNT162b2 that were unique or shared with other vaccines, we performed BTM-level GSEA on each vaccine dataset using genes ranked by pre-versus postvaccination t-statistic at each timepoint. We found that on day 1, the boost dose induced a robust innate response including upregulation of modules involved in antigen presentation, dendritic cell and monocyte activation, interferon signaling, and inflammatory responses, all of which were commonly induced by the Ebola and HIV live viral vector vaccines, seasonal influenza vaccine, and both doses of H5N1+AS03 (Fig.   4c). Conversely, the prime dose produced a much narrower response, only activating a limited set of interferon signaling modules. Instead, on day 7 the inverse trend occurred, with the boost dose having almost no commonly enriched modules but the prime dose sharing a cell cycle-related transcriptional signature with many vaccines (Extended Data Fig. 9a). However, in most other vaccines, the cell cycle signature is also associated with upregulation of B cell and plasma cell modules, reflecting the emergence and expansion of antibody-secreting cells 19,23 . This induction of B cell and plasma cell modules was absent in the BNT162b2 prime dose day 7 response (Extended Data Fig. 9b). Given that BNT162b2 successfully promoted a robust antibody response ( Fig. 1), the lack of any detectable plasma cell or B cell signature on day 7, particularly post-boost, was surprising and unlike other vaccine responses that we are aware of. I. It is possible that the kinetics of the plasma cell response to this vaccine are more delayed and this signature was therefore not captured at the day 7 timepoint.

Transcriptional correlates of antibody and T cell responses
As cellular and humoral immunity are the chief functional components mediating protection from infection, a key question is whether early transcriptional signatures exist that are associated with either the antibody or T cell responses following vaccination with mRNA vaccines. We therefore used GSEA to identify transcriptional modules whose expression at various timepoints post-vaccination was correlated with either the day 42 nAb or day 28 CD8+ IFN-g+ T cell response (Fig. 7a). In general, there was little overlap between signatures associated with the nAb and IFN-g CD8 T cell response, suggesting that there are distinct molecular pathways leading to cellular and antibody responses to BNT162b2. On day 22 (1-day post-boost), there was a striking divergence in signatures, with monocyte-related modules correlated with nAb responses while interferon and antiviral signatures highly associated with the later day 28 IFN-g CD8 T cell response (Fig. 7b). Surprisingly, plasma cell and cell cycle modules, which have previously been identified at day 7 as robust signatures of antibody response to other vaccines such as inactivated influenza vaccine 23 , were not associated at day 7 following either the prime or boost dose with the day 42 nAb titers.
The continued evolution of SARS-CoV-2 variants is a serious concern for the success of ongoing vaccination efforts. Therefore, we evaluated whether there are innate correlates of the cross-neutralization potential induced by vaccination. To this end, we defined a cross-neutralization index, using a ratio of variant to WA1 nAb titers, and used an enrichment approach to identify correlates of cross-neutralization.
Monocyte and neutrophil BTMs as well as TLR and innate immune pathways were highly associated with cross-neutralization index (Fig. 7c) suggesting a central role for myeloid cells in the effective immunity induced by mRNA vaccination. More interestingly, the frequency of classical monocytes at peak, 2 days post-secondary vaccination, as measured by CyTOF, strongly correlated with cross-neutralization index (Fig. 7d). To further evaluate this, we determined a gene score that defines C3, the classical monocyte cluster in the CITE-seq data, and asked if this gene score in the bulk RNAseq also correlates with crossneutralization. Clearly, there was a strong correlation (Fig. 7e). Collectively, these data demonstrate that while IFN signatures are associated with the CD8 T cell responses, monocyte and neutrophil gene signatures and TLR signaling BTMs strongly correlate with an nAb response against the WA1 as well as the B.1.351 strain.
In summary, our study describes a systems-based analysis that provides novel insights into the innate and adaptive immune responses to the Pfizer-BioNTech mRNA vaccine. We used a multiomics approach to define early transcriptional correlates of T cell and antibody responses. The data also demonstrate an enhanced innate immune response following secondary immunization indicating an innate memory-like response 25     inflammatory monocytes (right panel) and plasma IFN-g levels.
In b and f, the statistically significant differences between the peak and baseline time points were measured using two-sided Wilcoxon matchedpairs signed-rank test. The difference peak time points were measuring using two-sided Mann-Whitney rank-sum test (*p < 0.05, **p < 0.01, ***p < 0.001 and **** p < 0.0001). Blue and red dots indicate female and male participants, respectively.

Human subjects and experimentation
Fifty-six healthy volunteers were recruited for the study under informed consent. The study was approved by Stanford University Institutional Review Board (IRB 8269) and was conducted within full compliance of Good Clinical Practice as per the Code of Federal Regulations. The demographics of all participants were provided in Extended Data Table 1.

Focus Reduction Neutralization Titer assay
Neutralization assays with authentic SARS-CoV-2 virus were performed as previously described 27  plates were washed twice with PBS and foci were visualized on a fluorescence ELISPOT reader (CTL ImmunoSpot S6 Universal Analyzer) and enumerated using Viridot 29 . The neutralization titers were calculated as follows: 1 -(ratio of the mean number of foci in the presence of sera and foci at the highest dilution of respective sera sample). Each specimen was tested in two independent assays performed at different times. The FRNT-mNG50 titers were interpolated using a 4-parameter nonlinear regression in GraphPad Prism 8.4.3. Samples with an FRNT-mNG50 value that was below the limit of detection were plotted at 10. For these samples, this value was used in fold reduction calculations.

Focus Reduction Neutralization Titer assay against the variants of concern
The wildtype infectious clone SARS-CoV-2 (icSARS-CoV-2), derived from the 2019- and imaged on an ELISPOT reader (CTL).

Intracellular cytokine staining assay
Antigen-specific T cell responses were measured using the ICS assay as described previously 31 . Live frozen PBMCs were revived, counted and resuspended at a density of 2 million live cells/ml in complete RPMI (RPMI supplemented with 10% FBS and antibiotics). The cells were rested overnight at 37°C in CO2 incubator. Next morning, the cells were counted again, resuspended at a density of 15 million/ml in complete RPMI and 100 µl of cell suspension containing 1.

Bead-based antigen arrays.
We used an existing bead-based autoantigen array, and a cytokine array with expanded content that was based on our recent COVID-19 studies 16 . A complete list of all antigens, vendors, and catalogue numbers can be found in Supplementary Tables 2 and 3. The "COVID-19 Autoantigen Array" included 55 commercial protein antigens associated with connective tissue diseases (CTDs) (Supplementary Table 2).
The "COVID-19 Cytokine Array" comprised 58 proteins including cytokines, chemokines, growth factors, acute phase proteins, and cell surface proteins (Supplementary Table 3). Antigens were coupled to carboxylated magnetic beads (MagPlex-C, Luminex Corp.) such that each antigen was linked to beads with unique barcodes, as previously described 16,32,33 . Prototype human plasma samples derived from participants High-dimensional analysis of phospho-CyTOF data was performed using an R based pipeline described in 35 . Briefly, the raw fcs files were imported into R and the data were transformed to normalize marker intensities using arcsinh with a cofactor of 5. For visualization, another transformation was applied that scales the expression of all values between 0 and 1 using percentiles as the boundary. Cell clustering was performed with 4,000 cells randomly selected from each sample using FlowSom and ConsensusClusterPlus. The transformed matrix was used as an input for FlowSom and cells were separated into 20 clusters. To obtain reproducible results (avoid random start), a seed was set for each clustering. The 20 clusters were manually annotated based on the lineage marker expression and were merged to produce the final clusters. The clusters were visualized in two-dimensional space using UMAP. The abundance of cell populations was determined using Plotabundance function. In parallel, the data were manually gated to identify 25 immune cell subpopulations that were not well-distinguished in UMAP and used for all quantitation purposes.

Plasma protein profiling using multiplex Olink panel
We measured cytokines in plasma using Olink multiplex proximity extension assay (PEA) inflammation panel (Olink proteomics: www.olink.com) according to the manufacturer's instructions and as described before (41). The PEA is a dual-recognition immunoassay, where two matched antibodies labelled with unique DNA oligonucleotides simultaneously bind to a target protein in solution. This brings the two antibodies into proximity, allowing their DNA oligonucleotides to hybridize, serving as template for a DNA polymerase-dependent extension step. This creates a double-stranded DNA "barcode" which is unique for the specific antigen and quantitatively proportional to the initial concentration of target protein. The hybridization and extension are immediately followed by PCR amplification and the amplicon is then finally quantified by microfluidic qPCR using Fluidigm BioMark HD system (Fluidigm Corporation. South San Francisco, California).

Bulk Transcriptomics
Whole blood samples were collected in Paxgene tubes (BD Biosciences) and were frozen at -80°C until RNA isolation. RNA was isolated from each sample using the miRNeasy Mini kit (Qiagen) and 10 ng of total RNA was used as input for cDNA synthesis using the Clontech SMART-Seq v4 Ultra Low Input RNA kit (Takara Bio) according to the manufacturer's instructions. Amplified cDNA was fragmented and clustering and UMAP projections. Clusters were identified with Seurat SNN graph construction followed by Louvain community detection on the resultant graph with a resolution of 0.2, yielding 18 clusters.
Differential expression across timepoints was calculated with MAST 40 (v 1.12.0) to account for interparticipant heterogeneity.
Pseudobulk profiles were constructed by taking the average expression across all cells in each participant, per day. When computing fold changes across timepoints, each participant's pseudobulk profile was compared to their baseline profile to reduce participant-specific biases. To calculate the impact of removing a cluster, each cluster across all timepoints was iteratively removed and resulting fold-changes were recomputed.
C8 was re-embedded and reclustered with UMAP and Louvain community detection, respectively.
Distances from each sub-cluster to the other clusters was calculated as the Euclidean distance between the average expression of all genes of each cluster.
Complexheatmap (v. 2.2.0) was used for all heatmaps. All analysis was performed in R (v 3.6.3)

Data availability
CITE-seq and bulk RNA data are publicly accessible in the Gene Expression Omnibus under accession numbers GSE171964 and GSE169159, respectively.