Systems-Level Immunomonitoring from Acute to Recovery Phase of Severe COVID-19

Summary Severe disease of SARS-CoV-2 is characterized by vigorous inflammatory responses in the lung, often with a sudden onset after 5–7 days of stable disease. Efforts to modulate this hyperinflammation and the associated acute respiratory distress syndrome rely on the unraveling of the immune cell interactions and cytokines that drive such responses. Given that every patient is captured at different stages of infection, longitudinal monitoring of the immune response is critical and systems-level analyses are required to capture cellular interactions. Here, we report on a systems-level blood immunomonitoring study of 37 adult patients diagnosed with COVID-19 and followed with up to 14 blood samples from acute to recovery phases of the disease. We describe an IFNγ-eosinophil axis activated before lung hyperinflammation and changes in cell-cell co-regulation during different stages of the disease. We also map an immune trajectory during recovery that is shared among patients with severe COVID-19.


In Brief
Immune dysregulation plays a pivotal role during severe COVID-19. Rodriguez et al. present a systems-level, longitudinal study of 37 COVID-19 patients monitored from acute to recovery phases of disease. This study reveals cell-protein dependencies and co-regulated features, as well as a shared immunological trajectory during recovery.

INTRODUCTION
Since its emergence in December 2019, the severe acute respiratory syndrome-coronavirus 2 (SARS-CoV-2) causing coronavirus disease 2019 (COVID-19) has infected millions of individuals and caused hundreds of thousands of deaths worldwide. The betacoronavirus has a high degree of sequence homology with previous SARS-CoV and Middle East respiratory syndrome (MERS) coronaviruses and binds to the angiotensin-converting enzyme 2 (ACE2) receptor to enter cells in the respiratory and intestinal epithelium. 1 Cells recognize the presence of the virus through pathogenrecognition receptors (PRRs) and elicit antiviral response programs. 2 The two main components of such antiviral programs involve the production of type I and III interferons (IFNs) that induce downstream transcription of hundreds of IFN-stimulated genes (ISGs) that interfere with viral replication in the cell. 3 The second element of the antiviral response program is the secretion of chemokines that recruit specialized cells of the immune system to clear the virus. SARS-CoV-2, like other viruses, has evolved countermeasures to these defenses, and, in particular, the virus efficiently interferes with IFN signaling and the induction of ISGs in SARS-CoV-2-infected cells. 4,5 In contrast, pro-inflammatory cytokine and chemokine responses are induced normally, and this imbalance between antiviral and pro-inflammatory responses is a key feature of COVID-19. 6 Another observation during the COVID-19 pandemic is the different disease courses among different individuals infected by the SARS-CoV-2 virus. Most individuals present with very mild disease, often asymptomatic, and a few develop a lifethreatening disease requiring intensive care. The strongest determinant of disease severity is age, with children presenting almost exclusively with mild disease, 7 while the elderly, those older than 70 years of age, are much more likely to develop severe COVID-19. Males and females are infected at similar rates, but males are much more likely to develop severe disease requiring intensive care. 8 Obesity, smoking, and hypertension are other risk factors for severe COVID-19. 9 However, COVID-19 contrasts with other respiratory viral infections in that pregnant women do not seem to be more likely to develop severe disease, and this is also true for patients with various forms of immunodeficiency. One likely reason for these observations is that severe disease is associated with exuberant immune responses and a skewed immune regulation against proinflammatory responses in pregnancy and T cell deficiencies in transplant patients make such hyperinflammatory responses less likely. To treat hyperinflammation in severe COVID-19, we need to better understand what cells are involved, their interactions, and protein mediators used to orchestrate their responses. To this end, we performed systems-level analyses of the immune system in blood from 37 patients, from acute to recovery phases of COVID-19, with up to 14 blood samples collected from a given patient. These analyses reveal a sequence of responses involving many immune cell populations at different stages of the disease. A transient response involving IFN-g upregulating CD62L on eosinophils before lung hyperinflammation are exemplified when you look at coregulated cell populations, and immune correlates of productive antibody responses to SARS-CoV-2, as well as an integrated immune trajectory shared across patients recovering from severe COVID-19.

