key: cord-0753242-jj2baaph
authors: Mukund, Kavitha; Nayak, Priya; Ashokkumar, Chethan; Rao, Sohail; Almeda, Jose; Betancourt-Garcia, Monica M.; Sindhi, Rakesh; Subramaniam, Shankar
title: Immune Response in Severe and Non-Severe Coronavirus Disease 2019 (COVID-19) Infection: A Mechanistic Landscape
date: 2021-10-13
journal: Front Immunol
DOI: 10.3389/fimmu.2021.738073
sha: 6635de05a5e74f8620cfa9a476dd43f4d1d9ea22
doc_id: 753242
cord_uid: jj2baaph

The mechanisms underlying the immune remodeling and severity response in coronavirus disease 2019 (COVID-19) are yet to be fully elucidated. Our comprehensive integrative analyses of single-cell RNA sequencing (scRNAseq) data from four published studies, in patients with mild/moderate and severe infections, indicate a robust expansion and mobilization of the innate immune response and highlight mechanisms by which low-density neutrophils and megakaryocytes play a crucial role in the cross talk between lymphoid and myeloid lineages. We also document a marked reduction of several lymphoid cell types, particularly natural killer cells, mucosal-associated invariant T (MAIT) cells, and gamma-delta T (γδT) cells, and a robust expansion and extensive heterogeneity within plasmablasts, especially in severe COVID-19 patients. We confirm the changes in cellular abundances for certain immune cell types within a new patient cohort. While the cellular heterogeneity in COVID-19 extends across cells in both lineages, we consistently observe certain subsets respond more potently to interferon type I (IFN-I) and display increased cellular abundances across the spectrum of severity, as compared with healthy subjects. However, we identify these expanded subsets to have a more muted response to IFN-I within severe disease compared to non-severe disease. Our analyses further highlight an increased aggregation potential of the myeloid subsets, particularly monocytes, in COVID-19. Finally, we provide detailed mechanistic insights into the interaction between lymphoid and myeloid lineages, which contributes to the multisystemic phenotype of COVID-19, distinguishing severe from non-severe responses.

