key: cord-0686171-6wdd25xr
authors: Talla, Aarthi; Vasaikar, Suhas V.; Lemos, Maria P.; Moodie, Zoe; Lee Pebworth, Mark-Phillip; Henderson, Kathy E.; Cohen, Kristen W.; Czartoski, Julie L.; Lai, Lilin; Suthar, Mehul S.; Heubeck, Alexander T; Genge, Palak C.; Roll, Charles R.; Weiss, Morgan; Reading, Julian; Kondza, Nina; MacMillan, Hugh; Fong, Olivia C.; Thomson, Zachary James; Graybuck, Lucas T.; Okada, Lauren Y.; Newell, Evan W.; Coffey, Ernest M.; Meijer, Paul; Becker, Lynne A.; De Rosa, Stephen C.; Skene, Peter J.; Torgerson, Troy R.; Li, Xiao-jun; Szeto, Gregory Lee; McElrath, M. Juliana; Bumol, Thomas F.
title: Longitudinal immune dynamics of mild COVID-19 define signatures of recovery and persistence
date: 2021-08-19
journal: bioRxiv
DOI: 10.1101/2021.05.26.442666
sha: 8cba55e23bbfee4f4d0073dee64ddb5466318a5f
doc_id: 686171
cord_uid: 6wdd25xr

SARS-CoV-2 has infected over 200 million and caused more than 4 million deaths to date. Most individuals (>80%) have mild symptoms and recover in the outpatient setting, but detailed studies of immune responses have focused primarily on moderate to severe COVID-19. We deeply profiled the longitudinal immune response in individuals with mild COVID-19 beginning with early time points post-infection (1-15 days) and proceeding through convalescence to >100 days after symptom onset. We correlated data from single cell analyses of peripheral blood cells, serum proteomics, virus-specific cellular and humoral immune responses, and clinical metadata. Acute infection was characterized by vigorous coordinated innate and adaptive immune activation that differed in character by age (young vs. old). We then characterized signals associated with recovery and convalescence to define and validate a new signature of inflammatory cytokines, gene expression, and chromatin accessibility that persists in individuals with post-acute sequelae of SARS-CoV-2 infection (PASC).

CD16+ intermediate monocytes, and activation marker-positive B and T cells (Fig. 3B) . CD1c+ DC2s and CD4+ TCMs increased over time, consistent with a return from virus-infected tissues for DC2s (Bosteels et al., 2020; Gill et al., 2008) and expansion or differentiation of CD4+ TCM cells after acute infection. BAFF-R+ class-switched memory B cells and IgG+ B cells increased over time, which suggests recent memory differentiation and class-switching from activated B cells during acute infection . IgA+ plasmablasts also increased over time and may be a consequence of acute infection-induced plasmablasts contracting. Bulk proliferating lymphocytes identified in scRNAseq also significantly declined over time (Fig.  S3B) .

To assess coordination of the immune response, we correlated the change in cell type proportions measured by flow cytometry and scRNAseq with changes in the serum proteome over time (Fig. S3C) . Cell types that positively correlated with the activity of the most serum proteome pathways included CD14+ CD16+ monocytes, Tph, memory B cells, and proliferating lymphocytes (Fig. 3C) . Decreases in these cell type proportions over time were associated with decreased type I and II IFN responses, inflammatory cytokine and immune activation signaling, and PRR signaling. In contrast, these decreases were correlated with increases in BAFF-R+ B cells and serum proteome pathways for chemokines, wound healing, and tissue repair, suggesting temporal coordination between inflammatory resolution and homeostatic repair characterizes successful recovery from mild COVID-19.

We next used scRNAseq and scATACseq data to longitudinally map cell type-specific changes in function over time. Genes demonstrating significant changes over time were identified per cell type, and these changing genes were then mapped to their corresponding pathways. Pathways that showed the top enrichment scores and decreased the most over time, as weighted by number of impacted cell types and participants (see Methods), included TNF signaling, IFN responses, TLR signaling, and apoptosis (Fig. 3D,  Fig. S3D, Supplementary table 3B) . These pathways were consistent with inflammatory protein decreases in serum. Innate immune cells (monocytes, cDC2s, NK subsets) represented the majority of the top cell types expressing these longitudinally changing pathways. scATACseq showed that most transcriptional motifs enriched in acute infection rapidly declined in accessibility over time, as demonstrated by IRF family motif shifts (Fig S3E) . AP-1 motifs were among the most accessible features in innate immune cell types during early acute infection for nearly all COVID-19 participants (Fig. S2H) . Accessibility of AP-1 motifs peaked in CD14+ monocytes during acute infection, followed by a rapid decline to homeostatic baseline by 30 days PSO (Fig.  3E) . PASC participants, especially the sickest PASC, were an exception, maintaining higher accessibility >30 days PSO, which may represent persistent innate immune activation. Principal coordinates analysis (PCoA) on longitudinally changing immune features revealed the individual trajectories of COVID-19 participants from acute infection through convalescence (Fig. 3F) . Early acute infection samples (black circles) showed clear differences from uninfected controls (light green circles, Fig. 3F , top inset left), while most participants around 30 days PSO were closer to uninfected controls (Fig. 3F , bottom inset left). COVID-19 participant trajectories (lighter colored circles and arrows) converged closer to the uninfected control space by ~30 days PSO (Fig.  3F , middle panel). Distances for each sample were calculated from the centroid of the uninfected control group to visualize the timing of inflammatory resolution (Fig. 3F, right panel) . Most participants' trajectories resolved near uninfected controls around day 30, with notable exceptions for the two sickest PASC participants who remained >2.5 distances from the uninfected group. Overlaying comorbidity annotations on the PCoA revealed that the sickest PASC participant with the most comorbidities had the most distinct trajectory from other COVID-19 participants (Fig. S3F) . Altogether, longitudinal analyses across data types demonstrated that most mild COVID-19 participants resolved into a convalescent immune phenotype similar to uninfected controls, while PASC participants have persistent differences and unique trajectories.

We sought to perform deeper immune profiling of PASC participants compared to recovered COVID-19 participants. Recovered COVID-19 participants had no symptoms after ~2 weeks PSO (range: 0-17 days PSO) However, the five participants who were classified as PASC had continued symptoms lasting for > 60 days PSO. All PASC participants were female, similar to the sex bias reported in previous studies (Davis et al., 2020; Evans et al., 2021; Sigfrid et al., 2021) . Two participants who had more severe illness and persistent symptoms were categorized as sickest PASC due to the significant impact of ongoing symptoms on activities of daily living, quality of life and ability to continue employment. These symptoms included neurological symptoms, alterations in smell and taste, chest tightness and joint pain (PTID 795172; illness severity score 71, WHO score 3), and cardiovascular abnormalities, headache, fatigue, and change in taste (PTID 523731; illness severity score 141, WHO score 3) (Supplementary Table 1A) . PASC participants had qualitative differences in CoV-2-specific responses (Fig. 4A) , but none were statistically significant, likely due to the small sample size. There are however trends suggesting that PASC participants may have dysregulated adaptive immune responses. Most notably, the two sickest PASC participants had low or absent SARS-CoV-2-specific CD4+ and CD8+ T cell responses after 30 days PSO, high RBD IgA, and low RBD IgM. PASC participants also had SARS-CoV-2 response estimates below the lower quartile of recovered COVID-19 participants RBDand S-specific memory B cells, and N IgG (Fig. 4A) .

To identify immune pathways and cell types that define PASC participants, we evaluated scRNAseq data from PASC compared to those who had an uneventful recovery from acute SARS-CoV-2 infection. Differentially expressed genes (DEGs) were identified by comparing all PASC participants to recovered COVID-19 participants ≤ 15 days PSO by cell type. CD14+ monocytes, CD8+ TEMRAs and CD8+ effector memory T cells (CD8+ TEM) showed the highest number of DEGs (Fig. S4A) . Mapping DEGs to pathways by cell type revealed dysregulated signaling in early acute infection that persisted longitudinally in PASC. Consistent with evidence of dysregulated immune responses described above, early infection in PASC was characterized by significantly lower IFN responses compared to recovered participants in CD14+ monocytes, CD8+ TEMRAs and CD8+ TEM cells ( Fig. 4B upper panel, Fig. S4B ). These signatures remained persistently low throughout disease, while expression levels in recovered participants were higher during acute infection and declined over time. Recovered participants showed coordinated expression of genes from these pathways, which define an antiviral response (Fig. S4C, Supplementary Table 4A ). CD14+ monocytes from PASC participants also demonstrated significantly higher expression of proinflammatory pathways (IL-1β signaling, TLR signaling) in early acute infection compared to recovered participants, with the highest responses observed in the sickest PASC participants and persisting throughout infection as well as the postacute period (Fig 4B, lower panels) . The up-regulated gene signature included IL1B, CXCL8, HIF1A, S100 family genes (S100A8, S100A9, S100A12), TLRs (TLR2, TLR4, TLR8), and the AP-1 transcription factor subunit (FOS) (Fig. S4D, Supplementary Table 4A ). In contrast, IFN responses in CD14+ monocytes were lower in PASC participants during early acute infection compared to recovered, but these responses also persisted in PASC participants after 30 days PSO (Fig. S4E) . These results suggested PASC is characterized by persistent inflammation compared to timely resolution in successful recovery. In addition, both CD8+ TEMRAs and TEMs had significantly higher gene expression of TNFα signaling via NFκB in PASC compared to recovered participants during early acute infection, with highest elevation in sickest PASC, and this persisted throughout infection (Fig. S4F) .

The presence of an inflammatory gene signature led us to examine whether an inflammatory protein signature also defined PASC. We analyzed the longitudinal proteomes of an expanded PASC cohort including 60 PASC participants and 26 recovered COVID-19 participants. Outlier analysis on the longitudinal serum proteome identified protein signatures that distinguished PASC from recovered COVID-19 participants. Most recovered COVID-19 participants had a decreasing number of differential proteins after early acute infection compared to uninfected controls indicating their path back to normalcy (Fig. S4G) . In contrast, PASC participants had stable or increasing numbers of differentially expressed proteins over time. PASC participants were defined by a persistently elevated inflammatory signature throughout the course of infection. These included significantly up-regulated GO terms like NK cell mediated immunity, RIG-I-like receptor signaling, T cell activation, T cell proliferation, and regulation of innate immune repsonse among others, which were enriched for proteins such as IFNγ, IFNλ1, IL-1β, IL-12, IL-18, and TNF (Fig. 4C, Supplementary Table 4B ). Significantly upregulated proteins in PASC also included members of the TNF superfamily in PASC, including TNFSF13B/soluble BAFF receptor and TNFRSF1B/TNF receptor 2, consistent with immune activation (Supplementary Table 4B ). Overall, analysis of the serum proteome of PASC participants provides further evidence of a hyperinflammatory state.

In an effort to understand whether there are particular immune cells either driving or responding to the signals identified in the serum proteomic studies, we performed transcription factor (TF) motif analysis of scATACseq data. This revealed key TF motifs that correlate with aberrant cell phenotypes in PASC participants compared to COVID-19 recovered participants at >30 days PSO. A set of AP-1 family motifs were significantly enriched in dendritic cells (DCs) and CD14+ monocytes. Paired with evidence of increased expression of the genes encoding the AP-1 subunits FOS and JUN, this suggests a state of persistent immune activation that is most prominent in innate myeloid cells (Fig. 4D) . Other prominent enriched motifs included BACH, IRF, and STAT families, all associated with persistent inflammatory cytokine signaling. The same motifs were also enriched, albeit to a lesser extent, in CD8+ TEMRAs. NK cells specifically showed enrichment of NFκ B motifs. Motif enrichments were not observed across all innate immune cells or for most adaptive immune cell types in PASC. Using scATACseq, we identified genes nearest to differentially accessible PASC peaks with AP-1 motifs and analyzed their expression using scRNAseq data to compare PASC to recovered participants after 30 days PSO. Multiple genes were significantly up-regulated in PASC and these were enriched for inflammatory pathways including TNF signaling, IL1B signaling, TLR signaling, IFN responses, and hypoxia (Fig. S4H) . This agreement between scATACseq and scRNAseq suggests that inflammatory signals activating AP-1 contribute to persistent inflammation in PASC.

To identify potential signals driving the phenotype of CD14+ monocytes in female PASC participants, scRNAseq gene expression changes were used to predict ligand-receptor interactions by the NicheNet method (Browaeys et al., 2020) . Ligand activity was predicted and ranked by the correlation between knowledge-based predictions and experimentally observed levels of target gene expression (Fig. S4I) . Ligand expression was assessed in all cell types along with the fold increase in their expression per cell type compared to recovered participants to identify PBMC sources of ligand signals (Fig. S4J) . TNF, NAMPT and IFNG were among the top predicted ligands with highest target gene activity in CD14+ monocytes (Fig. S4I, S4J, 4E) . Diverse potential sources for these ligands among PBMCs included CD14+ monocytes, dendritic cells, proliferating CD8+/CD4+ T cells and proliferating NK cells (Fig. S4J) . Target genes of these predicted ligands included IL1B, FOS, NAMPT, BCL2A1, and CXCL8 (Fig. 4E, S4K) , which confirm observations in Fig 4B and Fig. S4D . Overall, these data suggested a key role for the cytokine milieu driving a persistently proinflammatory state in CD14+ monocytes, identifying biomarkers for diagnosis of PASC, and revealing multiple therapeutic targets.