Longitudinal Profiling of Patients with COVID-19
Given the enormous diversity among immune systems in humans, longitudinal monitoring of patients is required to appreciate the immunological changes occurring during the disease process. Also, systems-level analyses methods such as mass cytometry 10 enable all immune cell populations to be distinguished and analyzed in a given blood sample, allowing for coordinated changes across cell populations to be revealed. We have combined these cellular measurements with analyses of 180 unique plasma proteins using Olink analyses 11 ( Figure 1A). To understand systems-level immune responses during moderate to severe COVID-19, we monitored longitudinal samples from 17 patients, some treated in the intensive care unit (ICU) and some treated in regular hospital wards with oxygen support but no mechanical ventilation, in addition to 20 recovered patients comprised of 18 mild COVID-19 recovered patients and 2 hospitalized COVID-19 recovered patients ( Figure 1B). Patients did not receive immunomodulatory therapies in this cohort apart from one before ICU admission, and the immunological changes thus reflect the natural course of severe COVID-19 infection. All of the patients in this cohort survived their infection.

The Characteristics of Acute and Recovery Phases of COVID-19
Clinical measurements were taken from acute and recovered patients, including body temperature, white blood cell (WBC) counts, and lymphocyte counts. Milder cases of COVID-19 showed lower body temperatures as well as faster normalization of body temperatures compared to severe cases, which fluctuated more over time (Figure 2A). The WBC counts changed during the stages of infection. High WBC counts prototypically occur during acute inflammation and immune responses. In severe patients we observed fluctuating levels of WBC over time ( Figure 2B). More important, there were no signs of secondary bacterial infection in the patients in this cohort. Lymphopenia is one of the common features of COVID-19 and the degree of lymphopenia predicts disease severity. 9 Lymphocyte counts were measured and seemed to recover faster in milder as compared to severe cases, although this trend was not seen uniformly (Figure 2C). This is in line with other previous reports. 12 Plasma protein levels were measured and compared among acute and recovered phases and reflect the immune dynamics of severe COVID-19 ( Figures 2D-2G). Pro-inflammatory cytokines such as interleukin-6 (IL-6) and IFN-g predict disease severity. A decreasing trend was observed in IFN-g and IL-6 from early admission to the hospital through recovery during the weeks of the study (Figures 2D and 2F, respectively). Similarly, DDX58, the innate immune response receptor, also called RIG-I, and the monocyte chemoattractant protein MCP-3, were elevated during acute disease and decrease during recovery ( Figures 2E  and 2G, respectively).

The Immune Cell Changes from Acute to Recovery Phases of COVID-19
A defining feature of the acute immune response during COVID-19 is dramatic changes in immune cell composition. These changes can be informative about likely driving factors and triggers and can help us understand the disease process better. To understand severe COVID-19 better, we plotted relative proportions of 57 immune cell populations over time in the 37 patients (Figures 3, S1, and S2). These cell populations were defined using a recently developed tool for automated cell classification based on known immune cell phenotypes. 13 We confirmed the overrepresentation of neutrophils over lymphocytes during acute infection that is slowly reversed during the recovery phase ( Figure 3). This is in line with reports suggesting that the neutrophil:lymphocyte ratio (NLR) and degree of lymphopenia are predictive of disease severity in COVID-19. 12 The plasmablast response is clear and occurs during the first week after hospital admission in these patients ( Figure 3). The recovery of T cells after the initial lymphopenia occurs over the following 2-3 weeks and is dominated by CD127-expressing effector and central memory CD4 + T cells, as well as CD57-expressing and differentiated memory CD8 + T cells ( Figure 3). Also, all dendritic cell (DC) subsets increased from acute to recovery phases-CD1c + DCs, CD11c + DCs, and plasmacytoid DCs (pDC) (Figure 3). Despite a clear reduction in the relative abundance of neutrophils over time, the other granulocyte subsets, basophils and eosinophils, increased clearly from acute to recovery phases (Figure 3), and both of these were among the most dynamic cell populations during severe disease, which is suggestive of important contributions to antiviral defense and immunopathology.