On track to becoming endemic, the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) virus has triggered the global pandemic of coronavirus disease 2019 (COVID- 19) , a complex infection with multisystemic effects (1) . Most respiratory viruses such as the other beta-coronaviruses, respiratory syncytial virus (RSV), influenza, and rhinovirus, trigger a potent innate and adaptive immune response leading to a rapid resolution of infection and generation of circulating memory cells, which can combat a reinfection (2) (3) (4) . Current research however has indicated a more heterogeneous response to SARS-CoV-2 in humans, with symptoms ranging from mild to severe, resulting in hospitalization and mortality (4) (5) (6) (7) in some. Severe COVID-19 involves an extensive cross talk between the activated immune system and other physiological mechanisms, in many cases leading to multisystem comorbidities including acute respiratory distress, septic shock, seizure, renal failure, heart attack, and thromboembolism (8) . Recent research has extensively utilized high-throughput techniques such as single-cell RNA sequencing (scRNAseq) and high-throughput flow cytometry techniques to catalogue the immune cell-state changes contributing to SARS-CoV-2 response [e.g., (5, 7, (9) (10) (11) ]. Broadly, these studies have identified a dramatic remodeling of the major immune players. Specifically, monocytes and neutrophils (and their several novel subtypes) have been described to undergo expansion correlating with disease severity in COVID-19 subjects. COVID-19 research has also focused largely on deciphering the role of adaptive immunity, particularly CD4T, CD8T, and B cells (naïve and mature) in severe, non-severe, and convalescent/recovered COVID-19 [e.g., (10, (12) (13) (14) (15) (16) (17) (18) (19) ]. These studies and many others have identified impaired activation and reduced cytotoxicity arising from the adaptive arm including B and CD8T cells within severe disease. A higher percentage of activated and proliferative CD8T population has been documented within less severe infections (19) . Despite reduced CD4T frequencies, these studies have suggested normal activity for CD4T within COVID-19 (12) . T cells from COVID-19 patients have shown significantly higher levels of exhaustion with increasing severity (14) . Research has shown that though the plasma B and proliferative T-cell repertories correlate with severity, compositional differences of their precursors are influenced heavily by age and sex (16) and exhibit reasonably robust longterm memory against SARS-CoV-2 (18) .

Several of these studies, however, have been limited by cohort sizes, leading to interpretations that while contextually correct may be non-comprehensive. More importantly, these studies do not provide an integrative mechanistic understanding of the remodeling of the immune system and the concomitant alterations in immune response. To overcome this limitation, we analyzed scRNAseq data from four published studies from peripheral blood mononuclear cells (PBMCs) of COVID-19 patients: Lee et al. (4) , Wilk et al. (10) , Schulte-Schrepping et al. (5) , and Arunachalam et al. (20) . These data together comprise a substantial cohort size (111 patients across~350K cells) across 20 immune cell compartments of both lymphoid and myeloid lineage. Despite the larger cohort size, a simple linear analysis such as principal component analysis (PCA) on donors shows a heterogeneity among patient responses with no clean separation between non-severe and severe patients (results not shown). Such lower dimensional analyses demonstrate the need for more complex dimensionality reduction methods for delineating the differences between healthy, non-severe, and severe subjects, as demonstrated by our work. Our integrated analyses not only enhances the statistical power to enable functional and mechanistic insights into COVID-19, but allows for detection of low-frequency or transient cell types that may otherwise be poorly or superficially characterized in smaller cohorts. In contrast to two recently published integrated studies (16, 21) , we focus largely on these less characterized immune cell subsets, which have either low frequencies or are not well explored within COVID-19 literature. Additionally, we identify a set of consensus gene signatures across severe and nonsevere disease enabling the characterization of transcriptional signatures across all immune cell compartments. Taken together, our analyses allow us to put forth a mechanistic framework underlying the interconnected host immune responses, driving and distinguishing severe from non-severe immune responses within COVID-19.

Datasets for each of the published studies referenced above henceforth referred to as Lee (4), Wilk (10) , SS_C1, SS_C2 [two cohorts from SS, (5) ] and PA (20) were downloaded from Gene Expression Omnibus (GEO) or the sources as indicated in the original publications (see Materials and Methods). All donor relevant metadata information (from each study) used within our integrated dataset are provided in Supplementary Tables S1-S3. The original datasets were processed and integrated using Seurat v3.2 as outlined in the Materials and Methods, resulting in an integrated dataset of 375,438 single-cell transcriptomes from 111 donors ( Figure 1A and Supplementary Figures S1A, B) . The resulting parent Seurat object was scaled and clustered into 70 distinct clusters and visualized using UMAP embedding ( Figure 1B) . Grouping the 70 clusters based on their cellular identities revealed 20 distinct immune cell types ( Figure 1C and Supplementary Figure S1C , see Materials and Methods). Monocytes were the largest major cell type identified, particularly CD14 + monocytes at~31% (103,659 cells) , with lowdensity basophils being the smallest population of recognizable cells. A visual inspection of the cellular frequencies across severities revealed some notable differences between the 20 compartments ( Figure 1D and Supplementary Figure S1D) . To discern the trends in cell-type changes and between cell-type interactions, we grouped the mononuclear cell types by origin into myeloid and lymphoid cell types ( Figure 1E , see Section 4.6). Assuming that the cellular abundances correlated with levels of circulating cells, we observed significant gains (p < 0.05) in abundances within COVID-19 (compared with healthy) for monocytes and several low-frequency cell types including, low-density neutrophils (LDNs), megakaryocytes (MKs), and plasmablasts (PBs), consistent with extant research (5, 6, 21) . A progressive loss of several other cell types including dendritic cells (DCs), gdT cells, natural killer (NK) cells, and mucosal-associated invariant T (MAIT) was also observed ( Figure 1E and Supplementary Figure S1D ).

We sought to first understand if a common program of transcriptional dysregulation exists across the myeloid and lymphoid milieu and if it differs across severities. Towards this, we extracted differentially expressed genes (DEGs) by comparing cells from severe and non-severe with healthy cells, independently, for each of the 20 immune cell compartments (see Materials and Methods). A consensus gene set was defined as DEGs that were present in at least 50% of all cell types contained within each lineage for severe and non-severe disease (Supplementary Figure S2A , Supplementary Tables S5 and  S6 ). An inspection of the consensus gene signatures across lineages and severities highlighted the following four broad features ( Figure 1F ). First, significant dysregulation of genes associated with interferon type I (IFN-I) response was observed in both lymphoid and myeloid cells irrespective of COVID-19 severity (Supplementary Figure S2A ). An evaluation of the transcription factor (TF)-target enrichment for the consensus DEGs consistently revealed significant enrichment of TFs STAT1, STAT2, and IRF9 associated with IFN-I signaling (Supplementary Figure S2B) . Several RNA-binding proteins (RBPs) were also enriched within the consensus gene set (Supplementary Figure S2C , see Materials and Methods). Particularly ADAR (IFN-I induced) and YWHAZ were ubiquitously and significantly differentially expressed across severities. ADAR1 is multifunctional and has been extensively studied in the context of innate immunity (22) , while YWHAZ expression in SARS-CoV-2 infection has been hypothesized to contribute to the associated neurological deficits seen in COVID-19 (23) . Second, lymphoid cells showed an upregulation of cytokines and chemokines/receptors such as CXCR4 and interleukin (IL)-32, in both severe and non-severe samples.

Notably, IL-16 was particularly suppressed across multiple lymphoid cell types. Reduced IL-16 levels have been reported in plasma from convalescent COVID-19 patients (24) . Third, myeloid cells in both severe and non-severe disease showed an activation of genes influencing myeloid cell differentiation including HIF1A, TRIB1, HCLS1, and chemokine receptors such as CCR1. Several other chemokines and cytokines such as CXCL16, IL1RN, and IL17RA were significantly dysregulated among classical monocytes (CMs) and myeloid DCs (mDCs) but not commonly across all cell types (see Materials and Methods). Fourth, within severe disease alone, myeloid cells showed an increased activation of genes associated with leukocyte-mediated immune response and degranulation (25) including annexins and their receptors (ANXA2 and FPR1) and complement receptors such as C5AR1. These cells broadly suppressed expression of genes associated with mitochondrial respiration/ oxidative phosphorylation and antigen processing and presentation. Alarmins S100A8 and S100A9, broadly dysregulated across both lineages and severities, showed a consistent upregulation in myeloid cells. Both lymphoid and myeloid cells showed signifi cant dysregulation of immunoglobulin genes including IGHA1/IGHM/IGKC/IGLC2/ IGLC3 involved in complement activation, phagocytosis (recognition and engulfment), and regulation of humoral immune response. Interestingly, TRAFD1 and ETV7, negative regulators of IFN-I response, were enriched in severe disease across both cell lineages (Supplementary Figure S2B ). In the following sections, we assess several immune subtypes that exhibit significant cell abundance differences between healthy and COVID-19 subjects, including mDCs, MAITs, gdT, PBs, and MKs. We also investigate select subsets within monocytes and NK cells, which contribute to the abovementioned gene signatures and have been less explored (functionally and mechanistically) in COVID-19.

Monocytes represent a class of hematopoietic-derived innate immune cells whose function can range from inflammatory to anti-inflammatory. Traditionally classified into discrete subsets classical (CD14 high ), non-classical (CD16 high ), or intermediate monocytes (CD14 + CD16 + ), monocytes are increasingly acknowledged for their cellular, molecular, and functional plasticity (26, 27) . Concurrently, current research in COVID-19 has highlighted the overactivation of CD14 + monocytes and the emergence of novel monocytic subsets including myeloid-derived suppressor cell (MDSC)-like suppressive monocytes within COVID-19 patients (5, 12, 21, 28, 29) . Acknowledging this heterogeneity, we sought to functionally characterize the monocytes focusing on novel subsets and distinguishing them between severe and non-severe COVID-19 (5).

We subsampled and reclustered 22 monocyte clusters (CD14 + , CD16 + , and ITM) from the parent Seurat object, resulting in a total of 19 clusters (Figures 2A, B, Supplementary Figure S3 ). Using automatic annotation (SingleR), expression of known monocyte markers, and the cluster markers (Supplementary  Tables S7, S8) , we reannotated the subsampled monocyte space, identifying 11 distinct monocytic subsets including monocytic MDSCs (mMDSCs; Figure 2C ). Details of transcriptional signatures defining each subset are described in the Materials and Methods (Section 4.8). Eight of 11 (CM1-CM8) subsets were identified to be CD14 + ( Figure 2D ), of which two subsets were particularly interesting and are discussed below.

Subset CM1 (clusters 2, 6, 10, and 14) was characterized by an increased expression of SELL, alarmins (S100A8/9), CD11b/ ITGAM, and CD163 + . Notably, these clusters do not show a loss of HLA-DR expression but express at levels lower than mean (p.adj < 0.05), as compared with other clusters. This observation is consistent with published research, which has observed a distinct upward shift in the S100 high HLA-DR low monocytic population with increasing severity (5, 6, 13) . This subset particularly has been suggested to be influenced by and contribute to both the inflammatory proteomic and altered metabolomic plasma profiles seen in severe disease. An evaluation of the cluster markers within CM1 indicated increased expression of S100A8/A9/A12, CD163, PLAC8, ALOX5AP, CLU, SELL, CREG1, and VCAN, as compared with other clusters ( Figure 2E ). These markers have been previously identified to define a unique monocytic subset (MS1) among bacterial sepsis patients and in COVID-19 and may arise from hematopoietic progenitors via emergency myelopoiesis, with a potentially suppressive function in pathology (28) . TF-target enrichment of all cluster markers within CM1 identified SPI1 (PU.1) as a potent regulator of this subset. CM7 (cluster 9) expressed markers associated with an activated inflammatory state, including CD63, SELL/CD62L, HLA-DR, and ITGAX/CD11c ( Figure 2D) . Notably, this subset expressed MYD88, RELA, and EIF2AK2 ( Figure 2E ). Signaling through the adaptor protein MYD88 (an essential transducer for IL-1B and Toll-like receptor pathways) and the subsequent IFN-I stimulation has been implicated in conditioning of MDSC differentiation during sepsis (30) . The frequency of subsets CM1, CM3, and CM7 was increased in severe compared with non-severe COVID-19. Among CM subsets with increased abundances in non-severe disease, CM4 (clusters 0 and 8) and CM8 (cluster 15) exhibited upregulation of cluster markers associated with IL-1b response including HIF1A, NLRP3, EGR1, ICAM1, CCL3, RIPK2, and ANXA2 (p.adj < 0.05). Constitutively expressed HIF1A plays a role in functional reprogramming of monocytes from proinflammatory to an immune suppressive phenotype, via its regulation of IRAK3 (31) .

Differential gene expression analysis (DGEA) of CM1 and expanded subsets in severe disease (CM7, ITMs, and mMDSC) showed significant enrichment of proinflammatory markers and response to IFN-Is (p.adj < 0.05). CM1 and ITM showed a suppression of genes associated with the mitochondrial respirasome (K) mDCs identified in the integrated dataset represent a mix of conventional DCs (cDC2 and CD1C + DCs) (cluster 35) and pre-DCs (cluster 62). (L) A functional map that highlights the role of these various subsets. Relevant cluster markers are highlighted for each subset. We identified two distinct phenotypes associated with the monocytic subsets-immunosuppressive and proinflammatory. The * by CM1, a suppressive subset, indicates the differential interferon type I (IFN-I) responses between severe and non-severe disease, in particular, a more suppressed IFN response within severe subsets due to likely action of repressive factors including FOSL1 (seen in G above). While the theme of interferon response is shared by the subsets, the more nuanced analysis informs us that CM1 and mMDSC lead to immunosuppressive, specifically T cell-suppressive phenotype. CM1 also expressed genes indicative of emergency myelopoiesis in severe coronavirus disease 2019 (COVID-19) infection. Several CM classes share homotypic aggregation and express genes contributing to platelet-monocyte aggregation, thus leading to the hypercoagulability phenotype in severe infection. IL-1 dynamics within each of these subsets were particularly interesting, with increased expression within nonsevere disease for CM and mDCs indicative of a dysfunctional mDC state within severe COVID-19.

including complex I (such as NDUFA11, NDUFB1, NDUFB11, NDUFB7, NDUFC2, and NDUFV2), complex IV (COX4I1, COX7A2, and COX7C), and mitochondrial metabolism including TGFB1, GUK1, PGAM1, UQCR11, and UQCRB, suggestive of considerable metabolic remodeling associated with these subsets in both severe and non-severe disease (Supplementary Figures  S3C, D Table S8 ). Given the marked expansion of CM1 in severe compared with non-severe disease, we sought to delineate the transcriptional differences by performing DGEA between severe and non-severe cells in CM1 alone. Interestingly, CM1 elicited a more potent suppression of genes responding to IFN-I/II (IFI27, IFIT1, IFIT3, IRF1, ISG15, ISG20, MX1, OAS2, OAS3, STAT1, CCL4, CCL5, and FCGR1A) but upregulated gene targets associated with neutrophil degranulation (RNASE2, ANXA2, CD63, KLF4, KLF2, RETN, C5AR1, SELL, and CLU) in severe compared with non-severe disease ( Figure 2F , Supplementary Figure S3E ) Additionally, several members of the proteasomal degradation pathway including PSMB8/9/10, PSME1/2, UBC, RNF213, and UBE2LG were suppressed within severe disease for this subset. TF-target analysis revealed increased STAT/IRF activity in non-severe patients and an increased activity for FOSL1, FOSL2, and SMAD3 in severe patients ( Figure 2G ). Increased viremia has been previously reported in FOSL1 knockout chimeric mouse. In the presence of increased viremia and IFN-I production, FOSL1 has been shown to serve as negative feedback inhibiting expression of IFN-I (32). On the other hand, FOSL2 and SMAD3 have been previously shown to cooperatively regulate TGFb signaling (33) . TGFb can also suppress IFN-I responses by disrupting mitochondrial bioenergetics in alveolar macrophages during respiratory viral infections (34) . Taken together, these evidences emphasize the suppression of IFN-I responses in CM1 monocytes especially within severe subjects compared with non-severe subjects.

We observed significant enrichment of gene programs associated with homotypic cell adhesion among several subsets of CMs. CD63 (tetraspanin), which is associated with platelet/neutrophil degranulation and intracellular protein trafficking within monocytes, and L-selectin or SELL (CD62L), which is a major regulator of monocytic adhesion migration and signaling, were significantly upregulated in COVID-19, irrespective of severity ( Figure 2H ). Patients with certain autoimmune conditions (35) and viral infections (36) are known to exhibit increased circulating monocyte aggregates (in both the presence and absence of platelets), which are characterized by an overexpression of SELL and CD63, among other markers. We identified several of the other previously documented monocyte aggregate markers within our analysis including VCAN, and integrins such as ITGA4, ITGAM, and ITGB2 ( Figure 2I ). Additionally, the interaction between platelet specific isoforms of selectin (P-selectin) and CD63 has been recently implicated in the increased formation of platelet-monocyte aggregates contributing to the hypercoagulability phenotype seen in severe COVID-19 (37).

Two clusters within the primary Seurat object were identified as mDCs. Using previously published markers of human DCs (38), we identified the mDC subset to represent a mix of conventional DCs (cDC2, CD1C + DCs) (cluster 35 with transcriptional activation of ID2, ZEB2, IRF4, KLF4, CD1C, ITGAX, CLEC4A, CLEC10A, SIRPA, and FCER1A) and pre-DCs (cluster 62, AXL/ CLEC4C(CD303)/SIGLEC1), which are typically poised towards cDC2 ( Figures 2J, K) . Relative cDC2 cell abundances reduced with increasing severity (2,577 cells from healthy subjects, 1,426 cells from non-severe subjects, and 734 from severe subjects), consistent with recent observations in severe COVID-19 (39) ( Figure 1D ). DGEA with respect to healthy subjects indicated a robust activation of interferon-stimulated gene (ISG) signatures and response to IFN-I, as observed in monocytes, irrespective of severity (Supplementary Figure S3F) . Notably, we observed a suppression of key cDC2 mDC markers involved in DC signaling and response within severe disease compared with non-severe disease (Supplementary Table S9 ) including CD83 [marker for activation of mature mDCs (40) ], NR4A2 [expressed in immunogenic DCs and promotes anti-inflammatory cytokines (41) ], and FCER1A [inflammatory mediator likely promotes immune homeostasis (42) ].The mDC also exhibited increased expression of IL-1B in non-severe compared with severe disease (p.adj < 0.05, Supplementary Table S9 ). This suppression of IL-1B expression in mDCs has been suggested to lead to increased production by peripheral monocytes (43) ( Figure 2E ). IL-1b production in both monocytes and mDCs involves the canonical NLRP3 inflammasome-induced IL-1b cleavage and release, dependent on purinergic receptors such as P2RX4 (44) . These subsets exclusively expressed purinergic receptors such as P2RX4 (p.adj > 0.05). The impaired IL-1b dynamics arising in severe compared with non-severe disease in addition to the suppression of key markers involved in DC signaling are indicative of a dysfunctional role for mDCs within severe patients. An integrated perspective on the CM and mDC subsets and their phenotypic contributions is provided in Figure 2L .

LDN is an umbrella term often used to represent a heterogeneous group of highly adaptable, dynamic neutrophil-like cells composed of a mixture of immature and low-density mature neutrophils, progenitor cells, and granulocytic/polymorphonuclear MDSCs (PMN-MDSC). Given this cellular heterogeneity, their function is suggested to exist on a spectrum, ranging from immunosuppressive to proinflammatory (45, 46) . Subsampling and reclustering "LDNs" from the parent Seurat object resulted in nine clusters, grouped subsets into three distinct subsets based on the expression of a combination of established markers, namely, low-density granulocytes (LDGs), PMN-MDSCs, and progenitor-like ( Figures 3A-D) . Details of the transcriptional characterization of these subsets have been provided in Section 4.8). The progenitor (like) subset exhibited signatures reminiscent of a "proneutrophil" state (5) ( Figure 3E ). We identified this subset to also express AZU1, RNASE3, CSTA, CTSG, and RETN (p.adj < 0.05, Figure 3E ), which are crucial mediators of neutrophil communication via the cytoneme and contribute extensively to NETosis. Interestingly, previous studies have identified several additional markers expressed in this subset including MS4A3, PCNA, HMGN2, SRGN, PTMA, STMN1, DUT, and TUBA1B to be indicative of precursor neutrophils, which differentiate into neutrophils via alternate maturation mechanisms (47, 48) . This subset also uniquely expressed MKI67 (p.adj < 0.05) suggestive of active proliferation ( Figure 3F ). Taken together, these findings highlight a population of progenitor-like LDNs, which are highly proliferative and activated predominantly in severe COVID-19.

Characterizing the transcriptional control (Supplementary Figure S4A ) revealed increased activity of TFs such as LEF1, MYC, and MYCN within this subset. These TFs are crucial determinants for neutrophil granulopoiesis, proliferation, and differentiation (49), further emphasizing the progenitor-like state for cells contained within this subset.

We observed LDGs to be significantly proinflammatory (expressing chemokines) and IFN-I responsive (Supplementary Figures S4A-D) (50) . Transcriptional control of LDGs was likewise heavily influenced by TFs such as STAT1/3 and NFKB1/RELA. The top ranked cluster markers identified within LDGs have been previously implicated to play significant roles in neutrophil maturation, activation, and degranulation including NAMPT, IFITM2, IFNAR1, SOD2, G0S2, AQP9, and CXCR2 ( Figure 3E and Supplementary Table S10 ). Mature neutrophils in certain conditions, such as inflammation and cancer, can present at lower densities, at steady state without additional activation. Notably, activation of AQP9 is associated with changing densities of mature neutrophils due to water uptake, irrespective of activation/degranulation status and their subsequent shift to lower densities (45, 51) . A robust IFN response has been previously reported to originate from LDGs in vivo and in vitro, with distinct roles in the pathogenesis of diseases such as lupus and malaria (52, 53) . DGEA (with reference to healthy cells) indicated that LDGs exclusively had a statistically significant upregulation of genes associated with a potent IFN-I response and activation of ISGs in both severe and non-severe subjects (Supplementary Table  S11 ) . Dysregulated genes and TFs included BST2, IFI35, IFI6,  IFIT1, IFIT2, IFIT3, IFITM1, IFITM2, IFITM3, ISG15, ISG20,  MX1, OAS1, OAS3, OASL, RSAD2, XAF1, IRF1, and IRF7 . However, this response to interferon signaling seemed to be muted in samples from severe patients compared with nonsevere samples. In particular, genes associated with processes involved in FCg receptor signaling and IFN-I response including CYFIP2, HSP90AA1 HSP90AB1, WIPF1, WAS, ELMO1, CD3G, FYN, PIK3R1, CD247, VAV1, IFI6, OAS2, MX1, and IFIT1 were downregulated in severe disease (p.adj < 0.05). Interestingly, certain proinflammatory ISGs including IFITM1/2/3 and ISG20 were upregulated in severe disease. Alarmins (S100A8/S100A9/ S100A12) were upregulated in COVID-19, but significantly higher within severe disease (p.adj < 0.05). Additionally, LDGs within severe samples (compared with both healthy and nonsevere subjects) displayed a transcriptional repression of genes associated with neutrophil migration such as CXCR2, CXCR1 (Supplementary Figure S4E) , CD74, ITGB2, and RAC2 (p.adj < 0.05, Supplementary Table S12). CXCR2/CXCR1 serve as cognate neutrophil receptors driving neutrophil migration and are markers of mature neutrophils. In severe sepsis, nitric oxidemediated suppression of CXCR2 is associated with an impaired migration of neutrophils to the infection loci (54).

PMN-MDSCs comprise a functionally distinct phenotype of MDSCs with a neutrophil-like morphology that are known to play important roles in the immune dysregulation of several inflammatory states such as sepsis, cancer, and, as shown here, in severe COVID-19. We confirmed drastic expansion of PMN-MDSCs in an independent cohort of severe and hospitalized COVID-19 patients ( Figures 3D, G) . Given this expansion, we compared cells from severe and healthy subjects to establish the transcriptional landscape. Enrichment of DEGs (Supplementary Figure S4F ) revealed an expected suppression of genes, which regulate T-cell activation including CD47 (a neutrophil membrane protein), ARG1, CEACAM1, LILRB2, FYN, LYN, ITK, HLA-DPA1/B1, RAC2, PYCARD, RUNX3, ANXA1, ITGB2, and CCL5. The T cell-suppressive activity of PMN-MDSCs is suggested to be driven by a potent induction of reactive oxygen species (ROS), by the nicotinamide adenine dinucleotide phosphate (NADPH) oxidase system and neutrophil degranulation (55) . We consistently observed an increased expression of the NADPH generating superoxide genes such as NCF1 and CYBA among highly ranked DEGs (Supplementary Table S13 ) and increased regulation of electron transport chain genes COX5B, NDUFB7, UQCR11, COX6A1, and COX4I1, which participate in oxidative phosphorylation (OXPHOS). PMN-MDSCs in severe samples showed an upregulation of PGAM1 and GAPDH, suggesting increased oxygen consumption due to increased mitochondrial/metabolic energy metabolism via OXPHOS, and increased production of ROS (hydrogen peroxide (H 2 O 2 ), superoxide anions, and hydroxyl radicals through the NADPH system). Enrichment analysis also indicated a significant enrichment of Snitrosylation within PMN-MDSCs (Supplementary Figure  S4F ). S-nitrosylation has been suggested to have inflammatory consequences through a complex interplay of mechanisms especially within pathologies such as sepsis and cancer (56) . Interestingly, function of GAPDH, a well-studied target of S-nitrosylation, is mediated by alarmins S100A8/S100A9. Figure 3H captures the interplay of mechanisms highlighted in the sections above for proliferating progenitor-like, suppressive PMN-MDSCs and the proinflammatory LDGs, particularly within severe COVID-19.

MKs are the mature cells from which platelets are derived. Recent evidence suggests a direct role for MKs in viral infections (in part via IFITM3 upregulation) and in systemic inflammation, highlighting the importance of MKs and their interactions with other immune cells (57) . Two clusters from the parent Seurat object, 30 and 54, were annotated as MKs (see Materials and Methods) using previously published MK markers ( Figure 3I ). Extravasation to and expansion of MKs among PBMCs (6) and lungs (58) in severe COVID-19 has been recently reported. We concurrently observed a dramatic expansion of MKs within severe samples in the integrated dataset (1, 195 healthy cells, 1,153 non-severe cells, and 4,848 severe cells). DGEA in severe and non-severe disease compared with healthy subjects identified 434 and 304 DEGs in severe and non-severe disease, respectively (Supplementary Table S14) . Notably, both severe and non-severe samples upregulated ISGs including ISG15, IFITM3, and IFI6 and certain immunoglobulins including IGKC, IGLC2, and IGHA1. In severe disease, however, uniquely upregulated genes associated with increased mitochondrial energetics (genes such as COX4I1/5B/6A1/7A2/8A, ISCU, NDUFA1/A2/B1/B11/B3/UQCR11, SOD1/2, and CFL1), leukocyte degranulation (such as S100A8/A9, cathepsins including CTSD/W and SELP), and coagulation (such as FLNA, ITGA2B, ANXA5, MMRN1, and ANO6), alluding to their role in contributing to the thromboembolic phenotype of severe COVID-19 (Supplementary Figure S4G , Supplementary Table S15). Severe samples were also enriched for genes associated with an increased response to IL-12, a proinflammatory cytokine. JAK-STAT activation and subsequent IFN-g generation, after IL-12 stimulation, have been previously reported in patients with immune thrombocytopenia (59) .

Thrombocytopenia has been reported as a prominent feature of severe COVID-19 (60) . Immune thrombocytopenia is associated with increased platelet demand and additionally exhibit emperipolesis (61) . We hypothesized that thrombocytopenia seen in COVID-19 could be associated with the observed MK expansion and MK-mediated emperipolesis. Emperipolesis is a unique phenomenon of cell-cell interaction, which involves a bidirectional membrane transfer between cells such as neutrophils and MKs, where neutrophils enter the MK, fuse their membranes with the MK's demarcation membrane system (DMS), and then exit the MKs intact. Two receptor ligand pairs have been suggested to mediate emperipolesis including ICAM1/EZR and CD62P/PSGL1 (61) . DGEA identified significant upregulation of SELP(CD62P) in severe (compared to healthy) (Supplementary Figure S4G) . In pathological conditions such as idiopathic myelofibrosis, characterized by emperipolesis, CD62P (a granule protein) is increased and distributed abnormally to the DMS. The abnormal distribution of CD62P to the DMS is suggested to likely trap neutrophils on the DMS via its binding to PSGL-1 expressed on neutrophils. Notably, PSGL-1 expression was significantly upregulated in neutrophil subsets, especially LDGs, in severe COVID-19 samples within our study ( Figure 3F) . Emperipolesis is also a prominent feature of MKs in GATA1 low murine models, which also exhibit severe thrombocytopenia (62) . Consistent with these observations, samples from severe disease exhibited lower (than mean) levels of GATA1 expression (Supplementary Figure S4H) . Though much remains to be understood on whether emperipolesis is the cause or consequence of thrombocytopenia, current evidence suggests a potential role in the pathogenesis of COVID-19 ( Figure 3J ).

NK cells are a type of lymphocytes that respond rapidly to eliminate/control a host of pathological insults including tumor progression, microbes, and viruses. Recent research in COVID-19 has indicated a robust activation of NK cells, albeit with impaired cytolytic potential and reduced absolute cell counts, with increasing severity (9, 63) . Moreover, NK cells also show signs of exhaustion and increased interferon signaling with increased expression of the inhibitory surface proteins in moderate and severe patients (13) . The changing NK landscape in COVID-19 has been largely discussed in the context of CD56 dim /CD56 bright (NCAM1) NK cells. However, NK cell populations are increasingly acknowledged to be transcriptionally heterogeneous (64, 65) . We subset the NK cells to better characterize their heterogeneity and role within severe and nonsevere COVID-19.

The primary UMAP of the integrated dataset indicated a reduction in the NK cell population in severe compared with b o t h n o n -s e v e r e CO V I D -1 9 a n d h e a l t h y s u b je c t s (Supplementary Figure S1D ). Reclustering and re-embedding the subsampled NK cells after processing resulted in 18 clusters grouped into seven distinct subsets ( Figures 4A-C) . Subset annotation has been described in detail within Section 4.8 (Supplementary Figures S5A, B) . Majority of the subsets identified here, including NK3, NK4, NK5, NK7, and NK8, were CD56 low CD16 + and accounted for nearly 82% of the NK cells identified ( Figure 4D ). CD56 low NK cells represent the major circulatory subset of human PBMCs ( Figure 4E ). On the other hand, two subsets NK6 and NK2 were characterized by increased expression of CD56. Both the CD56 + subsets expressed PRF1, PRDM1, and S100A4 (a maturation marker) lower than mean ( Figure 4C ). CD56 + cells are less cytotoxic but actively produce cytokines such as IFNg in response to cytokinemediated stimulation with IL-2 and IL-18. Reduction of CD56 + cells, similar to what is observed here ( Figure 4E) , has been implicated in impaired IFN-g production and reduced cytotoxic activity in certain disease states (66) . Particularly interesting however was that the transcriptional profiles of NK2 differed significantly from those of the typical CD56 bright (NK6) population. This raised the question if this subset represented an intermediate stage in the maturation of NK cells or a novel subset of CD56 + cells with alternate roles. The top cluster markers for NK2 included genes such as NCAM1(CD56), IL18RAP, SPN (CD43), MACF1, AHNAK BTN3A1, SPOCK2, SORL1, PARP8, and ETS1 and functionally enriched for a potent response to IL2 (p.adj < 0.05).

Subset NK1 (cluster 13) on the other hand exhibited a low expression of CD16, along with a lack of CD56/CD44 expression (Supplementary Figure S5C) . Based on the increased mean expression of KLRF1 (Nkp80) and lack of CD56 expression, NK1 was determined to be a CD56 neg CD16 low population of NK cells (67) . There was a reduced representation of this subset within severe ( , with cluster markers represented within each subset. NK1 (CD56 − CD16 low ) represented a unique subset of cells, mostly seen in healthy and lost within COVID-19. These cells exhibit reduced NCR expression, subsequently implying a reduced cytotoxic potential and ability to communicate T and neutrophils. Within the CD16 + subsets, we identified an expanded subset that responded potently to interferon type I (IFN-I). The cluster markers of this subset were also significantly differentially expressed in severe and non-severe compared with healthy subjects. Red text indicates increased activity, while green text indicates reduced activity.

bacterial infections has indicated an expansion of the "dysfunctional and immature" CD56 neg /CD16 low NK population (68) . Cluster markers for this set included significant upregulation of alarmins (S100A8/A9), cytoplasmic ribosomal proteins, and OXPHOS genes (Supplementary Figure  S5E , Supplementary Table S16) . NK1 also showed reduced expression of cytotoxicity receptors of active and resting NK cells such as NCR3 (NKp30, which is also involved in interactions between NK cells and neutrophils) and CD247 (CD3z, a key subunit of natural cytotoxicity receptors suggestive of reduced cytotoxic potential) (Supplementary Figure S5F ). Regulon analysis for NK1 cells indicated increased activity of targets of repressive TFs such as FOXP1 [which inhibits JUN/MYC signaling and GZMB and IFNG production in T cells (69)], NR3C1 (glucocorticoid receptor involved in inhibiting IFNG production), and early NK cell TFs such as SOX9, suggestive of a precursor-like, immature state of cells ( Figure 4F ). Whether this subset represents a precursor-like immature population, which is subsequently lost due to maturation within COVID-19, warrants further research.

The top cluster markers and differential expression analysis comparing NK4 cells from severe and non-severe patients to healthy subjects indicated a strong activation of genes associated with IFN-I response including ISG15, IFI6, MX1, IFIT3, XAF1, RSAD2, OAS3, and STAT1 (enrichment shown in Supplementary Figure S5F and Supplementary Table S17 ). TF activity analysis indicated a robust regulation by IRF2/IRF9 and STAT1/2 for this cluster of cells, within both severe and nonsevere COVID-19 ( Figure 4G ). The NK4 subset also showed higher than mean expression of markers including CX3CR1 and CD38, previously attributed to terminally differentiated NK cells (64) . NK4 was the only subset to expand within COVID-19 (more prominently in non-severe) ( Figure 4B ). NK4 cells in COVID-19 exhibited cytolytic potential with increased expression of granzymes (GZMB and GZMA) and expressed CD56 (Supplementary Figure S5G) . In addition to a potent response to IFN-I, DGEA indicated a potent upregulation of alarmins such as S100A8 and TFs such as PLSCR1 in both severe and non-severe subjects, similar to what is seen within the interferon responding clusters of monocytes and LDNs. Notably, in the context of MCMV viral infection, a novel IFN-I dependent mechanism has been identified by which NK cells evade mechanisms of cell death via BCL2 signaling (70) . The observed expansion of this NK4 subset, in COVID-19, could occur via similar pro-survival mechanisms. We outline the abovementioned observations for all the NK cell subsets in Figure 4H .

Adaptive immunity is crucial for successful viral clearance and long-term immune memory, particularly T cells. As noted in the parent Seurat object, CD4T and CD8T subsets did not show very drastic changes between severe, non-severe, and healthy subjects ( Figures 5A, B, Supplementary Figures S1D and S6C) . To further delineate CD4/CD8T subsets, which might show differential signatures, we subset the CD4T/CD8T populations and reclustered them into eight distinct subsets as described within the Materials and Methods (Section 4.8). Abundance signatures of several CD4/CD8T compartments identified here indicated a downward shift, albeit not drastic, with increasing severity ( Figure 5B) , which is consistent with much of the published research (5, 7) . Specifically, we observed a higher proportion of CD8T/CD4T effector (Teff) cells in COVID-19 compared with healthy subjects, but lower in severe compared with non-severe subjects. A higher percentage of activated and proliferative CD8T effector population has been documented within less severe infections (12) . Likewise, a large proportion (~36.6%) of the cells identified as naïve CD4T cells were significantly increased in non-severe disease, compared with severe disease, consistent with published research (71) . In contrast to the abundance trends of CD4T/8T naïve and effector cells, we observed a progressive increase in the abundances of CD4 + Tregs, a low-frequency cell type, across severity ( Figure 5B ). Increase in the CD4 + Tregs has been previously reported through flow cytometry experiments in the blood of COVID-19 patients (71) . Tregs are crucial for regulating immune homeostasis and autoimmunity, controlling the quality and magnitude of immune responses in infections by modulating expression of key factors including FOXP3, CTLA4, and IL2RA (Supplementary Figure S6B ) (72) . Particularly interesting within the clusters identified were two clusters 16 and 18, which, in addition to CD8 Teff-like signatures (including expression of granzymes, IFNG, PRF1, and NKG1), showed expression of CD16. The CD16 + CD8T subset showed significant increases in proportion of cells within both severe and non-severe COVID-19 compared with healthy cells (Supplementary Figure S6C) . CD16 + CD8T cells have been reported previously in certain viral infections such as hepatitis (73) and in blood from smokers (74) and are suggested to emerge as a consequence of differentiation of T cells beyond terminally differentiated effector states, acquiring CD16 and NK-cell like functional properties. To further elucidate the functional relevance of this subset in COVID-19, we examined how the transcriptomes of this subset differed compared with those of healthy subjects. DGEA analysis identified 140 DEGs in severe and 201 DEGs in non-severe disease compared with healthy subjects (Supplementary Table S18 ). MSigDB analysis revealed both sets of DEGs to be significantly enriched for NK/T cell repertoires, consistent with the earlier characterization of NK-Tlike cells (Supplementary Figure S6D) . Functional analyses revealed that both severe and non-severe disease subsets exhibit significant response to IFN-I signaling (75) and an overexpression of alarmins including (S100A8/A9) (Supplementary Figure S6E) . Notably, however, only within severe disease do we see an upregulation of functions typically associated with platelet degranulation, homotypic adhesion, and coagulation including markers such as PPBP, ITGA2B/ITGB3, PF4, and TREML1. While most of these factors are typically associated with myeloid cells, emerging evidence suggests a role for these markers in modulating CD8T and NK differentiation and function in health and disease and is worth exploring in the context of this subsets expressing CD16 and NK-like CD8 + T cells.

A notable decrease in cellular abundances of low-frequency T subsets including MAIT (cluster 15) and gdT (cluster 29) was seen in COVID-19 samples compared with healthy subjects (Figures 1B, 5C , and Supplementary Figure S1D ). This change in cellular frequencies was further validated in an independent cohort ( Figure 5D ). Each cluster expressed markers consistent with cell types, MAITs expressed (KLRB1 (CD161), SLC4A10, NCR3, DPP4 (CD26), IL7R, and GZMK) (76) ; and gdT cells expressed transcriptional markers such as CD8A, CD8B, CD2, CD5, CD7, TRDC, and TRGV9. Interestingly, CD26 (DPP4), a suggested target for the SARS-CoV-2 spike proteins, was also found to be expressed only within subsets of the lymphoid compartment, especially MAIT (77) (Supplementary Figure S7A) .

MAITs are a class of non-conventional T cells, representing 1%-10% of the circulating T-cell population and preferentially respond to innate inflammatory signals including IL-12, IFN-g, and IL-18 with viral infections including COVID-19 (76, 78) . Likewise, gdT cells, which are a class of restricted T cells, are also activated preferentially by IL-12. Consistently, DGEA for MAIT and gdT cells identified an increased response to IL-12, particularly in severe disease (compared with healthy subjects), with an additional enrichment of genes associated with IL-2 production (more prominent in non-severe disease) (Supplementary Figures S7B, C) . Despite reduced abundances, MAIT cells in severe and non-severe disease showed increased activation [as detected through expression of CD69, an early activation marker (p.adj < 0.05), consistent with earlier reports (78)] (Supplementary Figure S7D) . Under inflammatory conditions, gdT recruits and activates neutrophils through the release of cytokines and chemokines. In the presence of H 2 O 2 , neutrophils can suppress gdT action (65) . Neutrophils are also thought to suppress MAIT cells in the presence of H 2 O 2 (66) . Increased expression of CD3E and a suppression of CD3z (CD247) was observed on MAIT and gdT cells in COVID-19 samples ( Figure 5E ). These expression patterns have been previously observed in MAIT and gdT in the presence of neutrophils. Suppression of CD3z chain further indicated compromised T-cell signaling without change to T-cell viability. These observations suggest that the reduced expansion of these cell types in COVID-19 is likely due to excessive activation of neutrophils/LDNs and increased oxidative stress. Whether loss of circulating MAIT and gdT to the airways and other tissues contributes to the observed reduction and subsequent turnover merits a more thorough investigation.

In the context of acute viral infections such as in COVID-19, immature B cells mature to naïve B cells and differentiate to memory cells or antibody-secreting cells (ASCs) upon antigen activation. Extant omics analyses of immune remodeling in COVID-19 have shown overall lymphoid dysfunction, with decreases in multiple cell types, including naïve B cells; however, PBs were a notable exception and showed expansion and heterogeneity with increasing disease severity (6, 10, 21) . A flow cytometric study of multiple B-cell subtypes showed overall decreases in most B-cell subtypes, with increasing disease severity, except for ASCs (17) . Likewise, other previous flow cytometric analyses report expansion of oligoclonal PBs in severe COVID-19 (11, 79) . To better understand how intracellular mechanisms, including transcriptional regulation, may contribute to B-cell dysfunction with increasing disease severity, we subset B cells as identified in the parent Seurat object (Figure 1 ) for further analysis (Supplementary Figure  S8 ). Of these, we analyzed in greater detail the populations that had the most marked changes between healthy, non-severe, and severe COVID-19 subjects: naïve B cells, which reduced with increasing severity (Supplementary Figure S8) , and PBs, which expanded in severe COVID-19 ( Figure 1D ).

From the initial UMAP ( Figure 6A and Supplementary Figure  S8) , we identified six groups of naïve B cells (A-F) by unsupervised hierarchical clustering (see Materials and Methods) ( Figure 6B) . Groups A-D comprised predominantly of healthy and non-severe cells, while in groups E and F, severe cells outnumbered healthy and non-severe cells. For clarity, groups A-D will hereafter be referred to as "B NS " (mostly nonsevere), and groups E and F will be "B S " (mostly severe). All the groups had expression of IGHD and IGHM, consistent with naïve B cells. Group F has additional expression of IGHA1/2, perhaps indicating priming for early stages of activation ( Figures 6C) . Expression of multiple markers was consistent with naïve B cells, and the top 15 significantly expressed cluster markers for each group of A-F are shown in Supplementary Figure S8 . Figure 6D summarizes the remodeling of naïve B cells in COVID-19. Signaling via the canonical B-cell receptor (BCR) is crucial for B-cell survival, development, and antibody production. Genes involved in BCR signaling had variable expression in non-severe and severe groups; for instance, CD22 and CD79A showed increased expression in B NS groups, while CD40 and CD79B were increased in B S . Genes downstream of BCR signaling involved in regulation of intracellular calcium showed increased expression within B NS . Several other cell surface markers, including MHC-II components (CD74 and CD83), TNFRSF13C, and IL2RG, showed increased expression in B NS relative to B S , while both groups expressed NF-kB inhibitor genes (NFKBIA/D), downstream of TNFRSF13C. Taken together, these indicate a severity-dependent difference in naïve B-cell activation, downstream intracellular signaling, and antigen presentation.

Although key TFs for early B-cell maturation (such as PAX5, BACH2, and FOXP1) were expressed in all groups, they were statistically significantly increased expression in B NS groups. In addition to subtle severity-dependent variations in maturity and activation, this could reflect multiple origins for naïve B cells in secondary lymphoid organs. For example, ZPF36L1, a TF required for the maintenance of marginal B cells, is increased B NS (80) . By contrast, multiple AP-1 family TFs including JUN, FOS, and FOSB had increased expression in the B S . Notably, group B within B NS showed a robust response to interferon stimulation, with increased expression of TFs such as IRF1, IRF7, and STAT1, as well as downstream genes including MX1, OAS1, and IFI6.

We observed a metabolic shift from oxidative phosphorylation to glycolysis with increasing disease severity, which likely reflects the primed state (for activation) of these cells. We additionally observed that B NS had expression of genes protective against oxidative stress compared with B S including HERPUD1, HIF1A, and GADD45B (may be protective against genotoxic stress in lymphocytes) (81), while B S had an increased expression of ubiquitin-proteasome genes (e.g., PSME1 and UBC). Differences in stress response could play a role in preparing B cells in COVID-19 for enhanced proliferation. Pro-proliferation genes such as MYC and BST2 were increased in Bs groups, likely suggestive of a state primed for proliferation in severe disease. Additionally, expression of cells markers related to adhesion and migration including a reduced expression for CXCR4 and SELL in B S with an increased expression of CD69 and IER2 could indicate an increased potential of B cells to differentiate into PBs (82, 83) . Taken together, these could indicate a state of naïve B cells primed for differentiation, perhaps more readily to the PB fate, within severe COVID-19.

PBs are short-lived, antibody-producing cells that are derived from antigen-activated memory B cells. We observed considerable expansion and increased heterogeneity of the PB cell population with increasing disease severity ( Figure 1D ). Previous flow cytometric analyses have reported expansion of oligoclonal PBs, notable even in the context of overall lymphoid cell dysfunction, in severe COVID-19 (11, 79) . While extant omics analyses of immune dysfunction in COVID-19 have shown PB expansion and heterogeneity (6, 10, 21) , the intracellular mechanisms that may contribute to cell dysfunction have not been explored in detail. To this extent, we subsampled and reclustered the PBs from the parent Seurat object resulting in 10 clusters, which were subsequently grouped into four major subsets (described in Section 4.8) (Figures 6E-G and Supplementary Figure S8 ). Subset-dependent expression of cell surface markers, including BCR components, cytokine/chemokine receptors, and HLA genes, reflected putative heterogeneity in cell growth, homing potential, and maturity ( Figure 6H ). Genes involved in promoting PB survival and growth TNFRSF17, BST2, and CD79B had increased expression in PB1, the dominant signature in non-severe patients. Chemokine and cytokine receptors, likewise, had varied expression across subsets including CCR10, ICAM2, and IFN-g in PB1; CCR2, IFN-a, and IL-6 receptor components in PB2; and CXCR4 and IL2RG in PB4. Increased expression of CCR10 in PB1 (more prevalent in non-severe disease) may be indicative of cells that home to mucosal surfaces (84) , while CCR2-and CXCR4-positive cells in PB2/4 (more prevalent in severe disease) may be indicative of cells that home primarily to the bone marrow (85) . Immunoglobulin chain expression was also subset-dependent, with IGHA1 more highly expressed in PB1, while PB2 had robust expression of the IGHG chains ( Figure 6I) . Finally, surface marker expression, together with TF expression (Figure 6J) , suggested varying maturation across PB subsets, with increased expression of IRF4, PRDM1, and XBP1 in PB2/4 indicating increased maturity (86) .

We next investigated the subset-dependent expression of key TFs and several of their targets that regulate PB commitment. All subsets expressed the core TFs involved in PB maturation-IRF4, PRDM1, XBP1, and POU2AF1-albeit with differential expression between the subsets: POU2AF1 was significantly expressed in PB4, while IRF4 and XBP1 expression was significantly increased in PB2. Concomitant with XBP1 upregulation in PB2, its downstream targets related to endoplasmic reticulum (ER) stress response, including genes involved in the unfolded protein response such as ATF4, were also significantly upregulated relative to the other subsets (87) . TFs JUN and FOS, though expressed in all subsets, were significantly upregulated in PB1. JUN targets involved in proteasome function were expressed in PB1 and significantly expressed in PB3. Together with the ER stress response, protein degradation pathways support increased immunoglobulin production in PBs. The differing significance of expression for JUN/FOS and their targets across PB1 and PB3, respectively, may reflect negative autoregulation of JUN and FOS within each group. We also noted active proliferation, especially in PB3 (more salient in severe samples) with significantly increased expression of MKI67, EZH2, and CDK1. JUN targets involved in cell cycle regulation (such as SRSF2 and NASP) and microtubule-related genes (such as TUBB) were significantly upregulated in PB3 and may also support increased PB proliferation. Finally, we observed a heterogeneous response to IFN-I stimulation across PB subsets. IRF1/2 and STAT1/3 were expressed in groups PB1 and PB2, respectively (albeit p.adj > 0.05). ISGs also had a subset-dependent expression with ISG15/ 20 significantly upregulated in PB1, and while MX1 and OAS2 were significantly upregulated in PB2. Taken together, our analyses indicated a more heterogeneous PB response, which may reflect more varied maturity and functional status, in severe compared with non-severe COVID-19 patients ( Figure 6J) . Further, while all subsets responded to IFN-I, PBs from severe COVID-19 patients were characterized by robust proliferation and ER stress response, which may support increased immunoglobulin production.

Current research utilizing high-throughput data from either limited cohorts or focused analyses of specific biological systems have generated valuable insights into the pathogenesis of COVID-19. Complementary to these reported observations, in this study, we present system level mechanistic insights into the pathogenesis of COVID-19, by integrating scRNAseq data from four sizable cohorts in both severe and non-severe COVID-19. Integration allowed for higher granularity in identifying and characterizing transcriptomic and cellular heterogeneity of immune response within COVID-19.

Modulation of innate immune cells manifested in several ways, distinguishing the responses between severe and nonsevere COVID-19 ( Figure 7) . Notably, the systemic and robust upregulation of gene signatures in response to IFN-I occurred across multiple cells within both myeloid and lymphoid lineages in both severe and non-severe disease as detected via consensus gene signatures ( Figure 1 ) and expression of IFN-I receptors IFNAR1/2 (Supplementary Figure S9) . Differing reports exist on the activation of a robust IFN-I response in COVID-19 within the PBMCs (88), while there is significant interferon response in lung epithelial cells (2) . We observed that specific subsets that expanded significantly in COVID-19 exhibited a more potent response to IFN-I including the NK4 subset within NK cells, CM1 within CMs, LDGs within LDNs, CD16 + CD8T subset within CD4T/8T, naïve B cells, and PBs. However, this response, specifically in CM1 and LDGs, was more muted in severe compared with non-severe disease. TF factor analysis within these subsets revealed several master regulators, which likely contribute to the observed dysregulation within severe and non-severe disease. Overall, while we cannot comment on the kinetics of IFN-I response, a robust response does exist in COVID-19 albeit muted in severe disease, compared with healthy subjects (12) . Increased expression of S100A8/A9 (alarmins) was seen arising from the abovementioned expanded subsets (such as CM1, LDGs, and CD16 + CD8T) in -19) . This figure encapsulates the main results identified within our manuscript. Rapid expansion of specific monocytic subsets, megakaryocyte, plasmablasts, and low-density neutrophils, along with reduction in cellular frequencies of mucosal-associated invariant T (MAIT), gamma-delta T (gdT), and natural killer cells, is seen within COVID-19, particularly severe disease. The observed reduction of myeloid dendritic cells could be a consequence of emergency monopoiesis, which results in the drastic expansion of suppressive subsets such as CM1 monocytes in severe disease. We observe an increased activation of interferon type I (IFN-I) response, interferon-stimulated genes (ISGs), and alarmins arising from both lymphoid and myeloid cells. However, low-density granulocytes (LDGs) and CM1 show a more muted response to IFN-I in severe disease, compared with nonsevere disease. The increased oxidative stress along with alarmins likely contributes to the suppression of T cells particularly MAIT and gdT. Neutrophils/LDNs have an increased tendency to spontaneously produce neutrophil extracellular traps (NETs), which has been observed in COVID-19 and suggested to contribute to the coagulopathy in COVID-19. Additionally, both monocytes and LDNs presented transcriptional signatures associated with aggregation, especially in severe COVID-19. The increase in MKs in circulation and likely emperipolesis within COVID-19 adds to the mounting evidence on the potent link between thromboembolic events mediated by platelets and their precursors and neutrophils/LDNs. A potential link between immune thrombocytopenia observed in patients and the expansion of MKs/emperipolesis, increased aggregation potential of various immune cell types, and generation of autoantibodies warrant further investigations especially within severe COVID-19 patients. The accompanying bottom panel highlights all differences observed in our manuscript between severe and non-severe disease.

severe disease. Alarmins have diverse roles and contribute extensively to neutrophil recruitment and degranulation and are increasingly acknowledged as crucial markers for COVID-19 infection (89) . The increased oxidative stress arising from LDNs as seen in our analysis [and likely other granulocytes (90) ], along with expression of alarmins, could significantly contribute to T-cell suppression including CD4T, CD8T, MAIT, and gdT, especially in severe disease.

Severe inflammation has been suggested to induce a state of emergency myelopoiesis, particularly monopoiesis and granulopoiesis (28) . We evidenced a drastic expansion of distinct subsets of MDSC (like) cells arising likely due to emergency myelopoiesis such as CM1 within CMs and PMN-MDSCs (through neutrophil progenitor-like cells) within the LDN subsets, especially in severe patients. In certain bacterial infections, experiments have revealed that infection-induced emergency monopoiesis occurs at the expense of DCs, particularly mDC development, and can be sustained for weeks after pathogen clearance (91) . It is therefore possible that the extent of emergency monopoiesis seen in severe COVID-19 (over non-severe COVID-19) contributes to the significant loss of mDCs (p < 0.01, Supplementary Figure S1D ). This loss of abundance taken together with their transcriptionally dysfunctional state in severe disease warrants further investigations of mDCs in the context of SARS-CoV-2 infections.

Nearly 35% of severe COVID-19 patients experience thromboembolic events at multiple sites including the microvasculature, brain, and lung, which can cause organ failure in these patients. Though the precise mechanisms are unclear, several functional links have been proposed, including an interaction between the plasminogen, complement, and platelet-activating systems in severe disease (2) . Both neutrophils and LDNs have an increased tendency to spontaneously produce neutrophil extracellular traps (NETs), which has been observed in COVID-19 and suggested to contribute to coagulopathy in COVID-19 (92, 93) . We identified that both monocytes and LDNs presented transcriptional signatures associated with aggregation as well as homotypic adhesion in severe COVID-19. LDGs and CD14 + monocytes express PTAFR, which indicates PAF activation (a phospholipid crucial for driving platelet aggregation). LDGs also express SELPLG (PSGL-1), a glycoprotein counter-receptor for the cell adhesion molecules P-, E-, and L-selectin, which are required for tethering of leukocytes to activated platelets or endothelia expressing selectins. MKs were found to express the ligand SELP, suggesting interactions between MKs and LDNs, and an active role for MKs in the pathogenesis of severe COVID-19. Our earlier observations on the likelihood of emperipolesis within COVID-19 adds to the mounting evidence for the potent link between thromboembolic events mediated by platelets (aggregation) and neutrophils/LDNs (2, 92) . Current research provides evidence for the role of MK cytoplasm and a-granules as soluble mediators with potential to impact neutrophil migration, behavior, and survival in vivo and in vitro (61) . It is hard to say whether the expansion of MK in severe disease is a cause or effect of increased neutrophil/LDN mobilization and the observed thrombocytopenia. Circulating LDN levels have been reported to correlate with disease severity, contributing to enhanced cytokine production and systemic increases d-dimer, IL-6, and TNF-a levels, in humans with COVID-19 (93, 94) . Ratio of circulating gdT to LDN cells has been suggested to correlate with COVID-19 severity and serves as an early marker of severe disease (94) . Neutrophils are known to mediate inhibition of NK development, function, and homeostasis (95) while also being influenced by MAIT population levels, for a controlled immune response (76, 96) , suggestive of similar roles for LDNs. Unlike other viral infections, we observed no significant changes to majority of the CD4T/CD8T milieu; however, there is evidence for functional exhaustion. Even though there is no significant remodeling of active B cells, naïve B cells are poised for differentiation, and PBs downstream of the B cell maturation spectrum expand drastically within severe COVID-19, in contrast to other respiratory infections.

In summary, our analyses provided novel functional insights into both systemic and cell type-specific changes that contribute significantly to the pathogenicity of COVID-19 in severe and non-severe patients, as shown in Figure 7 . Our analysis of the differential interferon response within immune cell subsets in severe and non-severe COVID-19 allowed for mechanistic insights that help explain prior work. We further highlighted that LDN responses could play a crucial role in COVID-19 pathogenesis and in distinguishing severe from non-severe responses. Particularly, they can serve as a crucial nexus between various myeloid and lymphoid cells affecting their interaction dynamics and contributing extensively to a proinflammatory and pro-thrombotic phenotype seen in severe COVID-19.

We acknowledge several limitations of the current study including limited control over patient cohort selection and initial sample processing techniques due to the use of published data. Factors such as definition of disease severity, sampling time after symptom onset, variability in patient treatment conditions, comorbidities, and sex could all affect clinical presentation and subsequent transcriptomic landscape. As our proposed mechanistic framework is based on integrated analysis of scRNAseq-derived transcriptomic signatures from the PBMCs and assessed based on the current understanding of immune cells, future experiments in an independent cohort are required for targeted insights and further validation.

In our current study, we integrated publicly available scRNAseq data derived from PBMCs of COVID-19 (and healthy) patients. Table S1 ). We utilized Seurat for all scRNAseq processes indicated in this study.

