key: cord-273587-nja58vxw
authors: Rendeiro, A. F.; Casano, J.; Vorkas, C. K.; Singh, H.; Morales, A.; DeSimone, R. A.; Ellsworth, G. B.; Soave, R.; Kapadia, S. N.; Saito, K.; Brown, C. D.; Hsu, J.; Kyriakides, C.; Chui, S.; Cappelli, L.; Cacciapuoti, M. T.; Tam, W.; Galluzzi, L.; Simonson, P. D.; Elemento, O.; Salvatore, M.; Inghirami, G.
title: Longitudinal immune profiling of mild and severe COVID-19 reveals innate and adaptive immune dysfunction and provides an early prediction tool for clinical progression
date: 2020-09-09
journal: medRxiv : the preprint server for health sciences
DOI: 10.1101/2020.09.08.20189092
sha: 
doc_id: 273587
cord_uid: nja58vxw

With a rising incidence of COVID-19-associated morbidity and mortality worldwide, it is critical to elucidate the innate and adaptive immune responses that drive disease severity. We performed longitudinal immune profiling of peripheral blood mononuclear cells from 45 patients and healthy donors. We observed a dynamic immune landscape of innate and adaptive immune cells in disease progression and absolute changes of lymphocyte and myeloid cells in severe versus mild cases or healthy controls. Intubation and death were coupled with selected natural killer cell KIR receptor usage and IgM+ B cells and associated with profound CD4 and CD8 T cell exhaustion. Pseudo-temporal reconstruction of the hierarchy of disease progression revealed dynamic time changes in the global population recapitulating individual patients and the development of an eight-marker classifier of disease severity. Estimating the effect of clinical progression on the immune response and early assessment of disease progression risks may allow implementation of tailored therapies.

Coronavirus disease-2019 (COVID- 19) , caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), is a global pandemic that (as of August 2020) has infected over 25 million people worldwide, caused over 840,000 deaths, and strains health systems on an unprecedented scale. COVID-19 has heterogeneous clinical manifestation, ranging from mild symptoms such as cough and low-grade fever to severe conditions including respiratory failure and death 1, 2 . While most patients with mild disease develop an appropriate immune response that culminates with viral clearance 3 , severe disease manifestations have been linked to lymphopenia and immune hyperresponsiveness leading to cytokine release syndrome (CRS) 1, 2, 4, 5 . The most effective therapeutic approaches developed so far for severe cases involve either general immunosuppression with glucocorticoids 6 or selective neutralization of interleukin 6 (IL-6) with tocilizumab 7 , a monoclonal antibody used to manage CRS in indications such as rheumatoid arthritis 8 . The efficacy of these therapies strongly support a key role for immune dysregulation in the pathogenesis of COVID-19. However, neither treatment has achieved high clinical remission rates in patients with severe COVID-19 9, 10 , suggesting that other immunological or immune-independent attributes may contribute to severity, treatment failure, and ultimately patient death. Thus, in-depth characterization of immune responses to SARS-CoV-2 infection is urgently needed.

Recent characterization efforts have uncovered broad dysregulation of the innate immune system 11 coupled with altered inflammatory responses 12 and impaired adaptive immunity 13 . Specifically, the adaptive immune compartment of COVID-19 patients exhibits marked lymphopenia 4, 14, 15 , polarization of T cells toward a memory phenotype 14 and functional exhaustion [16] [17] [18] , demonstrating that SARS-CoV-2 infection induces both cellular 19, 20 and humoral responses 21 . However, the molecular and cellular mechanisms through which SARS-CoV-2 infection induces these broad immunological derangements in only some patients remains to be elucidated. Further, little is known about the role of the innate immune responses that constitute the first defense against SARS-CoV-2 infection. Moreover, the degree of interaction between various immune compartments and demographic factors and medical comorbidities is unclear. The most prominent risk factors for severe disease and death by COVID-19 include age, cardiovascular or oncological comorbidities, and immunosuppression [22] [23] [24] . Additionally, men appear to be at significantly higher risk for severe COVID-19 than women 25 . While mortality rates are estimated at 4-6% in the general population, high-risk populations experience mortality rates >60% 26 .

Clarifying the early immunological alterations associated with mild versus severe COVID-19 may not only offer therapeutically actionable targets, but also enable the identification of cases at highest risk for clinical deterioration and death. The development of an effective clinical decision-making tool rooted in immunological monitoring has the potential to optimize patient care and resource utilization.

By profiling mild and severe COVID-19 patients and healthy donors with flow cytometry, we demonstrate that SARS-CoV-2 is associated with broad dysregulation of the circulating immune system, characterized by the relative loss of lymphoid cells coupled to expansion of myeloid cells. Severe cases demonstrated enrichment of natural killer (NK) cells expressing the immunosuppressive receptor killer cell immunoglobulin-like receptor, two Ig domains and short cytoplasmic tail 4 (KIR2DS4; CD158i), as well as alterations in the B cell compartment marked by reduced CD19, CD20, and IgM+ cells. These immune profiles enable reconstruction of a hierarchy of disease progression with pseudo-temporal modeling, which allows estimation of dynamic longitudinal changes within individual patients. Our approach also estimates the effect of clinical factors on immune dysregulation and thus establishes an immune-monitoring tool for disease progression.

We conducted an observational study of 45 individuals with COVID-19 that were treated at New Y Presbyterian Hospital and Lower Manhattan Hospitals, Weill Cornell Medicine (IRB 20-03021645) in-or outpatients between April and July, 2020. The disease was categorized as "mild" if the pati was not admitted or required <6 L non-invasive supplemental oxygen to maintain SpO 2 >92% (n = 2 Patients with "severe" disease required hospitalization and received >6 L supplemental oxygen mechanical ventilation (n = 15). Blood samples were collected at enrollment and, when permissib approximately every 7 days thereafter. Samples were also collected from non-hospitalized individu who had recovered from mild, laboratory-confirmed SARS-CoV-2 infection ("convalescent" group, n 9) and from healthy COVID-19 negative donors (n = 12) (Figure 1a) . The median age of COVID patients was 65 years, which was significantly higher than healthy donors (30 years) (Supplementa Tables 1-2; Supplementary Figure 1 ). Figure 1 : Immuno-profiling of COVID-19 patients reveals a disarrayed immune system. a) Composition of the study cohort. b) Description of immune panels and their target epitopes. c) Composition major immune compartments as a percentage of all live CD45 + cells. d) Abundance of major lymph compartments as a percentage of all lymphocytes. For (c) and (d), the upper panels divide patients by gene disease status and three lower panels further divide the study subjects by clinical intervention or outcom Significance was assessed using Mann-Whitney U tests and corrected for multiple testing with the Benjam Hochberg false discovery rate (FDR). **, FDR-adjusted p-value <0.01; *, FDR-adjusted p-value of 0.01-0.05.

