Mapping the recovery of critically ill COVID19 patients by high-dimensional proling identies longitudinal blood immunotypes

The COVID-19 pandemic poses a major burden on health-care and economic systems across the globe. Even though a majority of the population only develops minor symptoms upon SARS-CoV2 infection, a signicant proportion are hospitalized at intensive care units (ICU) requiring critical care. While insights into the early stages of the disease are gradually expanding, the dynamic immunological processes occurring in critically ill patients throughout their recovery at ICU are far less understood. Here, we have analysed longitudinally collected, whole blood samples of 40 surviving COVID-19 patients during their recovery at ICU using high-dimensional cytometry by time-of-ight (CyTOF) and cytokine multiplexing. Based on the neutrophil to lymphocyte ratio (NLR), we dened 4 sequential immunotypes during recovery that correlated to various clinical parameters, including the level of respiratory support at concomitant sampling times. We also identied classical monocytes as the rst immune cell type to recover by restoring HLA-DR-positivity and by reducing the immunosuppressive CD163+ monocyte population, followed by the recovery of CD8+ and CD4+ T cell, and mDC populations. The identied immunotypes also correlated to aberrant cytokine and acute-phase reactant levels. Finally, integrative analysis of cytokines and immune cell proles showed a shift from an initially dysregulated immune response to a more coordinated immunogenic interplay, highlighting the importance of longitudinal sampling to understand the pathophysiology underlying recovery from severe COVID-19.