For data from the PA and Lee studies, where the CellRanger outputs were available, count matrices were constructed using Read10X function, and a Seurat object was subsequently generated using Seurat v3.2 (97) . We excluded flu samples from the Lee study prior to Seurat object generation. Likewise, a Seurat object was generated from the Wilk study as indicated in the original publication and its Supplementary Material. However, at the time of data download for our current analysis (September 20, 2020), raw count files were as yet unavailable for download from Schulte-Schrepping et al. (EGAS00001004571) . Subsequently, we downloaded processed Seurat objects from fastgenomics.org (files seurat_COVID19_PBMC_cohort1_10x_jonas_FG_2020-08-15.rds and seurat_COVID19_PBMC_jonas_FG_2020-07-23.rds), as indicated in the original manuscript. We extracted the "counts" slot within the Seurat objects and utilized them for downstream integration. Given that the metadata for each study captured varying amounts of information (see Supplementary  Table S2 ), we restructured the metadata columns to facilitate downstream integration; specifically, we generated four new columns, "orig.ident", which captures the donor information from each study. Multiple samples from the same donor were aggregated for the purposes of this analysis: "Infection_status" (levels: Covid, Healthy) captures the condition of cells/samples, "Severity" (levels: Severe, Non-severe and Healthy) captures the severity of diseases among COVID-19 patients, and "Study" captures the dataset name (Wilk, Lee, SS_C1, SS_C2, and PA). All moderate, mild, and asymptomatic (from Lee et al, see Supplementary Table S2 ) donors across studies were reannotated to "Non-severe," while all "ventilated" donors were annotated to "Severe" (from the Wilk study, Supplementary  Table S2 ), for the sake of consistency. Filtering of the data was performed at cell and gene levels, to obtain reduced datasets, which were then used for downstream integration (see Data Filtering). Only the 14,063 common features across all five datasets were considered for all downstream analyses. Notably, some genes of interest such as IFNB1 are absent from 14,063 features and subsequently absent within our analysis here.