Eosinophil Activation and Homing during Acute COVID-19
Given the changes in eosinophil abundance described above, we decided to study eosinophils more carefully. There are reports of strong granulocyte-macrophage colony-stimulating factor (GM-CSF) responses in the lungs of COVID-19 patients, 14 and GM-CSF is known to stimulate eosinophils, particularly in interstitial pneumonia and allergic inflammation. 15 Taking advantage of the detailed longitudinal sample series, we used partition-based graph abstraction (PAGA), 16 to reconstruct single-cell phenotypic changes in blood eosinophils during acute COVID-19 ( Figure 4). Leiden clustering found 12  Figure 4C). This phenotype of eosinophils is reminiscent of lung-resident eosinophils, rather than induced inflammatory eosinophils in circulation, and such lung-homing cells have previously been reported to be important homeostatic regulators of inflammatory responses in the lung 18 (Figure 4D). It is tempting to think that this transient expansion of CD62L + eosinophils just before the development of severe lung hyperinflammation at $1 week after admission is related to this immunopathology of the lungs in COVID-19 patients.
To this end, further investigation into this eosinophil-IFN-g axis is required and may suggest novel therapies targeting this response to mitigate acute respiratory distress syndrome (ARDS) and lung inflammation.

Adaptive Immune Cell Dynamics during Severe COVID-19
Adaptive responses to SARS-CoV-2 are seen in most individuals, with one study reporting CD4 + T cell responses and CD8 + T cell responses in nearly all patients. 19 Similarly, the majority of symptomatic patients seroconvert within a few days and most developed high-titer antibody responses 20 ; however, one study has reported that a significant proportion of patients with COVID-19 do not develop neutralizing antibody responses. 21 To investigate the dynamics of adaptive immune cell responses in our cohort, we used the same PAGA approach as described above. We find a clear plasmablast response early after admission ( Figure 5A). The CD4 + T cell response was initially dominated by effector and central memory responses, followed by an increase in regulatory T cells (Tregs) $4 days after admission ( Figure 5B). The CD4 + T cells were split into two effector cell populations based on CD4 expression level, possibly reflecting an activation-induced downregulation in a subset of CD4 + T cells ( Figure 5B). 22 The CD8 + T cell responses are dominated by activated cells expressing high CD38 and also a subset of effector cells upregulating the CD147 receptor from $1week onward ( Figure 5C). Gamma-delta T cell receptor (TCR) T cells (gdT cells) and CD8 + T cells progressively upregulated the marker of terminal maturation CD57 from $1 week onward ( Figures 5C and 5D). These results are largely in agreement with other recent reports 23 and highlight the strong innate and adaptive immune activation during acute COVID-19.