We performed high-dimensional immune cell profiling of circulating blood by flow cytometry based seven independent fluorochrome-conjugated antibody panels, each targeting a specific surface prot 3/27 Supplementary Tables 3-7) . Longitudinal sampling was performed in eight patients, one in the "mild" and seven in the "severe" group, and included at least three samples per patient, making a complete dataset including 102 samples from 57 individuals.

Consistent with previous reports 12, 15, 27, 14 , we observed global loss of lymphocytes among CD45 + cells and enrichment of the myeloid cell compartment in the peripheral blood of COVID-19 patients compared with healthy donors (Figure 1c, top) . This was exacerbated in patients with severe disease, particularly among hospitalized patients who required mechanical ventilation or died, compared with individuals with mild disease (Figure 1c, bottom) . This lymphocyte depletion was primarily observed in the T and NK cell compartments (Figure 1d) . There was no difference in the abundance of B cells between mild and severe groups. These results highlight a major shift in peripheral immune cell absolute abundance from the lymphoid to myeloid lineage (Supplementary Tables 6-7) .

We next profiled CD4+ and CD8+ T-cells in COVID-19 patients and healthy donors (Supplementary Figure 2) . The CD4/CD8 ratio correlated positively with disease severity (Figure 2a ) 14, 19 . There was also an expansion of memory T cells (CD45RO + ) with reciprocal contraction of the naive compartment (CD45RA+) in severe cases relative to mild disease or healthy donors (Figure 2b) . We next quantified the abundance of populations expressing C-C motif chemokine receptor 7 (CCR7; CD197), selectin L (SELL; CD62L), and FAS cell surface death receptor (CD95). Within CD45RA + cells, effector CCR7 -(T EFF ) populations were increased in COVID-19 patients and those with severe disease, especially in the CD8 + compartment (Figure 2c) . Conversely, there was significant depletion of CD8 + CD45RO + CD95 -T cells in patients, which was exacerbated with severe disease.

To characterize these populations more objectively and independently of manual gating, we analyzed the expression of eight surface proteins at the single-cell level by jointly embedding CD3 + cells from all samples with the Uniform Manifold Approximation and Project (UMAP) method and clustering them (Figure 2d ). Not all clusters contained cells from all severity groups proportionally. Specifically, clusters 12, 18, and 21, which are characterized by reduced FAS expression, were enriched for cells from healthy donors (Figure 2e) . Moreover, there was increased expression of CD95 in samples from COVID-19 patients that correlated with disease severity (Figure 2f-g) , and FAScells were particularly depleted in all patients (Figure 2g) . Indeed, CD95 + CD25 + T cells were increased in severe cases, while no difference was observed between convalescent patients and healthy donors (Figure 2h ).

. CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

The copyright holder for this preprint this version posted September 9, 2020. . https://doi.org/10.1101/2020.09.08.20189092 doi: medRxiv preprint is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

The copyright holder for this preprint this version posted September 9, 2020. . https://doi.org/10.1101/2020.09.08.20189092 doi: medRxiv preprint 6/27

Next, we assessed the frequency of CD4 regulatory T cells (T REG , characterized by CD127 dim CD25 bright ). As markers of follicular helper T cells (T FH ), we also measured CD4 + , CXCR5 + , PD1 + , ICOS + T FH , which are critical to B cells in the initiation and maintenance of humoral immune responses 28 . We found a significant but modest increase in T FH in mild and severe COVID-19 cases, with their presence in convalescent patients similar to in healthy donors (Supplementary Figure  4a/b) . However, upon considering a broader spectrum of T FH cells regardless of ICOS expression (Supplementary Figure c) , CD4 + , CXCR5 + , PD-1 + T FH were most abundant in COVID-19 patients with mild disease (Figure 2i) . In T REG , severe COVID-19 patients showed significant increase compared with healthy donors (Supplementary Figure 4c) , while previous reports showed an increase in patients with mild course 5, 29 .

To investigate T cell functional phenotypes, we assessed the expression of co-inhibitory T cell receptors. We observed sustained increase of Programmed Cell Death 1 (PD-1) in COVID-19 patients compared with healthy donors in both CD4 and CD8 compartments. At the same time, V-set immunoregulatory receptor (VISTA) and lymphocyte-activating gene 3 (LAG3) were upregulated in mild cases (Figure 2j ). Exhausted T cell phenotypes, with high expression of VISTA and LAG-3, can be encountered in chronic viral diseases 30 , including chronic SARS-CoV-2 infection 17 . This phenotype suggests that these inhibitory receptors may operate at least partially via non-overlapping immunosuppressive signals that negatively regulate T cell responses during chronic viral infection 16 .

These results highlight a shift toward an activated T cell memory phenotype in COVID-19 patients, with a potential role for CD95-mediated cell death. By and large, convalescent patients and healthy donors displayed similar immunotypes in comparison with COVID-19 patients. However, we did identify populations such as CD45RA + , CCR7 + , CD62L -, FAS -CD8+ T EFF cells, which remained significantly different to healthy donors up to ~2 months into recovery (Figure 2c ). These cells may represent "T stem memory (T SM ) cells" with poor expansion potential 31 and/or aberrant terminally differentiated effector memory (T EM ) cells 32 .