For each of the individual datasets identified above, we performed cell-level filtering by excluding cells based on the following quality criteria: more than 25% mitochondrial reads, more than 25% hemoglobin genes (heme reads), less than 250 expressed genes or more 2 std. deviations from the mean, less than 500 detected transcripts, and a complexity (log10 genes/ UMI) of less than 80%. Additionally, we performed gene-level filtering by excluding genes that were expressed in less than 10 cells. To facilitate downstream integration, we only utilized the 14,063 features/genes common across all datasets for further analysis. Each sample cell contribution after filtering is captured within S2-meta (Supplementary Material).

Each filtered dataset was log normalized prior to integration. The original gene counts for each cell were normalized by total UMI counts, multiplied by 10,000 (TP10K), and then log transformed by log10(TP10k + 1). A total of 2,000 most variable features were detected within each dataset using the vst method implemented in Seurat. Subsequently, the five normalized datasets were integrated using the IntegrateData option of Seurat v3.2, on 2000 gene anchors identified via FindIntegrationAnchors (ndims = 30). The final integrated dataset considered for this analysis contained 14,063 features across a total of 375,438 cells from 111 donors (33 severe, 28 non-severe, and 50 healthy subjects, contributing a total of 116,234, 118,023, and 141,181 cells respectively) ( Figure 1A, Supplementary Table S3 ).

