key: cord-0753708-d66tcbtc
authors: Frishberg, Amit; Kooistra, Emma; Nuesch-Germano, Melanie; Pecht, Tal; Milman, Neta; Reusch, Nico; Warnat-Herresthal, Stefanie; Bruse, Niklas; Handler, Kristian; Theis, Heidi; Kraut, Michael; van Rijssen, Esther; van Cranenbroek, Bram; Koenen, Hans JPM.; Heesakkers, Hidde; van den Boogaard, Mark; Zegers, Marieke; Pickkers, Peter; Becker, Matthias; Aschenbrenner, Anna C.; Ulas, Thomas; Theis, Fabian J.; Shen-Orr, Shai S.; Schultze, Joachim L.; Kox, Matthijs
title: Mature neutrophils and a NFkB-to-IFN transition determine the unifying disease recovery dynamics in COVID-19
date: 2022-05-17
journal: Cell Rep Med
DOI: 10.1016/j.xcrm.2022.100652
sha: 0f6a93b8b5c9c1af81b714a645f3c05e109fabf0
doc_id: 753708
cord_uid: d66tcbtc

Disease recovery dynamics are often difficult to assess, as patients display heterogeneous recovery courses. To model recovery dynamics, exemplified by severe COVID-19, we apply a computational scheme on longitudinally sampled blood transcriptomes, generating recovery states, which we then link to cellular and molecular mechanisms, presenting a framework for studying the kinetics of recovery compared to non-recovery over time, and long-term effects of the disease. Specifically, a decrease in mature neutrophils is the strongest cellular effect during recovery, with direct implications on disease outcome. Further, we present strong indications for global regulatory changes in gene programs, decoupled from cell compositional changes, including an early rise in T cell activation and differentiation, resulting in immune rebalancing between interferon and NFkB activity and restoration of cell homeostasis. Overall, we present a clinically-relevant computational framework for modeling disease recovery, paving the way for future studies of the recovery dynamics in other diseases and tissues.

Recovery from many diseases is characterized by heterogeneous dynamics in individual patients, thus it is often difficult to identify generalizable cellular and molecular mechanisms characterizing and contributing to the recovery process. So far, many studies utilized transcriptomics data to investigate the biological alterations between mild and severe patients mostly at early time points describing disease progression, as well as discovering major changes in molecular mechanisms and immune cell compositional changes in both the myeloid and lymphoid compartment [1] [2] [3] [4] [5] [6] . Additionally, it was established that the recovery process in severe COVID-19 varies dramatically between patients with severe disease and can last between several days to many months [7] [8] [9] . However, whether common molecular mechanisms of recovery from severe COVID exist, despite the heterogeneity of disease recovery, is still unclear. We used data of severe COVID-19 patients admitted to the Intensive Care Unit (ICU) to study recovery processes at the cell, pathway and gene level.

Studying the molecular and cellular principles of ICU recovery over time requires longitudinal data, including multiple time points from different patients, during their recovery process. As the kinetics of recovery are patient-specific, and patients may enter the ICU at varying disease stages, analysis of fixed time points (e.g. using days after admission to the ICU), might not be adequate to determine the underlying molecular and cellular changes over time in a patient cohort. This is similarly true when addressing disease progression, which also differs between patients. To capture disease progression dynamics in a patient cohort as well as its molecular underpinnings, a method called TimeAx was recently developed 10 . Briefly, TimeAx assembles a disease consensus trajectory, shared across all patients, based on longitudinal molecular data. It then allows the projection of any new sample on top of this trajectory, predicting its corresponding disease state.

To study the disease dynamics during recovery from severe COVID-19, we utilized the TimeAx modeling approach, so far only used to describe disease progression, on whole blood transcriptomes in a longitudinal cohort of COVID-19 patients during their stay at the ICU. Despite very heterogeneous time courses of disease, we identified general principles of cellular and molecular changes across the recovery of these patients, which could not be deduced from studying the original time-course data. Finally, studying fatal outcomes using the recovery trajectory angle translated to the identification of putative disease outcome biomarkers.

To identify disease trajectories towards recovery in severe COVID-19, we conducted a longitudinal blood transcriptomics study in patients with severe disease (cohort 1) and validated the findings in patients derived from different hospitals across Europe (cohort 2 and 3) ( Figure  1A) . The longitudinal cohort 1, collected at Radboud university medical center in Nijmegen, the Netherlands, included 55 COVID-19 patients admitted to the ICU. Between 2 and 9 samples per patient were obtained after ICU admission for whole blood transcriptome analysis, resulting in 260 samples in total ( Figure 1B ). In addition, flow cytometry data covering all major immune cell subsets were generated in 41 patients for 1 to 5 time points, a total of 86 data points. Finally, plasma cytokine levels were measured across patients, ranging between 147 and 252 data points per cytokine. The first validation cohort (cohort 2) included 39 severe COVID-19 patients, from whom whole blood transcriptome samples were collected at a single time point. These patients were admitted to the ICUs of 3 medical centers in Germany, including the University Hospital RWTH in Aachen, Kiel University Medical Center and Saarland University Medical Center in Saarbrücken, as well as to the National and Kapodistrian University of Athens, Greece 3,11 (see Methods, Figure 1A) . The second validation cohort (cohort 3) included data from 45 severe COVID-19 patients admitted to the ICU of either the Addenbrooke's or the Royal Papworth hospital in Cambridge, UK 12 . For each patient, 1-2 consecutive whole blood transcriptome samples were collected during their stay in the ICU, resulting in 81 samples in total. In addition, for 55 of the samples, paired flow cytometry data was available. Finally, similar to cohort 1, plasma cytokine levels were measured for all patients in this cohort.

Specifically, median age within cohorts 1-3 was 65, 63 and 57 years, respectively, ranging from 21 to 86 years ( Figure 1C, Figure S1A ). In all three cohorts, the majority of the patients (76%, 82% and 84%) were males ( Figure 1D and Figure S1B ). Survival rate was lower in older patients ( Figure  1C and Figure S1A ) corroborating previous studies 13, 14 . In cohort 1, 45 of 55 patients recovered from severe COVID-19 and eventually were discharged from the ICU, while 10 patients deceased (survival rate of 82%). Time from onset of symptoms to discharge or death (15 to 97 days) and ICU length-of-stay (9 to 84 days) was highly variable irrespective of outcome ( Figure 1B) . Similar mortality rates were observed in cohort 3, with 8 deceased patients and 37 recovered patients (survival rate of 82%). Cohort 2 had a higher mortality rate, with 16 deceased patients and 23 recovered patients (survival rate of 59%) ( Figure  S1C ). In cohort 1, the patients' body mass index (BMI) ranged from 21 to 37 with 69% of the patients being overweight (Figure 1E ). During their stay in the ICU, all patients were ventilated for 165 to 1652 hours (median of 588 hours; Figure 1F ).