Cell-Cell Regulation Varies over Time during Severe COVID-19
Immune responses are always concerted efforts made by multiple, specialized cell populations communicating via direct interactions and secreted cytokines and other mediators. By studying such cell-cell relationships, a better understanding of the systems-level response can be obtained. We generated cell-cell correlation matrices using longitudinal cell population frequencies and binned the samples into four phases, from acute disease to recovery phase ( Figure 6A). We find that the first phase (days 0-4) is dominated by an inverse correlation between neutrophils and a number of myeloid and lymphoid cell types, as reflected in the elevated NLR, associated with severe disease 12 ( Figure 6A). The following phase (days 6-8) is characterized by a strong coordinated plasmablast, B cell, and abT cell module, and this is inversely correlated with a strong Treg and CD11c + DC module ( Figure 6A). From day 9 onward, a change is apparent, with a shift toward a co-regulated module involving eosinophils, pDCs, CD11c + DCs, and CD8 + T cells. This module is largely maintained in recovered patients, possibly reflecting a more normal cell-cell relationship ( Figure 6A). A prototype of a coordinated immune response to viruses is the appearance of virus-specific immunoglobulin G (IgG) antibodies, because such responses elicited by B cells require help from CD4 + T cells. Here, we investigated the seroconversion in this cohort and found a strong induction of IgG antibodies to the SARS-CoV-2 Spike protein (receptor-binding domain [RBD]) in the majority of patients ( Figure 6B). This is in line with similar analyses in other COVID-19 patients. 20,24 We were unable to test the neutralizing capacity of these antibodies at the time of the study, but another recent report has shown that a significant proportion of patients mount antibodies that lack such neutralizing capacity. 21 To understand the immunological correlates of IgG responses to SARS-CoV-2, we devised a mixed-effects model, using both plasma protein levels and cell frequencies as predictors, taking days after admission into account as a fixed effect. 25 It is important to note that this is not a simple correlation analysis since days from admission is taken into account as a fixed effect in the analysis. We found several features significantly associated with IgG responses, and in particular, strong proinflammatory cytokines IFN-g and IL-6 and chemokines CXCL10 and MCP-2 (CCL8) are negatively associated with anti-CoV-2 IgG responses ( Figure 6C). In contrast, the neutrophil-recruiting chemokine CXCL6 is positively associated with anti-CoV-2 IgG responses as was the fraction of circulating basophils ( Figure 6C). It is known that basophils are able to bind antigens on their surface and potentiate humoral immune responses 26 ; since basophils are depleted during acute and severe COVID-19 (Figure 3), our data collectively suggest that the degree of basophil depletion may influence the efficacy of IgG responses to SARS-CoV-2. It is believed that basophil-mediated enhancement of B cell responses occurs through the production of either IL-4 or IL-6, but levels of the latter were found to be inversely associated with antibody responses (Figure 6C), so it is more likely that another mechanism is responsible for the basophil enhancement of IgG responses in COVID-19. Collectively, these results indicate a coordinate adaptive immune response to SARS-CoV-2, enhanced by basophils and possibly suppressed by hyperinflammatory cytokine responses with high IL-6 levels during acute COVID-19.

A Shared, Integrated Trajectory of Recovery across Patients
Since none of the patients in this cohort were treated with immunomodulatory agents, apart from one patient who received Oxiklorin treatment before ICU admission, and have recovered with supportive care alone, we reasoned that a deeper investigation into the immunological changes during recovery from severe COVID-19 would be informative about the underlying immune processes involved. Given the strong interactions among immune cells and proteins in the immune system, we applied an integrative analysis method to search for a multiomic trajectory of immune recovery. We used multiomics factor analysis (MOFA). 27 This method allowed us to search for the latent factors (LFs) that best explain the variance across data types and use these to discern any possible relationship with the process of recovery from the disease. We found 10 LFs that explained the variance in the combined dataset ( Figure 7A), and of those, LF2 was associated with the transition from acute to recovery phases of the disease ( Figure 7B). There were no clear differences among ICU or non-ICU ward patients, and the levels of LF2 were highest in the samples taken from recovered patients ( Figure 7B). To understand the biology of immune recovery, we assessed the underlying features contributing to LF2. The plasma proteins changing the most decreased during recovery. The most prominent were IL-6, monocyte-chemotactic protein 3 MCP3, Keratin19 (KRT19), CXC chemokine ligand 10 (CXCL10), amphiregulin (AREG), and IFN-g ( Figure 7C). Conversely, the cells that changed the most during recovery were classical and non-classical monocytes, CD56 dim natural killer (NK) cells, eosinophils, and gdT cells, all increasing in relative proportions during recovery ( Figure 7D). This shared, integrated trajectory reveals the markers most indicative of recovery in patients with severe COVID-19, and if verified in independent sets of patients, these features could be valuable biomarkers to monitor during disease progression to detect a switch from acute to recovery phases in severe COVID-19.
In this article, we have performed an in-depth, longitudinal analysis of the immune system in patients with severe COVID-19 during acute disease and up until spontaneous recovery. The natural course of this process is mapped and found to be similar among patients. We find changes in cell populations, such as CD62L-expressing eosinophils, triggered by IFN-g and likely contributing to hyperinflammation and ARDS during acute disease. We show that basophils are depleted during acute dis-ease but recuperate during recovery, and the levels of basophils are significantly correlated with the titers of IgG antibodies to SARS-CoV-2 produced by B cells. In contrast, high levels of IL-6 and IFN-g are negatively associated with humoral responses. Finally, we uncover an immunological trajectory of disease recovery shared among patients. These results can be useful for the development of better immunomodulatory strategies to mitigate hyperinflammatory responses, optimize antiviral IgG responses, and monitor disease progression and recovery in patients with severe COVID-19.