Data scaling and dimensionality reduction (nPCs = 30) were performed on the integrated dataset. Linear regression was performed to remove unwanted variation due to % mitochondrial or heme reads. Optimal dimensions for clustering was identified using "ElbowPlot" function as 24. Subsequently, for two-dimensional data visualization, we performed UMAP based on the first 24 dimensions of the integrated dataset, with the cells clustered using the Louvain algorithm. Given the very large number of cells (~300K), we set the resolution at 2.5 to get meaningful clusters by using the "FindClusters" function.

Cluster identity of the integrated Seurat object was established and confirmed using a threefold strategy. We first ran SingleR, an automatic cluster annotation tool that compares the transcriptome of each single cell to reference datasets to determine cellular identity. We utilized a combination of SingleR's inbuilt references including DatabaseImmuneCellExpressionData, MonacoImmuneData, and NovershternHeamtopoieticData to improve its performance and domain knowledge on established markers of cell types. We additionally identified "cluster markers," which were defined as DEGs, which were either more or less expressed within a chosen subset of cells/clusters compared with all other subsets/clusters within the Seurat object. This test was performed using "FindAllMarkers" function in Seurat with Wilcoxon rank sum test. Genes with >0.25 log-fold changes, at least 25% expressed in tested groups, and Bonferroni-corrected pvalues <0.05 were regarded as cluster markers. The final cluster annotations were determined and verified based on the expression of cluster markers, legacy knowledge, and SingleR annotations. In the situations where multiple clusters were grouped, cluster markers were recomputed using the same methods, at the level of subsets.