Having observed myeloid expansion in COVID-19 patients (Figure 1c) , we next investigated the abundance of the MDSC subset. These elements are activated by interleukin 6 (IL-6) 33 and have immunomodulatory functions in cancer 34, 35 and viral infections 36 . Our flow cytometry panel considered CD3 -, CD56 -, CD19 -, HLA-DR -/dim , CD33 + , CD11b + cells and focused on distinguishing CD14 -, CD15 + granulocytic cells (G-MDSCS); CD14 + , CD15 -/dim monocytic-like cells (M-MDSC); and CD14 -, CD15 -/dim immature cells (I-MDSC) from each other. G-MDSCs were rarely detected in healthy donors but were prevalent in mild and severe COVID-19 patients (Figure 3a/b) . Convalescent patients showed numbers of G-and M-MDSCs closer to healthy donors, with a non-significant increase in I-MDSCs compared with healthy donors. Conversely, I-MDSC cells, while relatively rare as a fraction of all immune cells, were further reduced with disease (Figure 3a/b ). Since neutrophils are phenotypically similar to MDSCs, we compared their abundance with MDSCs. While there was a positive correlation between G-MDSCs and a high neutrophil count, neither could account entirely for the other (Supplementary Figure 3) . Next, we created a joint embedding of 2.4 million CD16 + cells from all samples using the UMAP method, deriving clusters based on similar cells (Figure 3c) . Clusters containing CD15 + cells were disproportionately enriched in samples from COVID-19 patients, while clusters with CD3 + , IL4R (CD124) were mostly composed of cells from healthy donors (Figure 3d ). In addition, CD15 expression was most prominent in COVID-19 patients, particularly in severe cases, but when selecting . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted September 9, 2020. . https://doi.org/10.1101/2020.09.08.20189092 doi: medRxiv preprint for CD3or CD3 -CD33+ cells, convalescent patients possessed a number of CD15 + cells m similar to patients with active disease than healthy donors. is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted September 9, 2020. . https://doi.org/10.1101/2020.09.08.20189092 doi: medRxiv preprint 8/27 severity (bottom). h) Expression levels of all four measured KIR receptors in each disease state. Significance was assessed using Mann-Whitney U tests and corrected for multiple testing with the Benjamini-Hochberg false discovery rate (FDR). **, FDR-adjusted p-value <0.01; *, FDR-adjusted p-value 0.01-0.05.

Next, we focused on innate lymphoid cells and determined the expression of KIR receptors in CD56 + , CD16 bright NK cells. While we observed no significant differences in the relative abundance of KIR receptors among COVID-19 patients with mild disease and healthy controls (Figure 3f) , a significantly higher proportion of cells expressed CD158i (NKG2A) in severe patients compared with mild or convalescent individuals. Moreover, we observed fewer CD158e (KIR3DL1) cells in patients with mild disease compared with severe patients and a lower proportion of cells not expressing any of the measured receptors (KIR -) in patients with severe disease. To further explore NK cell subsets independent of conventional gating, we harnessed single-cell analysis and integrated >500,000 cells in a UMAP representation, identifying cell clusters based on surface marker expression (Figure 3g ). Clusters significantly enriched in CD158i expressing cells were paucicellular in healthy donors compared with COVID-19 patients (Figure 3g, bottom) , and the relative frequency of CD158iexpressing cells was lower in healthy donors, regardless of the expression of other KIR receptors (Figure 3h ). Since the expression of KIR variants is stochastic, the apparent selection of KIRexpressing cells in severe COVID-19 patients could indicate that a viral antigen presented by major histocompatibility complex (MHC) class I molecules with higher affinity for CD158i could select for NK cells expressing this receptor.

Since B cells play a critical role in adaptive immunity, we investigated the expression levels of surface CD19, CD20, IgM, and IgG in circulating cells. Despite the backdrop of a relative decrease in B cell numbers as disease progresses, we observed only a mild, non-significant increase in plasmacytoid cells in patients with severe COVID-19 compared with healthy donors (Figure 4a) . However, the number of IgM + CD19 + CD20 + B cells was decreased in patients with severe disease compared with mild, while IgG + , CD19 + , CD20 + cells remained comparable across all patients (Figure 4b ). Next, we visualized single cells from all patients in a common UMAP plot and assigned clusters based on surface marker expression (Figure 4c ). This approach identified two distinct groups based on the expression of surface IgM, with the total number of IgM + cells within clusters increased in severe COVID-19 patients (Figure 4d) . Conversely, healthy donors displayed B cells with high expression of surface CD19 + and CD20 + antigens (Figure 4d ). Closer inspection of CD19 and CD20 expression identified two distinct populations that differ in CD20 levels (Figure 4e) . This approach also revealed that the relative abundance of circulating CD19 and CD20 bright B cells was lower in COVID-19 patients compared with healthy individuals regardless of disease severity. To shed light on the functional relevance of these different B cell subsets, we quantified the expression of IgG and IgM in each population identified based on CD19 and CD20 co-staining ( Figure  4f ). Circulating CD19 low B cells (populations A and B) were enriched for IgG + cells in patients with mild and severe COVID-19 and IgM + cells in severe COVID-19 patients, while convalescent patients resembled healthy donors. No such difference was observed with CD19 + and CD20 bright B cells (population C) and CD19 + CD20 + and CD19 + CD20 -B cells (populations D and E). Overall, despite dwindling numbers of B cells overall, specific subsets of B cells, especially those with lower CD19 expression, have distinct immunoglobulin expression patterns in COVID-19 patients, with severe patients more frequently bearing IgM + B cells. We speculate that these findings may be related to the . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted September 9, 2020. . https://doi.org/10.1101/2020.09.08.20189092 doi: medRxiv preprint plasmacytoid differentiation and immunoglobulin switching programs, which may be dysfunctional d to SARS-CoV-2 infection.

The abundance of total B cells, plasma, and IgG+ and IgG+ cells between disease states. c) Uniform Manifold Approximation and Projection (UMAP) projection of all cells colored by surface receptor expression, cluster assignment, or disease severity. d) Immunophenotype of each cluster (top) and its composition by disease severity (bottom). e) Identification and quantification of five populations of B cells dependent on CD20 and CD19 expression. f) Comparison of the abundance of the populations identified in e) between disease states. Significance was assessed using Mann-Whitney U tests and corrected for multiple testing with Benjamini-Hochberg FDR. **, FDR-adjusted p-value <0.01; *, FDR-adjusted p-value 0.01-0.05.