DISCUSSION
A number of researchers are studying the immune response to SARS-CoV-2, and we are learning about viral evasion of IFN-I/ III signals and prevention of the normal induction ISGs and the antiviral state. 4,5 At the same time, the proinflammatory response is strong. The secretion of chemokines and proinflammatory cytokines leads to the influx of neutrophils and myeloid cells into the lung, with strong local inflammatory responses and immunopathology. 6 Autopsy findings in patients who have succumbed to COVID-19 are characterized by perivascular T cell infiltration, microangiopathy, and widespread thrombosis in lung tissue. 28 The induction of IL-6 during severe COVID-19 has led to trials of blocking antibodies to the IL-6 receptor with mixed results. This is inspired by cytokine release syndrome (CRS), seen in cancer immunotherapy, which is also often treated with IL-6-blocking agents. However, there are a number of differences between severe COVID-19 and CRS, such as lower IL-6 levels and death caused by respiratory failure and thrombosis, rather than from circulatory failure and status epilepticus, as seen in CRS. 29 In this respect, the mechanisms of severe COVID-19 are incompletely understood, and better understanding is required for improved immunomodulatory therapies to be devised and immunopathology and mortality limited.
Human immune systems are highly variable, 30 and most of this variation is explained by environmental exposures, 31 particularly early in life. 32 The role of genetic variation in immune variation in general and in COVID-19 in particular is under investigation. 33 Systems-level analysis methods are useful in human immunology because they capture the many variable cell populations, proteins, and transcriptional programs involved in a complex immune response. Systems-level analyses also allow for the inference of relationships among such immune system components. 34 With this study, we add to the rapidly growing literature by providing a longitudinal, systems-level perspective on the immune system changes from acute to recovery phases of severe COVID-19 disease. Longitudinal analyses are important because cross-sectional analyses carry the risk of capturing snapshots of patients at different stages of the immune response and thereby misinterpret differences as qualitatively different. The longitudinal sampling presented herein is a strength of this study. Another Article ll OPEN ACCESS important aspect of this work is its use of whole blood, rather than peripheral blood mononuclear cells, allowing neutrophils and other granulocyte populations to be included in the analysis and also reduce the technical sources of variation caused by cell preparation and freezing. 35 By using this more holistic and longitudinal approach to analyze the immune response during severe COVID-19, we find previously unappreciated roles of eosinophils in the acute response. These cells play important roles in other respiratory infections, 36 but they have not been implicated much in COVID-19. The population of eosinophils that expand a few days after admission to the hospital were characterized by high CD62L expression, a previously reported marker of lung eosinophils, 18 and it is possible that such IFN-g-mediated upregulation of CD62L on eosinophils leads to the influx of these cells into the lung tissue. The development of ARDS and clinical deterioration is typically seen after $1 week in severe patients, and it is possible that this IFN-g-eosinophil axis contributes to this, but a causal role is beyond the scope of our study. The finding that basophil levels are positively associated with humoral responses to SARS-CoV-2 is intriguing and in line with previous studies in other viral infections. 26 Further investigation will be required to understand the mechanisms involved, but it likely would not involve the production of IL-6 by basophils, given that plasma levels of this cytokine were inversely associated with anti-RBD IgG titers. Another possible mechanism involves IL-4 production by basophils, known to potentiate B cell responses to infection in other settings. 37 It is worth noting that time after admission is taken into account in the mixed-effects model, and the reported associations are the significant ones after time from admission is taken into account.
There has been a lot of concern around antibody responses to SARS-CoV-2, and although nearly all of the patients with severe disease do produce antibodies in rather high titers, 19,24 the neutralizing capability of such antibodies are variable. 21 One hypothesis brought forward as a possible explanation of the severe disease occurring often after 1 week or so of stable disease is antibody-dependent enhancement (ADE). 38,39 This occurs when non-neutralizing antibodies bind a virus and via Fc receptors bring viruses into new cell types, not expressing the receptor required for viral entry-in this case, ACE2. Such responses are well known for dengue virus infection and could induce hyperinflammatory responses also in COVID-19. We have found that a significant proportion of CD4 + T cells in some patients showed CD4 downregulation as a sign of possible cell activation, but such downregulation can also occur if T cells are directly infected. 40 CD4 + T cells do not express ACE2, 41  Article ll OPEN ACCESS express Fc receptors and thus be subject to viral infection and replication via ADE. This is speculative at this time, but as more data surface on determinants of neutralizing antibody responses, the theory of ADE as a cause of severe COVID-19 will be testable and have important implications for vaccine development. 38 The influence of basophils in modulating humoral responses to SARS-CoV-2 uncovered herein should also be taken into account, as basophils are depleted during acute disease and the severity of such depletion may be an important determinant of the antibody response to the virus.