These results suggest that PASC is characterized by blunted innate and adaptive antiviral immune responses during acute infection that leads to a persistent inflammatory state including increased expression of proinflammatory genes and proteins. DCs and monocytes demonstrated the most significant increases in chromatin accessibility including motifs for AP-1, STAT, and IRF family transcription factors known to drive inflammatory gene expression. There is clear evidence of persistent inflammation in PASC, but It remains unknown whether these are the primary pathogenic mechanisms causing PASC or merely correlates.

Antibody responses are a key component to controlling viral infection, disease outcome, and protection from reinfection. Clinical trials have directly shown that neutralizing antibodies (nAb) to SARS-CoV-2 can reduce mortality and length of hospitalization . Immune responses in acute infection are critical for establishing the coordinated immune response required for development of antibody and memory B cell responses in convalescence. Serum protein expression and flow cytometry proportions were estimated at day 7 PSO. RBD IgG titers and S-specific memory B cell frequencies at day 30 PSO and neutralizing antibody (nAb) titers at day 42 PSO were estimated by linear mixed-effects models. Observed and estimated values were strongly correlated (Fig. S5A, S5B) .

We identified a signature of circulating proteins, present during early acute infection (day 7 PSO), that strongly correlated to peak antibody titers for RBD IgG and neutralizing antibody (nAb) at days 30 and 42 respectively (Fig. 5A, 5B) . A circulating protein signature was also identified that correlated with spike-specific memory B cell frequencies (Fig S5C) . Correlations were also found between specific immune cell subsets identified by flow cytometry and RBD IgG, nAb, and the S-specific memory B cell frequency (Fig. 5D, 5E; Fig.  S5D ). Analysis of the serum protein milieu at day 7 PSO revealed that higher expression of proteins in inflammatory pathways were significantly correlated with higher RBD IgG titers at day 30 PSO. Among the most prominent were ISGs (TRIM21, SAMD9L), complement (CD55, GPB1A), and innate sensor signaling proteins (TANK, MAVS, RIG-I) (Fig. 5A) . Neutralizing antibody titers at day 42 were also positively correlated with day 7 innate immune proteins (IFNλ1, DDX58), while negative correlates included a subset of inflammatory proteins (IL-34, IL-17A, CD209/DC-SIGN) and chemokines (CCL15, CCL28) (Fig. 5B) . These underscore the importance of innate immunity in early infection for encoding the magnitude of the convalescent antibody response. Sixteen serum proteins were common correlates for both RBD IgG and nAb (Fig. 5C) . Soluble correlates were complemented by positive correlations with flow cytometry populations, including selected activation marker-positive T cells (PD-1+/TIGIT+) and naive B cells (IgD+IgM-), as well as CD56 high NK cells (Fig. 5D) . This links the magnitude of humoral response to the degree of T and B cell activation, potentially through support of B cell maturation to plasma and memory cells. Negative correlates of RBD IgG included pre-switch memory B cells (CD27+IgD+), and CD11b+ B cells. Cell proportions that positively correlated with nAbs included circulating bulk CXCR5+ PD-1+ Tfh cells, and CD56 high and CD16-NK cells. Phenotypic markers of activation on T cells (CD38, CD71, PD-1) and plasmablasts (BCMA) were positively correlated with neutralizing titers. Negative correlates of neutralizing titers included naive CD4+ T cells and inhibitory receptor-positive (KLRG1, TIGIT) CD56 low NK cells (Fig. 5E ). CD56 high NK cells and CD38+ CD8+ T cells were shared positive correlates of both neutralizing antibodies and RBD IgG titers (Fig. 5F ). Coordination of multiple aspects of immunity (lymphocyte activation, innate immune signaling) in early infection are strong predictors of effective priming for IgG and neutralizing antibody responses against SARS-CoV-2.

A critical role for memory B cells has been posited to provide rapid protection upon rechallenge (Dan et al., 2020) . We identified serum proteins and cell frequencies from early infection that correlate with convalescent S-specific IgG+ memory B cell responses. Similar to IgG titers, we observed strong positive correlations with secreted proteins enriched in innate immunity, IFN responses, and RIG-I signaling (Fig. S5C) . In contrast, proteins negatively correlated with memory B cells included complement, chemokines (CCL16, CCL20), and interleukin (IL-1, IL-2) signaling. Cell frequencies that positively correlated with memory B cell response included DC2s and CD25-ILCs, and activated plasmablasts and memory B cells (BCMA) (Fig. S5D) . KLRG1 on non-naive CD4+ T cells and PD-L1 on CD14+CD16+ intermediate monocytes were also positively correlated with memory B cell responses. Overall, these findings identified diverse cellular and molecular features of early infection that correlate with antibody and memory B cell responses, identifying candidates for immunomonitoring and novel linkages between the acute antiviral response and the magnitude of humoral immune responses.

We sought to trace the longitudinal resolution of immune signals that characterized acute infection through longitudinal resolution and convalescence. Several serum inflammatory proteins that were elevated during acute infection compared to uninfected controls resolved over time in recovered participants but showed prolonged elevation in PASC participants, including IFNγ, IFNλ1, IL-12p40, IL6, and TNF among others (Fig.   6A, Fig. S6 ). To infer the origin and signaling activity of these proteins, we performed intercellular communication analysis of scRNAseq data using NicheNet (Methods). We first identified cell type-specific expression of proteins elevated in early acute infection (from Fig. 6A) , and subsequently determined candidate secreting/sender cell types for these signals (Fig. S7A) to define the ligand-cell type interaction network in PBMCs from acute infection (Fig. 6B , 46 predicted ligands in >10% of all cell types). Ten cell types were analyzed as potential receivers of signals based on changes observed in early acute infection (Figs. 2, 3) : plasmablasts, CD14+ monocytes, CD16+ monocytes, NK cells, proliferating cells (CD4, CD8, NK), and DCs (DC1, DC2, pDC). The top 10 high-confidence predicted ligands among these receiver cell types were determined (Fig. 6C) and ranked based on the percent of cell types with predicted ligand activity. These ligands demonstrated the prevalence of inflammatory signaling among cell types during COVID-19 infection, including IFNγ, IFNλ1, IL-12, and TNF (Supplementary Table 5A ; Fig 6B-D) . Multiple downstream target genes were shared between the predicted ligands, including STAT1 (Fig. 6D) , DDIT4 (regulator of mTOR, proliferation, autophagy), and other transcriptional regulators (NFKBIA, JUN, FOS).

Next, we analyzed LR interactions of differentially expressed target genes from 3 distinct comparisons: 1) early acute infection compared to uninfected, 2) longitudinal timepoints (>15 days PSO) compared to early acute infection, and 3) PASC participants (all timepoints) compared to recovered COVID-19. There was significant overlap of LR pairs (467 pairs, ~21%) identified in early acute infection, longitudinal timepoints, and PASC participants (Fig. 6E, Supplementary Table 5B ), indicating common signals driving each subset of mild COVID-19. PASC had the highest number of unique LR interactions (475 pairs, 21%) followed by early acute infection timepoints (404 pairs, ~18%), and longitudinal timepoints (216 pairs, ~10%). Notably, a subset of LR interactions (134, ~6%) were shared between early acute infection and PASC, which may represent signals from acute infection that persist abnormally in PASC. We dissected LR usage by participant subgroups to identify ligands and targets that were more abundantly used in PBMC networks for early acute infection, longitudinal resolution, and PASC (Figs. 6F, G). Multiple innate immune signals and inflammatory signaling were upregulated in early acute infection compared to uninfected controls, and persisted in PASC compared to recovered COVID-19 participants (Fig. 6A, Fig. S6 ). Increased activities of some predicted ligands in PASC participants were consistent with significantly elevated levels in the serum proteome, including IFNγ, IL-6, and IL-12 among others (Fig. 6F, Supplementary Table 5C, Fig. S7B) . A subset of downstream gene targets from predicted ligands was used in at least two subgroups, and primarily consisted of genes coding intracellular proteins: DDIT4, NFKBIA, and CDKN1A (Fig. 6G) .

Overall, data from multiple omics demonstrated that PASC participants showed enrichment for a subset of cytokines, chemokines, transcription factors (STAT1, AP-1 subunits FOS and JUNB), and genes proximal to these open transcription factor sites that were involved in inflammatory cascades such as TNF and IL-1β signaling pathways. The signaling landscape from PASC suggests neutralizing soluble cytokines, such as TNF or IFNγ, or blocking their downstream signaling, after acute infection could provide therapeutic benefit in PASC. Additionally, PASC may be driven by activation of stress response programs as evidenced by upregulated AP-1 and DDIT4 with elevated IL-1β signaling, which could result from an imbalance between autophagy and inflammasome activation. Overall, network analyses provide a platform for generating mechanistic hypotheses on key pathways in mild COVID-19 from acute infection through resolution that may prime adaptive immune responses, mediate successful resolution of infection, and enable diagnosis and treatment of complications such as PASC.

Our study provided an in-depth longitudinal analysis of the immune response to SARS-CoV-2 natural infection and mild COVID-19 by integrating serum proteomics, single-cell transcriptomics and epigenomics, and cellular immunophenotype by flow cytometry with clinical metadata and comprehensive analysis of the SARS-CoV-2-specific adaptive immune response in T cells, memory B cells, and antibodies. To our knowledge, this is the deepest longitudinal systems immunology study to-date in mild COVID-19, and reveals numerous new insights. First, we define the immune response to early acute infection including inflammatory cytokine and innate sensor signaling, stronger IFN responses in older participants, and a potential IFNplasmablast regulatory circuit. A subset of these changes were correlated with the humoral response to SARS-CoV-2 in convalescence. We then confirmed that the longitudinal resolution of these inflammatory pathways was coordinated with re-establishment of homeostasis in most participants. Five PASC participants were exceptions to this resolution, and could be distinguished by dampened IFN and antiviral responses during acute infection coupled with prolonged inflammation. Finally, we integrated these data to identify potential master regulatory nodes in early infection, longitudinal resolution, and PASC for hypothesis generation and validation of biomarkers and therapeutics.

A defining characteristic in our mild COVID-19 cohort was robust immune activation in the first 2 weeks of acute infection that resolved over time. This included inflammatory cytokine responses (IFNs, TNF), innate immune sensor signaling pathways, and activation in adaptive and innate immune cells. The key innate immune sensors triggered in natural SARS-CoV-2 infection are not confirmed, but multiple data types strongly implicate RIG-I. Serum proteomics identified increased extracellular levels of multiple PRR pathway members including RIG-I during acute infection, but their source and function are unclear. These may derive from recent cell death or extracellular vesicles, which have been reported to potentially transfer TLRs between cells (Zhang et al., 2019) . Innate danger sensors are key drivers of the IFN response, which was also robustly induced in our cohort along with other inflammatory signals such as TNF. A subset of upregulated proteins (SAMD9L, CXCL10, CXCL11) is selective to type I/II IFN responses, suggesting bias away from type III IFNs systemically (Allenspach et al., 2021; Forero et al., 2019) . As these pathways waned over time, activation marker-positive cells and inflammatory proteins largely returned to uninfected control levels around day 30 PSO. This temporal control is likely critical for successful resolution of mild COVID-19. This contrasts with persistent CRS reported in severe COVID-19, which includes mechanisms of inflammatory damage to tissue such as TNF/IFNγmediated cell death (Karki et al., 2021) . Proteins involved in homeostatic functions (EMT, coagulation, angiogenesis) increased from acute infection to convalescence. The longitudinal increase of multiple coagulation pathway proteins may contribute to reported increases in risk of immune thrombocytopenia, a complication associated with severe COVID-19 infection (Guan et al., 2020) . We found levels of THBD were significantly increased in convalescence of mild COVID-19, and has been reported to strongly correlate with duration of hospitalization and risk of mortality in hospitalized COVID-19 (Goshua et al., 2020) . These results suggest a link between the inflammatory response in acute infection, the kinetics of inflammatory resolution, and their dysregulation in long-term coagulopathy risk and severe COVID-19.