Having characterized the main circulating compartments of the immune system, we next sought leverage the high-dimensionality of the dataset and relate individual immune fingerprints of differ patients using hierarchical clustering (Figure 5a) . Not only did individuals with similar clini conditions tend to cluster together, but clinical factors such as hospitalization or intubation a appeared to be linked to immunotypes. We hypothesized that this underlying data structure would useful for reconstructing the clinical course of COVID-19. Thus, we employed pseudotime inference reconstruct an underlying latent space from a healthy state to a severe disease state (Figure 5b/c) . is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted September 9, 2020. . https://doi.org/10.1101/2020.09.08.20189092 doi: medRxiv preprint lymphocytes, gain of myeloid cells (G-MDSCs in particular), and a terminally activated/exhausted cell phenotype (Figure 5d) . Besides discovering immune signatures associated with each degree severity, this analysis allows the relative positioning of each time point in relation to the continuo changes characterized by pseudotime (Figure 5e ). Variable changes associated with the pseu temporal axis could be classified in three clusters (Figure 5f ). The first was composed of populations, with an increase toward higher disease severity with representatives such as the fract of myeloid cells, PD-1 + CD4 + T-cells, and CD62Lcells among CD45RA + , CD8 + T cells (Figure  left) . The second corresponded to a virtually stable cluster with 52 populations such as IgM + B ce with only mild fluctuation in the intermediate stage (Figure 5g, center) . The third included a clus with a steady decrease by disease severity, encompassing the overall lymphoid population as well B cells and CD45RA + CD4 + T cells (Figure 5g, right) . This effectively establishes a tempo hierarchy of changes as disease progresses in which populations such as CD62L + , CCR7 + , CD45R CD8 + T cells have a steady decline and others such as B cells have a stronger decline toward severe end of the pseudo-temporal timeline. Additionally, the dynamic character of changes raises possibility of using flow cytometry to improve COVID-19 patient stratification based on real-ti immunological monitoring. Although our observations do not indicate causality, immunologi variations in the pseudo-temporal dimension may offer testable hypotheses on COVID-19 progress mechanisms. is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted September 9, 2020. . https://doi.org/10.1101/2020.09.08.20189092 doi: medRxiv preprint 11/27 space that reconstructs the hierarchy of disease progression. The x-axis represents disease progression in the pseudo-temporal space. c) Distribution of samples grouped by disease state along the pseudo-temporal axis derived in (b). d) Immune populations associated with the pseudo-temporal axis represented by either the absolute change in percentage in their extremes (x-axis) or strength of linear association (y-axis). e) Heatmap of immunotypes and immune populations sorted by their order or relative abundance in the pseudo-temporal axis, respectively. f) Clusters of immune populations based on their abundance along the pseudo-temporal axis. g) Examples of immune populations from each cluster in (f).

Since various clinical and demographic factors influence disease incidence and mortality 1,2,37 , we investigated the interaction between SARS-CoV-2 infection, the circulating immune system, and various demographic and clinical factors. Thus, we fit regularized linear models to the proportional flow cytometry data with co-variates such as sex, race, age, disease severity, presence of comorbidities, hospitalization, intubation, and death (Supplementary Figure 5a) . We also estimated the interaction of sex with clinical variables such as disease severity, hospitalization, intubation, and death. The resulting network of significant effects identified several clinical factors associated with specific immune cell populations, highlighting how age, sex, and disease severity jointly influence the circulating immune systems in patients with COVID-19 (Figure 6a) . As a baseline, we could recover known effects independent of disease, such as a higher CD4:CD8 ratio in females than males and an overall decrease of the lymphoid population with age (Supplementary Figure 5b) . Lastly, we found associations between sex and clinical variables such as a significantly higher fraction of CD62L + , CCR7 + , CD45RO + , CD4 + T cells in males that died compared with females (Figure 6b, left) and much lower total lymphocyte levels in females that died compared with males (Figure 6b, right) . Regarding the effect of tocilizumab on the immune system, we compared post-treatment samples from eight treated severe patients to seven severe untreated patients. While we observed the largest effect in certain subsets of CD4 + T cells, there was also an increased relative abundance of B cells and a decrease in T cells expressing the co-inhibitory receptor hepatitis A virus cellular receptor 2 (HAVCR2; TIM3) (Supplementary Figure 5c) . Moreover, the signature associated with severe versus mild patients was broadly counteracted by tocilizumab (Figure 6c) . Associations between sex and clinical variables were found, such as a lower fraction of CD62L+, CCR7+, CD45RA+, CD8+ T cells in females treated with tocilizumab compared with males, contrary to the opposing trend in untreated individuals (Figure 6d) , or the lower frequency of CD158a NK cells in female intubated patients (Figure 6e) .

. CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted September 9, 2020. . https://doi.org/10.1101/2020.09.08.20189092 doi: medRxiv preprint CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted September 9, 2020. . https://doi.org/10.1101/2020.09.08.20189092 doi: medRxiv preprint

Since there is a need to stratify patients to provide better, more effective, and less costly care, particularly in the earlier stages of disease, we hypothesized that the high dimensionality of the immunotypes would make it possible to train a classifier to predict disease severity early on. A random forest classifier was trained to distinguish patients with mild from severe disease using only the earliest available sample of each patient in a cross-validated manner (Figure 6f) . We observed good performance of the classifier (median area under receiving operator curve [ROC AUC], 0.81) compared with one with randomized severity labels (median ROC AUC, 0.49) (Figure 6g) , providing good balance between true positive and false positive rates. Since our dataset is composed of immune populations from seven flow cytometry panels, we tested whether a smaller number of variables could discern patients with mild and severe disease courses. With only eight variables, the classifier could distinguish patients with different disease severities, albeit with lower performance (ROC AUC, 0.73 vs. 0.49 with randomized labels) (Figure 6h) . Furthermore, we hypothesized that our classifier could be used for real-time immuno-monitoring of COVID-19 patients. Thus, we applied it to subsequent samples of patients with more than three samples collected over the disease course, while withholding those samples from the training set (Figure 6i ). Patient 26, who had an overall mild disease course, had all samples classified as mild; severe patients often showed dynamic severity probabilities over time, with at least one timepoint classified as severe disease. To exemplify how this prediction relates back to flow cytometry data, we illustrate the aggregated expression of the activation marker CD25 and CD45RA/RO in single T cells over time (Figure 6j) . Patients with lower predicted severity toward the end of their course (e.g., patient 23) tended to have less CD25 expression and increased CD45RA expression, while the opposite was also true (e.g., patient 16). Patients with predictions that were either more stable or dynamic over time (patients 26 and 24, respectively) showed dynamics of expression in accordance to their overall predicted pattern over time. This proofof-principle work demonstrates our ability to leverage high-content immune profiling to predict overall disease course and provides the basis for real-time immune-monitoring of COVID-19 patients.

