Main

BNT162b2 has demonstrated 95% efficacy in preventing severe COVID-191. Although adaptive immunity to BNT162b2 has been characterized in humans2,3, little is known about the innate immune response to this vaccine (or to any mRNA vaccine). Systems immunology enables the comprehensive characterization of the cellular and molecular networks that drive innate and adaptive immunity to vaccines and infections4,5,6,7,8. Here we used systems tools to analyse immune responses in 56 volunteers (Extended Data Table 1) who received two doses of BNT162b2 (Extended Data Fig. 1a). More volunteers reported mild side effects after secondary than after primary vaccination (Extended Data Table 2).

Antibody and T cell responses

Primary vaccination induced binding antibody and neutralizing antibody responses in all but three individuals; these responses were boosted significantly after the secondary vaccination (Fig. 1a, b, Extended Data Fig. 1b–d). The neutralizing antibody response reduced about 2-fold by days 90–120, this being comparable to the response to the Moderna mRNA-1273 vaccine9 (Fig. 1b, Extended Data Fig. 1c). There was no gender-associated difference in antibody responses; however, the neutralizing antibody response inversely correlated with age (Extended Data Fig. 1e). Four participants with previous confirmed SARS-CoV-2 infection did not have high baseline titres, but after the first immunization three of these individuals had titres that were greater than 30-fold the geometric mean titres of uninfected individuals following the first dose, but did not increase further after the boost10 (filled black circles in Extended Data Fig. 1b, c). BNT162b2 vaccination also induced a neutralizing antibody response against the B.1.351 variant of concern, albeit at a tenfold-lower magnitude than against the wild-type WA1/2020 (WA1) strain (Fig. 1c, Extended Data Fig. 1f), consistent with previous studies11. The cross-neutralization potential, calculated as a ratio of neutralizing antibody responses between B.1.351 and WA1 strains, also showed a negative association with age (Extended Data Fig. 1g). Vaccination also stimulated spike-specific T cell responses, which were more readily detectable seven days after the secondary immunization (Fig. 1d, e, Extended Data Fig. 2a–e). Consistent with previous studies3, the CD4 T cell responses were primarily of the T-helper-1 type, although there was a low-level T-helper-2 (IL-4) response (Extended Data Fig. 2b). IFNγ and TNF were the dominant responses in CD8 T cells; three individuals with no known exposure to SARS-CoV-2 responded even at baseline, suggestive of CD8 T cells that are cross-reactive to related viruses, as has previously been shown12 (Extended Data Fig. 2d). There was no significant correlation between T cell responses and age or neutralizing antibodies (Extended Data Fig. 2f–k).

Fig. 1: BNT162b2 vaccination induces robust antibody and T cell responses.
figure 1

a, b, Serum anti-S IgG titres (a) (half-maximal effective concentration (EC50)) and live-virus neutralizing antibody titres (b) (half-maximal inhibitory concentration (IC50)), measured by enzyme-linked immunosorbent assay and focus reduction neutralization titre assays, respectively. c, Live-virus neutralizing antibody response against the wild-type (WA1) or B.1.351 variant measured in day-42 sera. d, e, Spike-specific CD4 (d) and CD8 (f) T cell IFNγ response in blood. Each dot represents a participant (n = 56 (a, b), 30 (c) and 38 (d, e)). Statistical differences were calculated using two-sided Wilcoxon matched-pairs signed-rank test. *P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001. Boxes show median and 25th–75th percentiles, and whiskers show the range.

Source data

Full size image

Lack of autoantibodies or anticytokine antibodies

Several studies have demonstrated the presence of serum autoantibodies13,14 and anticytokine antibodies15 in individuals infected with SARS-CoV-2, and the development of new-onset antibodies in a subset of patients who are hospitalized with COVID-1916. We screened sera of 31 vaccinated individuals for IgG autoantibodies and anticytokine antibodies on days 0, 21 and 42 using a 55-plex antigen and 58-plex cytokine arrays. We included sera of 17 patients with autoimmune and immunodeficiency disorders as positive controls. Five vaccinated individuals had pre-existing autoantibodies (suggestive of autoimmune thyroiditis, primary biliary cirrhosis or connective tissue disease) (Extended Data Fig. 3a). Anticytokine antibodies were largely absent or were observed at a low mean fluorescence intensity with no significant changes (Extended Data Fig. 3b). Two individuals had anti-IL-21 autoantibodies, and two additional individuals had anti-IL-1 antibodies (Extended Data Fig. 4). Importantly, none of the individuals with pre-existing autoantibodies or anticytokine antibodies experienced adverse events, nor did levels of pre-existing autoantibodies or anticytokine antibodies change in response to vaccination.

Innate immune responses

We first assessed whole-blood samples of 27 individuals using cytometry by time of flight (CyTOF). Unsupervised clustering identified 14 major cell types (Fig. 2a, Extended Data Fig. 5a, b), which we further subtyped manually (Extended Data Fig. 5c). The frequency of intermediate monocytes (CD14+CD16+ monocytes) increased significantly two days after the primary vaccination, and was substantially higher two days after the secondary vaccination (Fig. 2b, Extended Data Figs. 5d, 6). In addition, there were enhanced levels of phosphorylated (p)STAT3 and pSTAT1 in multiple cell types on day 1 after secondary vaccination, relative to day 1 after primary vaccination (Fig. 2c, d). These data suggested that BNT162b2 vaccination induced a heightened innate immune response after secondary immunization relative to primary immunization.