We observed a clear increase in immune responses with advancing age in mild COVID-19. Age is among the strongest risk factors for severe COVID-19 and mortality, but the mechanisms underlying these effects remain poorly defined (Williamson et al., 2020) . Many studies are confounded by age when comparing uninfected controls and mild COVID-19 cases, which are often younger, with typically older individuals in moderate and severe COVID-19. A male-specific age effect was correlated to poorer CD8+ T cell responses and disease severity with advanced age (Takahashi et al., 2020) . Our cohort was age-matched between COVID-19 participants and uninfected controls, and age was not significantly correlated with illness severity score. However, IFN responses showed a dramatic age-related effect: increased responsiveness of PBMCs from older COVID-19 participants to IFNs, as evidenced by both more pathway enrichment and higher enrichment scores from scRNAseq, and cell type-specific enrichment of IRF motifs in scATACseq. Older participants also showed increased cytokine and danger sensor signaling in innate immune cells. Enhanced IFN responses and TLR signaling in older participants are contrary to expectations from prior studies, where advanced age is associated with impaired IFN and TLR responses and age-associated decrease of TLR responses in DCs (Molony et al., 2017; Panda et al., 2010; Pillai et al., 2016) . These observations could be due to higher viral loads in older participants (Euser et al., 2021; Mahallawi et al., 2021) , cells from older participants being more intrinsically reactive to cytokines, the magnitude of inflammatory cytokine response being higher, or persistence of inflammatory cytokines longer in older participants. We also observed increased adaptive immune cell activation, which contrasts with expectations from prior studies showing immunosenescence in elderly healthy individuals. Differences between our results and prior studies on agerelated effects in immunity may also be due to our cohort including more young and middle-aged adults <55 years old compared to elderly adults typically >65 years old. Collectively, our results indicate that SARS-CoV-2 infection triggers enhanced inflammatory responses in older individuals in mild COVID-19, which may underlie the increased risk for severe COVID-19 in older populations.

Robust plasmablast expansion is a feature of viral infections (dengue, Ebola), vaccines, and chronic autoimmune diseases (Kim et al., 2016; McElroy et al., 2015; Wrammert et al., 2012) . Previous studies in moderate and severe COVID-19 reported robust plasmablast and extrafollicular B cell responses (Bernardes et al., 2020; Kaneko et al., 2020; Kuri-Cervantes et al., 2020; Mathew et al., 2020; Ren et al., 2021; Stephenson et al., 2021; Woodruff et al., 2020) . We also observed a robust plasmablast expansion in acute infection. A fraction of these plasmablasts were CoV-2 spike protein-specific, indicating these cells can in principle contribute to early control of viral replication and are at least partly CoV-2-specific. PD-1 high CXCR5-Tph cells were also positively correlated with plasmablasts, potentially providing help as observed in autoimmune diseases (Rao et al., 2017) . IFN responses are among the key drivers of extrafollicular B cell subsets in autoimmunity (Manni et al., 2018; Soni et al., 2020) , which was consistent with our scRNAseq data showing enhanced IFN signaling and our scATACseq data showing IRF motif enrichments in plasmablasts. The functional consequence of these plasmablasts in SARS-CoV-2 infection remains unclear, but may connect early IFN responses to antibody titers.

Priming of adaptive immunity is critical to successful resolution of acute infection and protection against reinfection. Correlates from acute infection were identified that explain interindividual heterogeneity in magnitude of humoral immune responses. Serum proteins involved in innate immune pathways, including IFN responses, chemokines, and PRR signaling, were positively correlated with RBD IgG titer, neutralizing antibody titer, and CoV-2-specific memory B cell frequency. These were also positively correlated with activation of T and B cells, indicating the importance of coordination in acute infection for an optimal humoral response. Complement proteins were also enriched, consistent with prior studies showing complement can facilitate antigen retention in follicular DCs (Phan et al., 2007) and enhancing BCR-mediated signaling (Fischer et al., 1998; Lyubchenko et al., 2005) . Plasmablasts were a positive correlate of antibodies and memory B cells, but it is unclear whether they play a functional role in clearing infection. Overall, these findings demonstrate that the immune response to acute infection is critical to an effective humoral response, emphasizing the importance of coordination between innate and adaptive arms of immunity.

PASC or long COVID is one of the most enigmatic consequences of the ongoing pandemic. The involvement of many organ systems coupled with the highly subjective nature of symptoms has made it difficult to define objective consensus criteria for diagnosis or clear therapeutic options. In our multi-omics discovery cohort, a subset of five COVID-19 participants progressed to PASC. All PASC participants in this discovery cohort were female, consistent with prior reports of female-biased presentation. Previous studies have suggested that females may have stronger inflammatory responses to vaccines and infections, and predisposition to autoimmune disease (Klein and Flanagan, 2016) ) but it is unclear whether this could contribute to the skewed incidence of PASC in females. Significant correlation was observed between PASC and number of initial symptoms as well as illness severity score during acute infection, as previously reported (Blomberg et al., 2021; Sudre et al., 2021) . Select SARS-CoV-2-specific humoral responses were negatively correlated with PASC, but previous studies report both positive and negative associations of SARS-CoV-2specific antibody titers and T cell responses with PASC (Blomberg et al., 2021; Files et al., 2021; García-Abellán et al., 2021) . Larger studies are likely required to clarify these correlations.

Our expanded PASC cohort was more diverse (60 female, 26 male) and validated the inflammatory protein signatures identified in serum that distinguished PASC participants from recovered COVID-19 participants after acute infection, particularly in those with the most severe symptoms. These signatures were coupled with evidence of persistent activation and inflammatory cytokine signaling in innate immune cells based on single-cell gene expression and chromatin accessibility after 30 days PSO. DCs and CD14+ monocytes in PASC showed the most differences in transcription factor motif accessibility enrichment. Among these, AP-1, STAT, IRF, BATF, and BACH suggest ongoing cellular stress, immune cell activation and differentiation, and inflammatory cytokine signaling during PASC. Gene expression signatures demonstrated increases in immunologically important transcripts from genes proximal to AP-1 sites that were differentially accessible in PASC monocytes. Furthermore, transcriptomics showed stronger TNF signaling in CD14+ monocytes from PASC participants. Early infection signaling and kinetics were also unique in PASC, including lower IFN responses in acute infection that did not wane longitudinally. This combination of changes mirrors severe COVID-19: dampened antiviral responses may fail to control viral replication in both, which can drive innate immune responses to persist beyond acute infection and cause ongoing pathology.

Integrative analysis defined a PASC-specific inflammatory cytokine signaling landscape that suggests multiple therapeutic targets. Single-cell gene expression and serum proteome data provided evidence of persistently elevated IFNγ, IFNλ1, IL-1β, IL-12, IL-6, and TNF in PASC. scATACseq provided evidence of active signaling from these cytokines, with increased accessibility of motifs for key transcription factors in PASC, including IRFs, STATs, and AP-1 in innate immune cell types, primarily DCs and CD14+ monocytes. Cumulatively, these findings suggest blockade of persistent inflammatory cytokine signaling may be therapeutic in PASC. JAK inhibitors are an appealing therapeutic option due to their ability to dampen signaling from several of the cytokines found to be elevated in PASC in addition to their short half-life that facilitates the ability to perform short therapeutic trials to assess efficacy. Similarly, the short half-life and safety of the IL-1 receptor antagonist anakinra could be a strong candidate for early therapeutic trials to probe the role of IL-1 in driving disease. Both JAK inhibitors and IL-1 receptor blockade have been tested in the context of severe COVID-19 immunopathology, with mixed results. Persistently elevated TNF may also be an appealing target given the potential for TNF-driven pathogenic cell death and correlations with COVID-19 disease severity and mortality (Del Valle et al., 2020; Karki et al., 2021) . TNF can also drive IL-1β expression, along with AP-1 and STAT/IRF enrichment observed in motif analyses, motivating therapeutic inhibition of TNF in PASC through TNF blockade. Another promising target is IFNγ, which was significantly elevated in serum from PASC participants compared to recovered participants, coexpressed among pathways in CD14+ monocytes by scRNAseq, and was a strongly predicted ligand signal from ligand-receptor interactions in the majority of cell types. Both TNF and IFNγ are also associated with severe COVID-19, suggesting similarities between severe disease and PASC. Broader anti-inflammatory therapy such as corticosteroids may be useful to attenuate pathogenic inflammation in PASC and has been tested in COVID-19 patients with persistent inflammatory lung disease (Myall et al., 2021) but chronic steroid use is associated with many side effects that make effective targeted therapies more desirable. These novel insights into PASC can focus future studies on pathogenic mechanisms and evaluating efficacy of therapeutic targets. Molecular classification of disease heterogeneity in PASC may identify disease subsets to further enhance clinical management and optimize therapeutic strategies.

There are several key limitations to our study, including 1) small sample size, especially for sickest PASC, 2) lack of geographic and racial diversity, 3) asynchronous sampling due to outpatient status, and 4) gaps in sampling during early infection (first ~7 days PSO) for some participants. Natural history studies are intrinsically limited to correlative associations. Confirming causality from our findings to mechanistically link early infection immune responses to convalescent CoV-2-specific immune responses will require preclinical models. Multi-omic assays were conducted independently from CoV-2-specific assays on parallel samples from each blood draw, so relationships between different -omics and CoV-2 responses were inferred. Direct multiplexed analysis, particularly adding TCR and BCR clonality, will better enable dissection of virus-specific vs. non-specific immune mechanisms. Despite these weaknesses, our study provides insightful data that advance our understanding of SARS-CoV-2 infection, mild COVID-19, convalescence, and PASC.

Overall, our study results provide a comprehensive longitudinal roadmap of immune activation and resolution in mild COVID-19, including a potential mechanism for an age-dependent effect on immune responses. We observed a robust plasmablast response that may be tightly regulated by early IFN responses, and identified key early correlates of antibody and B cell responses, both findings which should be broadly tested as potential shared features in diverse natural infections. A subset of participants who progressed to PASC revealed novel inflammatory and non-inflammatory signatures in serum proteins, and innate immunecentric hyperactivation. Multiple potential therapeutic targets in PASC are nominated by our analyses, and serum protein biomarkers may provide opportunities for objective diagnosis of inflammatory and noninflammatory PASC after validation in larger cohorts. A more personalized approach to immunomonitoring and therapy has the potential to improve outcomes across the spectrum of COVID-19 and PASC.

Cohort overview and participant demographics. Eighteen participants who had an initial blood draw within 15 days of developing symptoms (range 1-15 days) were followed longitudinally. B) Sample availability enumerated per assay type. PBMCs were analyzed by spectral flow cytometry, scRNAseq, scATACseq, and antigen-specific ICS assays. Serum was analyzed for SARS-CoV2 antibody serology or proteomics by Olink Explore 1536. Absent samples were due to limitations on material availability or timing. C) Longitudinal sampling timeline. PBMCs and serum were collected at 3-5 timepoints for each participant. Younger participants (top) are orange; older participants (bottom) are blue. Each Gantt is annotated with sex, documented comorbidities, and presentation of PASC. D) Self-reported symptoms occurring at any point during acute COVID-19 by participant. Illness severity score was calculated as described in Methods and shown vs. WHO ordinal severity. Each symbol represents one participant. Grey symbols indicate recovered COVID-19 participants, orange indicates PASC, and red indicates sickest PASC participants. (E-G) Cohort dynamics of SARS-CoV-2-specific adaptive immune responses: E) Aggregate frequencies of CD4+ and CD8+ T cells specific for E, M, N, S, and ORFs 3, 6a, 7a, 7b, and 8, F) titers of RBD IgG and neutralizing antibodies, and G) frequencies of S-specific of total plasmablasts and S-specific IgG+ of total memory B cells are shown over time for each participant with Loess smoothed curves overlaid and 95% confidence intervals shaded. Symbols for CD4+ and CD8+ T cell frequencies represent positive (filled) or negative (open) responses for each sample based on a positive response to any individual peptide pool stimulation determined using MIMOSA (Methods). H) Correlations between clinical metadata and estimates of SARS-CoV-2-specific immune responses. Day 30 estimates were used for all immune responses except neutralizing antibodies, which were estimated at day 42. Spearman's rank coefficients were calculated for each pairwise combination, multiple comparisons corrected by Benjamini-Hochberg, and only correlations with adjusted p-values < 0.05 are displayed. The size of symbols represents the magnitude of correlations, while color represents magnitude and direction of positive (red) or negative (blue) correlation.

A) Volcano plot showing differential expression of serum proteins between COVID-19 participants at ≤ 15 days PSO (n=18) compared to uninfected controls (n=22). Significantly differential proteins (FDR < 0.05) are colored as upregulated (red) or downregulated (green). B) Flow cytometry analysis shows cell type proportion differences between COVID-19 participants (n=18) compared to uninfected controls (n=23). Parent populations are plotted as % of CD45+ PBMCs, while phenotypic markers are plotted as % of parent. Differential populations were determined by fitting a linear model with change in proportions as a function of infection status adjusted for age and sex. Only significant proportions with adjusted p-values < 0.05 for infection status are shown. C) GSEA performed on scRNAseq using differentially expressed genes from comparison of older COVID-19 participants to younger COVID-19 participants in early acute infection within each cell type. Size of dots represents the percent of genes in each pathway that were enriched. Dot color represents the normalized enrichment score (NES) comparing older COVID-19 participants to younger COVID-19 participants. The top 10 most significant pathway enrichments are shown. D) Expression of genes enriched in IFN response pathways per cell type grouped by participant age and infection status. The column annotation bars represent the average gene expression across single cells per cell type and are scaled across samples. Gene expression was compared between groups using a Wilcoxon rank-sum test; asterisks indicate p < 0.05. E) Transcription factor motif analysis on scATACseq data showing differential accessibility of age-matched COVID-19 participants compared to uninfected controls in plasmablasts. Significantly enriched motifs were colored by motif family (Wilcoxon, FDR < 0.1), and the motifs with the largest shifts in each age subgroup (younger and older participants) are labeled. F) Predicted ligands with signaling activity in plasmablasts during early acute infection. Ligand-receptor interaction analysis was performed using DEGs from scRNAseq as input for NicheNet. The top 10 predicted ligands are shown ranked by Pearson's correlation coefficient (PCC) of target genes. Potential source cell types were identified among PBMCs based on expression of the predicted ligands.