To discern the trends in cell-type changes and between cell-type interactions, we grouped the mononuclear cell types by origin: myeloid (MKs, basophils, mDCs, CD14 + CMs, CD16 + non-CMs (NCM), intermediate monocytes (ITM), and LDNs) and lymphoid (CD4 + and CD8 + T-cell subsets, MAIT cells, gamma-delta T(gdT), PBs, B cells, and NK cells). As per the most recent lympho-myeloid model of hematopoiesis (38) , the DC phenotype exists on a spectrum from lymphoid (with plasmacytoid DCs most resembling lymphoid cells) to myelomonocytic (with mDCs most resembling myeloid/monocytic cells). We subsequently annotate the plasmacytoid DCs as belonging to the lymphoid lineage ( Figure 1E ).

The original parent Seurat object was subsampled, for specific cellular compartments (e.g., PBs, LDN, monocytes, and NK cells) to increase cluster annotation granularity. That is, we divided the parent Seurat object into individual Seurat objects comprising cells from each compartment of interest. The integrated space within these subsampled cells was reclustered similar to the processing of the parent Seurat object/primary UMAP. Resolution for each subsampled space was chosen between 0.8 and 1.2. The resolution chosen for a cell type is directly indicated within the results. The cluster identities of the clusters identified within the subsampled spaces (minor clusters) were established using the same strategy as for the parent Seurat object. Some clusters were renamed based on cluster markers identified from SingleR annotation.

Subsampling and reclustering the original 22 monocyte clusters (CD14 + , CD16 + , and ITM) resulted in a total of 19 clusters within the subsampled space, which was further grouped into 11 distinct monocytic subsets, as follows. Based on the expression of markers CD16, CD14, TNFRFA1, and TNFRFB2, clusters 3 and 13 were grouped into NCM ( Figure 2B ), while cluster 11 was identified as ITM based on the expression of markers CD14, CD16, HLA-DR, and TNFRFB2, and other known ITM gene markers including MARCO, APOBEC3A, MARCKSL1, and GBP4. ITMs showed significant enrichment (p.adj < 0.05) for members of the complement activation/signaling (C1QA, C1QB, C3AR1, C5AR1, and SERPINA1) and IFN-I response (Supplementary Table S5 ). Cluster 18 (expressing CD33 and ITGAM as cluster markers (p.adj < 0.05)) was annotated as CD33 + mMDSCs (27) . mMDSCs were markedly absent within healthy samples. Both ITM and mMDSC had the highest cell abundances within severe disease ( Figure 2C , Supplementary Table S8, Supplementary Figure S3A ). Abundance signatures for mMDSCs were verified within an independent patient cohort and are consistent with our observations (Supplementary Figure S3B) .