Here, we describe the circulating immune landscape of COVID-19 patients compared with healthy individuals. Consistent with previous reports 11, 14, 15, 29 , we demonstrate that disease progression is dominated by the progressive loss of circulating lymphocytes and gain of myeloid cells. We also detected selective expansion of NK populations and MDSCs, suggesting that the innate compartment may contribute to the immunological disarray of COVID-19 patients. We then harnessed this multidimensional dataset to generate a machine-learning classifier that could predict disease severity using a defined flow-cytometric signature. Our work provides a proof-of-concept that an immune-monitoring algorithm could provide a rapid and personalized approach to manage COVID-19.

While previous studies have focused on lymphocyte populations 12, 14, 15, 38 , to our knowledge the role of innate immune cells is less understood 39 . Our study highlights the expansion of MDSCs, especially G-MDSCs, in severe COVID-19 patients. Unlike their natural counterparts, these elements have suppressive function 40 that impairs immune responses in cancer 34 and derails effective responses against bacterial and viral infections by the adaptive immune system 41, 42 . Given the overall depletion of the immune system's lymphoid branch during COVID-19, an interesting hypothesis is that G-MDSCs and other myeloid cells represent uncontrolled negative feedback. These elements ultimately contribute to the establishment of pan-immunosuppression, leading to dysregulated responses from the adaptive immune system. It will be essential to establish whether they are actively recruited to . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted September 9, 2020. . https://doi.org/10.1101/2020.09.08.20189092 doi: medRxiv preprint 14/27 infected lungs and whether they are causally involved in disease pathogenesis or represent a systemic compensatory response to inflammation. Since MDSCs are virtually absent in healthy individuals, questions arise regarding the mechanisms of their genesis and tissue recruitment and how they interact with lung tissue. At the same time, our novel observation that NK cells expressing the CD158i variant are over-represented in patients with severe disease raises the question of whether this variant and other KIRs implicate NK cells in disease progression (Figure 3) .

Within the adaptive immune system, the mechanisms leading to severe immune depletion, a landmark seen with disease progression and markedly apparent in autopsy samples, are unknown 43 . To this end, we observed increased CD25 + T cells in COVID-19 patients, indicating a higher state of activation (Figure 2) , but also an increase in CD95 + with disease progression. This phenotype was significantly marked in severe patients, consistent with a recent study 38 . FAS has a crucial role in mediating cell death via FAS ligand engagement, as in activation-induced cell death (AICD), or by shifting cells to a more apoptotic-prone phenotype. While FAS is a natural regulatory checkpoint of T cells, it plays a role in autoimmunity 44 and cancer 45 , and AICD is involved in loss of CD4 + and CD8 + cells in HIV patients 46 . However, severely exhausted T cells can undergo apoptosis, and virus-specific T cell decline can favor viral escape [47] [48] [49] . Indeed, similar to previous reports 17, 18, 50 , T cells displayed an overall exhausted phenotype, with overexpression of VISTA, TIM3, LAG3, TIGIT, and PD-1 coinhibitory receptors in COVID-19 patient T cell populations. This likely results in inability of the adaptive immune system to keep viral proliferation in check. In the B cell compartment, we observed lower expression of CD19 in COVID-19 patients and higher expression of membrane-bound IgM and IgG in both mild and severe patients. These data suggest that under viral exposure, B cells undergo plasmacytoid maturation and immunoglobulin switching. Remarkably, several patients displayed higher IgM than IgG CD19 + CD20 +/cells, suggesting abnormal and delayed maturation of plasma cells (Figure 4) . While the implications remain speculative, they do warrant further investigation given the central role of B cells in the development of immunity by COVID-19 patients.

Taking advantage of our dataset's high-dimensional characteristics and pseudo-temporal modeling, we constructed a COVID-19 disease course landscape. This strategy reveals a continuum of disease progression between healthy state, mild disease, and severe disease. Remarkably, convalescent patients displayed immune phenotypes similar to healthy donors, suggesting a possible return to a largely healthy state, as previously suggested based on the exhaustion phenotype in adaptive responses 18 . Conversely, we could speculate that the immune landscape of mild/convalescent patients never achieved the level of disarray observed in severe patients. While there were marked differences between patients with prevalent mild or severe disease, their recognition remains a unique challenge. One interesting open question is whether the changes associated with mild vs. severe disease protect against disease progression or, conversely, which immune populations related to severe illness play a role in the progression to severe disease. While proof-of-principle, our classifier of severe disease shows robustness and overall value in predicting disease progression based on immune profiling and near-real-time disease monitoring. Thus, it may be valuable to inform clinical action like that proposed in chronic diseases with other high-dimensional assays [50] [51] [52] . Moreover, we demonstrated that a classifier with a limited number of markers retains good performance. If confirmed in large cohorts, it could provide a useful approach to stratify patients and predict clinical evolution using a rapid and economical assay.

