key: cord-0946815-6h9u2x4u authors: Pulendran, Bali; Arunachalam, Prabhu S title: Systems biological assessment of human immunity to BNT162b2 mRNA vaccination date: 2021-04-22 journal: Res Sq DOI: 10.21203/rs.3.rs-438662/v1 sha: ad93ec988df9e427cdf4e7a263e137f64da77626 doc_id: 946815 cord_uid: 6h9u2x4u 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 in 56 healthy volunteers vaccinated with the Pfizer-BioNTech mRNA vaccine. Vaccination resulted in robust production of neutralizing antibodies (nAbs) against the parent strain and the variant of concern, B.1.351, but no induction of autoantibodies, and significant increases in antigen-specific polyfunctional CD4 and CD8 T cells after the second dose. The innate response induced within the first 2 days of booster vaccination was profoundly increased, relative to the response at corresponding times after priming. Thus, there was a striking increase in the: (i) frequency of CD14+CD16+ inflammatory monocytes; (ii) concentration of IFN- y 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 vaccine responses against pandemic influenza, HIV and Ebola, within 2 days following booster vaccination compared to primary vaccination. Consistent with these observations, single-cell transcriptomics analysis of 242,479 leukocytes demonstrated a ~100-fold increase in the frequency of a myeloid cluster, enriched in a signature of interferon-response transcription factors (TFs) and reduced in AP-1 TFs, one day after secondary immunization, at day 21. Finally, we delineated distinct molecular pathways of innate activation that correlate with CD8 T cell and nAb responses and identified 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 immune responses induced by mRNA vaccines and demonstrate their capacity to stimulate an enhanced innate 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) . 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) . 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. 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. 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) . 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. 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. 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 . 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. The wildtype infectious clone SARS-CoV-2 (icSARS-CoV-2), derived from the 2019- and imaged on an ELISPOT reader (CTL). 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. 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. 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). 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) CITE-seq and bulk RNA data are publicly accessible in the Gene Expression Omnibus under accession numbers GSE171964 and GSE169159, respectively. Safety and Efficacy of the BNT162b2 mRNA Covid-19 Vaccine Phase I/II study of COVID-19 RNA vaccine BNT162b1 in adults COVID-19 vaccine BNT162b1 elicits human antibody and TH1 T-cell responses Systems biological assessment of immunity to mild versus severe COVID-19 infection in humans Systems biology of vaccination for seasonal influenza in humans Systems biology approach predicts immunogenicity of the yellow fever vaccine in humans Poor antigen-specific responses to the second BNT162b2 mRNA vaccine dose in SARS-CoV-2-experienced individuals. medRxiv Neutralizing Activity of BNT162b2-Elicited Serum Pre-existing immunity to SARS-CoV-2: the knowns and unknowns Anti-IFN-gamma autoantibodies in disseminated nontuberculous mycobacterial infections Critically ill SARS-CoV-2 patients display lupus-like hallmarks of extrafollicular B cell activation. medRxiv Broadly-targeted autoreactivity is common in severe SARS-CoV-2 Infection. medRxiv Prothrombotic autoantibodies in serum from patients hospitalized with COVID-19 Mapping Systemic Inflammation and Antibody Responses in Multisystem Inflammatory Syndrome in Children (MIS-C) Diverse Functional Autoantibodies in Patients with COVID-19. medRxiv New-Onset IgG Autoantibodies in Hospitalized Patients with COVID-19. medRxiv Safety and Immunogenicity of Two RNA-Based Covid-19 Vaccine Candidates Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles Molecular signatures of antibody responses derived from a systems biology study of five human vaccines Vaccines for the elderly Immunosenescence and vaccine failure in the elderly: strategies for improving response Singlecell analysis of the epigenomic and transcriptional landscape of innate immunity to seasonal and adjuvanted pandemic influenza vaccination in humans with autoimmune diseases with known reactivity patterns were purchased from ImmunoVision or were obtained from Stanford rheumatology clinics and had been characterized previously 16 . APS-1, PAP, or AMI serum samples were used for validation of ACA 16 Bound antibody was detected using R-phycoerythrin (R-PE) conjugated Fcγ-specific goat anti-human IgG F(ab')2 fragment (Jackson ImmunoResearch) prior to analysis using a FlexMap3D TM instrument (Luminex Corp.). A minimum of 100 events per bead ID were counted, and samples were studied in duplicate All data analysis and statistics were performed using R and various R packages 34 . For normalization, average MFI values for "bare bead" IDs were subtracted from average MFI values for antigen Fresh whole blood samples collected in sodium citrate cell preparation tubes (CPT) were fixed in proteomic stabilizer buffer. 270 μl of whole blood samples were mixed with 420 μl of PROT1 stabilizer After another 10 min at RT, they were centrifuged and washed in 1 ml CSM. They were then permeabilized, barcoded 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) The FCS files were bead-normalized before data export. The data were processed for debarcoding in Briefly, the bead-normalized file was used to gate single cells based on DNA content and event length using FlowJo. The single cells were reimported and debarcoded using dual-indexed bar codes using the NexteraXT DNA Library Preparation kit Libraries were validated by capillary electrophoresis on an Agilent 4200 TapeStation, pooled at equimolar concentrations, and sequenced on an Illumina NovaSeq6000 at 100SR, yielding 20 million reads per Three samples were more than 1.5 standard deviations away from the mean and were removed from the analysis. RLE plots were generated with EDAseq 36 ; samples with an RLE > 0.6 were removed from the analysis. Differential gene analysis was performed using DESeq2 37 (v 1.26.0), incorporating participant id 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 19 . 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). 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 analyzed with GSEA using the BTMs to obtain the BTM correlates with age. The same method was employed to obtain BTM correlates with IFN-g. IFN scores were computed by taking the per-participant mean fold change on day 22 of the unique set of genes 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 human IgG (Extended Data Table XX), 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 RNase Inhibitor (Sigma Aldrich) 10x Genomics scRNA-Seq and TotalSeq-A libraries were pooled and sequenced on an Illumina Cell Ranger v3.1.0 (10x Genomics) was used to quantitate transcript levels against the 10x Genomics GRCh38 reference (v3.0.0.) Raw count data was 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, or cells with fewer than 200 total reads. The filtered count matrix was used to create a Seurat 38 (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 center log normalization. Doublets were identified with scds 39 (v 1.2.0); cells with a doublet score in the top decile were removed. The remaining 242 Data curation and analysis All the authors read and accepted the manuscript Acknowledgements We thank all participants as well as the clinical staff who helped throughout the study. We particularly acknowledge the following group for their contribution to coordination, blood collection and processing of samples Collier from the Stanford Functional Genomics Facility at Stanford University for assistance with single-cell RNA-seq and the Yerkes Nonhuman Primate (NHP) Genomics Core [supported in part by National Institutes of Health (NIH) grant P51 OD011132 ); the Sean Parker Cancer Institute; the Soffer endowment (to B.P.); the Violetta Horton endowment (to B.P.); Next-generation sequencing services were provided by the Yerkes NHP Genomics Core, which is supported in part by NIH P51 OD 011132, and the data were acquired on a NovaSeq 6000 funded by NIH S10 OD 026799 was supported by National Institute of Allergy and Infectious Diseases of the National Institutes of Health We would like to thank the Stanford Biobank (PI, O'Hara) for their processing and storing of samples (funded through NIH CTSA)