Fig. 2: Innate immune responses induced by BNT162b2 vaccination.
figure 2

a, CyTOF-identified cell clusters from whole blood visualized by uniform manifold approximation and projection (UMAP). HPCs, haematopoietic progenitor cells; NK, natural killer. b, Frequency of inflammatory monocytes (CD14+CD16+ monocytes) as a proportion of live CD45+ cells. Boxes show median and 25th–75th percentiles, and whiskers show the range. c, Heat map of fold change of pSTAT3 and pSTAT1 levels in comparison to baseline in the cell types indicated on the y axis. Only statistically significant changes between the increase on day 1 after primary and day 1 after secondary immunizations, as measured using two-sided Mann–Whitney rank-sum test (P < 0.05), were plotted. mDC, myeloid dendritic cell; pDC, plasmacytoid dendritic cell; Treg, regulatory T. d, Fold change in pSTAT3 levels on days 1 and 22, compared to primary and secondary baselines, respectively. e, Volcano plots showing plasma cytokines significantly increased after primary (left) and secondary (right) vaccinations versus day 0 and day 21, respectively. f, Plasma IFNγ levels measured by Olink. NPX, normalized protein expression. g, Heat map of two-sided Spearman’s correlation between increase in plasma IFNγ and pSTAT3 or pSTAT1 levels in different cell types, shown on the y axis, on day 1 after secondary vaccination. The P values were corrected for multiple testing. h, Spearman’s correlation between pSTAT3 levels in CD4 T cells and plasma IFNγ levels. The error bands represent 95% confidence limits. In b, f, statistical differences between the peak and baseline time points were measured using two-sided Wilcoxon matched-pairs signed-rank test. The differences between peak time points were measured using two-sided Mann–Whitney rank-sum test. *P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001. Blue and red dots indicate female and male participants, respectively. n = 27 and 31 participants for CyTOF and Olink.

Source data

Full size image