The remaining 15 clusters were divided into eight subsets of annotated CD14 + CM via SingleR annotation ( Figure 2D) . Several of the CM subsets were reminiscent of subtypes explored in previous publications and are not explored in detail within this manuscript including CM2 (cluster 4, MHC-II high), CM5 (clusters 1 and 12), and CM6 (cluster 5) (5, 6, 10). These subsets in addition to NCM had reduced abundances within COVID-19. Subset CM3 (clusters 7, 16, and 17) was found to have reduced expression of HLA-DR CD14, alarmins, and/or SELL and increased expression of tetraspanins (e.g., CD63 and CD53). The remaining CM subsets are described and analyzed within the main manuscript.

Subsampling and reclustering "LDNs" from the parent Seurat object resulted in nine clusters (Figure 3 , Supplementary Figure  S4 ) grouped into three distinct subsets based on the expression patterns of established markers and cluster markers. Clusters 0, 1, 2, and 3 were identified as LDGs based on the mixed expression of activated mature neutrophil markers including CD11a/b/ITGAL/ITGAM, CD55, CD16, CD10/MME, and SELL and immature markers including CD16 high/low and CD10 low (33) . Clusters 4, 6, and 8 were identified as PMN-MDSCs, based on the significant expression (p.adj < 0.05) of several known markers including ITGAM, LCN2, CAMP, MMP8, ARG1, S100A8, and S100A9 (37, 38) and were subsequently grouped together ( Figure 3E ). OLR1 (LOX1), a recently validated PMN-MDSC marker, was exclusively increased within this subset (p.adj < 0.05) (38) (Figures 3E, F) . Clusters 5 and 7 expressed FUT4 (CD15), CD63, PRTN3, MPO, and ELANE reminiscent of a "proneutrophil" state (5).

Subsetting and reclustering NK cells from the primary UMAP resulted in 23 distinct clusters (Supplementary Figure S5A) . Three clusters (clusters 4, 6, and 19) showed significant enrichment of CD3 (CD3D/G) suggestive of an NKT (invariant NK-T) like population. Because most NK cells express CD7 + , one CD7 − cluster (cluster 16) was excluded from further analysis (Supplementary Figure S5B) (64, 65) . Reclustering and re-embedding of the Seurat object after processing from the original Seurat object resulted in 18 clusters ( Figure 4A ), grouped into seven distinct subsets based on previously published NK markers and considered for further analysis (Figures 4B, C, and Supplementary Figure S5C ). NK4 (cluster 8), NK5 (cluster 14), and NK8 (clusters 2, 9, and 12) were identified as CD56 low subsets. Subsets NK3 (clusters 4-7, 10-11, and 17) and NK7 (clusters 0 and 1) largely lacked expression of CD56 (CD56 neg ) (Supplementary Figure S5D) . Specifically, NK5, NK7, and NK8 represented activated/mature NK clusters based on the expression of cytotoxic markers including perforin (PRF1), granzymes (GZMA, GZMB, and GZMK), and remodeling markers including ACTB, ARPC3/4, CFL1, and CST7. NK3 showed an increased expression of inhibitory receptors including KIR2DL1/3, KIR3DL1/L2, KLRC2, IL32, and CD3E and cytolytic molecules including GZMH and CST7 and lower levels of FCER1G, in line with a profile of adaptive NK cells within increased cytotoxic potential (Supplementary Table  S16 ). Two subsets NK6 (cluster 3,~6.9% of total NK cells) and NK2 (clusters 16 and 15,~6.5% of total NK cells) were characterized by increased expression of CD56. NK6 expressed surface proteins commonly attributed to CD56 bright population including CD44, SELL, and GZMK chemokines XCL1, XCL2, DUSP1, FOS, JUND, IL7R, and IL2RB (64).

Subsampling and reclustering the original 19 CD4/CD8T (and naïve) clusters comprising 128,589 cells resulted in a total of 22 clusters within the subsampled space. Using the commonly used subtyping markers, we observed four clusters (clusters 5, 19, 20, and 21) to be non-T specific (with a small population of CD3 + T cells) and were removed. Additionally, we observed that cluster 5 was majorly contributed by healthy samples, from only SS_C1, and was subsequently eliminated. The remaining 18 clusters (119,683 cells) were further grouped into eight distinct subsets ( Figures 5A, B) using automatic annotation and verified based on the expression of commonly used T subtyping markers (Supplementary Figures S6A, B) (98) . Specifically, based on the expression of naïve markers CCR7, SELL, and TCF7, seven CD4 + clusters (clusters 1, 4, 7, 11, 12, 14, and 15) were grouped into the naïve CD4T subset, while a single CD8 + cluster (cluster 10) was identified as the naïve CD8 T-cell subset. Likewise, based on the expression of granzymes (GZMA/B), PRF1, and NKG7, clusters 2, 3, and 9 were identified to be terminally differentiated effector (Teff) CD8T cells, while CD4 + cluster 3 was identified as Teff CD4T. Clusters 6 and 13 were grouped into Treg CD4T subset based on the expression of CTLA4, FOXP3, IL2R, and ICOS (Supplementary Figure S6B) . Consistent with SingleR annotation, cluster 0 was identified as TFH (GATA3 + ), and cluster 17 was identified as memory T cells (Tmem, CD27 + , IL7R + , CD28 + , and FAS/CD95 + ). Studies have broadly identified increased exhaustion with increasing severity in the effector Tcell subset of COVID-19 patients (14) . Recent research has additionally demonstrated a presence of clonally distinct hybrid memory T-cell subpopulation with an exhausted phenotype (GZMK + and TOX + with markers of exhaustion including PDCD1, MAF, and NFATC2), which is in contrast to the classical understanding that memory and exhausted T cells arise from segregated pathways (99) . The CD8 + Tmem population identified in our study exhibits signatures consistent with this clonally distinct population, indicative of the transition of a subset of exhausted T cells to a memory stage in both severe and non-severe COVID-19 (cells in healthy = 361, non-severe = 656, and severe = 622).

We extracted the B-cell subsets (B lymphoblasts, naïve B cells, and non-switched memory B cells) from the original parent Seurat object (Supplementary Figure S8A) and identified the most significant remodeling occurring only within the naïve cell compartment. Subsampling and reclustering of naïve B cells resulted in 11 initial clusters (Supplementary Figure 8B) . Out of the 11 initial clusters, we excluded two (clusters 6 and 7), as their most significantly upregulated cluster markers indicated potential erythroid and myeloid lineage cells and were not entirely consistent with naïve B cells. Of the remaining nine clusters ( Figure 6A) , we verified expression of surface markers consistent for naïve B cells (Supplementary Figures S8C, D) . We further grouped these clusters into groups A-F, based on gross gene expression (via hierarchical clustering of scaled average gene expression values expressed in greater than 60% of all cells) in each cluster. In multiple instances, this unsupervised clustering method grouped together clusters that comprised similar proportions of healthy, non-severe, and severe cells ( Figure 6B ). In groups A (clusters 2 and 9) and B (clusters 5 and 8), non-severe cell numbers were much greater than in both healthy and severe cells. In group C (cluster 1), non-severe cells still dominate the cluster, but are followed by severe and then healthy cells. In group D (clusters 0 and 4), healthy and non-severe cells outnumber severe cells. Finally, in groups E (cluster 3) and F (cluster 10), severe cells are greater than non-severe and healthy cells.

Subsampling and reclustering PBs from the original parent object resulted in 10 clusters (Supplementary Figure S8F) , of which six showed surface marker expression most characteristic of PBs (CD19 − , MS4A1 − , CD27 + , and CD38 + ) and were retained for further analysis (Supplementary Figure S8G) (Figures 6E-G) . The six clusters were grouped into four subsets PB1 (clusters 0 and 2), PB2 (clusters 3 and 6), PB3 (cluster 1), and PB4 (cluster 5) based on expression patterns among the top 60 cluster markers. For clarity, the expression of the top 20 genes from each cluster is shown in Supplementary Figure S8H . The clusters comprising subsets PB1, PB2, PB3, and PB4 and the number of healthy, non-severe, and severe cells in each subset are detailed in Supplementary Figure S8I . Since PBs from healthy patients comprised only about 5% of all PBs, we considered cells only from severe and non-severe patients in our analysis. In addition to an expanded PB population relative to non-severe patients, the PB population in severe patients was also more heterogeneous: no single subset dominated the PB response in severe disease, as multiple subsets were expanded; by contrast, subset PB1 dominated the PB response in non-severe disease.

Dotplots used within this manuscript were all generated using the "DotPlot" function of Seurat. The average expression presented within each dotplot is the scaled/standardized expression values as defined by Seurat. DGEA was performed with respect to either healthy or non-severe samples, depending on the context of analysis, using the "FindMarker" function within Seurat. The thresholds for calling DEGs were as discussed above. Heatmaps of the DEGs highlighting the fold changes across conditions of interest were performed using "pheatmap" library in R/BioC. Gene Ontology (GO) enrichment on the DEGs identified was performed using ClusterProfiler v.3.10.1 (100) or Enrichr (101) . The top 10-15 categories of GO were utilized for visualization throughout the manuscript and Supplementary Material. Cluster profiler was also utilized for visualization of enrichment.

To characterize the transcriptional regulation of the altered genes programs within the subset, we queried TF activity via DoRothEA (Discriminate Regulon Expression Analysis) (102) , which utilizes the viper activity inference algorithm against a curated list of TFs and its target expression levels to predict TF activity (regulon). DoRothEA regulons have been generated from various sources and evidence types, and we subsequently included sets with "A" or "B" confidence ratings. We visualized the top 20 genes ranked by DoRothEA for visualization. We also utilized ChEA3 (103) to rank TFs. In contrast to DoRothEA, which utilizes the entire expression matrix to derive putative TFs involved, ChEA3 predicts TF enrichment based on the overlap between given lists of DEGs or genes of interest, and TF targets assembled from a compendium of resource. ChEA3 utilizes a Fisher's exact test, with a background size of 20,000, to compare the input gene set to the TF target-gene sets to establish significance. We utilized the top 10 highly ranked TFs defined via the "TopRank" metric for further visualization.

Immunophenotyping of Myeloid-Derived Suppressor Cells, Mucosal-Associated Invariant T Cells, and Gamma-Delta T Cells

Human Subjects: Blood samples from healthy subjects were enrolled and tested under Plexision IRB-approved protocol 6774. COVID-19 patients were enrolled under IRB-approved protocol #1551551 from Edinburg, TX. Patients with COVID-19 infection are categorized into three groups, COVID-19 patients who required no hospitalization (Mild), those who required hospitalization with no mechanical ventilation (Hospitalized), and those who were hospitalized and required ventilation for oxygen requirements (Severe). Patient cohort demographics and severity details are presented within Supplementary Table S19 .

To provide independent validation of expansion of certain immune cell populations in severe and non-severe COVID-19 and healthy controls, we performed flow cytometric analyses on PBMCs from the above independent patient cohort. Cells from PBMCs were labeled with fluorochrome-labeled antibodies to characterize mMDSC and PMN-MDSC subsets. The respective phenotypes were CD14 + HLA-DR − and CD14-CD15 + CD11b + . MAIT cells were gated based on CD3 + TCRgd-TCRVa7.2 + CD161 + cells, and gdT cells were identified based on CD3 + TCRgd + cells. Antibodies from BioLegend (San Diego, CA) and BD Biosciences (San Jose, CA) were utilized.

