key: cord-0898464-jcf0045k
authors: Krämer, Benjamin; Knoll, Rainer; Bonaguro, Lorenzo; ToVinh, Michael; Raabe, Jan; Astaburuaga-García, Rosario; Schulte-Schrepping, Jonas; Kaiser, Kim Melanie; Rieke, Gereon J.; Bischoff, Jenny; Monin, Malte B.; Hoffmeister, Christoph; Schlabe, Stefan; De Domenico, Elena; Reusch, Nico; Händler, Kristian; Reynolds, Gary; Blüthgen, Nils; Hack, Gudrun; Finnemann, Claudia; Nischalke, Hans D.; Strassburg, Christian P.; Stephenson, Emily; Su, Yapeng; Gardner, Louis; Yuan, Dan; Chen, Daniel; Goldman, Jason; Rosenstiel, Philipp; Schmidt, Susanne V.; Latz, Eicke; Hrusovsky, Kevin; Ball, Andrew J.; Johnson, Joe M.; Koenig, Paul-Albert; Schmidt, Florian I.; Haniffa, Muzlifah; Heath, James R.; Kümmerer, Beate M.; Keitel, Verena; Jensen, Björn; Stubbemann, Paula; Kurth, Florian; Sander, Leif E.; Sawitzki, Birgit; Aschenbrenner, Anna C.; Schultze, Joachim L.; Nattermann, Jacob
title: Early IFN-α signatures and persistent dysfunction are distinguishing features of NK cells in severe COVID-19
date: 2021-09-04
journal: Immunity
DOI: 10.1016/j.immuni.2021.09.002
sha: b6d02cdf33c062e3114bd88ca9aff671c4a58af8
doc_id: 898464
cord_uid: jcf0045k

Longitudinal analyses of the innate immune system including earliest time points are essential to understand the immunopathogenesis and clinical course of COVID-19. Here, we performed a detailed characterization of natural killer cells in 205 patients (403 samples, day 2-41 after symptom onset) from four independent cohorts using single-cell transcriptomics and proteomics together with functional studies. We found elevated IFN-α plasma levels in early severe COVD-19 alongside increased NK cell expression of ISGs and genes involved in IFN-α signaling, while upregulation of TNF-induced genes was observed in moderate disease. NK cells exert anti-SARS-CoV-2 activity but are functionally impaired in severe COVID-19. Further, NK cell dysfunction may be relevant for development of fibrotic lung disease in severe COVID-19, as NK cells exhibited impaired anti-fibrotic activity. Our study indicates preferential IFN-α and TNF responses in severe and moderate COVID-19, respectively, and associates prolonged IFN-α-induced NK cell response with poorer disease outcome.

The clinical presentation of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) infection is highly variable, ranging from asymptomatic to severe courses of coronavirus disease (COVID-19) Wang et al., 2020a; Zhou et al., 2020) . Besides epidemiological factors and certain comorbidities (Bennett et al., 2021; Williamson et al., 2020) , an imbalanced immune response underlies the clinical manifestation of COVID-19.

Patients with severe disease, in particular, present with elevated blood plasma levels of numerous cytokines and chemokines Giamarellos-Bourboulis et al., 2020a) , as well as a dysregulated type I interferon (IFN) response (Blanco-Melo et al., 2020; Hadjadj et al., 2020; Yao et al., 2021) . Further characteristics of severe COVID-19 are high frequencies of circulating CD14 hi CD16 hi monocytes, decreased CD14 lo CD16 hi monocytes (Hadjadj et al., 2020; Schulte-Schrepping et al., 2020; Su et al., 2020) , proliferating, type I IFN-activated HLA lo suppressive monocytes and emergency granulopoiesis. Metabolically hyperactive plasmablasts, IFN-activated circulating megakaryocytes, and erythropoiesis are increased in critically ill patients (Bernardes et al., 2020; Stephenson et al., 2021) . T and B cell compartments are also altered in severe COVID-19 Braun et al., 2020; Grifoni et al., 2020; Ni et al., 2020; Schulien et al., 2021) .

Despite many studies on important aspects of the immunopathology of COVID-19, our understanding of this disease is still incomplete. For example, the role of natural killer (NK) cells, a heterogeneous family of innate immune cells, has not been sufficiently studied.

Although there is clear evidence for their role in acute viral infections (Björkström et al., 2010; Blom et al., 2016; Kokordelis et al., 2014) , data on NK cells in SARS-CoV-2 infection are sparse (Maucourant et al., 2020; Rajaram et al., 2020) . A COVID-19 vaccine study demonstrated anti-Spike-dependent NK cell response in vaccinated macaques , suggesting that NK cells exert functions against SARS-CoV-2-infected cells.

Accordingly, a potential therapeutic benefit of NK cells in COVID-19 is currently being investigated in clinical trials (NCT04797975, NCT04634370, NCT04280224). On the other hand, NK cells can potentially exacerbate the extent of lung injury in viral respiratory infections (Rajaram et al., 2020) .

COVID-19 has been associated with NK cell activation, increased frequency of CD57 + adaptive NK cells (Maucourant et al., 2020; Varchetta et al., 2020) , impaired cytolytic activity (Osman et al., 2020) , reduced peripheral NK cells (Giamarellos-Bourboulis et al., 2020b; Jiang et al., 2020; Wang et al., 2020b; Wilk et al., 2020) and increased intrapulmonary NK cell frequencies (Chua et al., 2020; Liao et al., 2020; Xu et al., 2020) .

Here, we performed a detailed longitudinal characterization of NK cells in patients of different severities by combining single-cell transcriptomics and proteomics in four independent cohorts with comprehensive functional analyses, including studying NK cell activity against SARS-CoV-2-infected cells.

We demonstrate NK cells in early severe COVID-19 to display signs of a strong IFN-α response with increased expression of interferon stimulated genes (ISGs) and genes related to IFN-α signalling, whereas in early moderate disease NK cells were characterized by a TNF imprint. This differential gene expression pattern was specific for the first week after onset of symptoms and also enabled to discriminate between patients with fatal outcomes of COVID-19 and those who finally recovered. Moreover, we demonstrate an impaired anti-SARS-CoV2 NK cell activity which was particularly prominent and prolonged in severe COVID-19.

In summary, our data link persistent NK cell dysfunction, induced by an exaggerated IFN-α response, with an unfavourable disease course and thereby support a role for NK cells in the immunopathogenesis of COVID-19.

To assess the impact of SARS-CoV-2 infection on the function and composition of the NK cell pool, we analyzed longitudinally collected peripheral blood samples in a multi-center setting. Samples collected in Bonn (cohort 1) were analyzed on the cellular level by multicolor flow cytometry and a micro-well based scRNA-seq approach, while plasma levels of soluble factors were studied using a bead-based digital and planar-array ultrasensitive immunoassays. For cohorts 2 and 3, collected in Berlin and Kiel, respectively, as well as for the fourth cohort, assembled from two large cohorts carried out in the UK (Stephenson et al., 2021) and the US (Su et al., 2020) , scRNA-seq was performed using a droplet-based platform. Cohort 2 samples were additionally analyzed by mass cytometry (CyTOF) ( Fig. 1A and 1B, Table S1 ). Patients treated with dexamethasone were excluded from the analyses to avoid immunotherapy-induced biases in the results. A total of 205 patients and 81 controls, including 8 donors with flu-like symptoms, were studied.

In both cohorts 1 and 2, COVID-19 was associated with a decreased absolute number of circulating NK cells ( Fig. 1C and 1D ). Longitudinal analysis in cohort 1 demonstrated a similar loss of NK cells in both moderate and severe disease in week 1. In contrast, from week 2 onwards, patients with moderate disease showed normalization of NK cell counts, while severe COVID-19 was characterized by persistent NK cell depletion (Fig. S1A) .

Percentages of CD56 bright and CD56 dim NK cells did not differ between study groups ( Fig.   S1B and S1C). Correlation analysis demonstrated numbers of total and CD56 dim NK cells negatively correlated with C-reactive protein (CRP), an acute-phase protein reflecting the intensity of inflammation (Fig. 1E) . COVID-19 was associated with increased NK cell expression of the apoptosis marker active Caspase-3 and CD95 (Fig. 1F) . Nucleocapsid protein induced active Caspase-3 expression and a dose-dependent increase in CD95 (Fig. 1G) . In summary, our findings indicated COVID-19 to significantly affect the NK cell compartment.

In order to gain a more detailed insight into COVID-19-induced alterations of the NK cell pool, we assessed transcriptional changes of NK cells in the blood by scRNA-seq analysis. In cohort 1, NK cell transcriptomes were extracted from COVID-19 PBMC scRNA-seq data (Schulte-Schrepping et al., 2020) derived from 64 samples from 17 COVID-19 patients (8 moderate, 9 severe) collected between day 2-25 after symptom onset, and 13 sex-and agematched controls. UMAP visualization of the NK cells in cohort 1 revealed transcriptional alterations in diseased NK cells. Density coloring stratified for cells from controls, moderate or severe COVID-19 patients showed differential two-dimensional distribution ( Fig. 2A) . To investigate these disease-relevant differences, differentially expressed genes (DEGs) were calculated for severity groups (Fig. 2B) . Hierarchical clustering of the DEGs revealed 5 different gene modules with specific patterns according to disease groups. Gene enrichment analysis of the severe COVID-19-related module 3 and the moderate COVID-19-related modules 4 and 5 revealed enrichment in the Hallmark 'IFN-α response' and 'TNF signaling', respectively (Table S3) , indicating these pathways are discriminators for severe and moderate COVID-19 NK cells on the transcriptional level.

To further explore the transcriptional heterogeneity within the NK cell compartment, we performed clustering analysis of the single-cell transcriptomes, identifying 6 distinct subtypes J o u r n a l P r e -p r o o f (Fig. 2C) . Comparison to previously published NK scRNA-seq signatures and cluster marker expression, revealed these 6 subtypes comprised inflamed CD56 dim (high IFN-related genes), proliferating CD56 dim (MKI67), cytokine CD56 dim (CCL4, CCL3, IFNG), HLA hi CD56 dim (HLA-DP, HLA-DR), CD56 dim (FCGR3A) and CD56 bright (NCAM1) NK cells (Fig. 2D, Table S2 ). NK cell transcriptomes from the other 3 cohorts (cohorts 2-4) (Bernardes et al., 2020; Schulte-Schrepping et al., 2020) , comprising 49 samples from 18 COVID-19 patients and 22 control donors, 20 samples from 10 COVID-19 patients and 5 control donors as well as 201 samples from 110 COVID-19 patients and 39 control donors, respectively ( Fig. 1A and 1B , Table S1 and S6), resulted in 3 validation datasets of 6,964, 15,369, and 97,764 single-cell NK transcriptomes, respectively. Peripheral NK cell subtypes identified in cohort 1 were similarly found in cohorts 2-4, validating the subtype annotation ( Fig. S2A and S2B ).