Although our study confirmed some findings and provides new data on the innate immune landscape of COVID-19 patients, we recognize limitations that could be overcome with larger sample sizes and matched control populations. Indeed, large population studies may discover new immune populations . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted September 9, 2020. . https://doi.org/10.1101/2020.09.08.20189092 doi: medRxiv preprint 15/27 and find a relationship with disease progression or clinical factors. In this current pandemic, profiling a larger sample of the population and investigating multiple time points may help identify viral adaptation to the host (particularly when coupled to analysis of viral sequence) in patients with different outcomes. These programs may be achieved if an effective institutional organization, multicentric networking, and substantial financial support are available.

The targeted nature of flow cytometry interrogates limited sets of immune populations and implies that only certain molecules can be effectively profiled. In our study we used mainly proportional data when comparing the abundance of immune populations between patient groups. While this may not necessarily imply absolute changes in cell numbers, we observed good overall agreement between changes in proportions and absolute counts when comparing severe and mild disease status (Supplementary Figure 6 and Supplementary Table 7) . This highlights the importance of studies employing orthogonal modalities such as cytokine profiling 51 , single-cell RNA sequencing [52] [53] [54] [55] [56] , and their integration 57, 58 . Nevertheless, even without orthogonal studies, our machine learning approach for predicting disease severity demonstrates predictive potential, although it should be tested in a validation cohort before use in a clinical setting.

Lastly, we wish to note the work of others and their complementary findings. For example, Laing et al. 32 used peripheral blood flow cytometry and circulating cytokine measurements to demonstrate apparent immune dysregulation in COVID-19 patients. They highlighted additional interesting, complementary features, including increased IL-6, IL-10, and IP-10 and depletions of basophils, plasmacytoid dendritic cells, T H 1 cells, and T H 17 cells. Incorporating their most differentiating markers with ours could yield a more complete yet targeted panel of markers with more predictive power to determine which patients will rapidly progress to a severe disease state. This incorporation of additional differentiating markers should be pursued in future studies.

Collectively, our study highlights a profound imbalance in the COVID-19 immune landscape, characterized by G-MDSC expansion and T cell exhaustion that may open avenues for clinical translation. Further, our approach provides a powerful tool to predict clinical outcomes and tailor more effective and proactive therapies to COVID-19 patients.

The study was approved by the Institutional Review Board of Weill Cornell Medicine. Participants were recruited from patients hospitalized at New York Presbyterian Hospital from April to July 2020. Some participants in a COVID-19 convalescent plasma donor screening program with prior confirmed diagnosis (by RT-PCR or serology) were given the option to contribute a sample for this research. ARDS was categorized in accordance with the Berlin definition reflecting each subject's worst oxygenation level and with physicians adjudicating chest radiographs 59 . Informed consent was obtained from all participants.

For flow cytometric analysis of circulating leukocytes, peripheral blood was collected in Na-heparin. Except for the MDSC panel, in which PBMCs were prepared by density gradient centrifugation, erythrocytes were lysed with BD Pharm Lyse. Peripheral blood was washed in Dulbecco's PBS . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted September 9, 2020. . https://doi.org/10.1101/2020.09.08.20189092 doi: medRxiv preprint 16/27 (DPBS), lysed in 1× BD Pharm Lyse, and washed again in DPBS. PBMC cell suspensions were prepared with Ficoll-Paque following the manufacturer's protocol. Cells were stored briefly in storage medium (10% heat-inactivated fetal bovine serum/1% L-glutamine/1% pen-strep) before staining with antibody cocktails.

For each panel, one million cells were stained with specific cocktails of fluorochrome-conjugated antibodies (Supplementary Table 3-4) . Cells were washed with DPBS then stained with dead cell dye (BD Fixable Viability Stain 700) before washing with wash buffer (0.5% BSA/DPBS/NaN 3 ). Cells were then treated with 50 μ L of Fc-blocking solution (2% normal rabbit serum/10% BD Fc Block/DPBS) before application of a 100-μL antibody cocktail diluted in wash buffer. Samples were stained within 6 hours of sample collection and analyzed on a BD Biosciences FACSCanto flow cytometer within 2 hours of staining. The stopping gate was set to acquire 500,000 viable, nucleated single cells. The analysis of the T cell memory and checkpoint panels started with identifying T cells (CD3 vs. SSC), then the CD4 helper and CD8 cytotoxic subsets. To analyze the T cell checkpoint panel, individual exhaustion markers were gated on histogram plots. The T cell memory panel was further subdivided into CD45RA + /CD45ROand CD45RA -/CD45RO + subsets. Under these gates, two quadrant gates were placed on CD62L vs. CCR7 and CD62L vs. FAS.

Gating for the B-cell panel began with CD45 vs. CD19 then FSC-A and SSC-A to identify cells of the B lineage. The CD20 + and CD20subsets were gated (CD19 vs. CD20) and IgG and IgM were quantified within the CD20 + subset (IgM vs. IgG). NK cells were identified by sequential gates on CD56 vs. CD3, FSC-A vs. SSC-A, and CD56 vs. CD16. CD56 + , CD16 bright, mature NK cells were then interrogated for their reactivity with individual anti-KIR (CD158) and anti-NKG2A (CD159a) antibodies with gates on histogram plots. KIR-negative NK were identified by sequential gating on CD158a vs. CD158b double-negative, then CD158i and CD158e double-negative subsets.

. CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted September 9, 2020. . https://doi.org/10.1101/2020.09.08.20189092 doi: medRxiv preprint 17/27

Non-parametric Mann-Whitney U tests were employed to assess the significance of pairwise changes in the proportions of immune populations between severity groups using the Pingouin package, version 0.3.7. Multiple test correction was performed with the Benjamini-Hochberg FDR method.

To select cells from the events, single cells were gated using forward-side scatter height and area, CD45-positivity, viability dye-negativity, and the major marker of each panel (e.g.. CD3 for T cell memory panel). Compensation was applied using FlowKit 60 version 0.5.0, and an inverse hyperbolic transformation (AsinhTransform) was applied with parameters t = 10000, m = 4.5, a = 0. To construct a shared latent representation for all cells, dimensionality reduction was performed with Principal Component Analysis (PCA), a neighbor graph was computed using 15 neighbors per cell, UMAP 61 with default parameters, and Leiden clustering, all using the Scanpy package 62 version 1.5.1. For each discovered single-cell cluster, a proportion of cells was calculated in relation to a specific clinical factor after normalization by the frequency of the same factor in the cohort.