Introduction
The SARS-CoV2 coronavirus outbreak has infected more than 40 million people worldwide causing more than 1 million deaths and was o cially declared a pandemic by the World Health Organization (WHO) [1].
Due to the lack of immunity to SARS-CoV2 and in the absence of an effective vaccine, the unpredictable clinical outcome of COronaVirus-Induced Disease 2019 (COVID-19) has created signi cant strain on medical systems [2][3][4][5]. Until now, major research efforts have tried to unravel the pathobiology underlying the disease course [6,7]. As such, it has been shown that COVID-19 leads to the combination of pathologically elevated levels of pro-in ammatory cytokines, coagulopathy and a dysregulated immune response [6,8,9]. Deep immune pro ling in severe COVID-19 patients revealed excessive amounts of dysfunctional neutrophils, decreased levels of lymphocytes and low levels of antigenpresenting receptors on monocytes and dendritic cells, hindering e cient adaptive immune responses [8,[10][11][12]. In contrast to asymptomatic or mild/moderate symptomatic COVID-19 patients, critically ill patients typically show a biphasic disease course with an early viral stage of disease before ICU admission, followed by an hypoxic phase of disease, characterized by a cytokine storm and an acute respiratory distress syndrome, requiring (invasive) respiratory support, antivirals and immune-modulation to prevent multiple organ failure and death Considering the potentially severe and unpredictable impact of a SARS-CoV2 infection, research has also focused on predicting the disease course of COVID-19. Using a variety of multiplexing technologies, longitudinal blood sampling revealed a dichotomous pro le in cytokine expression levels that could early on differentiate patients prone to develop a mild/moderate versus a severe disease course within the rst two weeks following SARS-CoV2 infection [13][14][15]. In addition, severely diverging levels of immune cell populations, typically reported by aberrant neutrophil-over-lymphocyte ratios (NLR) and high levels of immature neutrophil frequencies, were also indicative in the early stages of the disease for a more severe disease course [16][17][18][19][20][21][22]. Based on these insights, various interventional trials were initiated, from which the use of corticosteroids showed promising results [23] although the optimal timing of these interventions remains to be determined and side effects (i.e. secondary infections) may occur. Therefore, until an effective vaccine is identi ed, patients exhibiting severe pathology will continue to be admitted to ICU requiring extensive medical care.
In contrast to the early stages of the disease, the biology from ICU admission to the recovery phase remains far less understood. A systems-biology study including the longitudinal follow-up of 10 ICU patients, suggested involvement of eosinophils early in the disease course [24]. Various immune cell populations were also identi ed to be increased in relative proportions during recovery while excessive cytokine levels returned to normal levels at the time of discharge from ICU [24,25] Also the restoration of T cell functionality, including the induction of effector-memory cells, has been described during the recovery phase of COVID-19 [26]. However, the exact sequence of changes occurring during the restitution of the peripheral immune system after severe infection remains elusive.
Within the transdisciplinary clinical COntAGIouS trial (COvid-19 Advanced Genetic and Immunologic Sampling, NCT04327570), we have collected and analysed longitudinal whole blood and serum samples of critically ill COVID-19 patients that survived a stay at ICU at the University Hospitals Leuven, Belgium.
By using high-dimensional analysis technologies, we identi ed 4 immunotypes (i.e. blood cell composition pro les) that were based on the NLR, de ning a common biological trajectory during patient recovery which correlated well to various clinical parameters, including the required level of respiratory support at concomitant sampling times. We identi ed classical monocytes as the rst cell type to recover, followed by CD8+ and CD4+ T cell, and mDC populations, while the NLR-based groups also correlated to cytokine and acute-phase reactant (APR) levels. Finally, integrative analysis of cytokines and immune cell pro les, allowed us to monitor the recovery from a dysregulated immune response, highlighting the importance of longitudinal sampling to better understand the pathophysiology underlying the recovery from severe COVID-19.

Patient cohort and sampling
Forty COVID-19 patients that were admitted to the ICU were enrolled in the Contagious trial (NCT04327570). Thirty eight patients were con rmed positive for SARS-CoV2 infection by nasopharyngeal swab and/or by bronchoalveolar lavage (BAL) sampling and subsequent PCR quanti cation. Two patients were con rmed positive for SARS-CoV2 via CT scan. Routine clinical laboratory tests at the time of sampling were performed in line with the 'COntAGIous' study design [12]. The patients enrolled in the study were sampled at several relevant time points: upon admission to the ICU (0-2 days post admission), at 6 +/-2 days post admission, and eventually at ICU discharge. Of note, if the medical condition of the patient was not improving, a bronchoscopy was performed during which whole blood samples were again collected for concomitant immune pro ling.
The healthy control group consisted of volunteers (recruited from UZ Hospital staff) with no prior diagnosis or recent symptoms compatible with COVID-19. This group was not tested by nasopharyngeal swab, but antibody assessment of plasma samples revealed absence of COVID-19 IgG antibodies.
Sample Processing and Staining Procedure for CyTOF Whole blood (WB) samples were collected into sterile anticoagulant lithium heparin blood collection tubes and processed for mass cytometry staining within 2-4 hours of isolation. WB was stained with the Maxpar® Direct™ Immune Pro ling Assay (DIPA) kit from Fluidigm© by following the work ow outlined for whole blood staining. Additional antibodies not included in the Maxpar DIPA kit were titrated on WB and RBC-free samples from donors to determine the optimal staining index. The antibodies are outlined in Suppl. Table 1. The additional markers were added after the red blood cell-removal step. The last step of the protocol is the 193-Ir intercalation, which was performed overnight at 4°C. The samples that could not be acquired on the instrument the very next day were immediately cryopreserved in Maxpar Fix/Perm buffer + 193-Ir at -80°C. The cryopreservation technique was validated in triplicate by dividing aliquots of donor samples stained on the same day and comparing cell viability and immune pro les between the fresh and cryopreserved samples. Batch effects were evaluated by running a reference sample derived from an aliquot of the same healthy donor over the period of the study.

Data Acquisition
Cells stained for mass cytometry were acquired the day after the staining procedure or within 1 week of cryopreservation. For CyTOF acquisition, the samples were pelleted in Maxpar Cell Acquisition Buffer (CAS) on the day of acquisition and transferred to the KU Leuven Flow and Mass Cytometry Facility to be acquired on a Helios® mass cytometer (Fluidigm). Prior to acquisition, 1 million cells/ml were resuspended in CAS containing EQ beads (1:10) and ltered through 35µm cell strainer cap tubes. The samples were acquired at a rate of 250-300 events per second on a Helios® Mass Cytometer. CyTOF software version 6.7.1016 and the Maxpar Direct Immune Pro ling Assay.tem template were used to acquire and normalize data from the stained samples.

Chemokine and cytokine assays
Chemokine and cytokine levels in plasma were assessed by Meso Scale Discovery using the V-plex human cytokine 30-plex kit, complemented with Human IL-1RA (V-plex) and human IL-18 (U-plex) kits.

Data Analysis
To analyse the data, we continued using the FlowSOM models built for [12]. To generate these models, normalized .fcs les were rst transferred to the Maxpar Pathsetter™ software (version 2.0.45), which performs a standardized, automatic and unsupervised quality check (bead removal and high-quality singlet selection) and gating of live single cells. Samples with fewer than 50000 cells were discarded, as were 4 samples which showed up as outliers in a PCA analysis of the 25, 50 and 75 percent quantiles of the marker values. The FlowSOM analysis was subsequently done in R, using 123 cleaned .fcs les, including 8 healthy controls and longitudinal samples of 40 COVID-19 patients, through various stages of the disease. Samples were rst pre-processed: margin events were ltered out, data was transformed with an arcsinh transformation with cofactor 5 and the PeacoQC algorithm (manuscript in preparation, tool available at https://github.com/saeyslab/peacoQC) was applied to remove any unstable signal regions during the measurement. On this cleaned data a rst FlowSOM model [27] was trained, using a random selection of cells for all samples, resulting in 3,000,093 cells to train on. The clustering made use of 11 markers (CD3, CD4, CD8a, CD11c, CD14, CD19, CD20, CD45, CD66b, TCR δ and NCAM), mapped the data onto a 10 by 10 SOM grid and resulted in 30 meta-clusters. 22 meta-clusters were selected as having CD66b values lower than 2 or CD45 values higher than 4, corresponding to non-granulocytes, while 8 meta-clusters were labelled as granulocytes. The full les were mapped onto this model, and for each of them new fcs les were generated corresponding to the two subsets of cells. These were then used to build two separate FlowSOM models, including either the non-granulocyte or the granulocyte cells, again using only a subset of 2,949,946 (granulocyte: 3,000,093) cells for training mapped onto a 10 by 10 SOM grid, this time using 33 markers (CD3, CD4, CD8a, CD11c, CD14, CD16, CD19, CD20, CD27, CD28, CD38, CD45, CD45RA, CD45RO, CD57, CD66b, CD69, CD294, CD161, CD163, CCR4, CCR6, CCR7, CXCR3, CXCR5, HLA-DR, IgD, IL-2Rα, IL-3R, IL-7Rα, NCAM, NKG2A and TCR δ). To be able to identify small populations, no meta-clustering was applied on these second models, and the 100 clusters were manually annotated by 3 independent experts according to their mean uorescence intensity (MFI) values. In the non-granulocyte model, 6 clusters were manually identi ed as still being mixtures of different cell types, and split into 2 or 3 clusters, resulting in 107 non-granulocyte clusters and 100 granulocyte clusters.
For this study, all samples were again fully mapped onto these models to identify their immune pro les. During this process, we noticed some les had a sub-optimal granulocyte / non-granulocyte split. Therefore, we rede ned this split using the le speci c MFI values of the metaclusters, selecting all metaclusters with CD66b values lower than 3 and CD45 values higher than 3 as non-granulocytes. For one le, 3 metaclusters were still wrongly assigned when inspecting the density distribution of the cells, and for this le the CD45 boundary was adapted to 4. Once all cells were assigned to the 207 clusters from the existing models, we used a manual annotation of these clusters to aggregate them in 24 cell populations, describing the immune pro le of each sample by the percentage of live cells each of these populations represented. Additionally, the total lymphocyte percentage (i.e. the sum of all populations which were annotated as lymphocyte subsets) was added, bringing the total to 25 parameters per sample. The samples were grouped in 4 equal sized groups depending on their NLR: 31 samples with a ratio higher than 7.1 (R4), 30 samples with a ratio between 7.1 and 4.8 (R3), 31 samples with a ratio between 4.8 and 2.7 (R2) and 31 samples with a ratio lower than 2.7 (R1). The four pro les were ordered by decreasing neutrophil/lymphocyte ratio, and are throughout the gures in this paper colored as black, green, blue and cyan respectively.
The differences between these four groups were then further characterized. To compute statistics, we made use of linear mixed models, and for the cytokines and ARP values we used mixed models with random effect for the patient, to correct for the fact that patients can have multiple longitudinal samples, spread over one or more groups. Correlations were calculated using the Spearman Rank test. All p-values were taken together to be corrected for multiple testing with the Benjamini-Hochberg (FDR) method. Curve tting analysis was done using the drc package in R from which the in ection points and the corresponding p-values were collected. Additionally, we also used the built-in immuno pro ling tools of the Maxpar Pathsetter™ software to analyse the overall immune cell pro le in a standardized and automated way. The Pathsetter-based blood pro les were grouped into 12 clusters using k-means clustering, and manually further grouped into 4 nal pro les. These pro les strongly corresponded with the NLR results and led to similar biological conclusions, thus giving additional trust in the results.
The code used to generate the results is available at https://github.com/saeyslab/CYTOF_covid19_study.

Results
Mapping the disease course of COVID-19 survivors at ICU To unravel the sequential recovery of the immune system throughout the course of a severe COVID-19 infection, we have analysed a cohort of 40 patients that survived a stay at the intensive care unit (ICU) at the University Hospitals Leuven, Belgium. These patients were enrolled in the COntAGIouS trial (NCT04327570) through which longitudinal whole blood samples were collected during their ICU stay (see methods and Table 1 for patient demographics and characteristics; a set of healthy controls was included for reference).
For each patient, a detailed timeline of their ICU stay was generated displaying the most clinically relevant events including onset of symptoms, hospital admission, ICU admission, start/stop of therapeutic interventions and nally discharge from ICU and hospital. If applicable, discharge from a specialized revalidation centre where the patient continued his/her recovery after their hospital stay, is also included.
Occasionally, blood was drawn at additional time points (such as during bronchoscopy), which is indicated on the timeline. Once individual time points were plotted ( Figure 1A), a large range in the duration of ICU stay (from 2 to 72 days) was found, indicating a highly variable recovery rate across patients. From the 40 enrolled patients, only 9 were women, corresponding to previous observations of males experiencing increased severity and hospitalization rates once infected with COVID-19 [28][29][30] ( Figure 1B). The age range and overall time spent at ICU of the admitted patients was however comparable for men and women ( Figure 1B and Table 1). Within this critically ill patient cohort, we did not nd a correlation between BMI and the overall duration of ICU stay (r=-0.15, p=0.35), despite >40% of the patients having a BMI above 30 (Table 1, Figure 1B) [31].
This timeline was further extended with the level of respiratory support that was required at concomitant sampling time points (Figure 1C, D). Here, respiratory support levels were classi ed from 0-5 where level 0 indicates no support, while level 4 and 5 indicate the need for mechanical ventilation (with level 5 indicating patients that were ventilated in prone position and/or requiring inhalation NO therapy or extracorporeal membrane oxygenation (ECMO)) ( Table 1). Subsequent correlation analysis showed that the maximal level of respiratory support correlated well to the overall duration of ICU stay (r = 0.67, p = 7e-17, Figure 1C).
Identifying peripheral blood pro les using longitudinal sampling Next, we performed deep immune pro ling on collected blood samples using high dimensional cytometry by time-of-ight (CyTOF) analysis. From this data, 33 surface markers were used to map changes in composition and phenotype of white blood cells during ICU stay. All data were mapped using unsupervised gating methods and manually annotated to identify the immune cell subtypes to reconstruct potential dynamic changes. However, given the variation in duration of hospitalization at the ICU in this patient cohort, mapping immune pro les on an absolute time scale did not produce coherent results. Alternatively, we ranked samples based on the neutrophil-to-lymphocyte ratio (NLR) and de ned 4 equal groups (Figure 2A-C) across all collected samples, regardless of when a sample was taken. As such, we de ned 4 "NLR" stages (R1-4) ranging from severe lymphopenia and neutrophilia (pro le R4, which showed the highest NLR values) to immune pro les resembling those of healthy controls, de ned by the lowest NLR values (R1) (Figure 2A, B and Suppl. Figure 1). In addition, a high NLR (R4) also correlated to increased levels of the acute phase reactants C-reactive protein (CRP) and ferritin ( Figure 2C, D), markers that are commonly monitored in ICU patients and de ne a critical state.
Integrating longitudinal immune-pro ling with clinical features Next, we assessed how the identi ed NLR groups correlated to a variety of clinical parameters. As indicated above, the level of respiratory support was measured at each sampling time point and correlation analysis of this feature with the 4 NLR pro les revealed that patients exhibiting R4 required increased respiratory support at that moment compared to lower level NLRs ( Figure 2E). Further correlation analysis also showed that the overall clinical status of each patient upon admission to ICU, as determined by calculating the Sequential Organ Failure Assessment Score (SOFA) [32,33], was signi cantly correlated to the NLR groups ( Figure 2F; r=-0.23, p=0.03) [33]. We did not nd any correlations of age or BMI with the NLR groups (Suppl. Figure 3A, B).
Using these groups, we reconstructed a detailed sequence of events for each individual patient and investigated how these evolved along their ICU stay ( Figure 1E, 2G, Suppl. Figure 4). For 36 out of 40 patients, we observed a recovery (overall downwards shift of the NLR score) or stabilisation of the immune system (and NLR score) along their stay at ICU (Suppl. Figure 4), with an enrichment of R1 towards ICU discharge ( Figure 2G). Strikingly, the R values at ICU discharge still ranged between R1-4 suggesting that the physical condition of these patients was still heterogeneous at that moment. In line with this observation, we found that patients exhibiting a higher NLR value upon discharge from ICU required signi cantly longer revalidation (either at the regular hospital ward and/or a specialized revalidation center) following their stay at ICU ( Figure 2H, r=-0.4, p=0.02). Additionally, Charlson comorbidity indices [34] were low for all patients ( Table 1), indicating that the majority of patients admitted to ICU were in good general health before their infection with SARS-CoV2. We also did not observe correlations of the Charlson comorbidity indices to either the NLRs or the required level of respiratory support (Suppl. Figure 3C, D).
Finally, we tracked alterations in patient's immune pro les who were treated with diverse immunomodulatory regimens. Commonly, the glucocorticoid methylprednisolone (MP) was given as a monotherapy (n=19). In some cases, MP was prescribed in a combination with either anti-IL-1 (anakinra; n=3) or anti-IL-6 (tocilizumab; n=5). One patient received anti-IL-6 exclusively, while the remaining 12 patients in our cohort relied on supportive care alone without immunomodulatory treatments (Table 1; Figure 1F). Unsurprisingly, our data analysis con rmed that patients who underwent steroid treatment, had a greater tendency towards a longer ICU stay. Correlation analysis of the NLR pro les with MP treatment, revealed that patients receiving MP exhibited a higher NLR upon ICU discharge, as opposed to patients that did not receive MP (Suppl. Figure 3E). Finally, while the cohort of patients that received anti-IL-1 and anti-IL-6 were small, the reconstitutional trajectory of their immune systems followed a similar pattern to those patients receiving standard-of-care treatment without immunomodulatory drugs and/or MP.

Reconstructing the cellular recovery of critically ill COVID-19 patients
Considering that the NLR levels correlated well with the clinical condition of the patients (see above), we next investigated the recovery of more speci c immune cell populations to eventually de ne a sequence of events. To do so, we used a dual approach: rst, we performed curve tting along the trajectory of the four NLR pro les, in which we determined the in ection points of each identi ed immune cell population to estimate the moment recovery would begin ( Figure 3A, B, Suppl. Table 2). The data from this approach were subsequently combined with the statistical comparison of the 4 NLR groups (Figure 2A, suppl. Table  3) to de ne the eventual order by which speci c cell populations were recovering. We found that, on average across this cohort of patients, classical monocytes were the rst immune cell population to recover followed by naive CD8+ T cells, naive and Th2-polarized CD4+ T cells. As indicated by our analysis, effector and memory T cell populations were only restored at later stages, including the antiviral Th1-polarized CD4+ T cells [35]. From the other professional antigen presenting cells, myeloid dendritic cells (mDC) were recovering immediately after the early T cell response, while plasmacytoid DCs (pDC) would only recover towards R1, similar to overall B and NK cell populations. This suggests that monocytes, as implied by other studies, could become important targets in both understanding COVID-19 disease progression, as well as, improving recovery at the early phase of ICU admission [36].

Functional recovery of monocytes
Antigen-presentation by monocytes, pDCs and/or mDCs plays a vital role in the initiation of an e cient adaptive immune response to viral infections [37][38][39][40][41]. As previously shown, declined levels of HLA-DR expression and hence a reduced antigen presentation seems to be an early hallmark of a trajectory towards a severe COVID-19 disease course compared to a mild/moderate one [10,12]. Considering that monocytes are among the rst immune cell types to recover, we further investigated the functionality of these cells across the 4 NLR pro les during recovery. As such, we observed a regain in both the numbers and the antigen presenting capabilities of the monocytes towards discharge ( gure 3C). The rst subset to re-establish were the classical HLA-DR+ monocytes (Cluster (Cl)3, Cl33 and Cl12; gure 3C), while it was only in later stages of disease progression that the non-classical HLA-DR+ monocytes were restored (see Cl11, Figure 3C; see also statistics in Figure 2A).
The early phase of severe COVID19 is characterized by a functional shift towards a more immunosuppressive spectrum of monocytes, as seen by a downregulation of HLA-DR and an enrichment of CD163+ monocytes in R4 (Figure 3) [10,12]. However, longitudinal follow-up shows that this shift is reversible, as seen by the recovery of HLA-DR expression and a relative reduction of the amount of immunosuppressive CD163+ monocytes ( see Cl53 and Cl32; Figure 3C), leading to a regain of the antigen-presenting phenotype of these cells.

Longitudinal cytokines pro ling
Next to immunophenotyping, we also performed multiplexed analysis of soluble serum proteins, which included 32 pro-in ammatory cyto-and chemokines in serum samples at concomitant sampling times.
The levels of these analytes were subsequently correlated to the NLR pro les. In line with previous observations, the initial R4 stage of patients arriving at ICU was characterized by an increase in IFN-, TNF-α, IL-2, IL-6, IL-7, IL-10, IL-15, IP-10, MCP-1, MIP-1α, MIP-1β levels and a reduction in TARC, MDC and IL1-α, suggestive of a pro-in ammatory cytokine signature, often referred to as a "cytokine storm" ( Figure  4A). This signature was steadily reversed as patients attained a normal R1 stage. Similarly, CRP and ferritin were also signi cantly higher in R4 and gradually decreased as patients reached R1 suggestive of a return to baseline following an acute phase induction of the immune system ( Figure 2D,E). Clustering analysis to uncover patterns in cytokine pro les across the various NLR groups further con rmed that pro-in ammatory cytokine levels collectively reduced along the duration of an ICU stay ( Figure 4B). On the other hand, MDC and TARC, two constitutive chemokines designated as CCL17 and CCL22 that are regulated at a post-translational level, increased to normal levels [42].
Even though the levels of cytokines within the R1 pro le remained comparable between COVID-19 patients and healthy individuals, we observed altered expression levels of TNF-α, IL-18, IL1-α, MCP-1 and MIP-1α, between the two study groups( Figure 4A; healthy controls are indicated as hollow black circles in the gure). It remains to be studied how long it would take for levels to normalise, however, in spite of their aberrant pro les, patients were still able to leave the ICU for their recovery.
Integrative mapping of the immune response during recovery Considering the important interplay between immune cells, cytokines and chemokines during severe COVID-19, we also performed an integrative, similarity matrix-based correlative statistical modelling analyses ( Figure 5), to uncover associative patterns along the duration of a patient's recovery at ICU. As recently determined, critically ill COVID-19 patients are characterized by a strong dysregulation of the immune reaction against the SARS-CoV2 virus, whereby normally highly concerted interplays of cyto/chemokines and speci c immune cell populations are now disentangled [8]. In line with above observations, our similarity matrix analyses also revealed an intense correlation between various immune cells and cyto/chemokines (a putative indication of ongoing immunological interactions) in the R4 critically-ill COVID-19 patients ( Figure 5A), such that these correlations progressively "normalized" or decreased in terms of number of intense clusters when traversing from R4 to R1 ( Figure 5A-D), thereby indicating that apart from above quantitative shifts, there was also a qualitative shift in possible immuneinteractions when going from critically-ill to recovering COVID-19 patients.
On the level of lymphocytes, the most striking phenotypes were observed in the R4 pro le ( Figure 5A) whereas patterns in R3/R2 were largely transitional and ultimately culminated into "contracted" lymphocytic correlative-compartment in R1 ( Figure 5B-D). More speci cally, in R4 patients ( Figure 5A), the typically infection-resolving lymphocytic compartment (e.g. Th1/Th2 cells, effector/memory CD4/CD8 T cells) had relatively less correlations with various effector-function cytokines including IFN-, TNF-α, IL-2, IL-15, thereby indicating a certain degree of immunological dysregulation. Most of these lymphocytes also had considerable negative correlation with neutrophils. Interestingly TNF-α, Th17 and Tregs formed a cluster together which might be a sign of pro-in ammatory signalling since Tregs and Th17 have been shown to reciprocally stimulate each other via TNF-signalling pathway [44]. We believe that in the current context, this crosstalk might play a disease-potentiating role in COVID-19 severity. These discords were largely ameliorated from R3 to R1 (e.g., Th17 gaining correlation with Th1/Th2/Tfh and activated-CD8 T cells in R1) thereby indicating better lymphocytic activity/regulation to be bene cial for recovery of COVID-19 patients.
In conclusion, above qualitative analyses revealed that neutrophil/lymphocyte-associated in ammation undergoes considerable changes in terms of co-associative immune-components between R4 to R1 such that a relatively contracted neutrophil cluster with a chemotactic/in ammasome-footprint and a contracted lymphocyte cluster de nes favourable recovery for COVID-19 patients.

Discussion
We have performed a longitudinal study in which consecutive peripheral blood samples of 40 surviving COVID-19 patients were analysed via high-dimensional mass cytometry analysis, allowing us to uncover detailed insights in the amount and functionality of the white blood cell (WBC) compartment. Importantly, by using fresh, non-cryopreserved whole blood samples, our study allowed the analysis of both granulocytes, from which neutrophils take the mainstay and play an important role in COVID-19 [45], in addition to the most important populations of lymphocytes.
Using those insights, we grouped all measured samples based on the NLR, as such identifying 4 blood pro les that de ned a biological trajectory through which patients evolved at a highly individualized pace during recovery. Importantly, these NLR pro les correlated well to various clinical variables, including the level of respiratory support needed at concomitant sampling times. Another striking observation was the large variety of NLR levels upon discharge from ICU, which indicates that patients upon discharge from ICU, still exhibit a broad range of physical tness, as suggested by the need for a longer revalidation following discharge from ICU when an increased NLR was measured.
Next to general clinical features, this cohort also included patients that received various therapeutic interventions, from which the use of corticosteroids (CS) was the most prominent. Our analysis showed that the likelihood of patients receiving CS therapy, correlated well with the duration of ICU stay, indicative of a good clinical selection of patients for this intervention. In addition, correlation analysis of the NLR pro les between patient groups that did and did not receive CS therapy, suggests that while both groups exhibited similar trajectories from high to lower ratio levels, the CS group showed a slower evolution towards a normal ratio. While it is known that CS therapy dampens the overall immune response [46], the excessive levels of cytokines and concurrent defective immune response as observed in COVID-19, warrants treatment with CS, in line with other studies [23]. Due to the small number of patients receiving anti-IL1 or anti-IL6, we cannot make conclusions towards changes in the recovery patterns of their peripheral immune system, even though we did not observe differences with those patients that did not receive these treatments. Overall, the timing of these therapeutic interventions remains critical and still requires further optimization [47,48]. Indeed, one of the major problems seems to be the delay of the innate immune system to e ciently respond to SARS-CoV2 infection, which further delays the subsequent adaptive immune response. While this cohort is small, it would be useful to track the phenotypes of severe "long-stay" vs "short-stay" patients as well as non-hospitalized COVID "long haulers" to identify early biomarkers in order to better predict immune responses to the virus [49].
The de ned biological trajectory also allowed us to further elucidate the sequence of events and study the order in which particular immune cell populations reappear in the peripheral blood including the analysis of certain functionalities. Using this approach, we identi ed classical monocytes to recover rst, followed by CD8+ and CD4+ T cells. While these observations are in line with other smaller studies [25,36], they will need further corroboration in larger cohorts. Overall, the underlying biology suggests that protecting and/or reactivating the innate immune system (from which the monocytes seem rst to recover), could be a potential therapeutic strategy. Indeed, based on these fundamental insights, therapies aimed at improving early monocyte functionality could be an interesting venue to improve the disease course, as suggested by the role for MCP1/CCR2 interference therapy in SARS-CoV1 [12,50] but also in COVID-19 as previously suggested [36]. We should however also look beyond monocytes: a reactivation of the other members of the innate immune system, including DCs and NK cells, could have the potential to ignite a stronger adaptive immune response and recovery. It also remains to be seen how speci c the observed patterns are for COVID-19 and how they diverge from other severe respiratory diseases to tailor future interventions. Moreover, the absence of a comparator cohort on non-survivors could make these observations stronger.
Finally, next to the immune cell pro les, we also assessed concurrent levels of various cytokines and chemokines and acute phase reactants (i.e. CRP and ferritin). Clustering analysis revealed the gradual release of the pathologically elevated cytokine and chemokine levels, often referred to as a "cytokine storm", even though the levels in COVID-19 patients are still lower compared to some auto-in ammatory disorders e.g. macrophage activation syndrome. When patients evolved toward lower NLR levels, the coordination of the immune system recovered with a better alignment of cytokines and in ammatory cell types. Indeed, qualitative, correlative matrix analyses further revealed that hyper-in ammatory neutrophil clusters and dysregulated lymphocyte clusters distinguished critically-ill COVID-19 patients. Accordingly, amelioration of these features (i.e., contraction of neutrophil-clusters down to neutrophil-speci c immunesignalling and consolidation of lymphocyte clusters to better crosstalk between different T cell subsets) marked the recovery phase in COVID-19.

Declarations
Ethics: This study was approved by the medical ethical committee of UZLeuven/KULeuven (study number S63881) and is also registered here: https://clinicaltrials.gov/ct2/show/NCT04327570 Competing Interests: The authors declare no competing interests.