To define whether blood transcriptomes harbor information associated with recovery of severe COVID-19 in the ICU and because the state of disease in an acute infection is difficult to determine a priori, we applied the TimeAx algorithm 10 , a method designed to discover shared dynamics of biological processes across multiple individuals based on time-series data. The TimeAx model predicts the specific pseudotime of each sample within the original time-series data and subsequently allows inferring the pseudotime of new samples, which were not part of the trajectory assembly, by projecting their profiles onto the model (Figure 2A) . To model recovery of severe COVID-19 patients, we assumed that survivors in cohort 1, who were discharged from the ICU, though individually different, also share a common underlying recovery process. We therefore trained the model using gene expression profiles of all non-deceased patients with 5 or more time points (n=20). TimeAx uses the temporal sample collection ordering of the samples from each patient to assemble a robust model of COVID-19 recovery over time, which spans a time course that is longer than that of any single patient (Robustness score = 0.86; Figure S2A ; See Methods). We next related the TimeAx model-based pseudotime assignment of each sample with clinical parameters, starting with disease chronological time. We observed a high correlation between time from onset of symptoms and the time from ICU admission (r=0.63; Figure S2B ), yet weak associations between predicted pseudotime and either of these chronological timelines (Cohort 1: r=0.21 and r=0.2; Figure S2C , Cohort 3: r=0.35; Figure S2D ), suggesting that the dynamics of the trajectories significantly differed between patients (from hereon we thus focus on 'time from onset of symptoms' as 'chronological time'). Furthermore, we identified similarities between disease trajectories of individual patients, with highly differing chronological times, indicating that chronological time was not capturing general disease trajectories, whereas pseudotime was ( Figure 2B) .

We next validated whether the TimeAx model-based pseudotime indeed represented a continuous description of the COVID-19 recovery process. Three lines of evidence strongly support our model: (1) . We observed a negative correlation of TimeAx model-based pseudotime with known severity markers [15] [16] [17] [18] such as CRP, creatinine, IL-6 and G-CSF (p<10 -5 , cohort 1, Figure  2C ; p<10 -6 , cohort 3, Figure S2E ), relationships that were much weaker or absent when chronological time was used instead. (2) We found that the Meta-Virus Signature (MVS), a specific gene signature which was shown to be associated with disease severity in different viral infections, including COVID-19 19,20 had a strong negative correlation with the pseudotime (r=-0.58; Figure 2D , left), a relationship that was lost when chronological time was used (r = -0.09; Figure 2D , right); (3) We demonstrate concordance between WHO COVID scores and trajectory positioning. Specifically, we used TimeAx to predict pseudotime positions of samples from the cross-sectional (i.e. single time point) Cohort 2 and observed that higher pseudotime positioning correlates with lower WHO scores (p=0.0002; Figure 2E ).

Next, we assessed gene expression trajectories as a function of pseudotime in cohort 1. We calculated the difference in correlations (pseudotime versus chronological time) and observed that while the agreement between association patterns were generally kept (r = 0.52, Figure S2F ), the associations of gene expression levels with pseudotime were much stronger, providing a significantly higher resolution of the regulatory processes governing the COVID-19 recovery ( Figure 2F) . Indeed, we observed significant associations for more than 23% of the genes only with the pseudotime, compared to merely 0.2% and 0.7% for associations with the chronological time or with both measures, respectively ( Figure 2G and Figure S2G ). Similar gene associations were obtained using additional time-series analysis tools (See Methods and Figure S2H ), supporting the pseudotime's ability in capturing recovery dynamics over time. We discovered many genes that were either positively or inversely-correlated with the pseudotime, but showed only minor associations with the chronological time ( Figure 2H ). For example, IL-17RA was identified as the top inversely correlated gene, a known marker of disease severity in COVID-19 21, 22 as well as a proposed drug target for COVID-19 23 . In addition, STAT5B, a gene from the STAT5 family which was shown to play a major role in the activation of the ACE2 receptor in COVID-19 24 , was negatively correlated with pseudotime while ITGA4 displayed a positive correlation, supporting previous findings concerning its elevation in mild compared to severe COVID-19 patients 3 ( Figure S2I ).

Overall, we devised a molecular model of COVID-19 recovery over time. Based on this model, we were able to predict a continuous pseudotime state for each sample, in each patient, which successfully simulated the recovery process and its associations with known disease covariates, showing consistently greater value over the use of the chronological time.

An unregulated recovery process is associated with clinical worsening In principle, though disease recovery is largely a unidirectional process, a molecular based trajectory describing this process may also capture 'molecular worsening' which may translate into different clinical phenotypes (ICU discharge versus fatal outcome). As we assumed that most patients in our study display recovery trajectories, and hence generally underwent a unidirectional transition from severe to recovered phenotypic states, and due to the observation that deceased patients tend to display lower pseudotime over time (Figure S3A) , we hypothesized that 'molecular worsening' over time can be captured by altering our viewpoint on the trajectory. To that end, we trained a second TimeAx model based on reverse ordering of the sampling time points of the original 20 patients of cohort 1 used for training of the recovery model. Next, we calculated a score for each cohort 1 patient for both the recovery model and the reverse-time generated 'worsening model' (model 2; Figure 3A ; see Methods for score). Contrasting the scores from the two TimeAx models yielded three major patient groups: patients with a high recovery score and a low worsening score (denoted as 'recovering'), patients with a high worsening score and a low recovery score (denoted as 'worsening') and patients with low recovery and worsening scores (denoted as 'stable'). Choosing a cutoff of 0.1 for maximizing the number of worsening and recovering patients, while keeping the number of out-group patients low ( Figure S3B) , resulted in nine patients which were classified as 'worsening' (17%), and 41.5% either as 'stable' or 'recovering' (Figure 3B ). Two patients that presented high inconsistencies in pseudotime across the two models, due to a sharp recovery after a worsening period, were omitted from the analysis ( Figure 3B and Figure S3C ). The 'worsening' group had the highest percentage of individuals with age above 70 (p=0.03; Figure 3C , left), slightly more males ( Figure  3C , middle), and higher BMI (Figure 3C, right) . Although not significant, we also observed more patients with comorbidities associated with unfavorable outcomes in the 'worsening' group, including diabetes, cardiovascular diseases and hypertension 25 ( Figure 3D) . Indeed, the mortality rate (56%) in the 'worsening' group was higher than in the 'stable' (18%) and 'recovering' (8%) groups, respectively (p=0.07; Figure 3E ).