Next, we analyzed the distribution of NK cell subtypes across different disease severities ( Fig. 2E, S2C) . In severe COVID-19 patients, both inflamed and proliferating CD56 dim NK cells were strongly overrepresented compared to moderate COVID-19. The fraction of cytokine CD56 dim NK cells was enriched in samples derived from patients with moderate disease. All these subtypes were rather low in controls, emphasizing their strong disease association. CD56 dim NK cells represented the main NK cell population in blood from control donors (Fig. 2E, S2C) . Taken together, inflamed and proliferating CD56 dim NK were associated with severe and cytokine CD56 dim NK with moderate disease, respectively.

In parallel, we applied flow cytometry in cohort 1 to study the peripheral NK cell compartment based on protein markers (Fig. 2F , S2D-E, Table S2 ). Analysis of NKp80 and CD94 excluded contamination with ILC1s within the NK cell gate (Fig. S2H) . We identified inflamed CD56 dim , proliferating CD56 dim , cytokine CD56 dim , HLA hi CD56 dim , CD56 dim , and CD56 bright subpopulations analogous to the transcriptome analysis (Fig. S2E) . Proportions of inflamed CD56 dim , proliferating CD56 dim , and HLA hi CD56 dim subpopulations were increased in COVID-19 patients (Fig. 2G) . Inflamed CD56 dim , cytokine CD56 dim , HLA hi CD56 dim , and CD56 bright NK subsets were also identified by CyTOF in cohort 2 (Fig. S2G) , where inflamed CD56 dim and HLA hi CD56 dim NK cells were elevated in COVID-19 (Fig. 2I) . The cytokine CD56 dim subset was particularly increased in patients with flu-like symptoms. The CITE-seq data from cohort 4 identified 6 NK cell subtypes as seen by transcriptome-based analysis further corroborating subpopulation structure at protein level ( Fig. 2D and 2E) . Singlemarker analysis confirmed elevated activation (both cohorts: CD69, HLA-DR; cohort 1: CD38) and proliferation (cohort 2: KI-67) in severe patients. Furthermore, an increase of NK J o u r n a l P r e -p r o o f cell-specific receptor expression was detected for severe (cohort 1: NKG2C) and flu-like illness (cohort 2: CD226), respectively ( Fig. S2H: gating; Fig. S2I and S2J) .

Together, NK cells stratified by disease severity revealed marked differences between severe and moderate COVID-19 in regard to gene expression and composition of NK cell subtypes.

To also include the aspect of disease dynamics in our analysis (Fig. 3A) , peripheral NK cells of cohort 1 were stratified by week after disease onset and DEGs were calculated comparing cells from the respective severity groups. UMAP representation revealed prominent timedependent changes (Fig. 3B) . Calculated DEGs between conditions were grouped into 15 modules by hierarchical clustering (Fig. 3C ) and used as input for functional enrichment analysis, transcription factor (TF) and upstream ligand prediction (Fig. 3D, Table S4 ).

Modules 1 and 2 were highly expressed in moderate COVID-19 NK cells and enriched for the terms 'TNF signaling via NF-κB' and 'response to interferon-gamma', indicating antiviral activity based on IFN and TNF signaling. The modules included IRF1, IFITM3, CCL3 and CCL4, which are induced by type I IFNs and genes such as TNFAIP, NFKBIA and FOSL2 relevant for TNF signaling. TF prediction further underlined IFN-induced response with STAT1/2 and the TNF impact with RELA among the top predicted TFs. RelA is a component of NF-kB that drives various transcriptional programs after TNF stimulation (Liu et al., 2017) . Module 3 comprised 46 genes characterized by a strong expression in the second week of severe disease. Functional enrichment analysis assigned the terms 'E2F targets' and 'DNA replication' to this module, indicating enhanced proliferative capacity. TF prediction pointed to members of the E2F-family as key TFs, further emphasizing the proliferative functionality of these genes. HMGB2, a factor related to cell proliferation in cancer , was predicted as the top potential ligand. Module 4 was enriched in genes specifically upregulated in week 1 in severe COVID-19 NK cells. Functional analysis of the 121 module genes revealed implications in 'IFN-α response' and 'negative regulation of viral processes'.

Correspondingly, the module contained numerous IFN-related genes (MX1, ISG15, ISG20, IFIH1). The top predicted ligands being members of the IFN-α family and the predicted TF including IFN-induced factors (STAT1, STAT2, IRF9) underlined the inflammatory character of this module. These results indicated the relevance of IFN-α signaling for severe COVID-J o u r n a l P r e -p r o o f 19 NK cells in early disease. In conclusion, early severe COVID-19 is dominated by IFN-α signaling (module 4) while, in contrast to early moderate COVID-19, showing lower enrichment for the TNF signaling pathway (module 1 and 2).

To assess the implication of the 6 NK cell subtypes (Fig. 2C ) in disease severity-and timespecific DEG modules ( Fig. 3C and 3D) , gene-set enrichment analysis of each module for each subtype was performed (Fig. 3E) . As expected by functional enrichment and TF prediction, modules 1 and 2, specific for early moderate COVID-19 NK cells, enriched especially in the cytokine-producing CD56 dim and partly in the inflamed CD56 dim NK cells, while module 4, enriched in week 1 after symptom onset in severe COVID-19 NK cells, was dominated by inflamed CD56 dim NK cells, stressing the early differences in severe and moderate COVID-19 and further highlighting the importance of TNF for a milder disease course. The proliferating CD56 dim NK cells contributed exclusively to module 3 including proliferation-related genes upregulated in the second week of severe COVID-19 ( Fig. 3C and 3D). Visualization of the proportions of different subtypes over time, showed an enrichment of inflamed CD56 dim NK cells in week 1 declining until mid-week 2 (Fig. 3F) . In contrast to severe COVID-19, moderate disease was characterized by a continuous presence of cytokine-producing CD56 dim NK cells. In severe patients, a strong increase of proliferating CD56 dim NK cells was observed starting approximately 11 days after symptom onset until the end of week 3 (Fig. 3F) .

Next, we defined time-dependent and severity-specific alterations of the 6 NK cell subsets in cohort 1 (Fig. 3G, 3H and S3A) and cohort 2 (Fig. S3B) . The proportion of inflamed CD56 dim NK cells was slightly higher in early severe COVID-19 compared to patients with moderate disease in cohort 1 (Fig. 3H ) and cohort 2 (Fig. S3C) . Consistent with the scRNAseq data, proliferating CD56 dim NK cells increased from week 2 to 3 in severe COVID-19 in cohort 1. Severe COVID-19 was also associated with increased protein expression of the activation markers CD38, CD69, and HLA-DR, especially in week 1 ( Fig. S3A and S3D) . In cohort 2, frequency of CD56 dim HLA-DR hi did not differ between moderate and severe disease but was increased in COVID-19 compared to controls and flu-like illness, respectively ( Fig. 3C and 3D) . Finally, the proliferation marker KI-67 was increased on NK cells in severe COVID-19.

Together, the predominant expression of activation markers was observed both on RNA and protein levels in the early phase of the disease course in severe COVID-19 patients.

To address the type I IFN system in more detail, we extracted the genes from the hallmark term 'IFN-α response' and visualized those that were DEGs (Fig. 4A) . Both moderate and severe COVID-19 patients showed elevated expression of these type I IFN signature genes at disease onset, which subsided in week 1 in moderate and in week 3 in severe patients. Several type I IFN target genes showed differential regulation between moderate and severe COVID-19, for example IFITM1 and IFITM3 were mainly elevated in moderate disease while GBP4, SELL, PSME2, CASP1, or TXNIP were only increased in severe COVID-19 (Fig. 4A) . Even when using all 'IFN-α response' genes for signature enrichment analysis, this response was elevated early after infection and persisted into the second week in severe disease (Fig. 4B) .

Examination of the data of cohorts 2-4 corroborated these findings as the IFN-α response was also enriched in COVID-19 NK cells in weeks 1 and 2 with a stronger signal in severe cases.

Investigating plasma levels of IFNs and proinflammatory cytokines (Fig. 4C) , we observed increased plasma concentrations of IFN-α together with other proinflammatory cytokines (TNF, IL-6, IFN-γ) in week 1, especially in severe disease (Fig. 4C) . In contrast to proinflammatory cytokines, plasma levels of IFN-α dropped after week 1. Plasmacytoid dendritic cells (pDCs), a main producer of IFN-α, were reduced in both moderate and severe COVID-19 in week 1 ( Fig. S4A and S4B) , which is in line with recent findings (Kuri-Cervantes et al., 2020) and argued against pDCs being the major source for elevated circulating IFN-α during this time.

When correlating IFNs and proinflammatory cytokines (week 1) with clinically determined disease severity, only IFN-α correlated with both WHO classification and SOFA score (Fig.   4D ). Hence, we used severe COVID-19 samples from all weeks after symptom onset and showed that the IFN response signature is elevated in patients from cohort 1 and 2 who succumbed to infection? (Fig. 4E) , which might therefore contribute to a fatal disease course.

For studying the role of TNF, we extracted the genes from the hallmark term 'TNF signaling via NF-κB' and visualized DE genes in cohort 1 (Fig. 4F) . These genes showed a distinct distribution from the 'IFN-α response' hallmark with very strong signals in moderate compared to severe COVID-19, particularly in week 1, with a prolonged expression for most genes. A few genes included in the TNF signaling hallmark (AREG, IL7R, CEBPD) were only induced in severe COVID-19. Enrichment analysis using the complete hallmark for 'TNF signaling via NF-κB' demonstrated a strong enrichment in NK cells from moderate J o u r n a l P r e -p r o o f COVID-19 patients that subsided over time with no enrichment in severe COVID-19 NK cells (Fig. 4G) . In cohorts 2-4, the TNF signature was enriched in moderate patients for the earliest time points available in the cohorts. In contrast to the IFN-α response signatures, the TNF signature was most elevated in NK cells from discharged patients both in cohorts 1 and 2 when analyzing severe samples from all weeks after symptom onset (Fig. 4H) .

The lack of enrichment of TNF signature genes in severe COVID-19 was discordant with the high level of TNF in plasma in these patients (Fig. 4C, S4C) . The interplay between the TNF and type I IFN pathways might be in part responsible for differential gene induction in NK cells (Schultze and Aschenbrenner, 2021) . We tested this possible interaction by incubating peripheral NK cells from control individuals with or without TNF in the presence of two different IFN-α concentrations and assessed IFN target genes MX-1, IFI6, and ISG15 ( Fig.   4I ). While addition of TNF in presence of high levels of IFN-α, reminiscent of severe COVID-19, led to a further increase of IFN target gene expression, this was not observed under low level IFN-α. In a second set of experiments, TNF target genes were assessed in presence of TNF with or without low level IFN-α (Fig. 4J) . Here, addition of IFN-α reduced expression of TNF target genes, mirroring the transcriptional signatures in severe COVID-19.

Collectively, we observe strong TNF signature gene induction in moderate but not severe COVID-19, while IFN-α response genes are predominant in NK cells from severe COVID-19 and linked to IFN signaling being associated with unfavorable outcome.