Limitations of Study
This study has several limitations. We performed longitudinal systems-level immunomonitoring of acute and recovered COVID-19 patients, but due to logistical limitations in the overwhelmed hospital wards, we were unable to collect longitudinal Article ll OPEN ACCESS samples for recovered patients and were also limited in the number of acute COVID-19 disease patients we could enroll. We were unable to include a healthy cohort, and we were not powered to robustly compare patients in ICU versus non-ICU with respect to their immunological trajectories.

STAR+METHODS
Detailed methods are provided in the online version of this paper and include the following:   Multiomics Factor Analysis (MOFA) was used to discover principal sources of variation in multi-omics datasets. 27 MOFA uses a set of data matrices as input formatted with features as columns and samples as rows, plasma protein expression and cell abundance datasets were used to build the MOFAobject with MultiAssayExperiment. The fitting step includes training the model with the multiomics data in order to be able to disentangle the heterogeneity into a small number of latent factors. The MOFAobject was trained in R through the reticulate package with 10 factors and a variance threshold of 0.01%. Both omics datasets were processed individually to remove any features resulting in zero or low variance before fitting the model. Convergence of the model was assessed using the change in ELBO (deltaELBO) to verify it fit the convergence threshold which is considered to be between 1 and 10. Multiple models were run under different initializations to validate that factors were consistent across trials for model selection. The fitted MOFA model could then be interrogated in R for downstream analysis to characterize these factors as technical or biological sources of variation.
Partition-based graph abstraction of single-cell data The CyTOF data were first preprocessed with arcsin h and scaled to unit variance and then partitioned into different subpopulations according to our in-house supervised learning algorithm. For each subpopulation, the phenotypic changes over different time points are inferred with a trajectory inference method called PAGA. 16 In brief, PCA was first applied to reduce the number of features to 20, and then an undirected kNN-like graph was constructed using the approximate nearest neighbor search within UMAP, while each node represents a single cell and each edge represents a neighborhood relationship. After the construction of graph, the highly connected clusters were detected with Leiden method. 45 Afterward, the clusters defined by Leiden were used by PAGA to infer a trajectory map. In the trajectory map, Leiden clusters are considered as connected if their number of inter-connections is larger than a fraction of the number of inter-connections expected under random assignment, and the threshold fraction is determined by a statistical model. Finally, the PAGA graph was taken as the initial position of a manifold learning method ForceAtlas2 (FA) 46 and produced topology-preserving single-cell embeddings for visualization.

Mixed effects modeling
A partially Bayesian method was applied with blme package on both datasets (plasma protein expression and cell abundance) to produce maximum a posteriori (MAP) estimates. 25 This provided the ability to nest the variables, and account for days from admission as well as RBD levels as fixed effects. Wald p values of covariates were extracted from models to assess significance.