We next tested whether pseudotime changes might be predictive markers of disease trajectories at early time points after ICU admission. By evaluating the difference between the first two pseudotime values in each patient, we observed a significant difference, including an increase and a decrease in recovering and worsening patients, respectively (p=0.003; Figure 3F ). We successfully validated the association between pseudotime dynamics and patient survival by calculating the difference in pseudotime between two consecutive time points in patients from cohort 3. We defined patients in cohort 3 as 'recovery-like' if their second pseudotime was higher than their first, or 'worsening-like' if we observed a decrease in their second pseudotime (see Methods). We found that the group of worsening-like patients had a higher mortality rate compared to the group of recovery-like patients ( Figure S3D) . Specifically, we show that a strong increase in pseudotime is related to recovery, while higher mortality rates are observed in patients with a less prominent increase or a decrease in pseudotime ( Figure S3E ).

We looked for potential pseudotime parameters that might explain the unfavorable outcomes within the 'worsening' group as opposed to the recovery group. Comparing the minimal and maximal pseudotime values of patients between the two groups yielded no significant differences, arguing against higher overall disease severity in the worsening group (Figure S3F -G). However, we observed a strong distinction between the two groups in the chronological time from ICU admission matching the lowest and highest pseudotime levels in each patient ( Figure  3G) , suggesting that not the pseudotime values themselves but their change over time can be utilized to predict patient trajectories. Specifically, we observed a significant difference with patients in the 'worsening' group showing minimum severity at earlier time points after ICU admission and maximum severity at later time points, while this was not the case for recovering patients ( Figure 3G and Figure S3H ).

Finally, we explored associations between pseudotime at the last sampled time point and selfreported outcomes three months after ICU admission. This analysis was performed on a subset of 24 patients from cohort 1 from whom these data were available (length of ICU stay of these patients ranged from 10 to 45 days; See Methods). Patients with high disease severity (i.e. lower pseudotime states) in their last sampled time point exhibited higher frailty scores at three months post-ICU and were more often categorized as being frail (r=-0.47, p=0.01 and p=0.006, respectively; Figure S3I ). Furthermore, lower pseudotime was associated with lower mental and physical component quality of life summary scores three months post-ICU (r=0.43, p=0.03 and r=0.45, p=0.03, respectively; Figure S3J -K).

Collectively, we extended the recovery model to also include disease trajectories that would be associated with unfavorable outcomes and determined that low level pseudotime at later time points after ICU admission are linked to this additional trajectory as well as to long-term outcomes of COVID-19.

Altered neutrophil dynamics distinguish between recovery and fatal outcome Next, we asked whether the pseudotime model would infer generalizable rules for immune cell composition changes during the recovery process directly from blood transcriptomes. As a clinical reference, we first overlayed routinely measured white blood cell counts onto the pseudotime model. Lowest pseudotime (reminiscent of most severe disease states) in cohort 1 was associated with lymphocytopenia, previously described as a key hallmark of severe COVID-19 ( Figure 4A ) 26, 27 . Indeed, a gradual increase in lymphocyte abundance across the pseudotime was also observed using flow cytometry data from cohort 3 ( Figure S4A ). Next, we assessed associations between pseudotime and the neutrophil-to-lymphocyte ratio, another welldescribed characteristic for COVID-19 disease severity 28, 29 . We observed a decrease in this ratio along the recovery pseudotime, which was linear, starting with high ratios in low states (low pseudotime) leading to low ratios in high states (high pseudotime), with a slight sharper decrease in low pseudotime (most severe states) ( Figure 4B ). With this clinical data framework in mind, we next deconvoluted the blood transcriptome data using CPM, a cellular deconvolution method that relies on single cell RNA-seq data to robustly predict cell state compositions in bulk RNA-seq samples 30 , using whole blood single cell RNA-seq data collected from COVID-19 patients as a reference for 29 different immune cell subsets 2 (Methods).