To further investigate this phenomenon, we measured plasma cytokines in 31 vaccinated individuals using Olink (https://www.olink.com/products/target/inflammation/). Of the 67 cytokines detected, the concentration of 2 cytokines (IFNγ and CXCL10) were increased significantly on days 1 and 2 after primary vaccination (Fig. 2e left). Similar to our observations on intermediate monocytes and pSTAT1 and pSTAT3, the concentrations of these cytokines were increased even further after the secondary immunization (Fig. 2e right). The concentration of plasma IFNγ was 11.3-fold higher at day 22 relative to day 1 (Fig. 2f). CXCL10 peaked on day 23 (Extended Data Fig. 7a). The anti-inflammatory cytokine IL-10 showed a similar trend with enhanced response following secondary immunization, although this did not reach statistical significance (Extended Data Fig. 7b). The concentration of IFNα2 (type I interferon) was not significantly increased on days 1 and 2 after vaccination, and the responses were just over the assay limit of the highly sensitive single-molecule array (SIMOA) (16 fg ml−1) (Extended Data Fig. 7c, d). Furthermore, there was a correlation between plasma IFNγ levels and pSTAT1 and pSTAT3 expression levels across several cell types (Fig. 2g, h, Extended Data Fig. 7e). Collectively, these data demonstrate that vaccination with BNT162b2 stimulates modest innate immune responses after primary immunization, which increase notably after the secondary immunization.

Transcriptional signatures of vaccination

Next, we performed bulk mRNA sequencing of whole blood from 32 participants. Four of 185 samples did not pass quality control and were removed from the analysis (Extended Data Fig. 8a, b). Secondary vaccination generated a greater transcriptional response (as has previously been seen in a recent study of adjuvanted hepatitis B vaccine17), with nearly a twelvefold increase of differentially expressed genes at day 22 as compared to day 1 (Fig. 3a); this is consistent with the increased markers of innate immunity detected by both CyTOF and Olink (Fig. 2). Gene-set enrichment analysis (GSEA) revealed that both doses of BNT162b2 stimulated antiviral and interferon response modules18 (Fig. 3b). However, the booster immunization led to a broader innate response. In addition to the induction of antiviral pathways, secondary vaccination led to increases in signatures of dendritic cell activation and the upregulation of Toll-like receptor signalling, monocyte and neutrophil modules on days 22 and 23 (Fig. 3b, Extended Data Fig. 8c, d). The results were consistent regardless of the baseline time point (that is, day 0 or 21) that we used (Extended Data Fig. 8e, f, Supplementary Table 1).

Fig. 3: Transcriptional signatures of BNT162b2 vaccination.
figure 3

a, Number of differentially expressed genes (DEGs) (absolute log2-transformed fold change > 0.2 and Wald P < 0.01) at each time point after prime (left) and boost (right). d, day. b, Interferon and innate BTMs that were significantly enriched (false discovery rate (FDR) < 0.05, absolute normalized enrichment score (NES) > 2) after prime (left) and boost (right) vaccination. Days 1, 2 and 7 were compared against day 0; days 22, 23 and 28 were compared against day 21 in a, b. c, BTMs on day 22 that were significantly associated with plasma IFNγ. GSEA was used to identify enrichment of BTMs within gene lists ranked by correlation with fold change in IFNγ between days 22 and 21.

Source data

Full size image

As mortality to COVID-19 is highest among individuals who are elderly and older populations are known to mount suboptimal responses to many vaccines19, we examined whether there were age-associated differences in response to mRNA vaccination. On day 22, younger participants tended to have greater changes in monocyte and inflammatory modules, whereas older individuals had increased expression of B and T cell modules (Extended Data Fig. 8g). Given that the plasma IFNγ concentration was significantly higher after secondary vaccination, we asked whether there was an association between IFNγ and the increased innate responses after the boost. Both interferon and inflammatory modules were significantly enriched by GSEA when using genes ranked by correlation with IFNγ on day 22 (Fig. 3c). The average fold changes of these modules also correlated with IFNγ (Extended Data Fig. 8h, i), which suggests that IFNγ may have a role in driving enhanced innate and antiviral responses after the boost.

Single-cell transcriptional response

Bulk transcriptomics signatures could reflect changes in cell composition as well as alterations in transcriptional activity within cells. We performed cellular indexing of transcriptomes and epitopes by sequencing (CITE-seq) of 45 peripheral blood mononuclear cell (PBMC) samples from 6 individuals (Extended Data Fig. 1a) to disentangle these two effects and to examine transcriptional changes at the single-cell level. We enriched dendritic cells and mixed them with total PBMCs at a ratio of 1:2 to represent minor subsets of dendritic cells sufficiently4. After quality control, we obtained 242,202 high-quality transcriptomes that segregated into 18 cell clusters (Fig. 4a, Extended Data Fig. 9a–c). Notably, cluster C8 (which expressed CD14, VCAN, CD1C, FCGR1A and CD274 mRNA or protein) emerged on day 22, one day after secondary vaccination (Fig. 4b). These cells uniquely expressed interferon-stimulated genes, including WARS (also known as WARS1), GBP1, IFI30 and IFITM3 (Extended Data Fig. 9d), and constituted about 0.01% of the Lin−HLA-DR+ population on day 1 after primary vaccination but increased almost 100-fold to about 1% one day after secondary vaccination (Fig. 4c). Iterative removal of each cluster from a pseudobulk score showed that C8 contributed to the increased interferon as well as monocyte blood transcriptional modules (BTMs) that we observed in the bulk transcriptomics data on day 22 (Extended Data Fig. 9e). To examine whether these cells uniquely emerge in response to mRNA vaccination or whether they are seen in natural SARS-CoV-2 infection (because they expressed CD14, CD1C and CD274, which are known to be expressed by myeloid-derived suppressor cells observed at elevated frequencies in the blood of patients infected with SARS-CoV-24,20,21,22), we combined and analysed innate immune cells (myeloid cells and plasmacytoid dendritic cells) from this study with data from two previous studies4,23 after batch correction using Harmony24. Well-annotated cell types such as monocytes and dendritic cells overlapped between the datasets, but the C8 cluster did not overlap (Extended Data Fig. 10a, b). The cluster closest to C8 was IFN-experienced CD14+ monocytes, previously defined in patients with COVID-194 (Extended Data Fig. 10b). Although the partial overlap was due to expression of interferon-stimulated genes in both clusters, C8 expressed higher levels of HLA-DR and other activation molecules and lower levels of S100 genes that are known to be expressed in myeloid-derived suppressor cells induced in patients with COVID-19 (Extended Data Fig. 10c). These data suggest that cluster C8 is not present in COVID-19 infection, and does not represent myeloid-derived suppressor cells.

Fig. 4: Single-cell transcriptional response to BNT162b2 vaccination.
figure 4

a, UMAP representation of cell types identified by single-cell transcriptional profiling. b, Feature plots across time points, showing cluster C8 in red. c, Frequency of C8 as a proportion of Lin−HLA-DR+ cell clusters (C3, C4, C7, C8, C11 and C13). Each dot or line represents an individual. Blue and red indicate female and male participants, respectively. d, UMAP representation of subclusters within C8 resolved using Louvain clustering. e, Heat map showing Euclidean distance between C8 subclusters and the rest of the cell types identified by single-cell profiling. f, DEGs determined between C8_2 and C8_1, and their closest cell clusters ranked by fold change between C8 and their parental clusters, are plotted. Genes in red, green and blue represent genes upregulated, unchanged and downregulated, respectively, in the C8 subcluster as compared to its closest parental cluster. g, Heat map showing scaled expression of key interferon-response and AP-1 transcription factors. h, Heat map showing an extended set of genes from ref. 25 that showed enhanced accessibility on day 42 after H5N1 + AS03 vaccination. i, Significantly enriched interferon BTMs (FDR < 0.05, absolute NES > 2) across clusters over time. Days 1 and 7 were compared against day 0; days 22 and 28 were compared against day 21.

Source data

Full size image

To further delineate the cellular composition of C8, we re-embedded C8 with uniform manifold approximation and projection using participant-corrected principal component analyses. Using Louvain community detection, we resolved seven distinct clusters within the original C8 cluster (Fig. 4d). C8 was a heterogeneous mix of classical monocytes (C8_0, C8_1 and C8_3), a classical dendritic cell subtype (cDC2; C8_2) and intermediate monocytes (C8_4), as evidenced by proximity to the original clusters calculated by Euclidean distance (Fig. 4e, Extended Data Fig. 10d–f). We analysed subclusters C8_1 and C8_2, which are closer to CD14+ monocytes and cDC2, respectively. In addition to expressing higher interferon-stimulated genes as compared to their parent clusters, they had a reduced expression of the AP-1 transcription factors FOS and JUN (Fig. 4f). This C8 population is similar to an epigenetically remodelled monocyte population in the blood of humans, 21 days after vaccination with two doses of an AS03-adjuvanted H5N1 pandemic influenza vaccine (H5N1 + AS03)25. These monocytes demonstrated an enhanced chromatin accessibility of interferon-stimulated genes and reduced accessibility of AP-1 transcription factors, and showed heightened resistance to infection with unrelated blood-borne viruses, such as dengue virus and Zika virus25. We asked whether C8 represents an analogous cell type at the transcriptional level. C8 had a relatively higher expression of IRF1, STAT1, STAT2, STAT3 and IRF8 and reduced levels of FOS, JUNB, JUND and ATF3—the same transcription factors that defined the monocyte population in the previous study25 (Fig. 4g). We confirmed this using an extended set of genes for which the chromatin accessibility profile was higher 21 days after H5N1 + AS03 vaccination (Fig. 4h). Notably, the emergence of C8 correlated with plasma IFNγ levels (Extended Data Fig. 10g, h). In vitro stimulation of purified healthy monocytes with IFNγ or day-22 plasma also induced a C8 signature, which suggests that IFNγ has a key role in inducing cluster C8, in response to mRNA vaccination (Extended Data Fig. 10i, j).

In addition to the emergence of C8, our CITE-seq analysis demonstrated that the interferon signature was broadly induced across cell types on day 1 and day 22, and the higher magnitude of response on day 22 was more evident (Fig. 4i). Furthermore, there was activation of natural killer cells on day 22, as observed by downregulation of natural-killer-cell-associated gene modules, and upregulation of AP-1 transcription factors26 on day 22 (Extended Data Fig. 10k, l).

Comparison 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 performed a comparative analysis of a set of published vaccine trials with BNT162b2 (Extended Data Table 3) by generating similarity matrices through pairwise correlations of mean gene fold changes between vaccines at days 1 and 7 after vaccination. The responses at day 1 after the prime and boost doses of BNT162b2 were broadly similar to the response induced by vaccination with adjuvanted vaccines (H5N1 + AS03), live viral vectors (Ebola and HIV vaccines), or inactivated influenza (which stimulates a recall response) (Extended Data Fig. 11a). At the BTM level, the shared signature consisted of innate immunity modules, including interferon signalling, dendritic cell activation and inflammatory responses (Fig. 5). Meanwhile, day-7 responses to both the prime and boost doses of BNT162b2 exhibited weak correlation both between themselves and with other vaccines (Extended Data Fig. 11b), but cell-cycle-related transcriptional modules after both doses were shared with many vaccines (Extended Data Fig. 11c). With most vaccines, this cell cycle signature is also associated with upregulation of B and plasma cell modules, reflecting the expansion of antibody-secreting cells18,27. However, this induction of B and plasma cell modules was absent in BNT162b2 (Extended Data Fig. 11d, e). Given that BNT162b2 successfully promoted a robust antibody response (Fig. 1), the lack of detectable plasma cell or B cell signature on day 7 (particularly after the boost) was surprising. Consistent with this, we observed less than a twofold increase in plasmablast responses by CyTOF (Extended Data Fig. 11f), in contrast to seasonal influenza or other vaccines that induce a higher frequency of plasmablasts27.

Fig. 5: Comparison of transcriptional responses with other vaccines.
figure 5

Circos plot of the overlap across vaccines in BTMs enriched on day 1. GSEA was performed on genes ranked by day 1 versus baseline t-statistic in each vaccine. Each segment of the circle represents one vaccine, and each point in a segment represents a single BTM. Bars in outer circle represent the NES of significantly enriched BTMs (FDR < 0.05). Lines connect BTMs with a significant positive enrichment shared between vaccines. Inner circle boxes and line colours represent the functional groups of the BTMs. ECM, extracellular matrix.

Source data

Full size image

Transcriptional correlates of adaptive immunity

As cellular and humoral immunity are the chief functional components that mediate protection from infection, we used GSEA to identify early transcriptional responses correlated with day-42 neutralizing antibody or day-28 CD8+IFNγ+ T cell responses. On day 22 (one day after the boost), monocyte-related modules correlated with neutralizing antibody responses, whereas interferon and antiviral signatures were associated with the CD8 T cell response (Extended Data Fig. 12a, b). The emergence of SARS-CoV-2 variants poses a serious challenge for the success of ongoing vaccination efforts. We asked whether there are early transcriptional correlates of the cross-neutralization potential induced by BNT162b2. We defined a cross-neutralization index (a ratio of variant-to-WA1 neutralizing antibody titres) and, using GSEA, found that a monocyte module was associated with this index (Extended Data Fig. 12c). In addition, the peak frequency of classical monocytes at two days after the boost correlated with cross-neutralization index (Extended Data Fig. 12d, e). Consistent with this, a gene score that defines C3 (the classical monocyte cluster in the CITE-seq data) also correlated with cross-neutralization (Extended Data Fig. 12f).

Discussion

Our study represents, to our knowledge, the first systems-level analysis of innate and adaptive immunity to an mRNA vaccine. In contrast to the dysregulated innate immune responses to SARS-CoV-2 infection4,20,28,29, BNT162b2 vaccination stimulated antiviral immunity with little type I IFN response after the first dose, but a notably enhanced innate response after the secondary immunization. Cluster C8, defined by single-cell transcriptional profiling, is a heterogenous mixture of myeloid cells that are uniquely induced by mRNA vaccination and are distinct from the IFN-experienced HLA-DRlow myeloid cells observed in natural infection4,20. A previous study has shown that vaccination with H5N1 + AS03 induces a similar transcriptional signature (enriched interferon-stimulated genes and diminished expression of AP-1 transcription factors) in monocytes and myeloid dendritic cells25. In that study, vaccination induced epigenetic reprogramming of these myeloid cells, leading to enhanced resistance against heterologous viruses such as dengue virus and Zika virus, even several weeks later25. In the case of BNT162b2, whether epigenetic reprogramming underlies the enhanced interferon-stimulated gene response in C8 after secondary immunization, and whether this confers enhanced resistance to viruses, remains an open question. Alternatively, it is conceivable that the enhanced myeloid cell response after secondary immunization reflects the response of these cells to a systemic cytokine response. Consistent with this, plasma IFNγ concentration was significantly higher one day after secondary immunization and associated with the emergence of C8. Thus, these data provide a model in which ‘cytokine feedback’ regulates the enhanced innate immune responses to secondary vaccination (Extended Data Fig. 10h–j). Whether or not T cells provide this feedback warrants further investigation, but we did not detect IFNγ within PBMCs (data not shown). Natural killer cells, tissue-resident T cells or innate lymphoid cells at the site of vaccination or draining lymphoid tissues could be potential sources of the rapidly induced circulating IFNγ. This model does not preclude complementary mechanisms, such as persistent epigenetic changes25,30. Finally, our analysis of transcriptional signatures of BNT162b2 vaccination relative to those induced by six other vaccines provides a useful benchmark to assess human immunity to mRNA vaccination in the broader context of immune responses to other vaccines.

Methods

No statistical methods were used to predetermine sample size. The experiments were not randomized, and investigators were not blinded to allocation during experiments and outcome assessment.

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 are provided in Extended Data Table 1.

Anti-spike binding enzyme-linked immunosorbent assay

SARS-CoV-2 spike protein was purchased from Sino Biologicals. Ninety-six-well high binding plates were coated with 100 ng of spike protein diluted at a concentration of 2 μg ml−1 in PBS. The next morning, the plates were washed once, blocked with 3% non-fat milk in PBS containing 0.1% Tween 20 (PBST) for 1 h at room temperature. Sera samples serially diluted in 1% non-fat milk containing PBST were added to the plates and incubated at 37 °C for 1 h. The plates were washed 3× with PBST, horseradish peroxidase conjugated goat anti-monkey IgG (γ-chain specific, Alpha Diagnostics, 1:4,000 dilution), in PBS-T containing 1% non-fat milk was added and incubated for 1 h at room temperature. Wells were washed 3× with PBST before addition of 3,3′,5,5′-tetramethylbenzidine (TMB) substrate solution. The reaction was stopped after 12 min by addition of 0.16 M sulfuric or 1 M hydrochloric acid. The optical density at 450 nanometres was measured with a Biorad microplate reader.

Focus reduction neutralization titre assay

Neutralization assays with authentic SARS-CoV-2 virus (infectious clone SARS-CoV-2 derived from 2019-nCoV/USA_WA1/2020 strain) were performed as previously described31. Sera samples were serially diluted (threefold) in serum-free Dulbecco’s modified Eagle’s medium (DMEM) in duplicate wells and incubated with 100–200 focus-forming units infectious clone derived SARS-CoV-2-mNG virus32 at 37 °C for 1 h. The antibody–virus mixture was added to Vero E6 cell (C1008, ATCC, no. CRL-1586) monolayers seeded in 96-well blackout plates and incubated at 37 °C for 1 h. After incubation, the inoculum was removed and replaced with pre-warmed complete DMEM containing 0.85% methylcellulose. Plates were incubated at 37 °C for 24 h. After 24 h, methylcellulose overlay was removed, cells were washed twice with PBS and fixed with 2% paraformaldehyde in PBS for 30 min at room temperature. Following fixation, plates were washed twice with PBS and foci were visualized on a fluorescence ELISPOT reader (CTL ImmunoSpot S6 Universal Analyzer) and enumerated using Viridot33. The neutralization titres 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 focus reduction neutralization titre (FRNT)–mNeonGreen 50% (mNG50) titres were interpolated using a four-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.

FRNT assay against the variants of concern

The wild-type infectious clone SARS-CoV-2, derived from the 2019-nCoV/USA_WA1/2020 strain, was propagated in Vero E6 cells (ATCC) and sequenced32. The B.1.351 variant was isolated as previously described34. Our laboratory plaque-isolated the virus on VeroE6 cells followed by a single round of propagation on Vero E6 cells (multiplicity of infection of 0.05), aliquoted to generate a working stock and sequenced. Viral titres were determined by focus-forming assay on Vero E6 cells. Viral stocks were stored at −80 °C until use.

FRNT assays were performed as described for the wild-type FRNT assay. The assay with each variant was performed simultaneously with wild-type controls. The samples were diluted at 3-fold in 8 serial dilutions using DMEM in duplicates with an initial dilution of 1:10 in a total volume of 60 μl. Serially diluted samples were incubated with an equal volume of SARS-CoV-2 (wild type or the variant) (100–200 foci per well) at 37 °C for 1 h in a round-bottomed 96-well culture plate. The antibody–virus mixture was then added to Vero cells and incubated at 37 °C for 1 h. After incubation, the antibody–virus mixture was removed and 100 μl of prewarmed 0.85% overlay was added to each well. Plates were incubated at 37 °C for 24 h. After 24 h, methylcellulose overlay was removed, and cells were washed 3 times with PBS. Cells were then fixed with 2% paraformaldehyde in PBS (Electron Microscopy Sciences) for 30 min. Following fixation, plates were washed twice with PBS and 100 μl of permeabilization buffer (0.1% bovine serum albumin (BSA), saponin in PBS), was added to the fixed Vero cells for 20 min. Cells were incubated with an anti-SARS-CoV spike primary antibody directly conjugated to biotin (CR3022–biotin) for 1 h at room temperature. Next, the cells were washed 3 times in PBS and avidin–HRP was added for 1 h at room temperature followed by 3 washes in PBS. Foci were visualized using TrueBlue HRP substrate (KPL, no. 5510-0050) and imaged on an ELISPOT reader (CTL).

Intracellular cytokine staining assay

Antigen-specific T cell responses were measured using the intracellular cytokine staining assay as previously described35. Live frozen PBMCs were revived, counted and resuspended at a density of 2 million live cells per ml in complete RPMI (RPMI supplemented with 10% FBS and antibiotics). The cells were rested overnight at 37 °C in a CO2 incubator. The next morning, the cells were counted again, resuspended at a density of 15 million cells per ml in complete RPMI and 100 μl of cell suspension containing 1.5 million cells was added to each well of a 96-well round-bottomed tissue culture plate. Each sample was treated with two conditions, no stimulation and a peptide pool spanning the spike protein at a concentration of 1 μg ml−1 of each peptide in the presence of 1 μg ml−1 of anti-CD28 (clone CD28.2, BD Biosciences) and anti-CD49d (clone 9F10, BD Biosciences) as well as anti-CXCR3 and anti-CXCR5. The peptides were custom-synthesized to 90% purity using GenScript, a commercial vendor. All samples contained 0.5% v/v DMSO in total volume of 200 μl per well. The samples were incubated at 37 °C in CO2 incubators for 2 h before addition of 10 μg ml−1 brefeldin-A. The cells were incubated for an additional 4 h. The cells were washed with PBS and stained with Zombie UV fixable viability dye (Biolegend). The cells were washed with PBS containing 5% FCS, before the addition of surface antibody cocktail. The cells were stained for 20 min at 4 °C in 100-μl volume. Subsequently, the cells were washed, fixed and permeabilized with cytofix/cytoperm buffer (BD Biosciences) for 20 min. The permeabilized cells were stained with intracellular cytokine staining antibodies for 20 min at room temperature in 1× perm/wash buffer (BD Biosciences). Cells were then washed twice with perm/wash buffer and once with staining buffer before acquisition using the BD Symphony Flow Cytometer and the associated BD FACS Diva software. All flow cytometry data were analysed using Flowjo software v.10 (TreeStar).

Bead-based antigen arrays

We used an existing bead-based autoantigen array, and a cytokine array with expanded content that was based on recent COVID-19 studies16. A complete list of all antigens, vendors and catalogue numbers can be found in Supplementary Tables 2, 3. The ‘COVID-19 Autoantigen Array’ included 55 commercial protein antigens associated with connective tissue diseases (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) such that each antigen was linked to beads with unique barcodes, as previously described16,36,37. Prototype human plasma samples derived from participants with autoimmune diseases with known reactivity patterns were purchased from ImmunoVision or were obtained from Stanford rheumatology clinics and had previously been characterized16. Serum samples from individuals with autoimmune polyendocrine syndrome type 1 (APS1), pulmonary alveolar proteinosis (PAP) or atypical mycobacterial infection (AMI) were used for validation of anticytokine antibodies16. Serum samples were tested at 1:100 dilution in 0.05% PBS-Tween supplemented with 1% (w/v) BSA. Bound antibody was detected using R-phycoerythrin (R-PE) conjugated Fcγ-specific goat anti-human IgG F(ab′)2 fragment (Jackson ImmunoResearch) before analysis using a FlexMap3D instrument (Luminex). A minimum of 100 events per bead identifier were counted, and samples were studied in duplicate. Binding events were displayed as mean fluorescence intensity (MFI). All data analysis and statistics were performed using R and various R packages38. For normalization, average MFI values for ‘bare bead’ identifier were subtracted from average MFI values for antigen-conjugated bead identifiers.

CyTOF analysis of whole-blood samples

Fresh whole-blood samples collected in sodium citrate cell preparation tubes (CPT) were fixed in proteomic stabilizer buffer. Two hundred and seventy μl of whole-blood samples were mixed with 420 μl of smart buffer, mixed and incubated at room temperature for 12 min and frozen at −80 °C until processing. Fixed frozen cells were thawed by gentle resuspension in CSM (PBS supplemented with 2% BSA, 2 mM EDTA and 0.1% sodium azide), washed twice with CSM and counted. Cells were permeabilized and barcoded using Cell-IDTM 20-Plex Pd Barcoding Kit (Fluidigm). The samples were washed with CSM, pooled and counted. One pooled sample containing a mix of all barcoded PBMC samples was stained for 30 min with surface antibody cocktail at room temperature. The sample was then fixed with 4% freshly prepared paraformaldehyde (Alfa Aesar) for 10 min at room temperature, washed with CSM, permeabilized with 100% methanol (Sigma) and kept at −80 °C overnight. The next day, the cells were washed with CSM, counted and stained with pre-titrated intracellular antibody cocktail for 30 min at room temperature. Cells were then washed with CSM, stained with iridium-containing DNA intercalator (Fluidigm), washed with MilliQ water and acquired on Helios mass cytometer (Fluidigm) in MilliQ water supplemented with 1× EQ four element calibration beads (Fluidigm).

The FCS files were bead-normalized before data export. The data were processed for debarcoding in Flowjo software v.10 (TreeStar). In brief, the bead-normalized file was used to gate single cells on the basis of DNA content and event length using FlowJo. The single cells were reimported and debarcoded using Helios software version 7.0.5189. The debarcoded samples were analysed using FlowJo or R version 1.2.1335 for analysis and visualization.

CyTOF data analysis

High-dimensional analysis of phospho-CyTOF data was performed using a previously described R-based pipeline39. In brief, 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 on the basis of 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 quantification 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. The PEA is a dual-recognition immunoassay, in which 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’ that 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).