Quantitative assessment of NK cell responses demonstrated a marked dysfunction of circulating NK cells after stimulation with K562 cells (Fig. 5A and S5A) . Disturbance of NK cell function was more pronounced in patients with severe COVID-19, who displayed reduction in percentages of IFN-γ + , and TNF + cells in both CD56 dim (Fig. 5B ) and CD56 bright NK cells (Fig. S5B) . However, cytotoxic degranulation was only impaired in the CD56 dim and not in the CD56 bright subgroup ( Fig. 5B and S5B ). Kinetic analysis demonstrated that cytokine production differed between moderate and severe COVID-19 over time ( Fig. 5C) with NK cells in severe COVID-19 showing a persistent functional disturbance after more than 2 weeks (Fig. 5C) . Similar observations were made when peripheral NK cells were stimulated with cytokines ( Fig. S5C and S5D ).

J o u r n a l P r e -p r o o f Next, we tested the anti-SARS-CoV-2 function of circulating COVID-19 NK cells. To this end, Caco-2 cells and Vero E6 cells were infected with SARS-CoV-2 and co-cultured with purified NK cells. Using SARS-CoV-2 Spike-specific nanobodies (Koenig et al., 2021) for quantification of virus protein levels in viable, active Caspase-3cells, we found peripheral NK cells from controls to reduce viral protein levels both in Caco-2 and Vero E6 cells (Fig.   5D ). In contrast, NK cells from both moderate and severe COVID-19 displayed impaired antiviral activity, independent of IL-2 prestimulation ( Fig. 5E and 5F ). Circulating NK cells increased the expression of active Caspase-3 in SARS-CoV-2-infected target cells, especially after pre-stimulation with IL-2. However, induction of active Caspase-3 did not differ between COVID-19 NK cells and controls ( Fig. 5G and 5H ). To test whether reduced NK cell IFN-γ/TNF production might be involved in impaired antiviral activity of COVID-19 NK cells, different concentrations of IFN-γ, TNF, or a combination of both cytokines were added to virus-infected cells. Both IFN-γ and TNF led to reduced viral RNA titers ( Fig. S5E and   S5F ) and decreased expression of Spike protein ( Fig. S5G and S5H ). In line with these findings, we found lower concentrations of IFN-γ and TNF in supernatants of COVID-19 NK cells compared to control cells after incubation with both cell lines ( Fig. 5I and 5J ).

Taken together, the antiviral activity of COVID-19 NK cells is markedly diminished and is associated with a decline in IFN-γ and TNF production.

Enhanced expression of immune checkpoint molecules on NK cells is suggested to be involved in ineffective antiviral immune responses (Hadjadj et al., 2020; Vabret et al., 2020; Wilk et al., 2020 , Kong et al., 2020 Schultheiß et al., 2020) . ScRNA-seq analysis revealed increased expression of several immune checkpoint genes in COVID-19, but no consistent differences were found between moderate and severe disease (data not shown). On the protein level, increased frequencies of PD-1 + , LAG3 + , and TIGIT + peripheral NK cells, especially in severe COVID-19, were observed in cohort 2 and higher proportions of TIM-3 + NK cells in cohort 1 ( Fig. S6A and S6B) . The proportion of TIM-3 + and PD-1 + NK cells was rather low and there was no correlation between IFN-γ production and frequency of TIM-3 + or PD-1 + NK cells (Fig. S6C) . Regarding TIGIT, we found increased IFN-γ production in TIGITthan in TIGIT + NK cells irrespective of COVID-19 severity. The severe COVID-19associated impairment of IFN-γ production was detected for both TIGIT + and TIGITsubpopulations (Fig. S6D) . In summary, little evidence was found for a definitive J o u r n a l P r e -p r o o f involvement of the checkpoint molecules, PD-1, TIGIT, LAG-3, or TIM-3, in functional NK cell impairment in COVID-19.

Given the increased concentrations of inflammatory and immunosuppressive cytokines observed in early severe COVID-19 (Fig. 4C) , we next incubated peripheral control NK cells with plasma from COVID-19 patients or controls. Incubation with severe COVID-19 plasma resulted in a marked functional impairment with decreased IFN-γ ( Fig. 6A and 6B ) and TNF ( Fig. 6C ) production, whereas plasma from patients with moderate disease had only minor effects. A more detailed analysis revealed these differences between moderate and severe COVID-19 mainly resulted from differences at weeks 2 and thereafter, resembling our ex vivo observations (Fig. 5C) . Ex vivo NK cell cytokine production per patient correlated with in vitro IFN-γ and TNF production of control NK cells after incubation with plasma from the respective COVID-19 patient (Fig. 6D) . These findings indicated soluble factors to be involved in COVID-19-associated NK cell dysfunction. In line with this hypothesis, we found IFN-α, TNF, and IL-6 suppressed NK production of IFN-γ (Fig. S6E) . However, neither blockade of individual cytokines nor simultaneous blockade of different cytokine combinations (data not shown) resulted in normalization of severe COVID-19 NK cell functions (Fig. 6E ). Yet, when culturing severe COVID-19 NK cells in the presence of plasma from controls, cytokine production and degranulation were almost completely restored ( Fig. 6F and 6G ). To test whether there is a direct effect of viral components, particularly the SARS-CoV-2 Nucleocapsid on NK cell function, we incubated NK cells from control donors with different concentrations of Nucleocapsid and analyzed IFN-γ production after co-culture with K562 cells. Neither with or without IL-2, control donor NK cells showed differences in IFN-γ production after incubation with Nucleocapsid ( Fig. S6F ).

While transcriptome analyses had illustrated altered transcriptional programming of NK cells in COVID-19, functional assays show that this is not an inherent cell-intrinsic characteristic, but a dysfunctional state triggered by severe COVID-19-associated soluble plasma components.

Severe COVID-19 beyond the second week is characterized by persisting clinical symptoms (Grasselli et al., 2020; Guan et al., 2020) . We therefore investigated molecular phenotype of NK cells in later stages of disease. Comparison of severe COVID-19 samples from week 3+ vs. all others in cohort 1 distinguished a group of differentially regulated genes, which were then assessed in the other cohorts (Fig. 7A) . While these genes also appeared to be J o u r n a l P r e -p r o o f differentially regulated in week 3+ for cohorts 2 and 3, the genes were already differentially regulated in week 2 for cohort 4. Fourteen genes showed similar average log fold changes in all cohorts (Fig. 7B) . Late phase NK cells from patients with severe disease were characterized by downregulated expression of IFN-related genes but higher expression of DUSP2 (a regulator of the ERK signaling pathway) (Jeffrey et al., 2006) as well as the glucocorticoid-inducible factor TSC22D3 and RNA-binding protein ZFP36L2, which are linked to immunosuppression (Salerno et al., 2018; Yang et al., 2019b) . In addition, we observed increased expression of the chemokine receptor CXCR4, and AREG (encoding for Amphiregulin (AR)), an epidermal growth factor receptor ligand, involved in pulmonary fibrosis (Ding et al., 2016) . Analysis of COVID-19 bronchoalveolar lavage fluid (BALF) samples (Liao et al., 2020) revealed the proportion of NK cells expressing higher levels of AREG and CXCR4 to be increased in severe COVID-19 (Fig. 7C) . MCFC confirmed CXCR4 upregulation on circulating CD56 dim NK cells (Fig. 7D , S7A and S7B), and AR expression on NK cells in late severe COVID-19 (Fig. 7E) .