We calculated the associations between the deconvoluted data with either the pseudotime or chronological time and revealed strong cell distribution changes across the pseudotime, which were not observed when using chronological time, further supporting the superior performance of the pseudotime model when describing recovery (20 of 31 cell types with p<0.05, Fisher's ztest; Figure S4B ). We plotted dynamics of lymphocyte deconvolved abundance over time ( Figure  S4C ), which illustrated that the same kinetics as determined by blood count analysis ( Figure 4A ) was captured in the blood transcriptomes. Similarly, the dynamics of the neutrophil-tolymphocyte ratio was also captured in the RNA-seq data following the pseudotime model ( Figure  S4D ), indicating that the combination of transcriptomes and pseudotime modeling can capture cell abundance dynamics during COVID-19 recovery. Signatures of different types of T cells including CD4+ T effector memory (TEM) and CD8+ central memory T (TCM) cells, but also B cells were correlated with recovery (high pseudotime values), while different subsets of mature neutrophils (MME+, IFIT+, CXCL8+ neutrophils), eosinophils and CD14+ monocytes were the cell types mainly linked to low pseudotime, reflecting the most severe disease states ( Figure 4C and Table S1), a finding that was corroborated by previously published single cell RNA-seq data 2 . The deconvoluted transcriptome-based pseudotime model data was significantly correlated with the flow cytometry data derived from a subset of patients and time points ( Figure S4E ; see Methods) further illustrating how pseudotime models applied to bulk blood RNA-seq data can be used to uncover cell subset dynamics over time in the recovery setting of COVID-19.

To determine if changes in cell composition are directly linked to outcome, we focused on the main differences between deceased and survived patients. Therefore, we extended our pseudotime model by including data from deceased patients in addition to all recovered patients (7 and 20 patients, respectively, model 3; Figure 4D ; See Methods). Since the inclusion of these data from diseased patients disrupts the recovery assumption, we termed it the 'disrupted' model. We observed a major difference compared to the original model, as the predicted pseudotime based on the 'disrupted' model, only partially overlapped with the recovery pseudotime (r = 0.5; Figure S4F ). Translating these differences into the cellular space, a group of cell subsets correlated differently with pseudotime within the recovery and the 'disrupted' models ( Figure 4E) . Specifically, mature neutrophils, which were strongly inversely-correlated with the recovery pseudotime, had much weaker correlations with the disrupted pseudotime (p<10 -5 ; Fisher's z-test). On the other hand, different lymphocytes, including a variety of T cells and B cell subsets, exhibited no significant change in their correlations with the recovery and the disrupted pseudotime, suggesting that these cell types are not strongly associated with unfavorable patient outcomes ( Figure 4E ). One hypothesis is that massive increase in mature neutrophil abundance in severe states (low pseudotime) is crucial for proper recovery, while in deceased patients, there seems to be a delay in the maturation of neutrophils, resulting in an unfavorable outcome 31, 32 . Indeed, we found that patients from the recovery group displayed higher levels of mature neutrophils compared to patients in the worsening group early on during their ICU admission (p=0.007; Figure S4G ). On the other hand, much weaker associations were found for lymphocytes and for the neutrophil-to-lymphocyte ratio (p=0.05 and p=0.02, respectively; Figure S4H ), suggesting that neutrophil abundance alone might be a better clinical predictor of patients' outcomes compared to the neutrophil-to-lymphocyte ratio. Moreover, we found even stronger discrimination between recovering and worsening patients, focusing at the change in abundance of mature neutrophils between the first and second time points (Figure 4F ). While the recovery group presented a decrease in the abundance of mature neutrophils over time, patients from the worsening group presented a significant increase in their abundance (p=2.2·10 -6 ; Figure 4F ), supporting the notion that there is a deviation in neutrophil dynamics. The discrimination between the groups was much weaker based on either the lymphocyte abundance or the neutrophil/lymphocyte ratio (p>0.05 and p=0.003, respectively; Figure 4G ). We also observed a superiority of changes over time in mature neutrophils over the neutrophil-tolymphocyte ratio, comparing the two time points in cohort 3 (p=9·10 -5 and p=0.0001, respectively; Figure S4I ), further strengthening that mature neutrophils have more value for patient outcome prediction.

Overall, we show that patients in severe states (low pseudotime) contain high abundance of mature neutrophils, which are reduced during the patients' recovery process. Dysregulation of this process, causing an increase in mature neutrophil quantities over time, is highly predictive for unfavorable disease outcomes and is superior to the neutrophil/lymphocyte ratio, currently proposed as a suitable biomarker for COVID-19 severity.

Downregulation of molecular activation pathways is associated with recovery progress and different outcomes, unrelated to changes in cell compositions

As the recovery pseudotime model (model 1) is molecularly driven, we reasoned that it may be used to identify underlying molecular phenotypes of recovery from severe COVID-19. We focused on differential gene expression patterns across the recovery process. Since a large part of the variance in gene expression comes from changes in cell compositions 33 , we decoupled cell type composition effects (denoted as 'compositional effect') from those that can not be explained by compositional changes (denoted as 'non-compositional effect'). We estimated the detected differential expression that is due to differences in cell abundance between samples (from hereon 'compositional-based data') by combining all cell subsets' gene expression profiles, weighted by their deconvolved composition within each sample. We contrasted compositionbased data with the total gene expression data (denoted as 'global data') ( Figure 5A , see Methods). We defined a gene to be significantly associated with recovery if it presented a significant p-value score (student's t-test) comparing its expression between low and high pseudotime samples (Methods). We observed a large difference in associations of the recovery pseudotime with the global data levels compared to the composition-based data ( Figure S5A ).

Overall, 21% (n=2262) and 17% (n=1786) of the genes showed either no significant effects or significant effects respectively in both the global and composition-based data ( Figure 5B) . Conversely, 54% (n=5747) of genes had solely compositional effects with no significant associations with the pseudotime using the global data, while an additional 8% (n=880) of the genes presented associations with the pseudotime that are not related to cell composition ( Figure 5B ). In addition, we discovered a large group of genes with positive compositional effects but negative associations with global data levels ( Figure S5B,C) . Overall, these results suggest that aside from molecular changes due to compositional effect, more general regulatory changes occur along the pseudotime across all cell types ( Figure 5B ).

Based on these results, we were able to study these global regulatory mechanisms of COVID-19 recovery, by focusing on genes with strong non-compositional effects and either non or opposite compositional effects in cohorts 1, 2 and 3 ( Figure 5C and Table S2 ). Among the genes with opposite associations between the recovery pseudotime and either the global data (negative association) and composition-based data (positive association) is the cell surface receptor ITGAL, which is involved in T cell regulation and indicative of activated memory/effector T cells 34 ( Figure  S5D ). While the increase in ITGAL over the pseudotime, based on the composition-based data, might be related to the increase in T cell abundance during the recovery process ( Figure 4C) , it is decoupled from global regulatory effects on T cell activation, resulting in a peak during severe states and decreases during recovery. Other genes associated with immunomodulatory functions such as GSDMD ( Figure S5D ) 35 and IL-16 ( Figure S5D ) 36 presented similar trends, further supporting the notion that T cell activation is decoupled from T cell abundance during the recovery process.

To determine major functional changes, we calculated enrichment scores based on gene sets taken from GO, MSigDB Hallmark gene sets, Reactome, and KEGG. Here, a high enrichment score accounts for enrichment of gene sets with strong associations with the recovery pseudotime (See Methods). Gene set enrichment scores highly overlapped between global and compositionbased data, suggesting that most gene expression changes along the pseudotime, including immune activation patterns as well as most of the cell metabolism, can be explained solely by cell type compositional changes ( Figure S5F ). One exception is heme metabolism, which presented high enrichment based on the global data, while only marginal enrichment was detected in the composition-based data ( Figure 5D) . Indeed, heme metabolism was shown before to be related with late COVID-19 12, 37 . Our results strongly support the importance of this program during COVID-19 recovery and suggest that regulatory changes and not cell compositional changes govern this effect. Additional pathways that were enriched in the global data but not in the composition-based data included programs linked to homeostasis and cell cycle, such as intracellular transport, centrosome and microtubule cytoskeleton, hinting at processes of cellular normalization along the pseudotime. On the other hand, the regulatory mechanisms of catabolic processes, apoptosis, and P53 activation were enriched, pointing towards a reduction of cell stress as a major component of COVID-19 recovery (Figure 5D , Figure  S5F and Table S3 ). These effects were normalized along the pseudotime, providing another indication for homeostasis in recovered patients ( Figure 5C and Figure S5E ). Alterations of NFkappa B and interferon signaling have been linked with acute COVID-19 [38] [39] [40] . A strong negative correlation of interferon activation with the pseudotime in the composition-based data, normalization during recovery, and weak association in the global data indicate that interferon activation is indeed impaired in severe patients ( Figure 5E ). Additionally, we observed strong negative correlations, based on the global data, with recovery pseudotime for four genes of the mammalian NF-κB family: NFKB1, NFKB2, REL, RELA, which was not explained by cell composition (Figure 5F ), suggesting that there is a change in balance between the interferon and NF-κB responses during the course of severe COVID-19, characterized by impaired interferon responses but enhanced NF-κB responses early on, which reverses during the recovery process.

Finally, we intended to define novel whole blood gene markers that relate to different disease outcomes. Therefore, we compared the correlations within the global data with the pseudotime of the recovery model or the disrupted model, respectively (model 3; Figure 4D ; See Methods). We detected 48 genes that were associated with the recovery but not the disrupted pseudotime. These were associated with innate immunity and especially TNF and interferon response pathways, including PTPRE, CEBPB, ZEP36, BCL3, CFLAR, MXD1, LITAF, FOSL2, IFITM2, IL4R and TRIM25 (Figure 5F , Figure S5G and Table S4 ). In addition, we identified genes associated with the regulation of T-cell activation and differentiation, such as THEMIS2 41 , ARBB2 42 , SEMA4A 43 , and IL-4 44 , emphasizing again the importance of this mechanism for the recovery process ( Figure 5F and Table S4 ). Conversely, we identified 12 gene markers (GPI, LAIR1, RFLNB, ABCA13, CEBPE, LCN2, BPI, LTF, POLE, ITGA9, CEACAM8, CEACAM6), which presented strong negative correlations with the disrupted, but not the recovery pseudotime (Figure 5F, Figure S5G and Table S4 ). CEACAM8 and CEACAM6 were previously pointed out as susceptibility and mortality markers of COVID-19 and linked to affecting the cross-talk between neutrophils and pneumocytes in patients' lungs 45 . These 12 gene markers may serve as blood-based severity markers for COVID-19 which are associated with patients' outcomes.

Collectively, we identified molecular mechanisms that govern COVID-19 recovery over time such as T-cell activation which was decoupled from changes in T cell composition, directly related with patients' outcomes. At the same time, recovery was accompanied by reduction in neutrophils, changes in the regulation of catabolic processes and apoptosis as well as regulatory processes within the interferon and NFkB response pathways.

The systematic identification of overarching molecular programs and cellular changes during the recovery from heterogeneous diseases is a major challenge in medicine and has been prominently exemplified by the difficulty to understand recovery processes in COVID-19. By focusing on the recovery phase in patients with severe acute COVID-19 admitted to the ICU, assuming that patients released from the ICU present a recovery trajectory, we provide strong evidence that whole blood transcriptomes combined with pseudotime modeling, utilized by the TimeAx method, cell type deconvolution, molecular pathway prediction and validation in additional cohorts, is well-suited to identify global principles of recovery kinetics compared to non-recovery despite the fact that individual disease trajectories differ dramatically over time. Specifically for severe COVID-19, we presented a gradual decrease in mature neutrophils during recovery as the major cellular determinant, which was also reflected in superior outcome prediction. Further, we identified numerous molecular changes ranging from normalization of catabolic processes to changes in the balance between NFkB-mediated and interferon driven gene programs during the recovery process.

We also assessed whether TimeAx-based pseudotime models could reveal molecular and cellular changes accounting for unfavorable outcomes rather than recovery. To do so, we generated a second model, the worsening model based on TimeAx. Here we reversed the order of sampling timepoints for each patient within our longitudinal cohort (cohort 1). We defined patients as 'recovering' or 'worsening', providing a clear association between a patient's temporal change in pseudotime and the clinical outcome of the disease. Specifically, increases in pseudotime values were directly associated with recovery, while reduced or low pseudotime values observed at later real time points correlated with fatal outcomes as well as with long-term outcomes in surviving patients, who reported to be more frail and experienced a lower quality of life. Importantly, our findings that the dynamics of changes between the two first measurements after ICU admission seem to be predictive for a recovering or worsening phenotype opens a window for future outcome prediction tools relying on patient pseudotime patterns. Furthermore, monitoring of pseudotime values in a later stage of admission could have predictive value for long-term outcomes and facilitate earlier and/or more aggressive interventions, such as rehabilitation programs.

Combining the TimeAx-based pseudotime model with cellular deconvolution of cell types within whole blood transcriptomes allowed us to define the general cellular changes occurring during recovery, but also those related to unfavorable outcomes. Low pseudotime was characterized by depletion of lymphocytes, as well as the increase in neutrophil/lymphocyte ratio, which has been previously shown to be associated with severe acute COVID-19 when using other methods (e.g. flow cytometry) 29 . These findings clearly indicate that blood transcriptomes combined with the TimeAx model and cellular deconvolution analysis reveal such findings without any additional measurements. However, in addition to these findings in the early and acute phase of the disease, our longitudinal study allowed us to generalize the major changes from the early severe disease stages towards the recovery trajectory. Here, we observed a continuous normalization in predicted frequencies of the cell types associated with severe acute COVID-19, discovering associations between the dynamics of cell compositions across the pseudotime and patient outcomes. Specifically, while gradual decrease in mature neutrophils was strongly related to recovery, higher mature neutrophil abundances at late time points were a good indicator for unfavorable disease outcomes, presenting higher predictive power compared to both lymphocyte counts and the neutrophil/lymphocyte ratio. Overall, our model strongly suggests that the course of mature neutrophils is more predictive for outcome than currently used parameters, in particular the neutrophil/lymphocyte ratio 29 .

Finally, we studied the molecular mechanisms that govern COVID-19 recovery over time. Importantly, by comparing the effects coming from cell composition changes only to the effects being independent of cell composition changes (global data), we identified molecular mechanisms changed irrespective of cellular changes. Specifically, we discovered a regulatory change in gene programs across the pseudotime that was associated with restoration of cell homeostasis. In addition, while T cell abundance increased over the pseudotime, we discovered initially elevated T cell activation in severe patients, which normalized as patients advanced in their pseudotime. This effect was disrupted by adding deceased patients into the model, suggesting the importance of this mechanism on determining patients' outcomes. Still, further experimental validation is required to fully support this new computationally-derived disease progression mechanism. Our analysis also revealed a normalization of the balance between the NFkB and the interferon system during the recovery process, which started out as impairment of the interferon system and enhanced NFkB activity in acute severe COVID-19, which may signify a maladaptive hyperinflammatory but poor antiviral state [38] [39] [40] .

By combining longitudinally sampled blood transcriptomes, TimeAx-based computational modeling generating pseudotime trajectories, overlaid cell deconvolution analysis and molecular pathway prediction, we provide a unique framework for the description of the recovery process including cellular and molecular changes, here exemplified for patients with severe COVID-19. Further, the pseudotime findings can be linked to clinical metadata and both short-term and longterm outcomes. As more data for other patient cohorts become available, it will be interesting to determine whether patients with other disease courses are characterized by similar cellular and molecular changes or whether additional recovery trajectories exist. Further, it will be interesting to study whether patients who develop persistent symptoms following mild acute COVID-19 are also distinct, when applying our analytical framework. Similarly, larger clinical validation trials for the genes and/or signatures, identified as related to the recovery process, as well as to unfavorable outcomes, are still required. Lastly, since many diseases are characterized by heterogeneous and patient-specific recovery trajectories, our model might also uncover a unifying recovery trajectory for these diseases.

As we have clear indications for general recovery mechanisms in COVID-19, experimental validation studies of these mechanisms will be needed to further support generalization of these processes. While we already included several cohorts, to better understand the trajectory towards fetal outcome, cohorts with a larger number of patients that died from the disease might uncover further mechanisms and biomarkers for this group of patients. The model underlying TimeAx currently presumes a unifying disease recovery dynamics. If there would be more complex trajectories with patient subsets including branching into different cellular and molecular trajectories, we would have to extend our current model to reflect such complexity. However, as for now, we have no indication for multiple trajectories of recovery.

FJT reports receiving consulting fees from ImmunAI and ownership interest in Dermagnostix. SSO holds equity and is a consultant of CytoReason.

We worked to ensure that the study questionnaires were prepared in an inclusive way. The author list of this paper includes contributors from the location where the research was conducted who participated in the data collection, design, analysis, and/or interpretation of the work. ratio (G, right) , between the first and second time point for each patient, for either the recovering or the worsening patient groups (x-axis) (n=55). In A, F and G, boxes represent the 25th, 50th, and 75th percentiles; whiskers show maxima and minima. 

Further information and requests for resources should be directed to and will be fulfilled by the lead contact, Prof. Joachim L. Schultze (Joachim.Schultze@dzne.de).

This study did not generate new unique reagents.

• Whole blood RNA-seq data for cohort 1 are available at European Genome-phenome Archive (EGA, ID:EGAS00001005735). • All original code has been deposited (https://github.com/amitfrish/COVIDRecoveryCode) and is publicly available as of the date of publication. • Additional Supplemental Items are available from Mendeley Data at http://dx.doi.org/10.17632/jsn2kpnmpg.1. • Any additional information required to reanalyze the data reported in this paper is available from the lead contact upon request.

Radboud university medical centre (Nijmegen, the Netherlands) cohort (Cohort 1): COVID-19 was diagnosed by a positive SARS-CoV-2 RT-PCR test in nasopharyngeal and throat swabs and/or by typical chest CT-scan findings. Patients with a pre-existing immunosuppressed status or other comorbidities that strongly influence prognosis were excluded. All blood samples were obtained from an arterial cannula already in place, so no venipunctures were required. Data were collected from the electronic patient files (EPIC, EPIC Systems Corporation, Verona, Wisconsin, USA) and recorded in the good clinical practice (GCP)-compliant data management system Castor (Castor EDC, Amsterdam, the Netherlands). The study was carried out in the Netherlands in accordance with the applicable rules concerning the review of research ethics committees and informed consent. All patients or legal representatives were informed about the study details and could decline to participate.

For 24 patients of cohort 1, self-reported outcome data at three months post-ICU admission were available. These data were collected in the MONITOR-IC study, an exploratory multicenter prospective cohort study in ICU survivors. The study protocol has been previously described 46 and approved by the local research ethics committee (CMO region Arnhem-Nijmegen, the Netherlands, No. 2016-2724). For the purpose of the current study, frailty and quality of life (QoL) outcomes were used. Frailty was measured using the Clinical Frailty Scale (CFS) 47, 48 , consisting of one item with a score ranging from 1 (very fit) to 9 (terminally ill). Patients were classified as "non-frail" (score of 1-4) or "frail" (5-9) 49 . QoL was assessed using the 12-Item Short Form Survey (SF-12) 50 , where scores were aggregated into two summary measures: a Physical Component Summary (PCS) and a Mental Component Summary (MCS). Scores of the PCS and MCS range from 0 to 100 and higher scores indicating a better QoL.

Cohort 2: public data of 39 severe COVID-19 patients admitted to the ICUs of the University Hospital RWTH in Aachen, Germany, Kiel University and University Medical Center, Germany, Saarland University Medical Center in Saarbrücken, Germany, and the National and Kapodistrian University of Athens, Greece, which were sampled once for whole blood transcriptome. Further information about this cohort can be found in 11 (Dataset E).

Cohort 3: public data of 45 severe COVID-19 patients admitted to the ICU of either the Addenbrooke's or the Royal Papworth hospitals in Cambridge, UK, which were sampled once or twice for whole blood transcriptome, during their stay in the ICU. For some of the samples, also paired flow cytometry data and plasma cytokines were measured. Further information about this cohort can be found in 12 .

J o u r n a l P r e -p r o o f RNA isolation, library construction, pre-processing, and read alignment Patient blood from cohort 1 was sampled into PAXgene Blood RNA Tubes (Qiagen) and RNA extraction was performed following manufacturer's instructions using either the PAXgene Blood RNA Kit or the QIAsymphony PAXgene Blood RNA Kit on a QIAsymphony SP automated liquid handling system (Qiagen). In both cases RNA integrity and quantity was determined via RNA assay on a Tapestation 4200 system (Agilent). 750ng total RNA were subsequently used as an input for NGS library prep according to the TruSeq Stranded Total RNA with Ribo-Zero Globin kit (Illumina) protocol using IDT for Illumina TruSeq UDI adapters (Illumina). Libraries were quantified using the Qubit HS dsDNA assay (Thermofisher) and library fragment size distribution was determined using the D1000 assay on a tapestation 4200 system (Agilent). Libraries were equimolarly pooled into pools of 96 samples and clustered lane wise at 200pM concentration on NovaSeq6000 S2 flow cells. Raw sequencing data was demultiplexed and processed into fastq files using bcl2fastq2 v2.20 and subsequently aligned to human genome build GRCh38 using STAR aligner. Reads were sequenced paired-ended for 50 bp each with an average depth of 2.4 10 7 reads per sample, with a mean alignment rate (unique and multiple mappings) of 89.4%.

In cohort 1, Ethylenediaminetetraacetic acid (EDTA)-anticoagulated blood was centrifuged (2000g, 10 min, 4 oC), after which plasma was stored at -80 o C until analysis. Concentrations of interleukin (IL)-6, and Granulocyte colony-stimulating factor (GCSF) were determined in one batch using a Luminex assay (Milliplex, Millipore, Billerica, USA). The lower detection limit was 3.2 pg/mL for both cytokines. In cohort 3, concentration of interleukin (IL)-6 was determined as well (see information in 12 ).

Whole blood cell counts in cohort 1 were obtained using a Coulter Ac-T Diff cell counter (Beckman Coulter, Inc., Brea, CA, USA) that was calibrated daily. 1 ~ 1.5 mL of whole blood was incubated in lysis buffer to lyse red blood cells. Remaining white blood cells were washed twice with PBS and resuspended in PBS + 0. All reagents were titrated and tested before they were used in the current study. Stained cells were measured on a 10-color Navios flow cytometer (Beckman Coulter, Inc., Brea, CA, USA) equipped with three solid-state lasers (488 nm, 638 nm, and 405 nm). Flow cytometry data were analyzed using Kaluza Analysis Software version 2.1 (Beckman Coulter, Inc., Brea, CA, USA). More information can be found in 51 .

Calculating recovery pseudotime using the TimeAx algorithm TimeAx 10 is a computational method that was introduced to align patient trajectories during disease progression. TimeAx relies on the order (or rank) of all samples of a dataset, following the similarities within the data over time and not the chronological time itself. Based on shared patterns within the complete data space across time and patients, and utilizing a feature selection pipeline to determine a seed of genes (here, using TimeAx default settings of 50 genes), it creates a longitudinal model consisting of a set of consensus trajectories, representing numerical measures of all possible disease states. TimeAx then allows the matching of new patient samples with the consensus trajectories, discovering their exact progression state (denoted as 'pseudotime'). Thus, pseudotime is a numerical continuous axis representing the biological dynamics across time, considering the molecular and physiological heterogeneity across patients. To make sure that the inferred pseudotime is indeed capturing patients' continuous dynamics over time, TimeAx produces a robustness score, measuring an additional pseudotime-like axis (robustness pseudotime), assuming that pseudotime can only be the same or higher for consecutive samples within the same individual. The robustness score is then calculated as the similarity between the pseudotime and the robustness pseudotime. Here, TimeAx and all its support functions were applied using the TimeAx R package. Recovery model (model 1): We used our longitudinal cohort (cohort 1), a time series RNA-seq data, to predict a specific pseudotime state for each single RNA-seq profile across multiple sampling time points in each patient. Based on the TimeAx's model, we were able to also predict the pseudotime states of additional RNA-seq samples from two validation cohorts (cohort 2 and 3). Based on the assumption that non-deceased patients, which eventually were released from the ICU, presented a general trajectory of recovery, we trained this model using all non-deceased patients with 5 or more samples each (total of 20 patients). Based on this model, we also predicted the pseudotime of samples from deceased patients. We show that pseudotime states based on this model, presented only a minor association with the patients' time from disease onset (Figure S2C) , and thus allowed us to discover many novel findings that could not be detected otherwise.

Worsening model (model 2): In this model, we trained TimeAx using the exact same patient selection as in the recovery model. By assuming that non-deceased patients in cohort 1 present recovery trajectories, we reversed the order of samples for each patient, resulting in a model for patient worsening instead of recovery. Therefore, the increase in pseudotime states within this model represents a worse disease state. This allowed us the calculation of an optimized fit score for each patient to either of the recovery or worsening models. For that purpose, we assumed that as patients advance in time, they can only go forward in pseudotime, thus, can only advance or stay at the same state. Therefore, patients that had a recovery phenotype, advanced their pseudotime mostly in the recovery model, but not in the worsening model. On the other hand, patients that had a worsening phenotype, advanced based on the worsening model but not the recovery model. We calculated a recovery and worsening scores as the standard deviation across the predicted pseudotime states in each of the models. To maximize the ratio between patients in the recovery and worsening groups to the number of out-group patients (Figure S3B) , we set the cutoff for high recovery / worsening to be above 0.1 and divided the patients into 3 groups: patients with a high recovery score and a low worsening score (denoted as "recovering"), patients with a high worsening score and a low recovery score (denoted as "worsening") and patients with low recovery and worsening scores (denoted as "stability"). To validate our findings using patients from cohort 3, which only included 2 time points per patient, we had to calculate an alternative score. In this analysis, we defined patients as 'advancing' or 'returning' if the pseudotime of their second time point was higher or lower compared to the pseudotime of the first time point, respectively. Disruption model (model 3): In this model, we trained TimeAx using patients with 4 or more samples, including both recovered and deceased patients (20 and 7 patients, respectively). By adding the deceased patients into the model, we disrupted the recovery assumption, resulting in pseudotime that does not represent the normal recovery process. We used this model to discover standard recovery mechanisms that are altered in deceased patients.

TimeAx pseudotime was validated using MEFISTO 52 , a recently reported method for factor analysis, using real-time as a covariate, for extracting time-smoothed biological factors. The values of the factors across patients are then aligned based on a selected patient as the backbone, to reduce the effect of patient heterogeneity. MEFISTO was applied using 2000 most highly variable genes for prediction of a single factor. This was shown superior over choosing a higher number of factors, each with a weaker association with time and less molecular interoperability.

Gene set enrichment scores were calculated based on the correlations of genes with the pseudotime. Gene sets were collected from GO, MSigDB Hallmark gene sets, Reactom, and KEGG. For each gene set, we compared the correlations of genes from the set and genes which are not part of the set, using a Kolmogorov-Smirnov test (ks-test). Finally, we accounted for multiple testing by applying a Benjamini-Hochberg false discovery rate (FDR) correction. To provide better comparability between gene set enrichment scores, the corrected p-values were then -log10 transformed.

We used the Cellular population mapping (CPM) deconvolution algorithm 30 ('scBio' R package, version 0.1.6) to cell composition abundances across in both the longitudinal and the validation cohorts. For the creation of a reference data, we used whole blood single cell RNA-seq data collected from COVID-19 patients and containing 31 different cell subsets from all major immune cell types, including Neutrophils, Monocytes, Dendritic cells, Eosinophils, T cells and B cells, as well as blood circulating cells, such as Erythrocytes and platelets 2 .

In brief, we used the whole blood data from cohort 2 from Schulte-Schrepping et al. and reannotated the cells in a two-step approach. First, we log-normalized the gene expression with a scaling factor of 10.000, selected the top 2000 variable genes using the vst algorithm and scaled the top variable genes regressing out "nCount_RNA" (number of transcripts per cell). After a principal component analysis based on the scaled variable genes, we selected the first 20 components for a UMAP dimensionality reduction. Cells were then clustered using louvain clustering at a resolution of 0.7. The granulocyte space was annotated into eosinophils, proneutrophils, pre-neutrophils, CXCL8+ neutrophils, MME+ neutrophils and IFIT+ neutrophils based on the top marker genes. The remaining cells were used to annotate them using the azimuth algorithm 53 according to the tutorial to increase resolution of the cell annotation in the lymphoid space.

Overall, we created our reference data based on 1463 cells, evenly distributed across the different subsets (~50 cells for each subset). To capture the variation across all cell subsets, we created the cell space for each cell subset as the first two principal components (PCs), obtained by applying a PCA analysis, using only genes with high expression compared to cells from other cell subsets (Fold change > 2). Cell spaces for all cell subsets were then combined into one big cell space. To robustly capture changes across this combined cell space, we used a model size of 31, selecting one cell from each subset in each CPM iteration. We calculate the associations between the pseudotime and the deconvolved cell compositions using a 'Pearson' correlation across all samples.

Studying gene effects by comparing global data with composition-based data

To observe molecular mechanisms in COVID-19 recovery which are not introduced by changes in cell compositions, we calculated bulk gene expression profiles corresponding to the contribution of cell composition and compared it to the measured gene expression levels. We calculated the cellular contributions for each sample by, first, weighting the gene expression profile of each cell subset (determined from single cell RNA-seq data; See cell type deconvolution in Methods), by its composition, inferred by cell deconvolution, and then inferring the sample's compositionalbased profile by summing all these weighted cell subset profiles. To determine significant gene associations with recovery, for both global and composition-based data, we defined a significance score by comparing gene values between low pseudotime samples (pseudotime < 0.5) and high pseudotime samples (pseudotime > 0.5) using Student's t-test. To account for multiple testing, we selected genes with significant p-values after applying a Benjamini-Hochberg false discovery rate (FDR) correction. Tables   Table S1 . Cell composition correlations with pseudotime. Relates to Figure 4 . Pearson correlations of deconvolved cell types (column 1) compositions with pseudotime, in cohorts 1-3 (columns 2-4). J o u r n a l P r e -p r o o f

Elevated Calprotectin and Abnormal Myeloid Cell Subsets Discriminate Severe from Mild COVID-19

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

Disease severity-specific neutrophil signatures in blood

Single-cell multi-omics analysis of the immune response in COVID-19

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

Complement activation induces excessive T cell cytotoxicity in severe COVID-19

Immune cell profiling of COVID-19 patients in the recovery stage by single-cell sequencing

Longitudinal transcriptome analyses show robust T cell immunity during recovery from COVID-19

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

Multiple trajectory alignment reconstructs disease dynamics for discovery and clinical benefit

Swarm Learning for decentralized and confidential clinical machine learning

Longitudinal analysis reveals that delayed bystander CD8+ T cell activation and early immune pathology distinguish severe COVID-19 from mild disease

Time course of risk factors associated with mortality of 1260 critically ill patients with COVID-19 admitted to 24 Italian intensive care units

Body mass index and mortality in coronavirus disease 2019 and other diseases: A cohort study in 35,506 ICU patients

Elevated level of C-reactive protein may be an early marker to predict risk for severity of COVID-19

The chronic kidney disease and acute kidney injury involvement in COVID-19 pandemic: A systematic review and meta-analysis

Role of interleukin 6 as a predictive factor for a severe course of Covid-19: retrospective data analysis of patients from a long-term care facility during Covid-19 outbreak

A neutrophil activation signature predicts critical illness and mortality in COVID-19

Multi-cohort analysis of host immune response identifies conserved protective and detrimental modules associated with severity across viruses

Integrated, Multi-cohort Analysis Identifies Conserved Transcriptional Signatures across Multiple Respiratory Viruses

Weighted Gene Co-Expression Network Analysis Combined with Machine Learning Validation to Identify Key Modules and Hub Genes Associated with SARS-CoV-2 Infection

The Role of Th17 Response in COVID-19

COVID-19: a case for inhibiting IL-17?

Activation of the SARS-CoV-2 Receptor Ace2 through JAK/STAT-Dependent Enhancers during Pregnancy

Coronavirus disease 2019 and cardiovascular complications: focused clinical review

Lymphopenia as a Biological Predictor of Outcomes in COVID-19 Patients: A Nationwide Cohort Study

Lymphopenia predicted illness severity and recovery in patients with COVID-19: A single-center, retrospective study

Early prediction of in-hospital death of COVID-19 patients: a machine-learning model based on age, blood analyses, and chest x-ray score

Longitudinal hematologic and immunologic variations associated with the progression of COVID-19 patients in China

Cell composition analysis of bulk genomics using single-cell data

Neutrophils in COVID-19

Neutrophils and COVID-19: Active Participants and Rational Therapeutic Targets

Computational deconvolution: extracting cell typespecific information from heterogeneous samples

COVID-19-Induced ARDS Is Associated with Decreased Frequency of Activated Memory/Effector T Cells Expressing CD11a+

GSDMD is required for effector CD8+ T cell responses to lung cancer cells

lnterleukin-16: the ins and outs of regulating T-cell activation

SARS-CoV-2 Proteins Bind to Hemoglobin and Its Metabolites

Early IFN-α signatures and persistent dysfunction are distinguishing features of NK cells in severe COVID-19

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

Functional Exhaustion of Type I and II Interferons Production in Severe COVID-19 Patients

THEMIS enhances TCR signaling and enables positive selection by selective inhibition of the phosphatase SHP-1

Loss of β-arrestin 2 exacerbates experimental autoimmune encephalomyelitis with reduced number of Foxp3+ CD4+ regulatory T cells

Human Semaphorin-4A drives Th2 responses by binding to receptor ILT-4

The IL-4 receptor: signaling mechanisms and biologic functions

The predicting roles of carcinoembryonic antigen and its underlying mechanism in the progression of coronavirus disease 2019

MONITOR-IC study, a mixed methods prospective multicentre controlled cohort study assessing 5-year outcomes of ICU survivors and related healthcare costs: a study protocol

A global clinical measure of fitness and frailty in elderly people

Changes in frailty among ICU survivors and associated factors: Results of a one-year prospective cohort study using the Dutch Clinical Frailty Scale

The impact of frailty on intensive care unit outcomes: a systematic review and meta-analysis

The MOS 36-item short-form health survey (SF-36). I. Conceptual framework and item selection

Differential effects of environmental and genetic factors on T and B cell immune traits

Identifying temporal and spatial patterns of variation from multimodal data using MEFISTO

Integrated analysis of multimodal single-cell data

A quantitative assessment of severe COVID-19 patients' recovery dynamics

Gradual reduction in mature neutrophils is a cellular hallmark of COVID-19 recovery

Early rise in T cell activation is a functional hallmark of COVID-19 recovery

Anti-human CD4-Pe-Cy5.5 (13B8.2) Beckman Coulter Cat#B16491; RRID: NA Anti-human CD56-APC (N901)

Anti-human CD8-APC-AF700 (B9.11)

Anti-human CD19-APC-AF750 (J3-119)

Anti-human CD25-Pe-Cy7 (M-A251) BD Biosciences RRID

Anti-human IgD-FITC (IADB6) Southern Biotech RRID

RRID: NA Anti-human CD3-ECD (UCHT1)

Coulter Cat#B21444; RRID: NA Anti-human CD38-PE-Cy7 (LS198-4-3) Beckman Coulter Cat#B49198; RRID: NA Anti-human CD24-APC (ALB9) Beckman Coulter Cat#A87785; RRID: NA Anti-human CD5-APC-AF700 (BL1a) Beckman Coulter Cat#A78835

We thank Y. Abraham for help with designing and creating figure illustrations. EK, NB, PP & MK were funded by a Radboudumc COVID-19 grant. This work was supported in part by the German Research Foundation (DFG) to JLS. under Germany's Excellence Strategy (DFG) -EXC2151 -390873048); the HGF grant sparse2big, the EU H2020 projects SYSCID (Grant Agreement No. 733100) and ImmunoSep (Grant Agreement No. 847422). JLS was further supported by the BMBF-funded excellence project Diet-Body-Brain (DietBB) (grant number 01EA1809A), iTREAT (FKZ: 01ZX1902A), and by NaFoUniMedCovid19" (FKZ: 01KX2021, project acronym "COVIM"). This study was funded in part by European Union's Horizon 2020 Research and Innovation