Heatmap of changes in serum proteins over time in mild COVID participants. Proteins that significantly changed over time were identified and enriched for pathways using Fisher's overlap test. Individual columns represent samples ordered by days PSO, with column annotations for days PSO, and recovery status as either PASC (orange), sickest PASC (red), or recovered (grey). Rows represent individual proteins. Features were considered significant with adjusted p-values < 0.05. B) Heatmap of longitudinally changing cell frequencies quantified by flow cytometry among COVID-19 participants. Columns represent individual samples ordered by days PSO, with column annotations for days PSO, and recovery status as either PASC (orange), sickest PASC (red), or recovered (grey). Rows represent significant cell frequencies shown as row-scaled % of total CD45+ cells or % of parent population. C) Pairwise correlations between functional pathways (circles, serum proteome) and cell frequencies (squares, flow cytometry; triangles, scRNAseq) over time. The "antibody response memory B cells" node combines plasmablasts (% of CD45+ cells), and BCMA+, CD71+, and CD86+ (% of post-switch memory B cells) as determined by flow cytometry. The "proliferating lymphocytes" node constitutes CD4, CD8, and NK proliferating cells (% of PBMCs) as obtained by scRNAseq. The "BAFF-R+ B cells" node constitutes BAFF-R+ (% B cells) and BAFF-R+ (% of post-switch memory B cells). The edges to these specific nodes represent the average Spearman correlation coefficient across these cell types to individual serum proteome features. Edges are colored by magnitude of positive (red) or negative (blue) Spearman correlation coefficients. Only significant correlations (adjusted p-values < 0.05) are shown. D) GSEA performed on genes from scRNAseq that were longitudinally decreasing over time per participant by cell type. Dot size represents the percent of COVID-19 participants that showed enrichment for the indicated pathway and cell type. Dot color represents the negative NES indicating the magnitude of decreasing pathway expression over time. The top 10 significant pathway enrichments are shown. E) Longitudinal transcription factor motif accessibility for AP-1 in CD14+ monocytes from COVID-19 participants. Uninfected participants are represented by black solid circles on the left. AP-1 significantly decreased over time from day 10 onwards in the overall cohort (p-value < 0.0001). F) Principal coordinates analysis on all longitudinally changing proteins in COVID-19 participants was used to visualize heterogeneity in the path to convalescence. Early acute COVID-19 samples (left top) and ~30 days PSO (left bottom) were plotted with uninfected controls. The middle panel shows the progression of COVID-19 participants from early (black symbols) to late (red symbols) days PSO. Most COVID-19 participants (solid circles) were well separated from uninfected participants (green triangles) during acute infection, but converged to the uninfected region as they recovered. The right panel shows the corresponding distances to the centroid of uninfected participants (solid black circles). The distances of uninfected participants were plotted with solid circles and randomly assigned days PSO between 0 and 5. The distances of recovered COVID-19 participants were plotted with black solid lines, while those of PASC and sickest PASC participants were plotted with orange and red solid lines, respectively. Loess smoothing (blue line with 95% confidence intervals shaded grey) was used to guide the eye.

SARS-CoV-2-specific adaptive immune responses estimated at day 30, except neutralizing antibodies at day 42, and comparing PASC and recovered COVID-19 participants. Each symbol represents 1 participant. Undetectable estimates are denoted by empty symbols. B) scRNAseq SLEA pathway scores of the significantly different pathways in CD14+ monocytes. GSEA was performed comparing all COVID-19 participants in early acute infection to uninfected controls. The p-values of GSEA were <0.05. The SLEA pathway score per sample was then calculated and tracked over time. C) Volcano plot (left panel) of the differentially expressed serum proteins over time comparing PASC to recovered COVID-19 participants. Each protein's expression was fitted to a linear mixed-effects model. Significantly differential proteins are shown with p-values < 0.05 are shown in red and purple dots. Results of the Gene Set Enrichment Analysis among all differential proteins are represented in the bar plot on the upper right panel, with the top five pathways (ranked by the NES) up-regulated (red) and down-regulated (purple) in PASC compared to recovered participants. Part of the proteins enriched among these pathways are also shown in the volcano plot. D) Transcription factor motif analysis of scATACseq data in samples >30 days PSO. PASC participants were compared to recovered COVID-19 participants and enrichments were calculated using ChromVar motif z-scores. The top 50 largest significant motif differences are shown (Wilcoxon FDR < 0.10). E) Ligand-receptor interaction analysis was used to identify predicted signals driving innate (CD14+ monocyte) cells in PASC participants. The top 5 predicted ligands and their target gene signatures are shown. Ligands were ranked by Pearson correlation coefficient between expected and observed target gene expression. shown. E) The overlap between the inferred ligand-receptor interactions from differential intercellular communication between three subgroup comparisons: 1) early acute COVID-19 participants ≤ 15 days PSO compared to uninfected, 2) longitudinal COVID-19 at >15 days PSO compared to early acute COVID-19, and 3) PASC participants compared to recovered COVID-19 participants at all timepoints. Longitudinal changes in F) ligand and G) target usage in early acute COVID-19, longitudinal and PASC participants respectively are shown. Changes were calculated as the difference in predicted ligand per cell type (F) or downstream target genes regulated by predicted ligands (G) in the subgroup specified compared to respective control conditions listed above. Columns represent each subgroup and differences >10% were shown. Fig. S1 : A) Principal component analysis (PCA) bi-plot of the serum proteome across all longitudinal COVID-19 participant samples (red points) and uninfected controls (green points). The loadings of the most varying proteins are represented by arrows; arrow length is proportional to the variance of the proteins contributing to sample distinction. Shapes of points represent sex and age represents age. A Wilcoxon rank-sum test was used to test if the distribution of principal components 1 and 2 was significantly different between COVID-19 participants and uninfected controls (p-value < 0.05). B) PCA of flow cytometry cell frequencies across all longitudinal COVID-19 participant samples and uninfected controls, including batch bridging controls. Points are colored by batch ID. Principal Variance Component Analysis (PVCA) was performed as indicated in the bar plot to assess batch effects. The proportion of variance (y-axis) in cell frequencies contributed by each factor (x-axis) is depicted. C) PCA bi-plot of the flow cytometry-based proportions across all longitudinal COVID-19 participant samples (red points) and uninfected controls (green points). The loadings of the top varying cell proportions are represented by arrows. A Wilcoxon rank-sum test was used to test if principal components 1 and 2 significantly (p-value < 0.05) differentiate infection status. D) UMAP density representation per batch as grids and PVCA analysis (bar plot) of the scRNAseq batch bridging control samples. The proportion of variance (y-axis) in gene expression contributed by each factor (x-axis) is depicted in the bar plot. E) UMAP density representation per batch as grids and PVCA analysis (bar plot) of the full scRNAseq cohort samples including all COVID-19 longitudinal samples and uninfected controls. The proportion of variance (y-axis) in gene expression contributed by each factor (x-axis) is depicted in the bar plot. F) UMAP of the cell type annotation on scRNAseq data using the Weighted Nearest Neighbors (WNN) approach implemented in Seurat v4 with a CITE-seq based reference dataset (Hao et al., 2020b) . Each point represents single cells color-coded by cell type label. G) Heatmaps of selected gene and predicted ADT (x-axes) expression by cell type (y-axes) of scRNAseq data. H) Correlation scatter plots between scRNA proportions (x-axis) and the equivalent population gating in flow cytometry (y-axis) across all samples. Significance of correlations was tested by a Spearman correlation test and only significant correlations are shown.

Box and jitter plots of scRNAseq-based cell type proportions that were significantly different (adjusted p-values < 0.05) comparing the early acute infection timepoints of COVID-19 participants (red) to uninfected controls (green) assessed by fitting a linear model adjusted for age and sex. B) Heatmap of the flow cytometry and scRNAseq cell type proportions (rows) that were significantly different (adjusted p-values < 0.05 assessed by Wilcoxon rank-sum test) comparing older (including middle aged and senior adults) to younger COVID-19 participants (columns) in early acute infection. C) Heatmap of serum proteins (rows) that were significantly different (adjusted p-values < 0.05 assessed by Wilcoxon rank-sum test) comparing older to younger COVID-19 participants (columns) in early acute infection. D) Bar plot of the number of DEGs (y-axis) comparing age-matched early acute COVID-19 participants to uninfected controls (per grid) within each cell type (x-axis) detected in scRNAseq. Number of DEGs that were significantly up-regulated (red) and downregulated (blue) at an adjusted p-value < 0.05 are reported. E) GSEA pathway enrichment (y-axis) among DEGs comparing age matched early acute COVID-19 participants to uninfected controls (per grid) in every cell type (x-axis) of scRNAseq data. The color gradient of points represents the NES per pathway indicating upregulation (red) or down-regulation (blue) of each pathway. The size of points represents the percent of genes enriched in a pathway per cell type and age group. F) GSEA pathway enrichment (y-axis) among DEGs comparing early acute infection timepoints of older COVID-19 participants to younger in every cell type (x-axis) of scRNAseq data. The color gradient of points represents the positive NES per pathway indicating upregulation of the pathways in older participants. The size of points represent the percent of genes enriched in a pathway per cell type and age group. Pathways are arranged as being expressed in most to the least number of cell types. G) Scatterplot of the -log10 of the lowest P value per cell type by the absolute value of the maximum z-score shift per cell type. CD16+ monocytes, CD14+ monocytes, DCs, and plasmablasts appear to be the most impacted during the early infection period. H) Heatmap of top 50 motif accessibility changes during early COVID19 infection. Z-scores representing motif accessibility as compared to background were calculated for each cell. Z-score shifts were tested for statistical significance by a Wilcoxon rank-sum test with an FDR < 0.1. Only significant shifts are colored.

participants that significantly (adjusted p-values < 0.05) decreased over time. Individual lines are color-coded by the participant ID. B) scRNA based cell type proportions (individual grids) among COVID-19 participants that significantly (adjusted p-values < 0.05) showed decrease over time. Individual lines are color coded by the participant ID. C) Pairwise correlations between functional pathways (circles, serum proteome) and cell type proportions (squares, flow cytometry; triangles, scRNAseq) overtime. The edge color gradient indicates significant (p-values < 0.05) positive (red) or negative (blue) Spearman correlations over time. D) GSEA pathway enrichment (y-axis) among longitudinal gene changes per COVID-19 participant in every cell type (xaxis) of scRNAseq data. The size of points represents the percent of COVID-19 participants the pathway is enriched in. Pathways are arranged as being expressed in most to the least number of cell types across participants. E) Longitudinal IRF motif accessibility per participant over disease course. Curves are colored for PASC and sickest PASC (orange and red respectively) and recovered (black) COVID-19 participants. Loess curves were added to visualize group trends. Uninfected participants are represented by black solid circles on the left. F) PCoA of serum proteome trajectories for individual COVID-19 participants and uninfected controls (black circles, no arrows). Symbols are colored by comorbidities and clinical features. Arrows are colored by visit interval over time.