Plasma from severe COVID-19 patients but not controls upregulated CXCR4 and AR ( Fig.   7F and 7G). We observed a positive correlation between post-culture expression of CXCR4 and AR (Fig. S7C ) resembling our findings on the transcriptome level (Fig. S7D) . Similar to severe COVID-19 (Fig. 7C, S7E) , upregulation of AREG, DUSP2, ZFP36L2, and TSC22D3 in pulmonary NK cells is also found in lung fibrosis (Habermann et al., 2020) (Fig. 7H) . To test the impact of COVID-19 NK cells on fibrotic activity of human lung fibroblasts, expression of the pro-fibrotic marker genes COL1A1 and ACTA2 were assessed ( Fig. 7I -J, S7F-H). Incubation with non-activated peripheral NK cells from COVID-19 patients had no effect on expression of pro-fibrotic genes in the fibroblasts (Fig. S7F) . However, after activation with IL-2, NK cells from control individuals reduced the expression of pro-fibrotic genes in fibroblasts, which was not the case after incubation with severe COVID-19 NK cells ( Fig. 7I ). Following activation with IL-2, NK cells from severe COVID-19 were impaired in inducing active Caspase-3 compared to controls and moderate COVID-19 in human lung fibroblasts (Fig. 7J, S7G ). Without activation, COVID-19 NK cells induced lower Caspase-3 than control NK cells, though no difference was observed regarding disease severity (Fig.

Natural killer cells are an essential part of the innate immune response and importantly involved in antiviral immune responses (Björkström et al., 2010; Blom et al., 2016; Kokordelis et al., 2014) . Increased intra-pulmonary NK cell frequencies (Chua et al., 2020; Liao et al., 2020; Xu et al., 2020) and anti-Spike-dependent NK cell responses observed in vaccinated macaques , suggest that NK cells also may play a role in SARS-CoV-2 infection. However, our understanding of the role of NK cells in COVID-19 is still limited. Here, we combined single-cell transcriptomics and proteomics together with comprehensive functional analyses for in-depth longitudinal characterization of natural killer cells during acute COVID-19. We analyzed a total of 205 patients (403 samples, day 2 to 41 after symptoms onset) from four independent cohorts which also allowed for cross-validation of our findings.

In line with earlier studies, we found COVID-19 to be associated with a decrease in circulating NK cells (Giamarellos-Bourboulis et al., 2020b; Jiang et al., 2020; Wang et al., 2020b; Wilk et al., 2020) and validated expression of NK cell activation markers, especially in severe COVID-19 (Maucourant et al., 2020; Varchetta et al., 2020) .

We found increased expression of ISGs and genes involved in IFN-α signaling to be characteristic for NK cells in severe COVID-19, whereas in moderate disease an upregulation of TNF-related genes was observed. Integrating our findings with earlier reports (Arunachalam et al., 2020; Blanco-Melo et al., 2020; Hadjadj et al., 2020; Lucas et al., 2020; Schulte-Schrepping et al., 2020; Stephenson et al., 2021; Su et al., 2020) , a picture emerges in which a type I IFN response is seen in early disease with subsequent decline of IFN-mediated signatures after week 1 in moderate COVID-19, while they stay elevated during week 2 in severe disease. Cross-regulation by different cytokines may play a role and explain our finding of downregulated expression of TNF-related genes in severe disease despite TNF plasma levels being similar or even higher than in moderate COVID-19. Indeed, TNF increased IFN-α-induced ISG expression, whereas IFN-α prevented upregulation of TNFrelated genes in NK cells, indicating a cross-regulatory interaction of these two cytokines (Cantaert et al., 2010; Karki et al., 2021) . Clinical relevance of the dysregulated IFN-α response in early COVID-19 was supported by our findings that plasma levels in week 1 were positively correlated to clinical parameters such as SOFA score and WHO severity grade and our observation that IFN-α imprint clearly discriminated between patients with J o u r n a l P r e -p r o o f fatal outcome and those that eventually recover. Thus, further studies are needed to fully address the role of differential IFN-α vs. TNF responses in COVID-19.

We further demonstrated that NK cells exert anti-SARS-CoV-2 activity but are functionally impaired in COVID-19. Type I IFNs have been shown to be of critical importance for IFN-γ production by NK cells in several viral infections (Baranek et al., 2012; Lee et al., 2017; ) .

Conversely, type I IFN can also suppress this NK cell function (Ahlenstiel et al., 2011; Lee et al., 2019) , depending on timing and extent of type I IFN produced (Marshall et al., 2006) . For instance, NK cells exert a basally high sensitivity to IFN-mediated STAT4 activation for IFN-γ production, but increase in IFN-α production during virus infection results in an increase in STAT1 thereby inhibiting IFN-γ production (Miyagi et al., 2007) . Such a scenario, in which a robust and punctual IFN-α response early after infection promotes effective antiviral immunity while a prolonged and excessive IFN-α production is detrimental, may also be relevant regarding the observed association of inborn errors in IFNα immunity or autoantibodies against type I IFNs with life threatening COVID-19. Here, dysregulation of IFN-α responses due to genetic defects or the pre-existence of autoantibodies may promote viral spread and propagation in the lung while longer lasting and excessive IFN-α production may finally result in impaired immune responses as observed in our study. However, impaired NK cell function was also observed after decline of IFN-α plasma levels and normalization of ISG expression.

Furthermore, blocking of IFN-α with a specific antibody was insufficient to prevent NK cell dysfunction induced by COVID-19 plasma, indicating additional factors being involved.

Our data also suggest that NK cell dysfunction not only affects antiviral immune responses but may also be relevant with respect to the development of fibrotic lung disease in severe COVID-19. NK cells have been shown to limit hepatic and cardiac fibrosis progression (Ong et al., 2015; Radaeva et al., 2006) and impaired antifibrotic NK cell activity has been associated with accelerated liver fibrosis (Glässner et al., 2013) . Here, we found that NK cells in late phase severe COVID-19 display a decreased antifibrotic capacity. Of note, NK cells in the later stage of severe COVID-19 expressed high levels of ZFP36L2 and TSC22D3 which have been linked with immunosuppressive effects in memory T cells (Salerno et al., 2018; Yang et al., 2019b) , and thus, may also interfere with NK cell activity. On the other hand, we found NK cells in late severe COVID-19 to display increased expression of DUSP2 and high surface expression of CD69 and CD38 indicating ongoing cell activation and inflammatory J o u r n a l P r e -p r o o f cell signaling (Jeffrey et al., 2006) , which have been shown to induce NK cell dysfunction (Alvarez et al., 2019; Merino et al., 2019) . Moreover, we observed elevated expression of AREG, encoding Amphiregulin both in circulating and lung NK cells. Data regarding the role of Amphiregulin in lung fibrosis are controversial (Branzk et al., 2014; Ding et al., 2016; Monticelli et al., 2011) , and little is known regarding the role of Amphiregulin-expressing NK cells. Interestingly, increased expression of AREG, DUSP2, ZFP36L2, and TSC22D3 was also shown in pulmonary NK cells in non-specific interstitial pneumonia (NSIP) and idiopathic pulmonary fibrosis, a fibrotic lung disease which resembles COVID-19 with respect to radiological and clinical findings.

Collectively, our study points to differential IFN-α vs. TNF responses as an important mechanism in the early phase of COVID-19 and describes a link between an exaggerated prolonged IFN-α-induced NK cell response and persistent NK cell dysfunction with an unfavorable course of the disease.

Small differences in scRNA-seq analysis between the 4 cohorts might be explained by different geographical location, local SARS-CoV-2 variants, different experimental setups or different sampling strategy, e.g. longitudinal (cohorts 1-3) vs. cross-sectional sampling (cohort 4). Still, longitudinal studies of COVID19 utilizing single cell omics are rare and it would be probably beneficial, if for cell types such as NK cells, but most likely even more important for other rare immune cell types, additional studies are being conducted that would allow to increase the number of patients and the number of cells to be analyzed.

Our study uncovered soluble factors to be responsible for NK cell dysfunction, as evident from experiments using patients' plasma. Yet, while we excluded many options, further studies are necessary to clarify which other components might account for this effect. 

(E) Pearson correlation between numbers of absolute and CD56 dim NK cells and serum CRP levels.

(F) Frequency of NK cells positive for active Caspase-3 or CD95 in cohort 1.

(G) Detection of CD95 and active Caspase-3 in control NK cells co-incubated without or with Nucleocapsid.

Kruskal-Wallis (KW) test corrected for multiple comparison by controlling FDR (Benjamini, Krieger, Yekutieli (BKY)); *p<0.05, **p<0.01, ***p<0.01.

For n see Table S6 . false discovery rate (FDR) cutoff of 5%. Hierarchical clustering of gene modules and functional enrichment using the KEGG and HALLMARK databases (Table S3 ).

(C) UMAP of NK cells from cohort 1 (scRNA-seq; 10,927 cells). NK subtypes defined by cluster marker expression and reference-based NK annotations (Table S2) For n see Table S6 . (Table S4 ).

(D) Selected results from functional enrichment analysis using the GO, KEGG and HALLMARK databases, transcription factor (TF) prediction and upstream ligand prediction for each identified heatmap module from (C) (entire list see Table S4 ).

(E) Heatmap of mean area under the curve (AUC) scores based on AUCell enrichment of heatmap gene modules from (C) for NK subtypes of cohort 1 (scRNA-seq).

(F) NK subtype occupancy over time in days after symptom onset as average of all samples stratified by severity.

(G) Density plot of cell frequency by disease severity and weeks after onset overlaid on the UMAP of cohort 1 (FC data).

(H) Heatmap divided by disease severity and weeks after onset showing the proportion of each severity group for the three identified NK subtypes of cohort 1 (FC data).

For n see Table S6 .

J o u r n a l P r e -p r o o f Table S6 . For n see Table S6 .

J o u r n a l P r e -p r o o f

Lead contact Further information and requests for resources and reagents should be addressed to and will be fulfilled by the Lead Contact Jacob Nattermann (Jacob.Nattermann@ukbonn.de).

This study did not generate unique reagents.

Single-cell RNA-seq data have been deposited at the European Genome-phenome Archive (EGA) and are publicly available as of the date of publication. Accession numbers are listed in the key resources Any additional information required to reanalyze the data reported in this paper is available from the lead contact upon request.

Samples from patients with COVID-19 were collected within four cohort studies designed to allow deep molecular and immunological transcriptomic and proteomic profiling of COVID-19 in blood. Patients were classified according to the highest score on the World Health Organization (WHO) Ordinal Scale for Clinical Improvement ever present ((WHO. R&D Blueprint -novel Coronavirus -COVID-19 Therapeutic Trial Synopsis. 2020.https://www.who.int/blueprint/priority-diseases/key-action/COVID-19_Treatment_Trial_Design_Master_Protocol_synopsis_Final_18022020.pdf.). Patients for which sufficient material was available for scRNA-seq, CyTOF or flow cytometry analysis, were included in this study. This study was designed to describe immunological deviations in COVID-19 patients without intention of the development of new treatments or new diagnostics, and therefore sample size estimation was not included in the original study design. Table S1 ) were included in the study. In patients who were not able to consent at the time of study enrollment, consent was obtained after recovery. Information on age, sex, medication, and comorbidities are listed in Table S1 . COVID-19 patients who tested positive for SARS-CoV-2 RNA in nasopharyngeal swabs were recruited at the Department of Internal Medicine I of the University Hospital Bonn or the Department of Gastroenterology, Hepatology and Infectious Diseases, University Hospital Düsseldorf, between March 30 and November 11, 2020 and allocated to moderate (WHO 2-4) or severe (5-7) disease according to the WHO clinical ordinal scale. Controls in cohort 1 were collected from healthy people or from otherwise hospitalized patients with a wide range of diseases and comorbidities including chronic inflammatory immune responses. These individuals were either tested negative for SARS-CoV-2, serologically negative or samples were collected before November 2019.

For validation of the findings from our prospective cohorts, data from three independent additional cohorts were analyzed:

This study includes a subset of patients enrolled between March 2 and July 02, 2020 in the Pa-COVID-19 study, a prospective observational cohort study assessing pathophysiology and clinical characteristics of patients with COVID-19 at Charité Universitätsmedizin Berlin . The study is approved by the Institutional Review board of Charité (EA2/066/20). Written informed consent was provided by all patients or legal representatives for participation in the study. The patient population included in all analyses of cohort 1 consists of 10 control donors (samples collected in 2019 before SARS-CoV-2 outbreak), 8 patients presenting with flu-like illness but tested SARS-CoV-2-negative, 12 moderate and 17 severe COVID-19 patients ( Figures 1A and 1B ; Table S1 ). Information on age, sex, medication, and co-morbidities is listed in Table S1. All COVID-19 patients were tested positive for SARS-CoV-2 RNA in nasopharyngeal swabs and allocated to mild (WHO 2-4) or severe (5-7) disease according to the WHO clinical ordinal scale. We also included publicly available single-cell transcriptome data derived from 22 control samples into the analysis; 3 samples were derived from 10x Genomics, San Francisco, CA 94111, USA (5k_pbmc_v3:

https://support.10xgenomics.com/single-cell-geneexpression/datasets/3.0.2/5k_pbmc_v3, pbmc_10k_v3: https://support.10xgenomics.com/single-cell-gene-expression/datasets/3.0.0/pbmc_10k_v3, pbmc_1k_v3:

https://support.10xgenomics.com/single-cell-geneexpression/datasets/3.0.0/pbmc_1k_v3), 19 samples derived from Reyes et al. (2020) .

In cohort 3 (Bernardes et al., 2020) COVID-19 patients were sampled in two independent University hospitals (Cologne, Kiel) between April 1, 2020, and May 20, 2020. From this study, patients were enrolled in our analyses if cell numbers were sufficient to enable identification and in-depth analysis of NK cells. In total, 8 COVID-19 patients and 2 controls were included. Information on age, sex, medication, and co-morbidities is listed in Table S1 .

In cohort 4, we combined datasets from the UK (Stephenson et al., 2021) For the UK cohort, each patient contributed to one sample while in the US, each patient was sampled exactly two times, therefore, the sampling strategy is rather cross-sectional. For comparison reasons, samples from COVID-19 patients that were sampled later than 4 weeks after symptom onset and samples from patients who received steroid treatment were removed from further analysis. After removal, a total of 30 controls with 40 samples and 110 COVID-19 patients with 161 samples were included in cohort 4. COVID-19 patients were allocated to moderate or severe disease by maximum WHO ordinal scale (US data) or "status_on_day_collection" (UK data). Information on age, sex, medication, and comorbidities are listed in Table S1 .

Vero E6 cells, a cell line originating from Chlorocebus aethiops, were continuously maintained in complete DMEM medium at 37°C and passaged upon reaching 80% confluence. Prior to infection experiments, cells were seeded in 96-well flat bottom plates and cultured for 48h at 37°C. Caco-2 cells, a cell line originating from a male human individual, were continuously maintained in EMEM medium at 37°C and cultured upon reaching 80% confluence. Prior to infection experiments, cells were seeded in 96-well flat bottom plates and maintained for 48h at 37°C. Primary human lung fibroblasts were commercially obtained (PromoCell) and cultured in fibroblast media (PromoCell) at 37°C according to the manufacturer's instructions. Cells were kept in culture no longer than to passage 4. For each passage, it was tested whether the cells could still be activated with recombinant TGFb1 after 2 days incubation at 37°C (10ng/ml; Miltenyi; readout PCR, see below). Prior to coincubation experiments, cells were seeded in 96-well flat bottom plates and maintained for 48 h at 37°C. Information on the sex of primary human lung fibroblasts remained undisclosed by provider.

All cells used are also listed in the key resource Table. No additional information on cell line authentication is provided.

Cell isolation (cohort 1) Peripheral blood mononuclear cells (PBMC) were isolated using Ficoll-Paque gradient centrifugation (Biochrom AG, Berlin, Germany), washed with DPBS, and directly cryopreserved in RPMI-1640 supplemented with 10% DMSO (Sigma). The processing of cells for scRNAseq analysis was described for the 3 cohorts here (Bernardes et al., 2020; Schulte-Schrepping et al., 2020).

Phenotypic analysis of cells was performed using an LSR-Fortessa Cytometer (BD Biosciences, USA). In brief, frozen cells were gently thawed at room temperature and transferred to a 14 ml tube. Then 2 ml of thawing medium (HBSS; 1% human serum albumin, CSL Behring; and 25 U/ml endonuclease, MoBiTec, Germany) was added dropwise, waited for 2 min and 8 ml of thawing medium was finally added gently. After 15 min at room temperature, cells were centrifuged at 300g for 10 min and incubated with viability dye (Zombie-Aqua, Biolegend, 1:500) for 10 min. After further washing with DPBS (10 min, 300g), cells were stained with appropriate antibody solutions in the Biolegend staining buffer. All antibodies were titrated, the panels were tested using FMO controls, and constant conditions were ensured by plate staining to guarantee an optimal staining result. To verify consistent fluorescence properties of the flow cytometer, calibration beads (ultra Rainbow beads, Spherotech) were applied before each measurement. The antibodies and panels used in this study are compiled for each corresponding Figure in Table S5 . NK cells were defined as CD45+CD56+Lin-lymphocytes (Lin: CD3, TCRab, TCRgd, CD34, CD20, CD19, CD14) with exclusion of CD94-NKp80-cells. For intracellular analyses of transcription factors the Foxp3 Transcription Factor Staining Kit (eBioscience, Germany) was used for permeabilization, fixation, and washing according to the manufacturer's specifications. FC raw data was analyzed by FlowJo software V10.6.1 (BD Bioscience, USA).

Purified NK cells from control donors were incubated with recombinant IFN-α (Immunotools, 1 or 10ng/ml) in combination with or without recombinant TNF (Immunotools, 10 or 25ng/ml) for 18h. Incubation was stopped with a lysis buffer from the RNA isolation kit (New England Biolabs, Monarch total RNA Mini Prep kit). After RNA isolation according to the manufacturer's protocol, quality control using NanoDrop and transcription into cDNA (QuantiTect RT Kit, Quiagen) were performed. ISG (MX-1, IFI6 and ISG15) and TNF hallmark transcripts (MAP3K, TNF1IP3, and LITAF) were analyzed by qPCR (96-well LightCycler; Roche). Relative gene expression (duplicates) was calculated by 2 -ΔCq method related to 2 housekeepers (geomean of RPL19 and EEF1). Primer sequences are listed in Supplementary Table S7 .

To assess ex vivo functionality of NK cells, PBMC from all study groups were thawed as described above and NK cells were isolated by negative magnetic separation according to the manufacturer's instructions (negative NK cell isolation kit, Miltenyi). After checking purity on the flow cytometer (at least 90% CD56 + CD3of CD45 + lymphocytes), NK cells (2x10 5 cells per ml) were seeded and incubated with and without IL-2 (10ng/ml, Miltenyi), for 18h in defined DMEM/F12 media containing DMEM (Gibco Life, USA) and F12 (2:1), 1% antibiotic and antimycotic (Gibco Life, USA), 20 mg/mL ascorbic acid (Sigma, USA), 24 mM 2-mercaptoethanol (Gibco Life, USA), 0.05 mg/mL sodium selenite (Sigma, USA), and 10% heat-inactivated human AB serum (Sigma, USA) based on previous protocols in 96-well round bottom wells. Alternatively, isolated NK cells (2x10 5 cells per ml) were incubated with plasma from healthy controls, moderate and severe COVID-19 patients for 18h at a ratio of 1 to 5 in defined DMEM/F12 media. Afterwards, NK cells were optionally co-incubated for 5h at a ratio of 1:2 with major histocompatibility complex-deficient K562 cells or co-incubated with a cytokine cocktail consisting of IL12 (10ng/ml, Immunotools), IL-15 (50ng/ml, Immunotools) and IL18 (50ng/ml, Immunotools) in defined DMEM/F12 media. After 1h stimulation, Brefeldin A (5μg/ml; Enzo, Germany) was added for the remainder of the incubation. After staining with the viability dye Zombie Aqua (Biolegend, USA) and surface antibodies (see key resource Table) , cells were permeabilized using the Cytofix/Cytoperm Kit according to the manufacturer (BD Biosciences, USA). NK cells were defined as CD56+Lin-lymphocytes (Lin: CD3, TCRab, CD34, CD20, CD19, CD14) with exclusion of CD94(-)Nkp80(-) cells.

IFN-γ and TNF were detected with specific antibodies by intracellular staining. FC raw data was analyzed by FlowJo software V10.6.1 (BD Bioscience, USA).

Further functional assays were performed to assess the impact of cytokine and SARS-CoV-2 Nucleocapsid stimulation on NK cell functionality. In this setting, purified NK cells (2x10 5 cells per ml) from controls were seeded and incubated with and without IL-2 (10ng/ml, Miltenyi), TNF (Immunotools, 10ng/ml), IFN-α (Immunotools, 10ng/ml), IL-6 (Immunotools, 10ng/ml), or IL-10 (Immunotools, 10ng/ml), or SARS-CoV-2 Nucleocapsid (Sinobiological, 20, 200, or 2000ng/ml) for 18h in defined DMEM/F12 media containing DMEM (Gibco Life, USA) and F12 (2:1), 1% antibiotic and antimycotic (Gibco Life, USA), 20 mg/mL ascorbic acid (Sigma, USA), 24 mM 2-mercaptoethanol (Gibco Life, USA), 0.05 mg/mL sodium selenite (Sigma, USA), and 10% heat-inactivated human AB serum (Sigma, USA) based on previous protocols in 96-well round bottom wells. After pre-stimulation, NK cells were optionally co-incubated for 5h at a ratio of 1:2 with major histocompatibility complex-deficient K562 cells or co-incubated with phorbol-12-myristate-13-acetate (PMA, 50ng/ml; Cell Signaling Technology Europe, Netherlands) and ionomycin (1µg/ml; Cell Signaling Technology Europe). Follow-up steps were performed as described above for detection of IFN-γ and TNF-α.

Furthermore, the impact of soluble plasma-derived factors on NK cell functionality was assessed. For this purpose, purified NK cells (2x10 5 cells per ml) from all study groups were seeded and optionally incubated with plasma from healthy controls, moderate or severe COVID-19 patients or alternatively combined with specific blocking antibodies (IL-10, IL-12, IL-6, IL1b, IL-4, IFN-AR, TNF, ISO) for 18h at a ratio of 1 to 5 in defined DMEM/F12 media. Afterwards, NK cells were co-incubated for 5h at a ratio of 1:2 with major histocompatibility complex-deficient K562 cells. Follow-up steps were performed as described above for detection of IFN-γ and TNF.

To test fibrotic factors in vitro, 1x10 4 primary human lung fibroblasts, seeded in 96-well flat bottom plates, were co-incubated with recombinant amphiregulin at 37°C and experiment was stopped with lysis buffer from the RNA isolation kit after 3 days (New England Biolabs, Monarch total RNA Mini Prep kit). After RNA isolation according to the manufacturer's protocol, performing quality control using NanoDrop, and transcription into cDNA (QuantiTect RT Kit, Quiagen), lung fibroblast activity was analyzed by qPCR (96-well LightCycler; Roche) with detection of genes for COL1A1 and ACTA-2. Relative gene expression (duplicates) was calculated by 2 -ΔCq method related to 2 housekeepers (geomean of RPL19 and EEF1). 2 -ΔCq values were further normalized by dividing the sample values by the mean of the control values without NK cells. Primer sequences are listed in Supplementary Table S7 .

To test the impact of NK cells on the activity of primary lung fibroblasts, 2 day seeded fibroblasts were co-incubated with 2x10 4 isolated NK cells (pre-stimulated with or without 10ng/ml IL-2 for 18h) from controls and patients with moderate or severe COVID-19 progression for 18 hours in defined DMEM/F12 media at 37°C. Then the supernatant was discarded and co-culture was carefully washed with PBS. Sufficient removal of NK cells was checked by light microscope (Zeiss, AxioVert 200M), the experiment was then stopped with a lysis buffer and qPCR was performed as described above.

To test the anti-fibrotic activity of NK cells, human lung fibroblasts were labeled with live dye e670 (1µM; ebioscience; to distinguish lung fibroblasts from NK cells) for 10 min at room temperature in the dark, washed twice with PBS, and seeded as described above in fibroblast media at 37°C. After 2 days, 2.5x10 4 isolated NK cells (pre-stimulated with or without 10ng/ml IL-2 for 18h) from the subject groups were co-incubated with the lung fibroblasts for 6h in defined DMEM/F12 media at 37°C. The supernatant was collected, remaining lung fibroblast were incubated with accutase (Sigma) for 5 min for detaching, sufficiently washed with PBS to remove all cells and reunited with the collected supernatant. After staining with the viability dye Zombie Aqua (Biolegend, USA), cells were permeabilized using the Cytofix/Cytoperm Kit according to the manufacturer (BD Biosciences, USA) and labeled with specific antibody for active Caspase-3 for detecting of apoptotic lung fibroblasts by flow cytometry (BD Canto 2). FC raw data was analyzed by FlowJo software V10.6.1 (BD Bioscience, USA).

For in vitro infection, 5x10 3 Vero E6 cells or 6x10 3 CaCo-2 cells were seeded in 96-well flat bottom plates and cultured at 37°C. After 48h cells were infected with SARS-CoV-2/human/Germany/Heinsberg-01/2020 virus at a MOI 0.1 (Vero E6 cells) or 1.0 (CaCo-2 cells). After one hour, the inoculum was removed, and cells were washed once with DPBS. Then, cells were cultured for additional 48h in the presence of increasing concentrations (0ng/ml, 1ng/ml, 10ng/ml) of recombinant human IFN-γ (Immunotools) or TNF-α (Immunotools) or a combination of both (10ng/ml each). Alternatively, purified NK cells (1x10 4 ) were incubated with or without recombinant IL-2 (10ng/ml) for 18h and were added 24h post-infection and co-cultured with Vero E6 or Caco-2 cells, respectively, in defined DMEM/F12 media for another 24h. Then, supernatant was removed and collected, the cells were washed twice with DPBS, controlled by microscope and tested for SARS-CoV-2 RNA replication or stained for SARS-CoV-2 Spike protein using specific nanobodies (Koenig et al., 2021) . In the supernatant, concentrations of the cytokines IFN-γ and TNF were detected using the Cytometric Bead Array (CBA) from Legendplex (Essential Immune Response Kit), and the beads were measured on the BD FACS Canto 2.

Detection of SARS-CoV-2 infection was performed via two approaches, assessing either expression of viral RNA by PCR or expression of viral proteins using fluorochromeconjugated nanobodies. For detection of SARS-CoV-2 specific genes, treatment of cells, RNA isolation and transcription into cDNA was performed as described above. Primers were purchased (N1/N2 target, SARS-CoV-2 (2019-nCoV) CDC RUO Kit, or sequences were obtained from other published work (M-gene, primer sequences in Supplementary Table  S7 ) (Toptan et al., 2020) . Relative gene expression was calculated with the use of duplicates by 2 -ΔCq method related to 2 housekeeping genes (RPL19, EEF1A, Supplementary Table  S7) .

For intracellular detection of SARS-CoV-2 Spike protein, in vitro infections of Vero E6 or CaCo-2 cells were performed as described above and cells were co-cultured with recombinant cytokines or isolated NK cells accordingly. arget cells were detached with accutase (Sigma) and washed with DPBS. After staining with viability dye Zombie Aqua (Biolegend, USA) cells were fixed/permeabilized using the Cytofix/Cytoperm Kit according to the manufacturer (BD Biosciences, USA) and subsequently incubated with the Spike specific nanobody VHH E (AF488 labeled) (Koenig et al., 2021) , or control nanobody LaM-4 (anti-mCherry, AF488 labeled) (Fridy et al., 2014) for 30 min in the dark by shaking at 4°C. Analysis was performed on a flow cytometer (BD Canto-2). FC raw data was analyzed by FlowJo software V10.6.1 (BD Bioscience, USA).

Plasma concentrations of the cytokines IL-5, IL-4, IL-12p70, IL-22, IL-8, IL-10, IL-6, IFN-γ, TNF-α, and IL-1b were measured with planar-array technology on the Quanterix® SP-X Imaging and Analysis System™ using the Simoa® CorPlex™ Human Cytokine Panel 1 assay (Item 85-0329). Plasma concentrations of IFN-α were also measured on the SP-X system using an assay in development, which will become commercially available in the future. Plasma concentrations of SARS-CoV-2 Nucleocapsid protein and SARS-CoV-2 anti-Spike IgG were measured using digital bead-based technology on the Quanterix® HD-X™ Analyzer with development versions of the assays Simoa® SARS CoV-2 N Protein Advantage Kit (Item 103806) and Simoa® SARS-CoV-2 Spike IgG Advantage Kit (Item 103769).

Processed and previously published PBMC scRNA-seq datasets from Schulte-Schrepping et al (cohort 1, Bonn data and cohort 2, Berlin data)(Schulte-Schrepping et al., 2020) and from Bernardes et al (cohort 3, Kiel data)(Bernardes et al., 2020) were downloaded from FastGenomics (https://www.fastgenomics.org) as Seurat objects and datasets from Su et al. (Su et al, 2020) and Stephenson et al.(Stephenson et al., 2021) (cohort 4, UK and US data) were received directly by the authors and downloaded from the COVID19 cell atlas (https://covid19cellatlas.org/), respectively and all datasets were imported to R 4.0.0. Subsequent gene expression data analysis analysis was performed using the R/Seurat package 3.2.0 (cohorts 1, 2 and 3) and 3.9.9 (cohort 4) (Butler et al., 2018; Stuart et al., 2019) .

Cohort 3 initially used a different disease severity group annotation than the other cohorts which had to be adjusted. Patients marked as "complicated (incremental)", "complicated (recovering)", "complicated with hyperinflammatory syndrome" and "critical" were binned into the group "severe", while "mild (recovering)" were annotated as "moderate". Patients annotated as "recovered" and "recovery/asymptomatic" were removed from further analysis. For cohort 4, the maximum WHO ordinal scale (US data) or "status_on_day_collection" (UK data) per patient was used as a discriminator between moderate and severe diseased patients. Patients who received steroid treatment were removed from further analysis. Moreover, samples with more than 4 weeks after onset of symptoms were taken out of cohort 3 and 4, since the main interests of this study are the changes in the first weeks. For all datasets, samples with 3 and 4 weeks after onset of symptoms were annotated as week "3+" after symptoms.

To analyze the data without having any influence of batch effects resulting from the two different studies (UK (Stephenson et al, 2021) and US (Su et al., 2020) and their locations, J o u r n a l P r e -p r o o f the Seurat implemented integration approach based on "anchors" across collection sites (Stuart et al., 2019) was used to harmonize and integrate the two datasets following the default settings but increasing the integration features to 10,000. Subsequently, the merged dataset was scaled, PCA was performed and UMAP was calculated based on the first 30 PCs.

Data integration approaches for cohorts 2 and 3 or 1 to 3 as an alternative to using a validation cohort approach revealed similar results (data not shown). Since the validation cohort approach did not require any correction for batch effects due to the use of different technologies (10x Genomics vs. BD Rhapsody) and since other parameters might also have been different between the three different clinical sites, we opted for the validation cohort approach to present the data within this manuscript.

NK cells present in each cohort were selected in a three-step process. This process is described as exemplary for cohort 1. First, the entire T and NK lymphocyte fraction was subsetted based on the cell type label provided by Schulte-Schrepping et al, including all T cells, NK cells and proliferating cells. This subset was subsequently normalized, scaled and dimensional reduction was calculated using the standard Seurat functions. For normalization, the gene expression values were normalized by total UMI counts per cell, multiplied by 10,000 (TP10K) and then log transformed by log10(TP10k+1). Subsequently, the data was scaled, centered and regressed against the number of detected transcripts per cell to correct for heterogeneity associated with differences in sequencing depth. For dimensionality reduction, PCA was performed on the top 2,000 variable genes identified using the vst method. For two-dimensional representation of the data structure, uniform manifold approximation and projection (UMAP) was calculated using the first 30 principal components (PCs). Next, the NK cells within this subset were identified: After UMAP calculation, cells were clustered using the Louvain algorithm based on the first 30 PCs and a resolution of 0.2. The cluster consisting of NK cells was identified using classical NK cell marker (KLRF1, GNLY, NKG7) and cluster-specific marker genes calculated with the Wilcoxon rank sum test using the FindAllMarkers Seurat function (parameters: min.pct=0.25, logFC.threshold=0.25). This NK cell cluster was selected and subsequently normalized, scaled and dimensionality reduction was calculated. Next, UMAP and clusters were calculated. All steps were performed with the same parameters as described above. Lastly, the NK cell subset was cleaned from non-NK cells. For this, the cells were overclustered using the Louvain algorithm with a resolution of 1. Cluster-specific marker genes calculated by the Wilcoxon rank sum test (same parameters as above) and the expression of KLRF1 were used to identify NK and non-NK cells. Clusters expressing classical marker genes related to other cell types, such as NKT cells (TRAC), CD8+ T cells (CD8A) and others, as well as clusters showing remarkably high expression of hemoglobin related genes (HBB, HBA1, HBA2) were removed from the dataset to yield clean NK cells (Table S2) . Hemoglobin-rich clusters may result from erythrocytes which contaminate the cells during sequencing preparation. To account for a donor-specific batch-effect in cohort 3, the first 30 PCs of the "harmony" algorithm (Korsunsky et al., 2019) were used instead of the PCs calculated by PCA, all other steps remained similar.

Differential expression (DE) tests were performed using FindMarkers or FindAllMarkers functions in Seurat with Wilcoxon Rank Sum test. Genes with a log-fold change greater than 0.2, at least 10% expressed in tested groups and with a bonferroni-corrected p-value < 0.05 were considered as DEGs. Group/subtype specific marker genes were identified by applying the DE tests for upregulated genes between cells in one group/subtype to all others in the dataset.

Gene set ontology enrichment analysis using the heatmap modules as input was performed on the gene sets from the Gene Ontology (GO) biological process (BP) database (Ashburner et al., 2000; The Gene Ontology Consortium, 2019) , the Kyoto Encyclopedia of Genes And Genomes (KEGG) database (Kanehisa, 2019) and the Hallmark gene sets (Liberzon et al., 2015) using the R package clusterProfiler (version 3.16.1) (Yu et al., 2012) . Ontologies with highest and statistically significant enrichment were used for presentation.

The R package RcisTarget (version 1.8.0) (Aibar et al., 2017) was used to predict the transcription factors potentially regulating heatmap module-specifically contained gene sets. The genomic regions of TF-motif search were limited to 10kb around the respective transcriptional start sites by using the RcisTarget-implemented "hg19-tss-centered-10kb-7species.mc9nr.feather" motifRanking file. Prediction was performed using the cisTarget function and the resulting top 3 predicted TF, according to their normalized enrichment scores (NES), were selected for each heatmap module.

Prediction of potential upstream ligands of each heatmap module gene set was performed using the R package NicheNetR (version 1.0.0) (Bonnardel et al., 2019; Browaeys et al., 2020) . For each heatmap module, the top 3 predicted ligands were selected according to their Pearson correlation coefficient (PCC).

Enrichment of gene signature sets was performed using the "AUCell" method (Aibar et al., 2017) implemented in the R package (version 1.10.0). We set the threshold for the calculation of the AUC to the top 3% of the respective ranked genes and normalized the maximum possible AUC to 1. Overlay analysis of donor origin of single cells onto the violin plots displaying IFN-α response and TNF signaling enrichments ensured that cells within NK cell clusters at defined time intervals and patient groups were not dominated by individual patients (data not shown).

For NK subtype definition of cohort 1, the combined information gained from cluster-specific DE gene expression and the enrichment of literature-based NK RNA signatures were used. Here, we evaluated sequencing-based NK annotations from the literature (Crinier et al., 2018; Smith et al., 2020; Yang et al., 2019a) and decided to use the signatures from Smith et al. . First, the entire NK population was clustered using the Louvain algorithm based on the first 30 PCs with a resolution of 0.7, resulting in a total of 8 clusters. Subsequently, DEGs for each cluster were calculated using the Wilcoxon rank sum test as described above. Due to a similar gene expression profile and close proximity, the first 3 clusters were united into one. Based on the upregulated DEGs of these 6 clusters, scRNA-seq NK signatures from Smith et al. were applied for gene set enrichment using the AUCell as described above. With the combined results, the 6 distinct NK clusters were annotated according to their subtype as inflamed CD56 dim (high IFN-related genes), CD56 bright (NCAM1), proliferating CD56 dim (MKI67), cytokine CD56 dim (CCL4, CCL3, IFN-γ), HLA hi CD56 dim (HLA-DP and HLA-DR related genes) and CD56 dim (FCGR3A) NK cells (Table  S2) . NK subtypes in cohorts 2, 3 and 4 were then annotated based on the markers identified in cohort 1.

To compare shifts in the NK subtypes stratified by disease group, the percentages of each subtype were quantified per sample and visualized together in boxplots. For determination of statistical significant differences in the distribution per disease group, a Kruskal-Wallis test with FDR correction was performed. Subtypes showing significant changes (FDR-corrected Kruskal-Wallis p value < 0.05) were further tested with a Dunn's Post-Hoc test using the "dunn.test" R package (version 1.3.5). Resulting p values were corrected for multiple testing using the Benjamini-Hochberg method.

For each NK subtype, the relative proportion across disease severity was visualized as a fraction of samples from the respective condition contributing to the NK subtype.

For each severity group, the proportional occupancy of the NK subtypes was calculated for all samples and their respective timepoints (in days after onset of symptoms) and the relative proportions were subsequently visualized as a function of time.

To analyse the correlation of IFN-related genes (IFIT2, IFIH1, IFI44L, XAF1, MX1, IFI6, ISG20, ISG15) with plasma cytokine levels, samples having both information (9 moderate, 13 severe) were selected and average gene expressions for the respective genes were calculated per sample. Subsequently, Spearman correlation was performed and the Spearman coefficients displayed using a heatmap. Significant correlations (p<0.05) were indicated.

Cellular Indexing of Transcriptomes and Epitopes by sequencing (CITE-seq) information of cohort 4 from US data (Su et al., 2020) and UK data (Stephenson et al., 2021) were used to strengthen the annotation of the clusters calculated for protein data based on the previously identified scRNA-seq clusters (see below). Since both studies of cohort 4 included different antibody markers, CITE-seq information was analyzed separately. In general, both CITE-seq datasets were first normalized using the CLR normalization method and subsequently scaled using the ScaleData function implemented in Seurat. For the US (Su et al.) dataset, a batchcorrection was performed by normalizing and scaling the major batches 6, 7, 9 and 10 individually and then merging the data again.

After pre-processing of the data (see section above), compensated fluorescence intensities were exported from FlowJo (BD, v. 10.7.1) for all the cells in the singlets/Lineage(-)/Living/CD56(+)/NKp80(+) gate. Exported .fcs files were imported in R (v. 3.6.2, Bioconductor v. 3.10) with the flowCore package (Bioconductor, v. 1.52.1) (Hahne et al., 2009 ). Fluorescence intensity values were then transformed with an auto-logicle transformation. Here we calculated the optimal width for the logicle transformation for each of the fluorescence parameters using the formula w=(m-log10(t/abs(r)))/2, where r is the most negative value to be included in the display (Chen et al., 2016) . Transformed fluorescence values from each experimental batch were used for dimensionality reduction using the UMAP algorithm (umap, CRAN, v. 0.2.6.0, n=15, mdist=0.2, metric="euclidean") (McInnes et al., 2018) . For each sample a maximum of 1.000 cells were randomly selected. To account for the small differences derived from the measurement of the samples in different experimental dates, the batch correction algorithm harmony (GitHub v. 1.0) (Korsunsky et al., 2019) was used to normalize the data. Normalized values were now used for UMAP dimensionality reduction (umap, CRAN v. 0.2.7.0) (McInnes et al., 2018) and Phenograph clustering (Rphenograph, GitHub v. 0.99.1, with number of nearest neighbours (k)=60) (Levine et al., 2015) . The expression of the markers CD69, CD38, CD94, KI67, CD161, CD95, FASL, NKp80, GrzB, HLA-DR, TIGIT, CD16, CD56 were used for the calculation of the UMAP and Phenograph clustering. After this step, data were visualized according to the Phenograph clustering after manual annotation of selected metaclusters to match the scRNA-seq annotation. Heatmaps of marker expression were calculated as scaled mean of the transformed fluorescence intensity of each marker for each Phenograph cluster. Confusion matrices were calculated normalizing first the cells from each condition to a total number of 1000 cells and later calculating the relative contribution of each condition on each cluster. Statistical testing for the difference in the frequency of each cluster was calculated with a Kruskal-Wallis test (with FDR correction, R v 3.6.2) to identify clusters where the grouping has an impact on the cluster frequency. Clusters having significant Kruskal-Wallis test (FDR corrected p value < 0.05) were tested with a Dunn's Post-Hoc test with Benjamini-Hochberg correction for multiple testing (CRAN rstatix package v 0.7.0). All heatmaps were calculated with the pheatmap package (CRAN, v. 1.0.12) and boxplots with the R package ggplot2 (CRAN, v. 3.3.3) . The analysis was performed in a dockerized environment (lorenzobonaguro/flowtools:R362_V2).

NK cells were identified based on exclusion of CD3 + , CD19 + , CD15 + and CD14 + cells and subsequent gating on all CD56 + cells. Because leukocyte counts were not available for all control samples, we use absolute NK cell numbers from our recently published control cohorts (Kverneland et al., 2016; Sawitzki et al., 2020) for comparison as also described in Schulte-Schrepping et al. (Schulte-Schrepping et al., 2020) . NK cells were clustered based on the expression of 16 markers: CD16, CD56, HLA-DR, CD38, CD69, Ki67, CXCR5, CXCR3, CCR6, PD1, TIGIT, CD226, CD62L, KLRB1, KLRG1, and KLRF1. The raw values obtained with CyTOF were first transformed with the inverse hyperbolic sine function (asinh) and then z-score normalised per marker. We clustered NK cells using Phenograph (Levine et al., 2015) with 100 nearest neighbours (k = 100). We found 15 clusters, which were annotated based on the average z-score transformed CyTOF expression of the markers in each cluster. Similarly, UMAPs were calculated with the selected markers mentioned above using the R package "uwot" with default parameters (Melville et al., 2020) . Statistical testing for the difference in the frequency of each cluster was calculated with the same approach as described below, i.e., Dunn's Post-Hoc test with Benjamini-Hochberg correction for clusters with significant Kruskal-Wallis test (FDR corrected p value < 0.05). For the non-weekly analysis, we considered the first sample per patient. For the weekly analysis, only the first sample per week was included and the repeated samples in the same week were excluded from the analysis.

The assignment of the FC/CyTOF sub-clusters to the 6 NK cell subsets (see above "NK subtype annotation") was done by comparing protein, transcriptome and CITE-seq expressions from the respective subset (cohort 1: Figures S2D, S2E ; cohort 2: S2F and S2G) for inflamed CD56 dim (cohort 1: CD69 high , CD38; cohort 2: CD38, CD69), proliferating CD56 dim (cohort 1: KI-67, HLA-DR), CD56 bright (both cohorts: high CD56 expression, no CD16), HLA hi CD56 dim (cohort 1: TIGIT high , CD95, HLA-DR ; cohort 2: HLA-DR), and cytokine CD56 dim cells (cohort 1: CD161 high , CD38, CD16 high cohort 2: KLRG1 high , CD161 high , CD16 high CD226). In cohort 2, the proliferating CD56 dim subset could not be assigned. Conversely, subclusters were found that did not match the subsets. This can be explained by missing markers, which would be necessary to uniquely assign all NK cell subsets.

Data variance was determined by controlling the False Discovery Rate for multiple comparisons following one-way ANOVA in GraphPad PRISM 9 (Graphpad software, Inc., La Jolla, CA). From a case number above 8 a normal distribution test according to D'Àgostino-Pearson was applied, below 8 individuals non-parametric tests were performed. Therefore, for unpaired/normally distributed data ordinary ANOVA and for unpaired/nonnormally distributed data Kruskal-Wallis test was performed. A p value below the limit of 0.05 was considered significant, and figures were produced using GraphPad Prism.

To find similarly DEGs in the later stage of severe COVID-19, a rank-rank analysis for the three cohorts was performed. First, DEGs for the comparison of severe week 3+ vs. all other conditions (cohort 1, cohort 2 and cohort 3), and severe week 2+ vs. all other conditions (cohort 4) were calculated using the Seurat implemented "FindMarkers" function defined by an expression in at least 10% of NK cells. Here, no fold-change cut-off was applied. Foldchanges of genes present in both cohort 1 and cohort 2, cohort 3 or cohort 4, respectively, were subsequently visualized in a rank-rank plot. Intersection of genes which were either upor downregulated in all three cohorts were highlighted and visualized in a heatmap representing their fold-change. For the definition of intersecting up-or downregulated genes the following log fold-change thresholds were used; cohort 1 = 0.15, cohort 2 = 0.2, cohort 3 = 0.15 and cohort 4 = 0.15.

To obtain scRNA-seq data of NK cells from the bronchoalveolar lavage fluid (BALF) in COVID-19, the sequencing data from Liao et al. (Liao et al., 2020) was downloaded from GEO with the accession number GSE145926. The authors used the droplet-based 10X Genomics technique and deposited the count matrices in h5 files. These files were imported into R 4.0.0 by using the Seurat-implemented "Read10X_h5'' function. Subsequent analysis was performed using Seurat (version 3.2.0). After removal of cells with less than 200 or more than 6,000 expressed genes and more than 10% mitochondrial reads, a total of 66,452 BALF cells remained. In addition, genes that were expressed in less than three cells were removed. Subsequently normalization was performed using the LogNormalization function. The gene counts for each cell were normalized by total UMI counts, multiplied by 10,000 (TP10K) and the log transformed by log10(TP10k+1). After normalization, the count data was scaled regressing for total UMI counts and mitochondrial read percentage, as described by Liao et al. and principal component analysis (PCA) was performed based on the 2,000 most variable features identified using the vst method implemented in Seurat. Since a batch-effect was observed for the different samples, the data was batch-corrected using the "harmony" algorithm (Korsunsky et al., 2019) based on the first 50 principal components. For twodimensional data visualization, UMAP was performed based on the first 50 principal components of the "harmony" data reduction. Subsequently, the cells were clustered using the Louvain algorithm based on the first 30 "harmony" dimensions with a resolution of 0.7, resulting in 19 clusters. For NK cell selection, the cluster expressing canonical NK cell markers (GNLY, NGK7, KLRF1) was sub-clustered. In the following analysis the NK cell population was cleaned from ambiguous cells as described for our cohorts. Cells originating from control donors had to be removed due to very low numbers of cells. The final BALF NK cell dataset contained 658 cells.

The processed scRNA-seq dataset (n=114,396 cells) from lung biopsies of patients with pulmonary fibrotic diseases and control by Habermann et al. (Habermann et al., 2020) was downloaded from GEO with the accession number GSE135893 and loaded into R 4.0.0 for analysis using Seurat. NK cell selection occured in a similar fashion as described for our cohorts. In brief, the entire T lymphocyte and NK cell fraction as defined by Habermann et al. was sub-clustered, next the cluster composed of NK cells was selected and finally, NK cells were cleaned from ambiguous cells. Furthermore, cells originating from patients with unclassifiable interstitial lung disease, sarcoidosis and chronic hypersensitivity pneumonitis were removed for further analysis. The resulting final NK cell subset comprised 1,550 cells from patients with nonspecific interstitial pneumonia (NSIP), idiopathic pulmonary fibrosis (IPF) and controls. For determination of statistical significant differences in the distribution per disease, a Kurskall Wallis test with FDR correction was performed. Subtypes showing significant changes (FDR-corrected Kruskal-Wallis p value < 0.05) were further tested with a Dunn's Post-Hoc test using the "dunn.test" R package (version 1.3.5). Resulting p values were corrected for multiple testing using the Benjamini-Hochberg method.

For data visualization the R packages Seurat, ggplot2 (version 3.3.2) (Wickham H (2016). ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York. ISBN 978-3-319-24277-4), (https://ggplot2.tidyverse.org), pheatmap (version 1.0.12) and ComplexHeatmap (version 2.4.3) (Gu et al., 2016) were used.

Supplementary Figure 2A+D+C IFITM1  HLA−C  B2M  LY6E  IFITM3  ELF1  IRF1  CD74  IRF2  LPAR6  UBE2L6  HERC6  LAP3  PARP14  SAMD9L  ISG20  IFIH1  STAT2  TAP1  TRIM14  PSMB9  IFIT2  PLSCR1  MX1  IFI44L  OAS1  BST2  RSAD2  EPSTI1  OASL  IFI44  EIF2AK2  CMPK2  SP110  SAMD9  PARP9  IRF7  ISG15  IFIT3  DDX60  ADAR  TRIM25  HELZ2  IFI35  GBP4  SELL  PSME2  CASP1 BTG1  DUSP1  CCL5  IER2  JUNB  PNRC1  EIF1  LITAF  ATP2B1  PPP1R15A  KLF6  FOS  NFKBIA  CCL4  TIPARP  PTGER4  BTG2  NR4A2  ZC3H12A  FOSL2  BHLHE40  IRS2  MAP3K8  ZFP36  CEBPB  TNFAIP3  CD69  IRF1  DUSP2  KLF2  IFIT2  FOSB  TAP1  DDX58  IFIH1  IL7R  CEBPD  AREG  TNFAIP8 AREG  ZFP36L2  CXCR4  DUSP2  TSC22D3  IFI6  IFI44L  XAF1  IFITM3  MX1  CX3CR1  ISG15 

Critical Role of Natural Killer Cells in Lung Immunopathology During Influenza Infection in Mice

Early changes in natural killer cell function indicate virologic response to interferon therapy for hepatitis C

SCENIC: single-cell regulatory network inference and clustering

Chemokines cooperate with TNF to provide protective anti-viral immunity and to enhance inflammation

Gene Ontology: tool for the unification of biology

Interferon-induced transmembrane protein 3 is a type II transmembrane protein

Differential responses of immune cells to type I interferon contribute to host resistance to viral infection

Autoantibodies against type I IFNs in patients with life-threatening COVID-19

The National COVID Cohort Collaborative: Clinical Characterization and Early Severity Prediction

Rapid expansion and long-term elevated NK cell numbers in humans infected with hantavirus

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

NK Cell Responses to Human Tick-Borne Encephalitis Virus Infection

Stellate Cells, Hepatocytes, and Endothelial Cells Imprint the Kupffer Cell Identity on Monocytes Colonizing the Liver Macrophage Niche

Neutrophils sense microbe size and selectively release neutrophil extracellular traps in response to large pathogens

SARS-CoV-2-reactive T cells in healthy donors and patients with COVID-19

NicheNet: modeling intercellular communication by linking ligands to target genes

Integrating singlecell transcriptomic data across different conditions, technologies, and species

Type I IFN and TNFα cross-regulation in immune-mediated inflammatory disease: basic concepts and clinical relevance

Clinical and immunological features of severe and moderate coronavirus disease 2019

Cytofkit: A Bioconductor Package for an Integrated Mass Cytometry Data Analysis Pipeline

COVID-19 severity correlates with airway epithelium-immune cell interactions identified by single-cell analysis

Global absence and targeting of protective immune states in severe COVID-19

Detection of 2019 novel coronavirus (2019-nCoV) by real-time RT-PCR

High-Dimensional Single-Cell Analysis Identifies Organ-Specific Signatures and Conserved NK Cell Subsets in Humans and Mice

Bone Marrow CD11c+ Cell-Derived Amphiregulin Promotes Pulmonary Fibrosis

A robust pipeline for rapid production of versatile nanobody repertoires

The Gene Ontology Resource: 20 years and still GOing strong

Complex Immune Dysregulation in COVID-19 Patients with Severe Respiratory Failure

Impaired CD4+ T cell stimulation of NK cell anti-fibrotic activity may contribute to accelerated liver fibrosis progression in HIV/HCV patients

Baseline Characteristics and Outcomes of 1591 Patients Infected With SARS-CoV-2 Admitted to ICUs of the Lombardy Region

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

Complex heatmaps reveal patterns and correlations in multidimensional genomic data

Single-cell RNA sequencing reveals profibrotic roles of distinct epithelial and mesenchymal lineages in pulmonary fibrosis

Clinical Characteristics of Coronavirus Disease 2019 in China

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

flowCore: a Bioconductor package for high throughput flow cytometry

Tumor necrosis factor (TNF)-alpha and TNF receptors in viral pathogenesis

Create Elegant Data Visualisations Using the Grammar of Graphics

Clinical features of patients infected with 2019 novel coronavirus in Wuhan

Positive regulation of immune cell function and inflammatory responses by phosphatase PAC-1

COVID-19 pneumonia: CD8+ T and NK cells are decreased in number but Immunology 218

Cytokines and Cell Surface Molecules Independently Induce CXCR4 Expression on CD4+ CCR7+ Human Memory T Cells

Toward understanding the origin and evolution of cellular organisms

Synergism of TNF-α and IFN-γ Triggers Inflammatory Cell Death, Tissue Damage, and Mortality in SARS-CoV-2 Infection and Cytokine Shock Syndromes

Structure-guided multivalent nanobodies block SARS-CoV-2 infection and suppress mutational escape

An effective interferon-gammamediated inhibition of hepatitis C virus replication by natural killer cells is associated with spontaneous clearance of acute hepatitis C in human immunodeficiency virus-positive patients

Storm of soluble immune checkpoints associated with disease severity of COVID-19

Fast, sensitive and accurate integration of single-cell data with Harmony

Comprehensive mapping of immune perturbations associated with severe COVID-19

Studying the pathophysiology of coronavirus disease 2019: a protocol for the Berlin prospective COVID-19 patient cohort (Pa-COVID-19)

Inflammatory monocytes require type I interferon receptor signaling to activate NK cells via IL-18 during a mucosal viral infection

Type I Interferon Receptor on NK Cells Negatively Regulates Interferon-γ Production

Data-Driven Phenotypic Dissection of AML Reveals Progenitor-like Cells that Correlate with Prognosis

Single-cell landscape of bronchoalveolar immune cells in patients with COVID-19

The Molecular Signatures Database Hallmark Gene Set Collection

NF-κB signaling in inflammation

Longitudinal analyses reveal immunological misfiring in severe COVID-19

Induction of interferon-γ from natural killer cells by immunostimulatory CpG DNA is mediated through plasmacytoid-dendritic-cell-produced interferon-α and tumour necrosis factor-α

Natural killer cell immunotypes related to COVID-19 disease severity

UMAP: Uniform Manifold Approximation and Projection

uwot: The Uniform Manifold Approximation and Projection (UMAP) Method for Dimensionality Reduction

Chronic stimulation drives human NK cell dysfunction and epigenetic reprograming

High basal STAT4 balanced by STAT1 induction to control type 1 interferon effects in natural killer cells

Innate lymphoid cells promote lung-tissue homeostasis after infection with influenza virus

Detection of SARS-CoV-2-Specific Humoral and Cellular Immunity in COVID-19 Convalescent Individuals

A Dynamic Immune Response Shapes COVID-19 Progression

Natural Killer Cells Limit Cardiac Inflammation and Fibrosis by Halting Eosinophil Infiltration

SARS-CoV-2 induces human plasmacytoid pre-dendritic cell diversification via UNC93B and IRAK4

Impaired natural killer cell counts and cytolytic activity in patients with severe COVID-19

Natural Killer Cells Ameliorate Liver Fibrosis by Killing Activated Stellate Cells in NKG2D-Dependent and Tumor Necrosis Factor-Related Apoptosis-Inducing Ligand-Dependent Manners

The Promise and Peril of Natural Killer Cell Therapies in Pulmonary Infection

Translational repression of pre-formed cytokine-encoding mRNA prevents chronic activation of memory T cells

Regulatory cell therapy in kidney transplantation (The ONE Study): a harmonised design and analysis of seven nonrandomised, single-arm

Characterization of pre-existing and induced SARS-CoV-2-specific CD8 + T cells

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

Next-Generation Sequencing of T and B Cell Receptor Repertoires from COVID-19 Patients Showed Signatures Associated with Severity of Disease

COVID-19 and the human innate immune system

Tumor Necrosis Factor Alpha Exerts Powerful Anti-Influenza Virus Effects in Lung Epithelial Cells

Diversity of peripheral blood human NK cells identified by single-cell RNA sequencing

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

Comprehensive Integration of Single-Cell Data

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

Optimized qRT-PCR Approach for the Detection of Intra-and Extra-Cellular SARS-CoV-2 RNAs

Immunology of COVID-19: Current State of the Science

Unique immunological profile in patients with COVID-19

Clinical Characteristics of 138 Hospitalized Patients With 2019 Novel Coronavirus-Infected Pneumonia in

Characteristics of Peripheral Lymphocyte Subset Alteration in COVID-19 Pneumonia

IL-4 and a glucocorticoid up-regulate CXCR4 expression on human CD4+ T lymphocytes and enhance HIV-1 replication

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

Factors associated with COVID-19-related death using OpenSAFELY

The differential immune responses to COVID-19 in peripheral and lung revealed by single-cell RNA sequencing

Heterogeneity of human bone marrow and blood natural killer cells defined by single-cell transcriptome

Stress-glucocorticoid-TSC22D3 axis compromises therapy-induced antitumor immunity

Cell-Type-Specific Immune Dysregulation in Severely Ill COVID-19 Patients

clusterProfiler: an R Package for Comparing Biological Themes Among Gene Clusters

DNA vaccine protection against SARS-CoV-2 in rhesus macaques

High-mobility group box 2 promoted proliferation of cervical cancer cells by activating AKT signaling pathway

Inborn errors of type I IFN immunity in patients with life-threatening COVID-19

Clinical course and risk factors for mortality of adult inpatients with COVID-19 in Wuhan, China: a retrospective cohort study

KEGG: Ribosome (RPL35

GO: Response to interferon-gamma (IRF1, IFITM1, IFITM3

GO: DNA replication

GO: Response to interferon-gamma (CCL5, IFITM1, CCL3, ...) KEGG: Ribosome (RPL11, RPS15, RPL21

KEGG: Antigen processing and presentation

KEGG: Natural killer cell mediated cytotoxicity

GO: Hydrogen peroxide metabolic process (PRDX3, HBB, HBA1) KEGG: Glycolysis / Gluconeogenesis (GAPDH, PKM, LDHB) GO: production of molecular mediator of immune response (CD226

We would like to acknowledge the assistance of the Flow Cytometry Core Facility at the Institute of Experimental Immunology, Medical Faculty at the University of Bonn. We are grateful to the patients and donors volunteering to participate in this study making this research possible. Berthold Gottgens, John Marioni, Menna Clatworthy, Sarah Teichmann, Paul Lyons, Ken Smith, Fernando Calero-Nieto and the Cambridge Institute of Therapeutic Immunology and