Bulk transcriptomics

RNA was isolated from blood samples stored in Paxgene tubes at the Yerkes Genomics Core (http://www.yerkes.emory.edu/nhp_genomics_core/). RNA quality was assessed using an Agilent 4200 TapeStation and concentration via the RNA HS assay on the Qubit. Globin transcripts in blood RNA were blocked with the FastSelect Globin Reagent (Qiagen) before library preparation. Libraries were prepared using the Clontech SMART-Seq v.4 Ultra Low Input RNA kit (Takara Bio) in combination with the NexteraXT DNA Library Preparation kit to append dual-indexed adaptor sequences (Illumina). Libraries were validated by capillary electrophoresis on an Agilent 4200 TapeStation, pooled at equimolar concentrations, and sequenced on an Illumina NovaSeq6000 at 100SR, yielding 25–30 million reads per sample.

Bulk transcriptomics analysis

Gene-level counts were filtered to remove those with a median expression less than 32. Principal component analysis (PCA) was performed on baseline samples to identify outliers. Three samples were more than 1.5 s.d. away from the mean and were removed from the analysis. Relative log expression (RLE) plots were generated with EDAseq40; samples with an RLE > 0.6 were removed from the analysis. Differential gene analysis was performed using DESeq241 (v.1.26.0), incorporating participant identifier into the model to account for inter-participant bias. Genes were ranked by the Wald statistic as reported by DESeq2 for GSEA using the BTMs18. Per-participant fold changes were computed by dividing the DESeq2 normalized expression data for the day of interest by either day 0 (for day 1, day2 and day 7) or day 21 (for day 22, 28 and 42). To obtain BTM correlates with age, the age of each participant was compared against the per-participant fold changes for day 22. The resulting correlation values were ranked by t-statistic and analysed with GSEA. The same method was employed to obtain BTM correlates with IFNγ. IFN scores were computed by taking the per-participant mean fold change on day 22 of the unique set of genes present in the 5 interferon BTMs (M75, M111.1, M150, M127 and M68) that significantly correlated with day-22 IFNγ fold change. Similarly, the per-participant M16 gene score was computed using average fold change on day 22 of the genes present in M16.

Vaccine dataset meta-analysis

Datasets were obtained from Gene Expression Omnibus (GEO) via the accession identifiers in Supplementary Table 4. CEL files of all the samples belonging to the same trial were grouped and normalized in Bioconductor by RMA42, which includes global background adjustment and quantile normalization. Probes mapping to multiple genes were discarded, and the remaining probes were collapsed to gene level in each dataset by selecting the probe for each gene with the highest mean expression across all subjects. The only non-microarray dataset was GSE97590, for which the normalized count matrix from GEO was used. Genes not present in all datasets were removed. Baseline normalized log2-transformed fold changes were then computed per subject for all genes. GSEA was then performed to identify enriched BTMs using gene lists for each dataset ranked by t-statistic from two-sided Student’s t-tests on the post-vaccination log2-transformed fold changes.

CITE-seq

CITE-seq analysis of PBMCs were assayed exactly as previously described4. In brief, live frozen PBMCs were thawed and 2× washed with RPMI supplemented with 10% FBS and 20 μg ml−1 DNase I (Sigma Aldrich). Dendritic cells were enriched using the Dynabeads DC Enrichment Kit (Invitrogen, 11308D) according to manufacturer’s instructions with 3–4 million PBMCs as starting material. The enriched cells were mixed with total PBMCs at a ratio of 1:2 and mixed cells were stained with a cocktail of TotalSeq-A antibodies in PBS supplemented with 5% FBS, 2 mM EDTA and 5 mg ml−1 human IgG, washed twice with PBS supplemented with 5% FBS, and 2 mM EDTA, and resuspended in PBS supplemented with 1% BSA (Miltenyi), and 0.5 U μl−1 RNase Inhibitor (Sigma Aldrich). About 9,000 cells were targeted for each experiment.

Cells were mixed with the reverse transcription mix and subjected to partitioning along with the Chromium gel-beads using the 10X Chromium system to generate the gel-bead in emulsions (GEMs) using the 3′ V3 chemistry (10X Genomics). The reverse transcription reaction was conducted in the C1000 touch PCR instrument (BioRad). Barcoded cDNA was extracted from the GEMs by post-GEM reverse transcription cleanup and amplified for 12 cycles. Before amplification, the cDNA amplification mix was spiked in with ADT additive primer (0.2 μM stock) to amplify the antibody barcodes. Amplified cDNA was subjected to 0.6× SPRI beads cleanup (Beckman, B23318). Amplified antibody barcodes were recovered from the supernatant and were processed to generate TotalSeq-A libraries as instructed by the manufacturer (BioLegend, TotalSeq-A antibodies with 10x Single Cell 3′ Reagent Kit v.3 3.1 protocol). The rest of the amplified cDNA was subjected to enzymatic fragmentation, end repair, A tailing, adaptor ligation and 10X-specific sample indexing as per manufacturer’s protocol. Libraries were quantified using Bioanalyzer (Agilent) analysis.

10x Genomics scRNA-seq and TotalSeq-A libraries were pooled and sequenced on an Illumina HiSeq 4000 using the recommended sequencing read lengths of 28 bp (read 1), 8 bp (i7IndexRead) and 91 bp (read 2). CellRanger v.3.1.0 (10xGenomics) was used to demultiplex raw sequencing data and quantify transcript levels against the 10x Genomics GRCh38 reference v.3.0.0.

CITE-seq analysis

10x Genomics scRNA-seq and TotalSeq-A libraries were pooled and sequenced on a Novaseq S4. Cell Ranger v.3.1.0 (10x Genomics) was used to quantify transcript levels against the 10x Genomics GRCh38 reference (v.3.0.0.) Raw count data were filtered to remove cells with a mitochondrial RNA fraction greater than 20% of total RNA counts per cell, cells with fewer than 100 unique features and cells with fewer than 200 total reads. The filtered count matrix was used to create a Seurat43 (v.3.1.4) object. Filtered read counts were scaled by a factor of 10,000 and log-transformed. The antibody-derived tag matrix was normalized per feature using centre log normalization. Doublets were identified with scds44 (v.1.2.0); cells with a doublet score in the top decile were removed. The remaining 242,479 cells were processed with the default Seurat pipeline. Specifically, the most variable 2,000 RNA features were used to perform PCA on the log-transformed counts. The first 25 principal components were used further downstream analyses, including 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 time points was calculated with MAST45 (v.1.12.0) to account for inter-participant heterogeneity.

Pseudobulk profiles were constructed by taking the average expression across all cells in each participant, per day. When computing fold changes across time points, the pseudobulk profile of each participant was compared to their baseline profile to reduce participant-specific biases. To calculate the effect of removing a cluster, each cluster across all time points 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 subcluster to the other clusters was calculated as the Euclidean distance between the average expression of all genes of each cluster. The Euclidean distances were calculated in the original data space. Specifically, the Euclidean distance was calculated using all genes as input to the dist function in R. The dist function calculates Euclidean distance d(x, y) as:

$$d(x,y)=\sqrt{\mathop{\sum }\limits_{i=1}^{n}{({y}_{i}-{x}_{i})}^{2}}$$

in which x is the value for gene i in cluster A and y is the value for gene i in cluster B.

Complexheatmap (v.2.2.0) was used for all heat maps. All analysis was performed in R (v.3.6.3).

Combined analysis of single-cell RNA sequencing

Data from ref. 23 were downloaded from https://covid19cellatlas.org/ as an .h5ad file and converted to a Seurat object in R. Both the resulting Seurat object and the vaccine data were subset to include only myeloid cells and combined using Harmony. Similarly, the data from ref. 4 were integrated with the myeloid cells from the vaccine study using Harmony24. Lymphoid cells from ref. 4 were removed after integration. For both integrations, UMAP was performed on the Harmony-corrected embeddings.

Monocyte purification and stimulation

Monocytes were negatively enriched from healthy PBMCs using Dynabeads Untouched Human Monocyte kit (Invitrogen, cat. no. 11350D) following manufacturer’s instructions. In brief, 50 million live PBMCs were stained with the antibody cocktail for 20 min at 4 °C. The cells were washed and mixed with 0.5 ml of premixed Dynabeads. The samples were incubated in a hulamixer for 15 min at 4 °C. The tubes were placed on the magnet and the unbound fraction containing purified CD14+ monocytes was aspirated using a pipette. The purified monocytes were washed thoroughly and resuspended at a density of 5 million per ml for stimulation. The purity of monocytes was estimated by flow cytometry and was over 95% in all the samples.

Monocytes (0.5 million) were stimulated per condition in 96-well round-bottomed plate for 24 h in 100 μl complete RPMI. Different concentrations of IFNγ, as shown in Extended Data Fig. 7j, were added in 100 μl complete RPMI. Day 0, 1 or 22 plasma samples were from the participant 2055 with <10 pg ml−1, <10 pg ml−1 and 300 pg ml−1 IFNγ, respectively, as measured by enzyme-linked immunosorbent assay. Fifty μl of plasma samples was added to appropriate wells. Fifty μl of complete medium was added to make up the volume to 0.2 ml in total. The plates were incubated at 37 °C, 5% CO2 cell culture incubators.

RNA isolation and qPCR

RNA was isolated using Aurum Total RNA minikit (Biorad, cat. no. 7326820) following the manufacturer’s protocol. cDNA was synthesized using iScript Advanced cDNA synthesis kit (Biorad, cat. no. 1725038) using 150 ng total RNA in 20 μl volume. The cDNA samples were diluted 5-fold by adding 80 μl sterile nuclease-free water and 5 μl of cDNA was used for PCR reaction. The PCRs were carried out using Biorad Prime PCR reagents and SYBR green chemistry (SsoAdvanced Universal SYBR Green Supermix (cat. no. 1725272)) in Biorad CFX384 real-time PCR.

Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this paper.