Systems vaccinology of the BNT162b2 mRNA vaccine in humans

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 vaccine 9 ( 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 boost 10 (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 studies 11 . 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 studies 3 , 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 shown 12 ( Extended Data Fig. 2d ). There was no significant correlation between T cell responses and age or neutralizing antibodies ( Extended Data Fig. 2f – k ).

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. 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.

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 sufficiently 4 . After quality control, we obtained 242,479 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-2 4 , 20 – 22 ), we combined and analysed innate immune cells (myeloid cells and plasmacytoid dendritic cells) from this study with data from two previous studies 4 , 23 after batch correction using Harmony 24 . 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-19 4 ( 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. 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 virus 25 . 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 study 25 ( 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

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 monocyte and inflammatory modules were highly 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 ).

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 .

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 described 31 . 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 virus 32 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 Viridot 33 . 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% (mNG 50 ) titres were interpolated using a four-parameter nonlinear regression in GraphPad Prism 8.4.3. Samples with an FRNT–mNG 50 value that was below the limit of detection were plotted at 10. For these samples, this value was used in fold reduction calculations.

Antigen-specific T cell responses were measured using the intracellular cytokine staining assay as previously described 35 . 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 CO 2 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 CO 2 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).

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.

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).

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 EDAseq 40 ; samples with an RLE > 0.6 were removed from the analysis. Differential gene analysis was performed using DESeq2 41 (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 BTMs 18 . 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.

CITE-seq analysis of PBMCs were assayed exactly as previously described 4 . 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.

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 Harmony 24 . Lymphoid cells from ref. 4 were removed after integration. For both integrations, UMAP was performed on the Harmony-corrected embeddings.

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.

CITE-seq and bulk RNA data are publicly accessible in the GEO under accession numbers GSE171964 and GSE169159 , respectively. Any other relevant data are available from the corresponding authors upon reasonable request. Source data are provided with this paper.

Extended Data Fig. 1 | Antibody responses to BNT162b2 vaccination. a , Schematic of the study design and the number of participants used in various assays. b , c , Binding ( b ) and neutralizing ( c ) antibody responses to BNT162b2 vaccination in female and male participants ( n = 56). d , e , Correlation between binding antibody and neutralizing antibody titres ( d ) and neutralizing antibody responses and age ( e ). f , Neutralizing antibody response to B.1.351 variant of concern in female and male participants ( n = 30). g , Correlation between age and cross-neutralization index, defined as the ratio of neutralizing antibody response to B.1.351 versus WA1 strains. Each dot represents an individual in all plots. Blue and red colour denote female and male participants, respectively. Boxes show median and 25th–75th percentiles, and whiskers show the range in all box plots. All correlations are two-sided Spearman’s correlations. The error bands represent 95% confidence limits. The statistical difference between time points is calculated using two-sided Wilcoxon matched-pairs signed-rank test and the statistical differences between groups were calculated using two-sided Mann–Whitney rank-sum test. *** P < 0.001, **** P < 0.0001. Extended Data Fig. 2 | T cell responses to BNT162b2 vaccination. a , Frequency of spike-specific CD4 T cell responses measured in blood at time points indicated on the x axis. b , Polyfunctional profiles of CD4 T cells. c , Frequency of spike-specific CD8 T cell responses measured in blood at time points indicated on the x axis. d , Polyfunctional profiles of CD8 T cells. e , Frequency of spike-specific CD4 T cells secreting IL-21 and CD154 at time points indicated on x -axis. f , Correlation between spike-specific CD4 (left) and CD8 (right) T cell frequencies and neutralizing antibody responses. g , Correlation between cross-neutralization index, ratio between neutralizing antibody responses against B.1.351 to WA1 strains, and spike-specific CD4 T cell frequencies, IFNγ + (left) or polyfunctional CD4 T cells expressing IL-2, IFNγ and TNF (right). h , Correlation of spike-specific CD4 (left) and CD8 (right) T cell responses with age. i , Correlation of spike-specific IL-21 + CD154 + T follicular helper-like cells on day 28 and neutralizing antibody response on day 42. j , Frequency of CXCR5 + CD4 T cells in PBMCs in DMSO-stimulation condition. k , Correlation of peak (day 7) CXCR5 + CD4 T cells and neutralizing antibody response on day 42. Each dot represents an individual in all plots. Blue and red colour denote female and male participants, respectively. Boxes show median and 25th–75th percentiles, and whiskers show the range in all box plots. The IFNγ response plots in CD4 ( a ) and CD8 ( c ) T cells are from Fig. 1d , e repeated here for completeness. The pie charts in b , d represent the proportion of T cells expressing one, two or three cytokines as shown in the legend. All correlations are two-sided Spearman’s correlations. The error bands represent 95% confidence limits. The statistical difference between time points is calculated using two-sided Wilcoxon matched-pairs signed-rank test. * P < 0.05, ** P < 0.01, *** P < 0.001, **** P < 0.0001. n = 38. Number of samples differ between time points as shown in Extended Data Fig. 1a . Extended Data Fig. 3 | Autoantibodies and anticytokine antibodies in vaccinated individuals. a , Heat map depicting serum IgG antibodies discovered using a 55-plex bead-based protein array containing the indicated autoantigens ( y axis). Autoantigens are grouped on the basis of disease (for example, scleroderma, myositis and overlap syndromes such as mixed connective tissue disease, systemic lupus erythematosus (SLE) and Sjögren’s, and gastrointestinal and endocrine disorders), DNA-associated antigens and antigens associated with tissue inflammation or stress responses. Vaccinated individuals are shown on the left ( n = 30 individuals, on day 0, day 21, and day 42 and n = 1 on day 0 and day 21), and eight prototype autoimmune disorders are shown on the right. Colours correspond to the MFI values shown at far right. b , Heat map using a 58-plex array of cytokines, chemokines, growth factors and receptors. Cytokines are grouped on the y axis by category (interferons, interleukins and other cytokines, growth factors and receptors), and serum samples are shown on the x axis. Vaccinated individuals are shown on the left ( n = 30 individuals, on day 0, day 21 and day 42 and n = 1 on day 0 and day 21). Data are displayed as fold change over baseline. Prototype samples from patients with immunodeficiency disorders are shown on the right, and include three patients with AMI, three patients with PAP and three patients with APS1. Colours for the prototype samples correspond to the MFI values shown at far right. Extended Data Fig. 4 | Pre-existing autoantibodies and autocytokine antibodies do not change in vaccinated individuals. a – f , Bar plo