To learn a latent manifold of the data, the non-linear method Laplacian Eingenmaps 63 was used as implemented in the "SpectralEmbedding" method of the scikit-learn framework 64 (version 0.23.0) with default parameters. A z-scored matrix of proportional data was input for all immune cell populations (variables) and patient samples (observations). To rank the features by their association with the learned space, Pearson's correlation was calculated between the first component and each variable, in addition to the fold and absolute change in the variable between the top and bottom 10% of the samples in each extreme of the embedding. The same procedure applied to a Uniform Manifold Approximation and Projection (UMAP) latent representation of the same data yielded similar results, with the exception that the spread of samples according to disease progression was parallel to multiple learned axes rather than single.

To rank variables by the amount of change in both real time since the reported start of symptoms for a single patient or over the learned latent space across all patients, GPy package 65 was used to fit Gaussian Process regression models on the learned pseudotime axis (independent variable) and the abundance of each immune cell population (dependent variables). A variable radial basis function (RBF) kernel and a constant kernel (both with an added noise kernel) were fitted and the log-likelihood and standard deviation of the posterior probability of the two were compared as described previously 66 . To cluster the abundance of immune populations based on their dynamics over the pseudotime axis, the same kernels were used to fit a Mixture of Hierarchical Gaussian Processes (MOHGP) as implemented in the GPClust package 67,68 using 8 as an initial guess of number of clusters.

Due to the proportional nature of the dataset, generalized linear models were fit using a Gammadistributed noise model with a log-link function. Ridge regularization was employed to ensure robust coefficients given the low abundance of some populations, and the model was fit with ordinary least squares optimization using the statsmodels package 69 version 0.11.1. Categorical variables were onehot encoded and numeric ones such as age or days since symptoms started were kept as years or . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted September 9, 2020. . https://doi.org/10.1101/2020.09.08.20189092 doi: medRxiv preprint 18/27 days, respectively; the date of acquisition was transformed into days and scaled to the unit interval. Since values for clinical categorical variables and comorbidities were only available to COVID-19 patients, various models were employed that aimed to explore different aspects of immune system change during COVID-19:

1. Comparison of healthy donors to COVID-19 patients: sex + race + age + batch + COVID19 2. Effect of clinical/demographic factors on COVID-19 patients: sex + race + batch + COVID19 + severity group + hospitalization + intubation + death + diabetes + obesity + hypertension + age in years + days since symptoms start 3. Effect of tocilizumab treatment on severe patients only: sex + age + batch + tocilizumab

To generate a graph of interactions between factors and immune populations, significant coefficients (FDR-adjusted p-value <0.05) were used as undirected edges between factors and immune populations. For edges between factors, the Pearson correlation between factors across immune populations was used. Exclusively for visualization, coefficients for the continuous variables "age" and "time since symptoms started" were multiplied by half of the median of the values of that variable (33.0 and 10.8, respectively) to make the range of coefficients comparable with the categorical variables. Visualizations were produced using Gephi version 0.9.2 with the Force Atlas2 layout with parameters "LinLog mode", "scaling factor" 8.0, and "gravity" 11.0.

A Random Forest Classifier was trained as implemented in scikit-learn framework 64 (version 0.23.0) to distinguish between cases with "mild" and "severe" disease using 10-fold cross validation. The crossvalidation loop was repeated 100 times and models were fit with real or randomized labels. Test set performance was assessed with the ROC-AUC. To investigate the performance of the classifier, feature importance was averaged across cross validation folds and iterations and the log fold importance of the real models over the randomized labels was calculated. A sign was added to the feature importance depending on the sign of the Pearson correlation of each variable with each class. Only the earliest temporal sample of each patient was used to ensure lack of data leakage (avoid training/testing on samples from the same patient without stratified cross validation) and to maximize the utility of the model. The same cross-validation scheme was used to develop a classifier using a subset of features but including feature selection using mutual information inside the cross-validation loop. To predict severity longitudinally for single patients, a model was trained on the initial samples from all other patients and tested on the samples of the patient in question.