Bar plot of the number of DEGs (y-axis) comparing early acute infection timepoints of PASC participants to recovered COVID-19 participants in every cell type (x-axis) detected in scRNAseq. Number of DEGs that were significantly up-regulated (red) and down-regulated (blue) at an adjusted p-value < 0.05 are reported. B) Box plot of the SLEA score (y-axis) of the interferon response pathways in scRNAseq significantly downregulated (adjusted p-values < 0.05) in CD8+ TEMRAs and TEM from PASC compared to recovered COVID-19 participants in early acute infection. The SLEA score of the pathways per sample was calculated and tracked over the course of disease (x-axis). C) Arc co-expression network of genes enriched in the interferon responses pathway in CD14+ monocytes, from early acute timepoints of recovered participants. Edges represent correlation between genes at greater than 0.8 Spearman correlation coefficient, while size of each gene node represents the degree or number of gene correlations above mentioned coefficient. D) Arc coexpression network of genes enriched in the IL1B and TLR signaling pathways in CD14+ monocytes, from early acute timepoints of PASC participants. Edges represent correlation between genes at greater than 0.8 Spearman correlation coefficient, while size of each gene node represents the degree or number of gene correlations above mentioned coefficient. E) Bar plot of the pathways (y-axis) enriched in CD14+ monocytes comparing PASC to recovered participants in early acute infection (red bars) and beyond 30 days PSO (black bars). X-axis represents the GSEA NES, with a positive NES indicating up-regulation in PASC, while negative NES indicating down-regulation in PASC compared to recovered. F) Box plot of the SLEA score (y-axis) of TNF signaling via NFKb pathway significantly upregulated (adjusted p-values < 0.05) in CD8+ TEM and TEMRAs from PASC compared to recovered COVID-19 participants in early acute infection. The SLEA score of the pathways per sample was calculated based on scRNAseq gene expression per cell type and tracked over the course of disease (x-axis). G) Number of differentially expressed serum proteins (y-axis) over time (xaxis) compared to uninfected controls for each COVID-19 participant. Outlier analysis was performed comparing COVID-19 participants to uninfected background and selecting features >2 standard deviations from the mean. Recovered, PASC and the sickest PASC COVID-19 participants are colored by grey, orange and red respectively. H) Heatmap of the Gene Set Enrichment Analysis results in scRNAseq beyond 30 days PSO among a curated AP-1 target gene set. These genes were identified as those nearest to differentially accessible scATAC-seq peaks with AP-1 motifs. Columns represent pathways while rows represent the genes enriched among pathways. The inset on the right shows genes unique to the "innate immune system" pathway. I) Ligand activity observed in CD14+ monocytes and their expression in sender cell types of PASC participants, J) dot plot showing expression of these ligands among PBMCs cell types, where the size of circles indicates the percentage of cells expressing the gene for each ligand and the color scale indicates the average expression level, and K) target genes of predicted ligands to CD14+ monocytes. ≥ 31 days. Wilcoxon rank-sum tests were used to compare the PASC to recovered participants across all timepoints, and a representative subset of significantly differential proteins are shown. P-values < 0.05 were considered statistically significant.

Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contacts: Gregory Szeto (greg.szeto@alleninstitute.org); M. Juliana McElrath (jmcelrat@fredhutch.org); and Thomas F. Bumol (tomb@alleninstitute.org).

This study did not generate new unique reagents.

The RNAseq and ATACseq data generated during this study are available at GEO under accession number GSE173590. All analytical code and figure generation will be available on Github.

Eighteen immunocompetent study participants who were diagnosed with SARS-CoV-2 were selected for these studies, as they had donated at least one sample of serum at or earlier than 15 days post-symptom onset. Both studies were recruited at the Seattle Vaccine Trials Unit (Seattle, Washington, USA). The cohort demographics are described in the following 

Blood collected in SST tubes was allowed to clot for at least 30 min at room temperature. Following cotting, the samples were centrifuged in the primary tubes at 1600-1800 X for 10 min. Serum was transferred to a 15ml tube, aliquoted into 300ul aliquots and stored at -80°C until use in the assays.

Blood collected in acid citrate dextrose tubes was transferred to Leucosep tubes (Greiner Bio One). The tube was centrifuged at 800-1000g for 15 min and the PBMC layer recovered above the frit. Peripheral blood mononuclear cells (PBMCs) were washed twice with Hanks Balanced Solution without Ca+ or Mg+ (Gibco) at 200-400xg for 10 min, counted, and aliquoted in heat-inactivated fetal bovine serum with 10% dimethylsulfoxide (DMSO, Sigma) for cryopreservation. PBMCs were cryopreserved at -80°C in Stratacooler (Nalgene) and transferred to liquid nitrogen for long-term storage.