The current knowledgebase of all RBPs was downloaded from RBPBase (104) containing 4,257 RBPs within humans. Of these, we considered only those RBPs that are present in the top 25% of the hits_Hs (>6), bringing down the list to 1,031 high-confidence RBPs. Consensus genes that overlapped with this list are provided in Supplementary Figure S2C .

The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding author.

The studies involving human participants were reviewed and approved by IRB-approved protocol #6774 and IRB-approved protocol #1551551. The patients/participants provided their written informed consent to participate in this study.

Immunological Characteristics Govern the Transition of COVID-19 to Endemicity

Plasmin Cascade Mediates Thrombotic Events in SARS-CoV-2 Infection via Complement and Platelet-Activating Systems

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

Immunophenotyping of COVID-19 and Influenza Highlights the Role of Type I Interferons in Development of Severe COVID-19

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

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

Deep Immune Profiling of COVID-19 Patients Reveals Distinct Immunotypes With Therapeutic Implications

Comorbidities and Multi-Organ Injuries in the Treatment of COVID-19

Impaired Natural Killer Cell Counts and Cytolytic Activity in Patients With Severe COVID-19

A Single-Cell Atlas of the Peripheral Immune Response in Patients With Severe COVID-19

Comprehensive Mapping of Immune Perturbations Associated With Severe COVID-19

Cell-Type-Specific Immune Dysregulation in Severely Ill

Multi-Omics Resolves a Sharp Disease-State Shift Between Mild and Moderate COVID-19

Reduction and Functional Exhaustion of T Cells in Patients With Coronavirus Disease 2019 (COVID-19)

Impaired Type I Interferon Activity and Inflammatory Responses in Severe COVID-19 Patients

COVID-19 Immune Features Revealed by a Large-Scale Single-Cell Transcriptome Atlas

Cell Subsets as Severity-Associated Signatures in COVID-19 Patients

Immunological Memory to SARS-CoV-2 Assessed for Up to 8 Months After Infection

Broad and Strong Memory CD4+ and CD8+ T Cells Induced by SARS-CoV-2 in UK Convalescent Individuals Following COVID-19

Systems Biological Assessment of Immunity to Mild Versus Severe COVID-19 Infection in Humans

Single-Cell Multi-Omics Analysis of the Immune Response in COVID-19

Editor-in-Chief

SARS-CoV-2 Dysregulation of PTBP1 and YWHAE/Z Gene Expression: A Primer of Neurodegeneration

Cytokine and Chemokine Levels in Coronavirus Disease 2019 Convalescent Plasma. Open Forum Infect Dis

Proteins Associated With Neutrophil Degranulation Are Upregulated in Nasopharyngeal Swabs From SARS-CoV-2 Patients

Identification of Novel Human Monocyte Subsets and Evidence for Phenotypic Groups Defined by Interindividual Variations of Expression of Adhesion Molecules

The Endless Saga of Monocyte Diversity

Plasma From Patients With Bacterial Sepsis or Severe COVID-19 Induces Production of Suppressive Myeloid Cells From Human Hematopoietic Progenitor Cells In Vitro

Single-Cell Analysis of Two Severe COVID-19 Patients Reveals a Monocyte-Associated and Tocilizumab-Responding Cytokine Storm

MyD88-Dependent Expansion of an Immature GR-1+ CD11b+ Population Induces T Cell Suppression and Th2 Polarization in Sepsis

Human Monocytes Undergo Functional Re-Programming During Sepsis Mediated by Hypoxia-Inducible Factor-1a

FOSL1 Inhibits Type I Interferon Responses to Malaria and Viral Infections by Blocking TBK1 and TRAF3/ TRIF Interactions

FOSL2 Positively Regulates TGF-b1 Signalling in Non-Small Cell Lung Cancer

TGF-b1 Suppresses the Type I IFN Response and Induces Mitochondrial Dysfunction in Alveolar Macrophages

Chronic Psoriatic Skin Inflammation Leads to Increased Monocyte Adhesion and Aggregation

Enhanced T-Cell Maturation and Monocyte Aggregation Are Features of Cellular Inflammation in Human T-Lymphotropic Virus Type 1-Associated Myelopathy

Platelet Activation and Platelet-Monocyte Aggregate Formation Trigger Tissue Factor Expression in Patients With Severe COVID-19

Human Dendritic Cell Subsets: An Update

COVID-19 Severity Associates With Pulmonary Redistribution of CD1c+ DCs and Inflammatory Transitional and Nonclassical Monocytes

CD83: Activation Marker for Antigen Presenting Cells and Its Therapeutic Potential

Nuclear Receptor Expression Atlas in BMDCs: Nr4a2 Restricts Immunogenicity of BMDCs and Impedes EAE

The Role of Fcϵri Expressed in Dendritic Cells and Monocytes

Monocytes and Macrophages as Viral Targets and Reservoirs

Ilr1b Production is Dependent on the Activation of Purinergic Receptors and NLRP3 Pathway in Human Macrophages

On the Origin of Low-Density Neutrophils

Neutrophil Diversity in Health and Disease

Identification of an Early Unipotent Neutrophil Progenitor With Pro-Tumoral Activity in Mouse and Human Bone Marrow

Defining the Emergence of Myeloid-Derived Suppressor Cells in Breast Cancer Using Single-Cell Transcriptomics

LEF-1 is a Decisive Transcription Factor in Neutrophil Granulopoiesis

Human Neutrophils in the Saga of Cellular Heterogeneity: Insights and Open Questions

Aquaporin 9 Phosphorylation Mediates Membrane Localization and Neutrophil Polarization

Type I Interferon Transcriptional Signature in Neutrophils and Low-Density Granulocytes Are Associated With Tissue Damage in Malaria

Low-Density Granulocytes Activate T Cells and Demonstrate a non-Suppressive Role in Systemic Lupus Erythematosus

Down-Regulation of CXCR2 on Neutrophils in Severe Sepsis is Mediated by Inducible Nitric Oxide Synthase-Derived Nitric Oxide

Activated Neutrophils Exert Myeloid-Derived Suppressor Cell Activity Damaging T Cells Beyond Repair

Human Megakaryocytes Possess Intrinsic Antiviral Immunity Through Regulated Induction of IFITM3

Increased Number of Pulmonary Megakaryocytes in COVID-19 Patients With Diffuse Alveolar Damage: An Autopsy Study With Clinical Correlation and Review of the Literature

Elevated Expression of IL-12 and IL-23 in Patients With Primary Immune Thrombocytopenia

Thrombocytopenia Is Associated With Severe Coronavirus Disease 2019 (COVID-19) Infections: A Meta-Analysis

Megakaryocyte Emperipolesis: A New Frontier in Cell-in-Cell Interaction

Increased and Pathologic Emperipolesis of Neutrophils Within Megakaryocytes Associated With Marrow Fibrosis in GATA-1low Mice

Natural Killer Cell Immunotypes Related to COVID-19 Disease Severity

Heterogeneity of Human Bone Marrow and Blood Natural Killer Cells Defined by Single-Cell Transcriptome

Diversity of Peripheral Blood Human NK Cells Identified by Single-Cell RNA Sequencing

CD56bright Natural Killer (NK) Cells: An Important NK Cell Subset

A NKp80-Based Identification Strategy Reveals That CD56neg NK Cells are Not Completely Dysfunctional in Health and Disease

A New Hope for CD56negCD16pos NK Cells as Unconventional Cytotoxic Mediators: An Adaptation to Chronic Diseases

Targeting Foxp1 for Reinstating Anticancer Immunosurveillance

Type I IFN Promotes NK Cell Expansion During Viral Infection by Protecting NK Cells Against Fratricide

Marked T Cell Activation, Senescence, Exhaustion and Skewing Towards TH17 in Patients With COVID-19 Pneumonia

Regulatory T Cells and Immune Tolerance

Elevated Numbers of Fcgriiia+ (CD16+) Effector CD8 T Cells With NK Cell-Like Function in Chronic Hepatitis C Virus Infection

Single-Cell Analyses Identify Dysfunctional CD16+ CD8 T Cells in Smokers

Paracrine Costimulation of IFN-g Signaling by Integrins Modulates CD8 T Cell Differentiation

MAIT Cells in Health and Disease

The MERS-CoV Receptor DPP4 as a Candidate Binding Target of the SARS-CoV-2 Spike

Expansion of Plasmablasts and Loss of Memory B Cells in Peripheral Blood From COVID-19 Patients With Pneumonia

Maintenance of the Marginal-Zone B Cell Compartment Specifically Requires the RNA-Binding Protein ZFP36L1

Gadd45b is a Pro-Survival Factor Associated With Stress-Resistant Tumors

Antibody-Secreting Cell Destiny Emerges During the Initial Stages of B-Cell Activation

CCR10 Expression is a Common Feature of Circulating and Mucosal Epithelial Tissue IgA Ab-Secreting Cells

CXCR4-CXCL12 Interaction is Important for Plasma Cell Homing and Survival in NZB/W Mice

Gene Expression Profiling of Plasma Cells and Plasmablasts: Toward a Better Understanding of the Late Stages of B-Cell Differentiation

Ufbp1 Promotes Plasma Cell Development and ER Expansion by Modulating Distinct Branches of UPR

The Type I Interferon Response in COVID-19: Implications for Treatment

Induction of Alarmin S100A8/ A9 Mediates Activation of Aberrant Neutrophils in the Pathogenesis of COVID-19

Devilishly Radical NETwork in COVID-19: Oxidative Stress, Neutrophil Extracellular Traps (NETs), and T Cell Suppression

Innate Immune System Favors Emergency Monopoiesis at the Expense of DC-Differentiation to Control Systemic Bacterial Infection in Mice

Neutrophil Extracellular Traps and Thrombosis in COVID-19

A Specific Low-Density Neutrophil Population Correlates With Hypercoagulation and Disease Severity in Hospitalized COVID-19 Patients

Whole Blood Immunophenotyping Uncovers Immature Neutrophil-to-VD2 T-Cell Ratio as an Early Marker for Severe COVID-19

The Defensive Alliance Between Neutrophils and NK Cells as a Novel Arm of Innate Immunity

Neutrophils Suppress Mucosal-Associated Invariant T Cells in Humans

Single-Cell Transcriptomics of Human T Cells Reveals Tissue and Activation Signatures in Health and Disease

Two Subsets of Stem-Like CD8+ Memory T Cell Progenitors With Distinct Fate Commitments in Humans

Clusterprofiler: An R Package for Comparing Biological Themes Among Gene Clusters

Enrichr: A Comprehensive Gene Set Enrichment Analysis Web Server

Robustness and Applicability of Transcription Factor and Pathway Analysis Tools on Single-Cell RNA-Seq Data

Transcription Factor Enrichment Analysis by Orthogonal Omics Integration

RNA-Binding Proteins in Human Genetic Disease

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2021.738073/ full#supplementary-material Conflict of Interest: Phenotyping of the independent cohort was performed at Plexision, a University of Pittsburgh spinoff in which RS and the University hold equity. RS and CA are currently employed at Plexision.The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.Publisher's Note: All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.Copyright © 2021 Mukund, Nayak, Ashokkumar, Rao, Almeda, Betancourt-Garcia, Sindhi and Subramaniam. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.