Quantification of immune cell populations is available as a Supplementary Table file. Hierarchical data format files with single cell data (h5ad) are available as indicated in the repository with source code for the study (https://github.com/ElementoLab/covid-flowcyto) and raw FCS files will be made available upon post-peer review publication.

The source code for the analysis is available on GitHub at the following URL: https://github.com/ElementoLab/covid-flowcyto . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

The copyright holder for this preprint this version posted September 9, 2020. is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

The copyright holder for this preprint this version posted September 9, 2020. . https://doi.org/10.1101/2020.09.08.20189092 doi: medRxiv preprint 21/27 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

The copyright holder for this preprint this version posted September 9, 2020. . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

The copyright holder for this preprint this version posted September 9, 2020. CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

The copyright holder for this preprint this version posted September 9, 2020. . https://doi.org/10.1101/2020.09.08.20189092 doi: medRxiv preprint 24/27 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

The copyright holder for this preprint this version posted September 9, 2020. 

Supplementary 

Clinical Characteristics of Coronavirus Disease 2019 in China

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

Breadth of concomitant immune responses prior to patient recovery: a case report of non-severe COVID-19

Clinical features of patients infected with 2019 novel coronavirus in Wuhan

Immunopathological characteristics of coronavirus disease 2019 cases in Guangzhou

Effectiveness and safety of glucocorticoids to treat COVID-19: a rapid review and metaanalysis

Tocilizumab in patients with severe COVID-19: a retrospective cohort study

Interleukin-6 inhibitors in the treatment of rheumatoid arthritis

Dexamethasone in Hospitalized Patients with Covid-19 -Preliminary Report

Tocilizumab for treatment of patients with severe COVID-19: A retrospective cohort study

Severe COVID-19 is marked by a dysregulated myeloid cell compartment

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

Acute SARS-CoV-2 infection impairs dendritic cell and T cell responses

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

Comprehensive mapping of immune perturbations associated with severe COVID-19

Coregulation of CD8+ T cell exhaustion by multiple inhibitory receptors during chronic viral infection

Marked T cell activation, senescence, exhaustion and skewing towards TH17 in patients with COVID-19 pneumonia

Functional exhaustion of antiviral lymphocytes in COVID-19 patients

Phenotype and kinetics of SARS-CoV-2-specific T cells in COVID-19 patients with acute respiratory distress syndrome

Targets of T Cell Responses to SARS-CoV-2 Coronavirus in Humans with COVID-19 Disease and Unexposed Individuals

Analysis of a SARS-CoV-2-Infected Individual Reveals Development of Potent Neutralizing Antibodies with Limited Somatic Mutation

Prevalence of comorbidities and its effects in patients infected with SARS-CoV-2: a systematic review and meta-analysis

Cardiovascular Implications of Fatal Outcomes of Patients With Coronavirus Disease 2019 (COVID-19)

SARS-CoV-2 Transmission in Patients With Cancer at a Tertiary Care Hospital in Wuhan, China

Sex Differences in Case Fatality Rate of COVID-19: Insights From a Multinational Registry

Clinical course and outcomes of critically ill patients with SARS-CoV-2 pneumonia in Wuhan, China: a single-centered, retrospective, observational study

Disease severity-specific neutrophil signatures in blood transcriptomes stratify COVID-19 patients. Infectious Diseases (except HIV/AIDS

Follicular Helper T Cells

High-dimensional immune profiling by mass cytometry revealed immunosuppression and dysfunction of immunity in COVID-19 patients

Increasing LAG-3 expression suppresses T-cell function in chronic hepatitis B: A balance between immunity strength and liver injury extent

Adoptive transfer of effector CD8+ T cells derived from central memory cells establishes persistent T cell memory in primates

A dynamic COVID-19 immune signature includes associations with poor prognosis

IL-6 and IL-8 Are Linked With Myeloid-Derived Suppressor Cell Accumulation and Correlate With Poor Clinical Outcomes in Melanoma Patients

The Nature of Myeloid-Derived Suppressor Cells in the Tumor Microenvironment

Myeloid-derived suppressor cells as regulators of the immune system

The Role of Myeloid-Derived Suppressor Cells in Viral Infection

COVID-19 and gender-specific difference: Analysis of public surveillance data in Hong Kong and Shenzhen

Increased CD95 (Fas) and PD-1 expression in peripheral blood T lymphocytes in COVID-19 patients

Expansion of myeloid-derived suppressor cells in patients with severe coronavirus disease (COVID-19)

Neutrophils and PMN-MDSCs: their biological role and interaction with stromal cells

Early Activation of Myeloid-Derived Suppressor Cells Participate in Sepsis-Induced Immune Suppression via PD-L1/PD-1 Axis

Identification of an Immunosuppressive Cell Population during Classical Swine Fever Virus Infection and Its Role in Viral Persistence in the Host

Histopathology and ultrastructural findings of fatal COVID-19 infections in Washington State: a case series

Why do defects in the Fas-Fas ligand system cause autoimmunity?

CD95 promotes tumour growth

Activation-induced CD4+ T cell death in HIV-positive individuals correlates with Fas susceptibility, CD4+ T cell count, and HIV plasma viral copy number

Virus persistence in acutely infected immunocompetent mice by exhaustion of antiviral cytotoxic effector T cells

Memory CD8 T-cell differentiation during viral infection

Effector and memory CTL differentiation

Elevated exhaustion levels and reduced functional diversity of T cells in peripheral blood may predict severe progression in COVID-19 patients

Longitudinal analyses reveal immunological misfiring in severe COVID-19

Single-Cell Omics Reveals Dyssynchrony of the Innate and Adaptive Immune System in Progressive COVID-19

A single-cell atlas of the peripheral immune response in patients with severe COVID-19

Tocilizumab treatment in severe COVID-19 patients attenuates the inflammatory storm incited by monocyte centric immune interactions revealed by single-cell analysis

Immunophenotyping of COVID-19 and influenza highlights the role of type I interferons in development of severe COVID-19

Single-cell sequencing of peripheral blood mononuclear cells reveals distinct immune response landscapes of COVID-19 and influenza patients

Multiomic Immunophenotyping of COVID-19 Patients Reveals Early Infection Trajectories

Suppressive myeloid cells are a hallmark of severe COVID-19

ARDS Definition Task Force et al. Acute respiratory distress syndrome: the Berlin Definition

Intuitive Python framework for flow cytometry analysis and visualization

Dimensionality reduction for visualizing single-cell data using UMAP

SCANPY: large-scale single-cell gene expression data analysis

Laplacian Eigenmaps for Dimensionality Reduction and Data Representation

Scikit-learn: Machine Learning in Python

A Gaussian process framework in python

Chromatin mapping and single-cell immune profiling define the temporal dynamics of ibrutinib response in CLL

Fast variational inference in the conjugate exponential family

Hierarchical Bayesian modelling of gene expression time series across irregularly sampled replicates and clusters

Econometric and Statistical Modeling with Python

This project was supported by a Translational Pathology Research COVID-19 grant to G.I., by a ISTM grant to M.S. and by the National Center for Advancing Translational Science of the National Institute of Health Under Award Number UL1TR002384. C.K.V. is supported by NIH K08 AI132739; A.M. is supported by grant KL2TR002385 of the Clinical and Translational Science Center at Weill Cornell Medical College. K.S. is supported by NIH K08 AI139360; C.D.B. is supported by NIH T32 AI07613-19 (PI: Gulick) and by the Kellen Foundation. L.G. is supported by from the Leukemia and Lymphoma Society (LLS), a startup grant from the Dept. of Radiation Oncology at Weill Cornell Medicine (New York, US), a Rapid Response Grant from the Functional Genomics Initiative (New York, US). We thank Andrew Marderstein, Fayzan Chaudhry, and Liron Yoffe for helpful discussions on the machine learning classifier for disease severity. We are grateful for the support of members of the Immunopathology laboratory at New York Presbyterian Hospital, Weill Cornell Medicine, whose dedication and contribution have been instrumental for the execution of this project. We are grateful to the patients and their family who agreed to be part of the study and all the medical staff who cared for them.

The authors declare no competing financial interests.