To assess cell type proportions, PBMCs were analyzed with four 25-color immunophenotyping flow cytometry panels (P1, P2, P3 and P4). 1×10 6 thawed PBMCs were centrifuged (750×g for 5 minutes at 4°C) using a swinging bucket rotor (Beckman Coulter Avanti J-15RIVD with JS4.750 swinging bucket, B99516), the supernatant was removed using a vacuum aspirator pipette, and the cell pellet resuspended in DPBS without calcium and magnesium (Corning 21-031-CM 

Flow cytometry was used to examine SARS-CoV-2-specific CD4+ and CD8+ T-cell responses using a validated ICS assay. The assay was similar to a published report (Dintwe et al., 2019; Horton et al., 2007) and the details of the staining panel are included in Reagents Table 2 . Peptide pools covering the structural proteins of SARS-CoV-2 were used for the six-hour stimulation. Peptides matching the SARS-CoV-2 spike sequence (316 peptides, plus 4 peptides covering the G614 variant) were synthesized as 15 amino acids long with 11 amino acids overlap and pooled in 2 pools (S1 and S2) for testing (BioSynthesis). All other peptides were 13 amino acids overlapping by 11 amino acids and were synthesized by GenScript. The peptides covering the envelope (E), membrane (M) and nucleocapsid (N) were initially combined into one peptide pool, but the majority of the assays were performed using a separate pool for N and one that combined only E and M. Several of the open reading frame (ORF) peptides were combined into two pools, ORF 3a and 6, and ORF 7a, 7b and 8. All peptide pools were used at a final concentration of 1 microgram/ml for each peptide. As a negative control, cells were not stimulated, only the peptide diluent (DMSO) was included. As a positive control, cells were stimulated with a polyclonal stimulant, staphylococcal enterotoxin B (SEB). Cells expressing IFNγ and/or IL-2 and/or CD154 were the primary immunogenicity endpoint for CD4+ T cells and cells expressing IFNγ were the primary immunogenicity endpoint for CD8+ T cells. The overall response to SARS-CoV-2 was defined as the sum of the background-subtracted responses to each of the individual pools. A sample was considered positive for CD4+ or CD8+ T cell responses to SARS-CoV-2 if any of the CD4+ or CD8+ T cell responses to the individual peptide pool stimulations was positive. Positive responses to a given peptide pool stimulation were determined using the MIMOSA (Mixture Models for Single-Cell Assays) method (Finak et al., 2014a) . The MIMOSA method uses Bayesian hierarchical mixture models that incorporate information on cell count and cell proportion to define a positive response by comparing peptide-stimulated cells and unstimulated negative controls. MIMOSA estimates the probabilities that peptide-stimulated responses are responders and applies a false-discovery rate multiplicity adjustment procedure (Newton et al 2004) . Responses with false-discovery rate q-values < 0.05 were considered positive.The total number of CD4+ T cells must have exceeded 10,000 and the total number of CD8+ T cells must have exceeded 5,000 for the assay data to be included in the analysis.

Half-well area plates (Greiner) were coated with purified RBD protein at 16.25ng/well in PBS (Gibco) for 14-24h at room temperature. After 4 150ul washes with 1X PBS, 0.02% Tween-2 (Sigma) using the BioTek ELx405 plate washer, the IgA and IgG plates were blocked at 37°C for 1-2 hours with 1X PBS, 10% non-fat milk (Lab Scientific), 0.02% Tween-20 (Sigma); IgM plates were blocked with 1X PBS, 10% non-fat milk, 0.05% Tween-20.

Serum samples were heat inactivated by incubating at 56°C for 30 minutes, then centrifuged at 10,000 x g / 5 minutes, and stored at 4°C previous to use in the assay. For IgG ELISAs, serum was diluted into blocking buffer in 7-12 1:4 serial dilutions starting at 1:50. For IgM and IgA ELISAs, serum was diluted into 7 1:4 serial dilutions starting at 1:12.5 to account for their lower concentration. A qualified pre-pandemic sample (negative control) and a standardized mix of seropositive serums (positive control) was run in each plate and using to define passing criteria for each plate. All controls and test serums at multiple dilutions were plated in duplicate and incubated at 37°C for 1 hour, followed by 4 washes in the automated washer. 8 wells in each plate did not receive any serum and served as blocking controls.

Plates then were plated with secondary antibodies (all from Jackson ImmunoResearch) diluted in blocking buffer for 1h at 37C. IgG plates used donkey anti-human IgG HRP diluted at 1:7500; IgM plates used goat anti-human IgM HRP diluted at 1:10,000; IgA plates used goat anti-human IgA HRP at 1:5000. After 4 washes, plates were developed with 25ul of SureBlock Reserve TMB Microwell Peroxide Substrate (Seracare) for 4 min, and the reaction stopped by the addition of 50ml 1N sulfuric acid (Fisher) to all wells. Plates were read at OD 450nm on SpectraMax i3X ELISA plate reader within 20 min of adding the stop solution.

OD 450nm measurements for each dilution of each sample were used to extrapolate RBD endpoint titers when CVs were less than 20%. Using Excel, endpoint titers were determined by calculating the point in the curve at which the dilution of the sample surpassed that of 5 times the average OD 450nm of blocking controls + 1 standard deviation of blocking controls.

N IgG was measured using SARS CoV-2 IgG Architect (Abbot), and index values used as a quantitative measure of the N binding activity.

VeroE6 cells were obtained from ATCC (clone E6, ATCC, #CRL-1586) and cultured in complete DMEM medium consisting of 1x DMEM (VWR, #45000-304), 10% FBS, 25mM HEPES Buffer (Corning Cellgro), 2mM L-glutamine, 1mM sodium pyruvate, 1x Non-essential Amino Acids, and 1x antibiotics. The infectious clone SARS-CoV-2 (icSARS-CoV-2-mNG), derived from the 2019-nCoV/USA_WA1/2020 strain, was propagated in VeroE6 cells (ATCC) and sequenced (Xie et al., 2020) .

Neutralization assays with SARS-CoV-2 virus were performed as previously described (Vanderheiden et al., 2020b) . Plasma/serum were serially diluted (three-fold) in serum-free Dulbecco's modified Eagle's medium (DMEM) in duplicate wells and incubated with 100-200 FFU infectious clone derived SARS-CoV-2-mNG virus at 37°C for 1 hr. The antibody-virus mixture was added to VeroE6 cell (C1008, ATCC, #CRL-1586) monolayers seeded in 96-well blackout plates and incubated at 37°C for 1 h. Post-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 (Katzelnick et al., 2018) . 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-mNG 50 titers were interpolated using a 4-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 20.

Sample preparation, hashing, and pooling: Single-cell RNA-seq libraries were generated using the 10x Genomics Chromium 3' Single Cell Gene Expression assay (#1000121) and Chromium Controller Instrument according to the manufacturer's published protocol with modifications for cell hashing (Stoeckius et al., 2018) .

To block off-target antibody binding, Blocking Solution (5 µL of Human Trustain FcX (BioLegend #422302) , and13.7 µL of a 10% Bovine Serum Albumin (BSA)) was added to 500,000 cells suspended in 50 µL Dulbecco's Phosphate Buffered Saline (DPBS; Corning Life Sciences #21-031-CM) and incubated for 10 minutes on ice. To stain samples, 0.5 µg (1 µL) of a TotalSeq™-A anti-human Hashtag Antibody was suspended in 31.3 µL DPBS/2% BSA, then added to each sample with . For each batch of samples, 100,000 cells from 8 hashedsamples with a distinct Hashtag Antibody were pooled into Pool 1; 8 additional samples were pooled using the same method into Pool 2. Roughly 20,000 cells from a Leukopak healthy control were also labeled with a distinct TotalSeq™-A Hashtag Antibody, and were spiked into each pool to serve as a batch control.

Droplet encapsulation and reverse transcription: From each pool, 64,000 cells were loaded into each well of a Chromium Single Cell Chip G (10x Genomics #1000073) (8 wells per chip), targeting a recovery of 20,000 singlets from each well. Gel Beads-in-emulsion (GEMs) were then generated using the 10x Chromium Controller. The resulting GEM generation products were then transferred to semi-skirted 96-well plates and reverse transcribed on a C1000 Touch Thermal Cycler programmed at 53°C for 45 minutes, 85°C for 5 minutes, and a hold at 4°C. Following reverse transcription, GEMs were broken and the pooled single-stranded cDNA and Hashtag Oligo fractions were recovered using Silane magnetic beads (Dynabeads MyOne SILANE #37002D).

Barcoded, full-length cDNA including the Hashtag Oligos (HTOs) from the TotalSeq™-A Hashtag Antibodies were then amplified with a C1000 Touch Thermal Cycler programmed at 98°C for 3 minutes, 11 cycles of (98°C for 15 seconds, 63°C for 20 seconds, 72°C for 1 minute), 72°C for 1 minute, and a hold at 4°C. Amplified cDNA was purified and separated from amplified HTOs using a 0.6x size selection via SPRIselect magnetic bead (Beckman Coulter #22667) and a 1:10 dilution of the resulting cDNA was run on a Fragment Analyzer (Agilent Technologies #5067-4626) to assess cDNA quality and yield. HTO libraries were purified further with SPRIselect magnetic bead (Beckman Coulter #22667) and amplified and indexed with a custom HTO i7 index on a C1000 Touch Thermal Cycler programmed at 95°C for 3 minutes, 10 cycles of (95°C for 20 seconds, 64°C for 30 seconds, 72°C for 20 seconds), 72°C for 1 minute, and a hold at 4°C. The resulting HTO libraries were purified with SPRIselect magnetic bead (Beckman Coulter #22667) postamplification and a 1:10 dilution of the resulting HTO libraries were run on a Fragment Analyzer (Agilent Technologies #5067-4626) to assess HTO quality and yield. A quarter of the cDNA sample (10 ul) was used as input for library preparation. Amplified cDNA was fragmented, end-repaired, and A-tailed is a single incubation protocol on a C1000 Touch Thermal Cycler programmed at 4°C start, 32°C for minutes, 65°C for 30 minutes, and a 4°C hold. Fragmented and A-tailed cDNA was purified by performing a dual-sided size-selection using SPRIselect magnetic beads (Beckman Coulter #22667). A partial TruSeq Read 2 primer sequence was ligated to the fragmented and A-tailed end of cDNA molecules via an incubation of 20°C for 15 minutes on a C1000 Touch Thermal Cycler. The ligation reaction was then cleaned using SPRIselect magnetic beads (Beckman Coulter #22667). PCR was then performed to amplify the library and add the P5 and indexed P7 ends (10x Genomics #1000084) on a C1000 Touch Thermal Cycler programmed at 98°C for 45 seconds, 13 cycles of (98°C for 20 seconds, 54°C for 30 seconds, 72°C for 20 seconds), 72°C for 1 minute, and a hold at 4°C. PCR products were purified by performing a dual-sided size-selection using SPRIselect magnetic beads (Beckman Coulter #22667) to produce final, sequencing-ready libraries.

Quantification and sequencing: Final libraries were quantified using Picogreen and their quality was assessed via capillary electrophoresis using the Agilent Fragment Analyzer HS DNA fragment kit and/or Agilent Bioanalyzer High Sensitivity chips. Libraries were sequenced on the Illumina NovaSeq platform using S4 flow cells. Read lengths were 28bp read1, 8bp i7 index read, 91bp read2.

To remove dead cells, debris, and neutrophils prior to scATAC-seq as described previously (Swanson et al., 2021) , PBMC samples were sorted by fluorescence activated cell sorting (FACS) prior to cell permeabilization. Cells were incubated with Fixable Viability Stain 510 (BD, 564406) for 15 minutes at room temperature and washed with AIM V medium (Gibco, 12055091) plus 25 mM HEPES before incubating with TruStain FcX (BioLegend, 422302) for 5 minutes on ice, followed by staining with mouse antihuman CD45 FITC (BioLegend, 304038) and mouse anti-human CD15 PE (BD, 562371) antibodies for 20 minutes on ice. Cells were washed with AIM V medium plus 25 mM HEPES and sorted on a BD FACSAria Fusion. A standard viable CD45+ cell gating scheme was employed; FSC-A x SSC-A (to exclude sub-cellular debris), two FSC-A doublet exclusion gates (FSC-W followed by FSC-H), dead cell exclusion gate (BV510 LIVE/DEAD negative), followed by CD45+ inclusion gate. Neutrophils (defined as SSC high , CD15 + ) were then excluded in the final sort gate. An aliquot of each post-sort population was used to collect 50,000 events to assess post-sort purity.

Sample preparation: Permeabilized-cell scATAC-seq was performed as described previously (Swanson et al., 2021) . A 5% w/v digitonin stock was prepared by diluting powdered digitonin (MP Biomedicals, 0215948082) in DMSO (Fisher Scientific, D12345), which was stored in 20 µL aliquots at −20°C until use. To permeabilize, 1×10 6 cells were added to a 1.5 mL low binding tube (Eppendorf, 022431021) and centrifuged (400×g for 5 min at 4°C) using a swinging bucket rotor (Beckman Coulter Avanti J-15RIVD with JS4.750 swinging bucket, B99516). Cells were resuspended in 100 µL cold isotonic Permeabilization Buffer (20 mM Tris-HCl pH 7.4, 150 mM NaCl, 3 mM MgCl2, 0.01% digitonin) by pipette-mixing 10 times, then incubated on ice for 5 min, after which they were diluted with 1 mL of isotonic Wash Buffer (20 mM Tris-HCl pH 7.4, 150 mM NaCl, 3 mM MgCl2) by pipette-mixing five times. Cells were centrifuged (400×g for 5 min at 4°C) using a swinging bucket rotor, and the supernatant was slowly removed using a vacuum aspirator pipette. Cells were resuspended in chilled TD1 buffer (Illumina, 15027866) by pipette-mixing to a target concentration of 2,300-10,000 cells per µL. Cells were filtered through 35 µm Falcon Cell Strainers (Corning, 352235) before counting on a Cellometer Spectrum Cell Counter (Nexcelom) using ViaStain acridine orange/propidium iodide solution (Nexcelom, C52-0106-5). 

To confirm some of the inflammatory signals observed in the cohort of 19 individuals (including 5 PASC), we selected additional PASC individuals from the longitudinal cohort regardless of their availability of serum collections within 15 days of diagnosis. We also included additional follow-up visits for the original PASC participants. All pre-vaccine visits from 55 PASC who had confirmed symptoms for more than 60 days PSO were included in the extended analysis. As controls, a group of 12 non-PASC recovered individuals were included in the analysis and they were selected for having: 1) at least a blood draw collected, 2) follow-up for at least 60 days to confirm their recovery, and 3) median age and age ranges comparable to the PASC individuals.

When combined with the cohort of 19 individuals described in Figure 1 , PASC and Recovered individuals had the following characteristics: 

Serum samples were inactivated with 1% Triton X-100 for 2h at room temperature according to the Olink COVID-19 inactivation protocol. Inactivated samples were then run on the Olink Explore 1536 platform, which uses paired antibody proximity extension assays (PEA) and a next generation sequencing (NGS) readout to measure the relative expression of 1472 protein analytes per sample. Analytes from the inflammation, oncology, cardiometabolic, and neurology panels were measured.

For plate setup, samples were randomized across plates to achieve a balanced distribution of age and gender. Longitudinal samples from the same participant were run on the same plate. To facilitate comparisons with future batches, sera from 15 donors was commercially purchased (BioIVT) and randomly interspersed amongst the above study samples. Commercial samples included serum from COVID-19 serology-negative, serology-positive, PCR-positive, and recovered (no longer symptomatic) participants.

Data were first normalized to an extension control that was included in each sample well. Plates were then standardized by normalizing to inter-plate controls run in triplicate on each plate. Data were then intensity normalized across all samples. Final normalized relative protein quantities were reported as log2 normalized protein expression (NPX) values.

Olink results and QC flags were reviewed for overall quality. Results for TNF, IL6 and CXCL8, which were measured on all 4 Olink panels, were reviewed prior to averaging to a single NPX value for analysis. Two samples had discrepant cross-panel measurements on these proteins. The results that trended most consistently with the participant's longitudinal measurements were kept and averaged. Serum samples were analyzed in two batches. Following the method recommended by Olink, results of the later batch were bridged to those of the earlier batch using a set of 42 cohort samples that were tested in both batches. A batch offset for each analyte was calculated as the median difference on the 42 samples as measured between the two batches, excluding samples with QC warning flags. The analyte-specific offsets were then added to the raw NPX values of the later batch.

For acute COVID-19, categorization of illness severity was classified by participant report of impact on Activities of Daily Living (ADLs) for each day of illness (U S Department et al., 2017) . Severe days hospitalized were also recorded as were any treatment or therapies received. Participants were scored according to their maximum severity for each day: 0, no symptoms; 1, mild impact on ADLs reported; 2, moderate impact on ADLs reported; 3, severe illness without hospitalization; 4, severe illness with hospitalization; 5, life threatening illness hospitalized with ICU care. Durations were assigned for days spent at each illness severity. A cumulative illness severity score was calculated for each participant by multiplying each severity score by the number of days spent at each level, then summing all values. Participants were classified as post-acute sequelae of SARS-CoV-2 infection (PASC) if any symptoms continued from acute illness or related to COVID-19 beyond 60 days. Sickest PASC participants had ongoing symptoms that significantly impacted activities of daily living, quality of life, and inability to continue employment. 

Following sample acquisition on the Cytek Aurora flow cytometer, samples were unmixed using library reference controls pre-recorded on the instrument using the batch control sample, and the data files were exported in FCS 3.1 format. Unmixing performance was visually checked and unmixing errors were adjusted with compensation on the batch control sample in each individual batch using FlowJo v10.7.1. Compensation adjustments from the batch control were applied to the entire batch. The fluorescent channels in each compensated FCS file were transformed using the logicle transform. The compensated and transformed data was analyzed using bivariate hierarchical gating. Initial placement of gates was implemented automatically using the openCyto R package using customized gating templates for each panel. Following the preliminary automated gate placement, all gates were visually checked for errors and corrected manually per sample in a Flowjo v10.7.1 workspace file generated with the CytoML R package . Gating plots for all samples from each participant were generated with the ggCyto R package (Van et al., 2018) , and gates were further reviewed for consistency and batch effects across longitudinal timepoints. For each population gate, a cutoff of 50 events was also set, below which populations were not reported. Proportions for each cell subtype were calculated by dividing the counts for the gated subpopulation by the total number of CD45+ leukocytes for each sample. Proportions of activated cell populations were calculated by dividing the counts for each positive activation marker population by the counts of the corresponding parent population. Samples were run in multiple batches with bridging controls and we determined that the proportion of variability due to batch was lower than biological variability across samples (Fig. S1b) .

Linear mixed effects models were used to estimate binding antibody RBD titer, are the between-person variation in the intercept and slope of log RBD titer responses respectively, Cov(b, c) is the covariance between the intercept and slope, and

The random effects, ܾ and ܿ , are each assumed to be independent for different individuals and the within-individual errors ߳ are assumed to be independent for different i, j and to be independent of the random effects. A similar model was fit to the untransformed IgG N indices, which are linear on the original scale. The function lme from the R package nlme was used to fit the models. Individual-level estimates of the intercepts, slopes, day 30 and day 180 were obtained from the models.

Data preprocessing: Binary Base Call (BCL) files were demultiplexed using cellranger mkfastq (10x Genomics v3.0.2) to produce FASTQ files for scRNA-seq and HTO barcode libraries. scRNA-seq FASTQ files for each hashed well were aligned to the 10x Genomics GRCh38 reference transcriptome (10x Genomics vGRCh38-3.0.0) using cellranger (10x Genomics v3.1.1) using default settings. HTO barcode FASTQ files were processed using BarCounter (Swanson et al., 2021) to quantify the count of each HTO barcode for each cell barcode. A custom R pipeline, BarcodeTender (to be released at https://github.com/AllenInstitute/BarcodeTender-pipeline), was used to deconstruct gene-by-cell matrices from cellranger outputs and assign cells to their originating participant sample based on HTO counts. For each well, a threshold was determined for each HTO by removing barcodes with low counts, then performing 1 dimensional k-means clustering to group cells into a "high" or "low" group. The minimum count for the "high" group was set as a cutoff, and all cells above that cutoff were considered positive for that HTO barcode. If the mean of the "high" group was not greater than 4-fold of the mean of the "low" group, the minimum cutoff was updated to the mean of the "low" group, and the process was iterated. If a separation could not be established, the cutoff was set to the maximum count, and all cells were considered negative for the HTO barcode. Cutoffs were then used to binarize cell barcodes as positive or negative for each HTO, and the number of positive HTOs was calculated for each cell barcode. Cell barcodes positive for only one HTO barcode were considered "singlets", and any cells with two or more were considered "multiplets". For each well, count matrices generated by cellranger were split to select singlet cell barcodes from each PBMC sample based on HTO barcodes, then cells from each PBMC sample were concatenated into separate matrices for each participant. Throughout the process, metadata tracking the originating Batch and 10x Chromium Chip and Well were retained to enable batch effect analysis. During the splitting and merging process, quality control reports were generated to allow review prior to downstream analysis.

We utilized the supervised PCA projection and anchor-based transfer approach implemented in Seurat v4.0.0 to map cells to a reference projection generated using the Weighted Nearest Neighbor (WNN) graph algorithm, as described in (Hao et al., 2020b) . A CITE-seq dataset from 8 healthy HIV vaccine volunteers during the course of vaccination was used as the reference dataset to label cells. The CITE-seq reference consisted of 228 surface proteins (antibody-derived tags, ADT) and approximately 160k single cells. The reference consisted of three levels of cell label granularity (Levels 1, 2 and 3). We constructed a level 2.5 label structure that consists of all cell type annotations from level 2 except merging the "CD8 TEM_4" and "CD8 TEM_5" from level 2, into a "CD8 TEMRA" label (based on marker gene expression) and replacing the "Tregs" label in level 2, with "Treg naive" and "Treg memory" labels from level 3, thus spanning a complete spectrum of cell type annotations. We assessed the distribution of label transfer scores per cell and retained cells with a label score ≥ 0.5. Batch variability was lower than biological variation contributed by cell type and participant-specific differences (Fig. S1d, S1e) .

The hurdle model implemented in the MAST package (Finak et al., 2015) was used to identify differentially expressed genes between the early acute infection COVID-19 participants compared to uninfected controls. An age matched comparison between infected and uninfected controls was performed, i.e older (> 40 years) COVID-19 participants versus older uninfected controls and younger (< 40 years) COVID-19 participants versus younger uninfected controls. Genes that were expressed in at least 10% of all cells were considered for downstream analysis. Gene expression was normalized by scaling raw counts by the total number of reads per cell multiplied by a scaling factor of 10,000, then using a log 2 transformation. A hurdle model was fit on the filtered and normalized data, modeling the infection status and adjusting for the batch for every cell type. A likelihood ratio test was then performed to assess if the coefficients are different from zero. The p-values were adjusted for multiple comparisons using the Benjamini and Hochberg (BH) method (Benjamini and Hochberg, 1995) . Adjusted p-values < 0.05 were considered significant.

Longitudinal changes in gene expression were also identified using the hurdle model implemented in the MAST package. A hurdle model was fit to each COVID-19 participant independently in order to identify participant-specific longitudinal transcriptomic changes. Genes that were expressed in at least 10% of cells per participant were considered for this analysis. The models were fit on the filtered and normalized data, modeling the days since symptom onset as a continuous variable within each cell type and adjusting for the batch only if any timepoints from the same participant were run across multiple batches. A likelihood ratio test was then performed to assess if the coefficients are different from zero. Obtained p-values are adjusted for multiple comparisons using the BH method. Adjusted p-values < 0.05 were considered significant.

Pathway enrichment analysis: Gene Set Enrichment Analysis (GSEA) (Subramanian et al., 2005) was performed among genes that defined early acute infection status and genes that defined longitudinal changes. A custom collection of genesets that included the Hallmark v7.2 genesets, KEGG v7.2 and Reactome v7.2 from the Molecular Signatures Database (MSigDB, v4.0) was used as the pathway database. The "Type III interferon signaling" gene set was manually curated from the Interferome database (Rusinova et al., 2013) . Genes were pre-ranked by the decreasing order of their log fold changes or coefficients. The running sum statistics and Normalized Enrichment Scores (NES) were calculated for each comparison. The pathway enrichment p-values were adjusted using the BH method and pathways with p-values < 0.05 were considered significantly enriched.

Sample-level enrichment analysis (SLEA, (Gundem and Lopez-Bigas, 2012 ) was used to represent the GSEA pathway expression results on a per-sample basis. The SLEA score was calculated by first calculating the mean expression value of genes (averaged across single cells) enriched in a pathway, then comparing it to the mean expression of random sets of genes (averaged across single cells) of the same size for 1,000 permutations per sample. The difference between the observed and expected mean expression values for each pathway was determined as the SLEA pathway score per sample.

Data preprocessing: scATAC-seq libraries were processed as described previously (Swanson et al., 2021) . In brief, cellranger-atac mkfastq (10x Genomics v1.1.0) was used to demultiplex BCL files to FASTQ. FASTQ files were aligned to the human genome (10x Genomics refdata-cellranger-atac-GRCh38-1.1.0) using cellranger-atac count (10x Genomics v1.1.0) with default settings. Fragment positions were used to quantify reads overlapping a reference peak set (GSE123577_pbmc_peaks.bed.gz from GEO accession GSE123577; (Lareau et al., 2019) ) which was converted from hg19 to hg38 using the liftOver package for R (Lawrence et al., 2009) , ENCODE reference accessible regions (ENCODE file ID ENCFF503GCK; (Vierstra et al., 2020) ), and TSS regions (TSS ±2kb from Ensembl v93; (Yates et al., 2020) ) for each cell barcode using a bedtools (v2.29.1; (Quinlan and Hall, 2010) ) analytical pipeline.

Quality Control: Custom R scripts were used to remove cells with less than 1,000 uniquely aligned fragments, less than 20% of fragments overlapping reference peak regions, less than 20% of fragments overlappingENCODE TSS regions, and less than 50% of peaks overlapping ENCODE reference regions. The ArchR package (Granja et al., 2021) was used to assess doublets in scATAC data. Doublets were identified using the ScoreDoublets function using a filter ratio of 8, and cells with DoubletEnrichment scores from 0-1.16 were considered "singlets" and retained for further analysis. Samples with particularly high doublet scores across all cells (>70% of cells with DoubletEnrichment scores > 1.5) were not considered for downstream analysis.

We used the ArchR package to generate a count matrix for the PBMC reference peak set described above (Lareau, et al. 2019 Nat Biotech) . Dimensionality reduction was performed using the ArchR addIterativeLSI function (parameters varFeatures = 10000, iterations = 2), and the addClusters function was used to identify clusters in LSI dimensions using the Louvain community detection algorithm. For visualization, UMAP was performed using ArchR's addUMAP function with default settings. The ArchR addGeneIntegrationMatrix function (parameters transferParams = list(dims = 1:10, k.weight = 20) was used to label our scATAC cells using the Seurat level 1 cell types from the Seurat v4.0 PBMC reference dataset (Hao et al., 2020b) . We observed that Louvain clusters contained cells with mixed level 1 identity assignments from label transfer, and cluster labels were often spread across the UMAP space. In comparison, cell type assignments in UMAP coordinates seemed to cleanly separate cell-type identities. To generate clusters that more closely matched label transfer results, we performed K-means clustering on the UMAP coordinates using a range of number of cluster centers from 3 to 50, and identified a set of K-means clusters that each had > 80% of cells sharing a single cell type identity. Almost all such clusters contained >= 98% cells from a single major cell type class (T Cells, B Cells, NKs, or Monocytes/DCs/other), with the exception of a single cluster with 88% purity. We used K-means clusters that shared cell class identities to subset the data into T Cells, B Cells, NKs, or monocytes/DC/other classes for downstream analyses. For each broad type, we performed dimensionality reduction by Iterative LSI using 500 bp genomic tiles. We then performed a second round of label transfer using the ArchR addGeneIntegrationMatrix function (parameters as described for level 1, above) using the higher resolution level 2.5 cell identities (described for scRNA-seq label transfer, above) from the Seurat PBMC reference dataset.

Peak calling and motif enrichment analysis: For each of the broad types above, we grouped cells within each level 2.5 cell type label by participant, and used the ArchR addGroupCoverages and addReproduciblePeakSets functions to perform de novo peak calling to identify putative regulatory regions throughout the genome. We ensured that pseudo-bulk replicates derived from participant and cell type intersections were robust by requiring a minimum of 100 cells per pseudo-bulk replicate, with a minimum of 2 replicates for each group.. After identifying de novo peak sets, we annotated transcription factor binding sites (TFBS) in all peaks using the ArchR addArchRAnnotations function (parameter collection = "EncodeTFBS") to label binding sites that have been previously observed in ENCODE datasets. We then used the ArchR addBgdPeaks and addDeviationsMatrix functions (parameter peakAnnotation = "EncodeTFBS") , based on ChromVar (Schep et al., 2017) , to generate a measure of increased or decreased binding site accessibility compared to a random GC-matched peakset. Differential TFBS usage was computed using the ArchR getMarkerFeatures function (parameters bias = c("TSSEnrichment", "log10(nFrags)"), testMethod = "wilcoxon") to identify differentially accessible genomic regions, followed by the peakAnnoEnrichment function (parameters peakAnnotation = "EncodeTFBS", cutOff = "FDR <= 0.1") to score differentially enriched TFBS deviations. Longitudinal TFBS usage scores were calculated by taking the mean TFBS deviation Z-score for each TF, participant, and time point among COVID19-positive samples.

Linear regression models were fit to assess changes in normalized protein expression (NPX) for each protein from Olink, as well as flow cytometry or scRNA-seq cell proportions as a function of infection status using the lm function from the package stats in R. Participant age and biological sex were used as fixed effects in the model in order to control for these potentially confounding variables. P-values were adjusted using the BH method to control false discovery rate (FDR). Adjusted p-values < 0.05 were considered significant.

Linear mixed-effects models (LME) were fit to assess the longitudinal changes in the serum proteome, flow cytometry and scRNA based cell type proportions as a function of days since symptom onset using the lme function (lme4 v1.1-26) implemented in R. Individual NPX protein expression and cell type proportion values were treated as dependent variables. However, COVID-19 participants appear to vary in the slopes and intercepts of days since symptoms. Thus, random effects for both the slope (days since symptom onset) and intercept (participant ID), and fixed effects for age and sex were used in the mixed effects model. P-values were calculated using Wald chi-square tests and adjusted using the BH method to control FDR. Adjusted pvalues < 0.05 were considered significant.

Pathways enriched among significant longitudinally changing proteins were identified using Fisher's overrepresentation analysis. Pathway p-values were adjusted for multiple comparisons using the BH method and pathways significant at adjusted p-value < 0.05 were reported

We performed outlier analysis to identify differential proteins in each COVID19 participant. To identify outliers, we calculated mean expression and standard deviation (SD) of each protein in uninfected participants. The normalized protein expression in COVID-19 participants was compared with the average expression from uninfected participants. Outliers were defined as proteins with expression greater than mean ± 2SD (uninfected participants).

Supervised principal coordinate analysis (PCoA) was performed on Olink data using the R function pcoa from the ape package (v5.5, Paradis and Schliep, 2019, Bioinformatics) . Euclidean distances between proteins that were significantly different between COVID-19 participants at Visit 1 and uninfected controls and proteins that significantly changed over time among COVID-19 participants based on the linear regression models and linear mixed effect models were used to define the distance matrix between samples. The distance of a particular sample to the centroid of uninfected controls was then defined as 

To identify the cell type specific ligand-receptor pairs significantly enriched in COVID-19 participants we selected scRNA-seq data from healthy participants and COVID-19 participant first visit samples for analysis using the nichnetr package (Browaeys et al., 2020) . The ligand-target model retrieved from the nichnetr package includes 688 ligands and 25,345 potential downstream target genes, with values denoting the prior potential that a particular ligand might regulate the expression of a specific target gene. Receivers and senders were identified per cell type for plasmablasts, proliferating CD4+ T cells, proliferating CD8+ T cells, CD14+ monocytes, CD16+ monocytes, NK cells, proliferating NK cells, cDC1, cDC2, and pDC types. Genes were considered for downstream analysis if expressed in at least 10% of the cells from the corresponding cell type. Background gene expression data was obtained from samples of uninfected control participants. Ligandreceptor scores were ranked based on Pearson correlation coefficients (PCC). The top 50 identified ligands were then used to infer putatively active ligand-target links using the nichenetr get_weighted_ligand_target_links function with default parameters. The top 10 ligands and their targets were used to visualize the cell type-specific interaction network. CASP8  IRF1  IRF8  CASP4  STAT3  TRIM21  TRIM5  DDX58  TRIM14  STAT2  IFIT3  IFIH1  IFITM1  IFI27  IFITM3  IFIT5  OAS2  OASL  MX1  MX2  IFI44  IRF7  IFI44L  IFI6  IFIT1  DDX60  OAS3  OAS1  ISG15  ISG20  IFITM2  IFIT2  IRF9  CASP1  IFI35  STAT1  IRF4  STAT4  NFKB1  IRF5  IRF2  IFNAR2  TRIM26 Target   IL10RB  IFNGR1  IFNGR2  FCN1  ITGB1  TGFBR1  TGFBR2  VDR  TNFRSF13B  TNFRSF13C  TNFRSF17  TRAF2  IL6R  TNFRSF1B  CD40  TNFRSF14  IL2RG  IL7R  CSF2RB IFNL1  IFNG  TGFB1  TNFSF13B  ADAM17  CD40LG  LTA  IL7  TNF  TNFSF13  MDK  IL23A  ANGPT2  GMFB  OSM  GZMB   IFNL1  IFNG  TGFB1  TNFSF13B  VEGFB  ADAM17  CD40LG  LTA  IL7 LAMP3  LAP3  SAMD9L  DDX58  BST2  CDC27  EREG  SIGLEC1  TNFRSF8  PTPN1  GRN  CRACR2A  BID  TNFAIP8  ATOX1  CASP2  SNAP23  C4BPB  IRAK1  BAX  TNFRSF14  CCL8  IL18  MAP2K6  CD48  CEACAM1  LILRA2  LAG3  TNFRSF1B  CD274  PSME2  IL2RA  PSME1  TNF  CD300C  CXCL10  LPO  BTN3A2  CD83  FABP5  CCL7  MARCO  MSR1  NUB1  ARSA  PDCD1  NT5E  SH2D1A  CSTB  LILRB4  GZMA  PILRA  LRIG1  CD200R1  LILRB2  LILRB1  LTA  TNFRSF4  CCL19  CLEC7A  SIGLEC6  KLRB1  IL12B  LAIR1  CDCP1  KLRD1  TREM2  KIR3DL1  SIGLEC5  FOLR3   CCN1  CDH6  FAP  THBD  MFAP5  TGFBR3  NRP1  COL6A3  CCL16  COMP  LRP1  CTSV  CXCL5  ITGAV  PDGFRB  SDC4 CCL4 CCN2

Adaptive immune system Innate immune system RET  DDX58  MAP4K5  TIMP1  IFNL1  GLRX  TNNI3  REN  SIGLEC7  LYPD8  TYMP  KIFBP  NUDT5  PPP1R9B  TXLNA  PSPN  FCRLB  FUS  SERPINB6  IVD  ENPP5  DAPP1  ATXN10  PDCD1  ENO1   CGREF1  BAMBI  APOH  BPIFB1  CD209  AGR2  B4GAT1  TFF3  FGF23  NELL2  AOC1  ULBP2  GHRHR  SCG2  PDGFRA  PNLIPRP2  ARTN  CCL28  SLC16A1  SPINK1  CCL15  CTSB  ADAM15  VSTM2L  GDNF  SLC39A14  XG  PI3  SCGN  BGN  SLITRK6  NPY  PLA2G1B  RASA1  GDF2  DCTN6  KIRREL2  TACSTD2  IL34  IL17A   RET  DDX58  MAP4K5  TIMP1  IFNL1  GLRX  TNNI3  REN  SIGLEC7  LYPD8  TYMP  KIFBP  NUDT5  PPP1R9B  TXLNA  PSPN  FCRLB  FUS  SERPINB6  IVD  ENPP5  DAPP1  ATXN10  PDCD1  ENO1   CGREF1  BAMBI  APOH  BPIFB1  CD209  AGR2  B4GAT1  TFF3  FGF23  NELL2  AOC1  ULBP2  GHRHR  SCG2  PDGFRA  PNLIPRP2  ARTN  CCL28  SLC16A1  SPINK1  CCL15  CTSB  ADAM15  VSTM2L  GDNF  SLC39A14  XG  PI3  SCGN  BGN  SLITRK6  NPY  PLA2G1B  RASA1  GDF2  DCTN6  KIRREL2  TACSTD2  IL34 IFNG  LTA  CXCL5  TNFSF14  HMGB2  IL10  PIK3CB  HBEGF  SELL  TNFSF13  HLA E  COL18A1  ANXA1  DSC2  ITGAM  XCL1  HLA DRA  MDK  TGFB1  PTPRC  TNFSF13B  CCL28  EFNA4  CD320  SEMA4D  ADAM17  OSM  VEGFB  ICAM1  ASIP  CD28  GZMB  IL12A  SELP  CXCL16  CCL3  CD40LG  HSPG2  IL18  ITGB1  IL7  IL23A  CCL5  IFNL1  HEBP1  AREG   JUNB  SOCS3  ZFP36L2  TNFAIP3  FKBP5  HIST1H2AG  HIST1H3B  MALAT1  HIST1H2AH  ETFB  LMNA  GAPDH  JUN  DUSP1  FOS  UBC  CXCR4  STAT1  PDE4D  EGR1  ID2  ZFP36  VEGFA  NFKBIA  DDIT4 DDX58  SAMD9L  LAP3  CXCL10  IRAK4  TNFRSF8  STAT5B  LAMP3  IL18  IL6  CXCL11  LGALS9  CCL7  IL1RN  C4BPB  LAT2  IFNL1  IRAK1  TNFSF14  CCL5  LAG3  CCL3  CCL8  CXCL13  CD300E  TNFRSF1B  LILRB4  LCN2  IL16  IL7  CD14  IL18R1  IL2RA  CCL23  CD274  CD300C  LTA  IL1RL1  LYN  CD48  LILRB2  TNF  CD163  TNFSF13B  TNFRSF14  LYAR  CD200R1  TNFRSF1A  CD74  IL18BP  CD38  CXCL12  TNFRSF10A  IL4R  LILRA2  TNFRSF21  C2  IFNLR1  CD83  TNFSF13  TNFRSF13B  CD59  CD34  CD1C  CD5  CCL20  CD97  CCL18  IL12B  LILRA5  CCL19  CCL22  LEP  IFNG  CD4  IL24   0  50  100   Days PSO  Status  Visit  Age  Sex  PSO 

Germline SAMD9L truncation variants trigger global translational repression

Systems biological assessment of immunity to mild versus severe COVID-19 infection in humans

Autoantibodies against type I IFNs in patients with life-threatening COVID-19

Controlling the false discovery rate: A practical and powerful approach to multiple testing

Longitudinal Multi-omics Analyses Identify Responses of Megakaryocytes, Erythroid Cells, and Plasmablasts as Hallmarks of Severe COVID-19

Imbalanced Host Response to SARS-CoV-2 Drives Development of COVID-19

Long COVID in a prospective cohort of homeisolated patients

Inflammatory Type 2 cDCs Acquire Features of cDC1s and Macrophages to Orchestrate Immunity to Respiratory Virus Infection

NicheNet: modeling intercellular communication by linking ligands to target genes

The first 12 months of COVID-19: a timeline of immunological insights

Age-related susceptibility to coronavirus infections: role of impaired and dysregulated host immunity

SARS-CoV-2 Neutralizing Antibody LY-CoV555 in Outpatients with Covid-19

Longitudinal analysis shows durable and broad immune memory after SARS-CoV-2 infection with persisting antibody responses and memory B and T cells

The Immunology of Multisystem Inflammatory Syndrome in Children with COVID-19

Immunological memory to SARS-CoV-2 assessed for up to eight months after infection

Characterizing Long COVID in an International Cohort: 7 Months of Symptoms and Their Impact

An inflammatory cytokine signature predicts COVID-19 severity and survival

SARS-CoV-2 viral load distribution reveals that viral loads increase with age: a retrospective cross-sectional cohort study (medRxiv)

Physical, cognitive and mental health impacts of COVID-19 following hospitalisation--a multi-centre prospective cohort study

Cytokine Storm

Plasma proteomics reveals tissue-specific cell death and mediators of cell-cell

Duration of post-COVID-19 symptoms are associated with sustained SARS-CoV-2 specific immune responses

Mixture models for single-cell assays with applications to vaccine studies

OpenCyto: an open source infrastructure for scalable, robust, reproducible, and automated, end-to-end flow cytometry data analysis

MAST: a flexible statistical framework for assessing transcriptional changes and characterizing heterogeneity in single-cell RNA sequencing data

CytoML for cross-platform cytometry data sharing

Dependence of germinal center B cells on expression of CD21/CD35 for survival

Differential Activation of the Transcription Factor IRF1 Underlies the Distinct Immune Responses Elicited by Type I and Type III Interferons

Untuned antiviral immunity in COVID-19 revealed by temporal type I/III interferon patterns and flu comparison

Antibody response to SARS-CoV-2 is associated with long-term clinical outcome in patients with COVID-19: A longitudinal study

Differential recruitment of dendritic cells and monocytes to respiratory mucosal sites in children with influenza virus or respiratory syncytial virus infection

Endotheliopathy in COVID-19-associated coagulopathy: evidence from a singlecentre, cross-sectional study

Author Correction: ArchR is a scalable software package for integrative single-cell chromatin accessibility analysis

IFN-λ is able to augment TLRmediated activation and subsequent function of primary human B cells

Clinical Characteristics of Coronavirus Disease 2019 in China

Sample-level enrichment analysis unravels shared stress phenotypes among multiple cancer types

Impaired type I interferon activity and inflammatory responses in severe COVID-19 patients

flowCore: a Bioconductor package for high throughput flow cytometry

FKBP5 Regulates RIG-I-Mediated NF-κB Activation and Influenza A Virus Infection

Integrated analysis of multimodal singlecell data

COVID Symptoms, Symptom Clusters, and Predictors for Becoming a Long-Hauler: Looking for Clarity in the Haze of the Pandemic

Loss of Bcl-6-Expressing T Follicular Helper Cells and Germinal Centers in COVID-19

Synergism of TNF-α and IFN-γ Triggers Inflammatory Cell Death, Tissue Damage, and Mortality in SARS-CoV-2 Infection and Cytokine Shock Syndromes

Viridot: An automated virus plaque (immunofocus) counter for the measurement of serological neutralizing responses with application to dengue virus

High-dose influenza vaccine favors acute plasmablast responses rather than long-term cellular responses

Sex differences in immune responses

Comprehensive mapping of immune perturbations associated with severe COVID-19

Droplet-based combinatorial indexing for massive-scale single-cell chromatin accessibility

BAFFR controls early memory B cell responses but is dispensable for germinal center function

rtracklayer: an R package for interfacing with genome browsers

Time-resolved systems immunology reveals a late juncture linked to fatal COVID-19

Sequelae in Adults at 6 Months After COVID-19 Infection

Longitudinal analyses reveal immunological misfiring in severe COVID-19

Coligation of the B Cell Receptor with Complement Receptor Type 2 (CR2/CD21) Using Its Natural Ligand C3dg: Activation without Engagement of an Inhibitory Signaling Pathway

Association of viral load in SARS-CoV-2 patients with age and gender

Regulation of age-associated B cells by IRF5 in systemic autoimmunity

A minimal common outcome measure set for COVID-19 clinical research

Deep immune profiling of COVID-19 patients reveals distinct immunotypes with therapeutic implications

Human Ebola virus infection results in substantial immune activation

Aging impairs both primary and secondary RIG-I signaling for interferon induction in human monocytes

Case Series of Multisystem Inflammatory Syndrome in Adults Associated with SARS-CoV-2 Infection -United Kingdom and United States

Persistent Post-COVID-19 Interstitial Lung Disease. An Observational Study of Corticosteroid Treatment

Age-associated decrease in TLR function in primary human dendritic cells predicts influenza vaccine response

Longitudinal system-based analysis of transcriptional responses to type I interferons

Subcapsular encounter and complementdependent transport of immune complexes by lymph node B cells

Mx1 reveals innate pathways to antiviral resistance and lethal influenza disease

BEDTools: a flexible suite of utilities for comparing genomic features

Pathologically expanded peripheral T helper cell subset drives B cells in rheumatoid arthritis

COVID-19 immune features revealed by a large-scale single-cell transcriptome atlas

Presenting Characteristics, Comorbidities, and Outcomes Among 5700 Patients Hospitalized With COVID-19 in the

Interferome v2.0: an updated database of annotated interferon-regulated genes

chromVAR: inferring transcriptionfactor-associated accessibility from single-cell epigenomic data

Severe COVID-19 Is Marked by a Dysregulated Myeloid Cell Compartment

COVID-19 and the human innate immune system

Adaptive immunity to SARS-CoV-2 and COVID-19

An adult with Kawasaki-like multisystem inflammatory syndrome associated with COVID-19

Long Covid in adults discharged from UK hospitals after Covid-19: A prospective, multicentre cohort study using the ISARIC WHO

Plasmacytoid Dendritic Cells and Type I Interferon Promote Extrafollicular B Cell Responses to Extracellular Self-DNA

mRNA vaccination boosts cross-variant neutralizing antibodies elicited by SARS-CoV-2 infection

The cellular immune response to COVID-19 deciphered by single cell multi-omics across three UK centres (medRxiv)

Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles

Attributes and predictors of long COVID

Simultaneous trimodal single-cell measurement of transcripts, epitopes, and chromatin accessibility using TEA-seq

Interferon-λ Enhances the Differentiation of Naive B Cells into Plasmablasts via the mTORC1 Pathway

Sex differences in immune responses that underlie COVID-19 disease outcomes

Division of AIDS (DAIDS) table for grading the severity of adult and pediatric adverse events, corrected version 2.1 (US Department of Health and Human Services

ggCyto: next generation open-source visualization software for cytometry

Type I and Type III Interferons Restrict SARS-CoV-2 Infection of Human Airway Epithelial Cultures

Development of a Rapid Focus Reduction Neutralization Test Assay for Measuring SARS-CoV-2 Neutralizing Antibodies

Deep Immune Profiling of MIS-C demonstrates marked but transient immune activation compared to adult and pediatric COVID-19

Global reference mapping of human transcription factor footprints

Diverse Functional Autoantibodies in Patients with COVID-19

Retrospective multicenter cohort study shows early interferon therapy is associated with favorable clinical responses in COVID-19 patients

The tidyverse. R Package Ver 1, 1

OpenSAFELY: factors associated with COVID-19 death in 17 million patients

Extrafollicular B cell responses correlate with neutralizing antibodies and morbidity in COVID-19

Rapid and massive virus-specific plasmablast responses during acute dengue virus infection in humans

Characteristics of and important lessons from the coronavirus disease 2019 (COVID-19) outbreak in China: summary of a report of 72 314 cases from the Chinese Center for Disease Control and Prevention

An Infectious cDNA Clone of SARS-CoV-2

RIG-I triggers a signaling-abortive anti-SARS-CoV-2 defense in human lung cells

Inborn errors of type I IFN immunity in patients with life-threatening COVID-19

Extracellular Vesicles with Exosomelike Features Transfer TLRs between Dendritic Cells

We thank the study participants for their dedication to this project; the Allen Institute founder, Paul G. Allen, for his vision, encouragement, and support; Ken Stuart, for constructive feedback on early drafts; the Human Immune System Explorer (HISE) software development team at the Allen Institute for Immunology for their support and dedication. This paper and the research behind it would not have been possible without the collaborative computational data analysis environment provided by HISE.