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

Understanding the mechanism that leads to immune dysfunction in severe coronavirus disease 2019 (COVID-19) is crucial for the development of effective treatment. Here, using single-cell RNA sequencing, we characterized the peripheral blood mononuclear cells (PBMCs) from uninfected controls and COVID-19 patients and cells in paired broncho-alveolar lavage fluid (BALF). We found a close association of decreased dendritic cells (DCs) and increased monocytes resembling myeloid-derived suppressor cells (MDSCs), which correlated with lymphopenia and inflammation in the blood of severe COVID-19 patients. Those MDSC-like monocytes were immune-paralyzed. In contrast, monocyte-macrophages in BALFs of COVID-19 patients produced massive amounts of cytokines and chemokines, but secreted little interferons. The frequencies of peripheral T cells and NK cells were significantly decreased in severe COVID-19 patients, especially for innate-like T and various CD8+ T cell subsets, compared to healthy controls. In contrast, the proportions of various activated CD4+ T cell subsets among the T cell compartment, including Th1, Th2, and Th17-like cells were increased and more clonally expanded in severe COVID-19 patients. Patients’ peripheral T cells showed no sign of exhaustion or augmented cell death, whereas T cells in BALFs produced higher levels of IFNG, TNF, CCL4, CCL5, etc. Paired TCR tracking indicated abundant recruitment of peripheral T cells to the severe patients’ lung. Together, this study comprehensively depicts how the immune cell landscape is perturbed in severe COVID-19.


Introduction
Coronavirus disease 2019 (COVID-19) caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) has been spreading rapidly worldwide, causing serious public health crisis. Although most SARS-CoV-2-infected cases have asymptomatic or mild-to-moderate diseases, around 10% of those infected may develop severe pneumonia and other associated organ malfunctions 1 . Old age, male sex, and underlying comorbidities are risk factors for causing severe COVID-19 2 ; however, the pathogenesis mechanisms remains unclear. Immune perturbations were thought to play crucial roles in the COVID-19 pathogenesis. Indeed, many reports showed that lymphopenia and increased blood cytokine levels were closely associated with the development or recovery of severe COVID-19 and proposed to treat the patients by inhibiting cytokine storms [3][4][5][6] .
Recent studies have revealed more aspects of different immune players in COVID-19 infections. Although SARS-CoV-2 infection in host cells elicits robust secretion of chemokines and cytokines, it only weakly produces interferons (IFNs) 7 . In our earlier studies, we showed increased recruitment of highly inflammatory FCN1 + monocyte-derived macrophages in patients' bronchoalveolar lavage fluids (BALFs), suggesting their implication in the cytokine storm 8,9 . Intriguingly, T and B cell responses and neutralizing antibodies recognizing SARS-CoV-2 were detected among most of the infected patients, with higher levels in those with old age and severe diseases [10][11][12] . The dysregulation of other immune cell compartments, such as monocytes, dendritic cells (DCs), and innate-like T cells is also likely present in severe COVID-19 [13][14][15][16] .
To unravel the barely known immune mechanisms underlying COVID-19 pathogenesis in an unbiased and comprehensive manner, we and others have applied the single-cell RNA sequencing (scRNA-seq) to profile the immune cell heterogeneity and dynamics in BALFs, blood, and respiratory tract samples from the COVID-19 patients. Collectively, those studies revealed a stunted IFN response, depletion of natural killer (NK) and T lymphocytes, loss of major histocompatibility class (MHC) class II molecules, and significantly elevated levels of chemokine production from monocytes in the patients 8,9,13,14,[17][18][19] . However, a complete picture of the COVID-19-induced immune perturbation has not been generated. Here we conducted scRNA-seq analysis of paired blood and BALF samples from the same COVID-19 patients. Our data revealed profound alterations of various immune compartments and depicted a dichotomy of peripheral immune paralysis and broncho-alveolar immune hyperactivation in COVID-19 patients.

Results
The overview of dysregulated peripheral immune landscape in severe COVID-19 patients A high-quality scRNA-seq dataset composed of 200,059 cells were generated that characterized peripheral immune cells from three healthy controls (HC), five mild, and eight severe COVID-19 patients (Fig. 1a). The metadata of these patients is listed in Supplementary  Table S1. Patients with mild diseases were all cured and discharged after 11-18 days of hospitalization, 2 of the 8 patients with severe diseases succumbed, whereas other severe cases recovered after 12-58 days of hospitalization. The clinical course of this patient cohort is similar to earlier reports, where elderly patients with underlying diseases were prone to develop severe symptoms and showed higher mortality 1,2,20 . In addition, the plasma levels of interleukin (IL)-6 and C-reactive proteins (CRP) in severe patients were higher, while lymphocyte counts were reduced, indications of both cytokine storm and lymphopenia.
The clustering analysis showed 29 clusters and 10 major cell types annotated by marker genes, including T cell (CD3D), NK cell (KLRF1), B cell (CD79A), monocyte (CD14, FCGR3A), myeloid DCs (mDCs) (CD1C), plasmacytoid DCs (pDCs) (IL3RA), and plasma cells (PCs) (IGKC), megakaryocyte (MYL9), cycling cells (MKI67), and erythrocyte (HBB) (Fig. 1b, c and Supplementary Fig.  S1a, b). By visual inspection, the data integration was efficient and showed no significant batch effect (Supplementary Fig. S1c). Erythrocytes were not included in subsequent analysis and cycling cells were reclustered into cycling T, cycling PC, and cycling NK cells based on specific markers ( Supplementary Fig. S1d, e). This dataset indicated significantly dysregulated peripheral immune landscapes in COVID-19 patients compared to HC, especially among severe cases (Fig. 1d). The most prominent changes included an expansion of monocyte and cycling T cells and a reduction of NK, T, and mDC populations, thus resulting in largely increased monocyte/ T cell ratios in COVID-19 patients (Fig. 1e, f). The frequency of pDCs was also decreased in severe COVID-19, although the difference was not statistically significant. Together, these data show that SARS-CoV-2 infection greatly perturbs the blood immune cell compartments, particularly in those with severe diseases.

Remodeling of circulating myeloid cell populations in patients with COVID-19
We observed that the proportions of monocytes were increased in COVID-19 patients, especially in those with severe diseases. To further understand the remodeling of myeloid cell compartment, we re-clustered myeloid cells and identified five distinct cell types including CD14 + classic monocyte, CD14 + CD16 + intermediate monocytes, CD16 + non-classic monocytes, DC1, and DC2 ( Fig. 2a and Supplementary Fig. S2a). The composition of myeloid cells in severe COVID-19 patients differed significantly from that of mild cases and controls. Proportion of CD14 + monocyte increased significantly in severe COVID-19 compared to that in the mild COVID-19 and control group, whereas those of CD16 + non-classical monocyte (vs controls), CD14 + CD16 + monocytes (vs mild COVID- 19), and DC2 significantly (vs mild COVID-19 and controls) decreased in severe COVID-19 (Fig. 2b, c).
CD14 + monocytes represent the major peripheral myeloid cell type, and differential Uniform Manifold Approximation and Projection (UMAP) projection patterns of CD14 + monocytes between COVID-19 and controls (Fig. 2b) indicated perturbed transcriptome features. Among the differentially expressed genes (DEGs) in CD14 + monocytes, we found 116 and 134 upregulated genes in mild COVID-19 and severe COVID-19 cases compared to controls, vs only 74 upregulated genes in comparison between the two COVID-19 groups. In contrast, we found 217 and 160 downregulated genes in severe COVID-19 compared to controls or mild cases, respectively, vs only 104 downregulated genes in comparison between mild COVID-19 and controls (Supplementary Fig. S2b and Table S2). Gene Ontology (GO) terms of upregulated DEGs included response to virus, type I IFN, and IFN-γ in both the COVID-19 groups, and neutrophil activation and energy metabolism pathways in severe COVID-19 (Fig. 2d). Unexpectedly, GO analysis of downregulated DEGs indicated deficient  c Comparisons of percentages of each peripheral myeloid cell types between the two COVID-19 groups and controls (two-sided Student's t-test, *P < 0.05, **P < 0.01). d Enrichment of GO biological process (BP) terms for upregulated genes (left) and downregulated genes (right) in blood CD14 + monocyte comparisons between mild and HC (M vs H), severe and HC (S vs H), and severe and mild groups (S vs M) (representative terms are shown, adjusted P < 0.01 as indicated by the colored bar). e The heatmaps show the selected differentially expressed genes and their associated GO terms as indicated in d (logFC > 0.41 or < -0.41, adjusted P < 0.01). Mean.Exp, average expression. f Density plots show the composite MHC II signature scores and calprotectin signature scores of peripheral CD14 + monocytes in 2D maps. The horizontal and vertical lines separating the four quadrants represent the median scores of all CD14 + monocytes. The percentages of cells in each quadrant are indicated. g Left panel shows the representative flow cytometric data of HLA-DR expression on CD14 + and CD14 -PBMCs. Right plot shows the summarized data from more subjects (two-sided Student's t-test). h The Pearson correlation of "MDSC-like signature score" and plasma CRP, IL-6 levels, blood neutrophil, CD3 + , CD4 + , and CD8 + T cell counts. monocyte functions mainly in severe COVID-19 cases, such as decreased type I IFN production, cytokine secretion, chemokine production, and antigen processing and presentation (Fig. 2d). We further examined the DEGs associated with those GO terms. Indeed, many canonical IFN-stimulating genes (ISGs), including ISG15, IFITM1, IFITM3, MX1, IRF7, IFI27, etc., were expressed at higher levels in COVID-19 patients than controls, while the genes associated with neutrophil activation, including S100A8, S100A9, S100A12, CLU, RNASE2, etc., were expressed at higher levels in severe COVID-19 than mild COVID-19 and controls (Fig. 2e). Despite of upregulated ISGs, we failed to detect higher production of type I or type III IFNs in COVID-19 than controls. Regarding the downregulated geneset, MHC II molecules, including HLA-DQA1, HLD-DRA, HLA-DRB1, HLA-DMB, HLA-DMA, etc., and cytokine/chemokine genes, including IL1B, TNF, CCL3, CCL4, and CXCL8 were expressed at lower levels in COVID-19 patients, especially those with severe diseases (Fig. 2e). Thus those upregulated DEGs in CD14 + monocytes from COVID-19 patients reflect the immune response to SARS-CoV-2 infection, while the downregulated DEGs in CD14 + monocytes from patients with severe COVID-19 suggest an immune paralyzed status of those cells.
Myeloid-derived suppressor cells (MDSCs) are a population of heterogeneous immature myeloid cells expanded during inflammatory conditions and could suppress T cell responses 21,22 . In peripheral blood, monocytic MDSCs have the phenotype CD14 + HLA-DR -/lo , whereas monocytes are HLA-DR positive 23,24 . Downregulation of MHC II molecules, increased calprotectin (S100A8 and S100A9), and immune-suppressive functions are reported features of MDSCs. Indeed, by their unique composite scores of MHC II molecules (lower levels) and calprotectin (higher levels) vs those scores in mild COVID-19 and controls, we identified that the monocytes in severe COVID-19 highly resembled MDSCs (Fig. 2f). Decreased levels of HLA-DR in CD14 + monocyte from severe patients were further validated by flow cytometry (Fig. 2g) and also reported by other studies 14,18 . Intriguingly, MDSC-like scores in our study positively correlate with serum CRP, IL-6 levels, and neutrophil-to-lymphocyte ratio and negatively correlate with decreased blood CD3 + , CD4 + , and CD8 + T cell counts (Fig. 2h).
Together, our scRNA-seq characterization revealed a multifaceted remodeling of peripheral myeloid compartment in COVID-19 patients. While the circulating monocytes in COVID-19 patients are featured by heightened ISG responses, they produce little IFNs, cytokines, and chemokines. Furthermore, the loss of DCs and emergence of MDSC-like monocyte suggest their involvements in the immune paralysis of severe COVID-19 patients.

Abnormally activated lung monocyte-macrophages in severe COVID-19
Recently, we discovered the aberrant activation of BALF monocyte-macrophages in severe COVID-19 8,9 . Here, to further understand the connection between the lung monocyte-macrophages and their blood counterparts and assess their differential roles in COVID-19, we studied paired BALF and blood samples from two mild and five severe patients. The integration analysis of BALF and circulating myeloid cells showed clusters of neutrophil (FCGR3B), mDC (CD1C), monocyte-macrophages (CD14, FCGR3A, and CD68) ( Fig. 3a and Supplementary Fig.  S3a). Macrophage subset classification markers, including FCN1, SPP1, and FABP4, were differentially expressed by circulating and BALF monocyte-macrophages from patients with mild or severe COVID-19 ( Fig. 3b). Analysis of differentiation trajectory of circulating and BALF monocyte-macrophages from the same patient revealed a consensus blood-toward-BALF course ( Fig. 3c and Supplementary Fig. S3b), consistent with the recruitment of peripheral monocytes into inflammatory tissues as expected, although the effect of lung immune microenvironment on the exact differentiation trajectory still requires investigation and validation.
Next, we performed transcriptome analysis of circulating and BALF monocyte-macrophages to understand their functional status. Among the DEGs, there were 524 shared upregulated genes and 501 downregulated genes in BALF monocyte-macrophages vs those in blood, identified from both mild and severe COVID-19 patients ( Fig. 3d and Supplementary Table S3). Such a large number of DEGs suggested that significant difference existed between the peripheral and lung monocytemacrophages. Indeed, GO analysis revealed broad activation of multiple immune pathways in BALF monocytemacrophages, including response to IFNs and cytokines, neutrophil activation, and leukocyte migration, while the pathway involved in myeloid cell differentiation, ATP metabolism, etc. were enriched in blood monocytes (Fig.  3e). In addition, these comparisons revealed perturbed pathways in BALF monocyte-macrophages relevant to severe COVID-19, including responses to hypoxia, high temperature, metal ion, wounding, and Fc receptor signaling pathways were especially upregulated (Fig. 3e), while pathways related to alveoli macrophage functions were downregulated, including lipid metabolism, apoptotic cell clearance, and antigen presentation (Fig. 3e). The representative DEGs involved in those pathways are shown in Supplementary Fig. S3c.
Monocyte-macrophages were thought to play key roles in driving the cytokine storm underlying the development of severe COVID-19 25 . Therefore, we examined the cytokine and chemokine levels in monocyte-macrophages in paired blood and BALF samples from the same patient.
We found that all types of IFNs (IFNAs, IFNB, IFNG, and IFNLs) were minimally expressed by monocyte-macrophages, whereas cytokines (IL1A, IL1B, IL1R2, IL1RN, IL18, IL6, TNF, IL10, and TGFB1) and multiple chemokines were highly expressed in monocyte-macrophages from BALFs but not in those from paired blood samples (Fig. 3f). These data suggest that, although monocyte-macrophages in BALF are activated, they are immune silent in periphery, possibly due to continuous engagement of viral stimulation or the pro-inflammatory milieu in the lung.
Intriguingly, we observed relatively higher levels of antiinflammatory cytokines (IL1R2, IL1RN, and TGFB1) and The expression of monocyte-macrophage markers FCN1, SPP1, and FABP4 were projected to UMAP from a. PM, peripheral cells of mild cases; PS, peripheral cells of severe cases; BM, BALF of mild cases; BS, BALF of severe cases. c Differentiation trajectory of the blood monocytes and BALF monocyte-macrophages from two representative COVID-19 patients, analyzed independently. d Venn diagram shows the number of upregulated and downregulated DEGs in monocyte-macrophage comparisons as indicated. logFC > 0.41 or < -0.41, adjusted P < 0.01. e Enrichment of GO biological process (BP) terms for upregulated genes (left) and downregulated genes (right) in monocyte-macrophage comparisons as indicated. Selected terms are shown, adjusted P value is indicated by the colored bar. f Heatmaps show the expression of selected interferon, cytokine, and chemokine genes in paired blood and BALF monocytemacrophages derived from the same patient. Stars indicate that the genes are differentially expressed in BALF monocyte-macrophages between mild and severe COVID-19. Purple and green stars show that the gene expression are significantly upregulated in severe COVID-19 and mild COVID-19 groups, respectively (MAST; P < 0.01). B, BALF samples; P, PBMC samples. g The levels of selected cytokines and chemokines in paired BALF and plasma samples were measured by CBA (two-sided Wilcoxon test between BALF and PBMC of severe patients).
lower levels of IL18 in BALF monocyte-macrophages from severe COVID-19 than mild cases, whereas classical pro-inflammatory cytokines (IL1A, IL1B, IL6, and TNF) were comparable between the two groups (Fig. 3f). In contrast, as shown in our earlier studies, monocyte-and neutrophil-recruiting chemokines (CCL2, CCL3, CCL4, CCL7, CCL8, CXCL1, CXCL2, CXCL3, and CXCL8) were highly expressed, whereas chemokines (CXCL9 and CXCL16) recruiting T cells were less expressed by monocyte-macrophages in BALFs of severe COVID-19 than those in mild cases (Fig. 3f). The higher levels of cytokines (IL-1β, IL-6, etc.) and IL-8 in BALFs than paired plasma was further validated at the protein levels, particularly exemplified by the extremely high levels of IL-8 in BALFs (Fig. 3g). Thus these paired analyses revealed the involvement of tissue monocyte-macrophages in cytokine storms during severe COVID-19, particularly through producing chemokines and recruiting more monocytes and neutrophils, but less likely attributed to excessive production of pro-inflammatory cytokines.
Cell density UMAP projections revealed an obviously perturbed T cell landscape in severe COVID-19 compared to the mild COVID and control groups (Fig. 4c). Within the T cell compartment, percentages of innate-like T cells, including MAIT and NKT cells, were significantly lower in severe COVID-19 than those in mild COVID-19. Percentages of CD8-Naive, CD8-GZMK, and CD8-GZMB subsets were also lower in severe COVID-19 patients than those in mild cases, although the difference in CD8-GZMB comparison was not statistically significant ( Fig.  4d and Supplementary Fig. S4b). In contrast, the percentages of several CD4 + T cell subsets, including CD4-Naive, CD4-LTB, CD4-ICOS, Treg-CTLA4, as well as cycling T cells, were significantly increased in severe COVID-19 than those in mild COVID-19. The percentages of CD4-GATA3 and CD4-CCR6 also showed an increasing trend in severe COVID-19 (Fig. 4d). In addition, single-cell TCR (sc-TCR) analysis revealed increased clonal expansion levels in several CD4 + but not CD8 + T cell subsets in severe COVID-19 than those in mild cases (Fig. 4e, f). Consistently, TCR sharing analysis among different T cell subsets revealed much more active interchange between different CD4 + T cell subsets and cycling T cells but not among CD8 + T cell subsets during severe COVID-19 (Fig. 4g). Together, our data revealed the preferential activation of CD4 + T cell responses but significant depletion of multiple innate-like T cell and CD8 + T cell subsets in peripheral blood as the featured T cell perturbations in severe COVID-19. To further explore the clues of CD8 + T cell lymphopenia, we conducted transcriptome comparisons of MAIT, CD8-GZMK, and CD8-GZMB subsets between the patients and control groups (Supplementary Fig. S4c and Table S4). There was no evidence of T cell exhaustion, activation of cell death pathways, and cytokine productions in those cells from COVID-19 patients, although the pathways related to virus infection and IFN responses were identified (Supplementary Fig. S4d). Similar findings were also noticed by other groups, 16 thus exhaustion and cell death are unlikely the major causes for T cell loss during COVID-19.

Tracking T cell function and migration across peripheral blood and BALFs
We sought to study T cell movement from blood to BALFs in paired samples from COVID-19 patients. First, we integrated data of NK and T cells from peripheral blood mononuclear cells (PBMCs) and BALF. Cells were re-clustered into nine major types, including NK cells, MAIT, CD4-Naive, CD4-Tm, Treg, CD8-Naive, CD8-Tm, CD8-IL7R, and cycling T cells (Fig. 5a and Supplementary  Fig. S5a). PBMCs included plenty of naive CD4 + and CD8 + T cells, while there were fewer naive T cells in BALFs, which were mainly composed of NK cells, CD4-Tm, CD8-Tm, and cycling T cells (Fig. 5b). We performed gene expression analysis to determine the functional divergence of NK and T cells across the peripheral blood and BALFs (Supplementary Table S5). Compared with peripheral counterparts, we observed that response to virus, type I IFN, and IFN-γ were commonly activated in NK, CD4-Tm, and CD8-Tm cells from BALFs (Supplementary Fig. S5b). We also noticed higher levels of cytokines including IFNG, TNF, CSF1, TNFSF10, and TNFSF13B; chemokines including CCL3, CCL4, and CCL5; and IL-15 signaling modules including IL15RA, IL2RB, IL2RG, JAK1, and STAT3 in BALF cells (Fig. 5c). However, the levels of other T cell relevant cytokines, including IL4, IL5, IL10, IL13, IL17A, IL17F, IL21, IL25, and IL33, were not detected in either blood or BALFs (Supplementary Fig. S5c).
Next, we utilized the TCR clonotype information to track the migrating T cells in blood and BALF compartment. The T cell clonal expansion status and clonotype sharing were assessed in patients with paired blood and BALF samples (Fig. 5d, e). Interestingly, considerable proportion of BALF T cells belonging to CD8-Tm, CD8 CTL, and cycling T cells subsets could trace their clonotypes back to the paired blood counterparts, especially in those patients with severe COVID-19 (Fig. 5f). Here only two mild cases (M1 and M2) had paired blood and BALF samples for the analysis, and it showed a lower degree of clonotype shared across the two compartments (Fig. 5f). The higher levels of chemokines from the lung of c Density plots show the UMAP projection of peripheral NK and T cells from COVID-19 patients and controls. d Comparisons of percentages of each peripheral NK and T cell types between the two COVID-19 groups and controls (two-sided Student's t-test, *P < 0.05, **P < 0.01, ***P < 0.001). e Clonally expanded T cells from COVID-19 patients and controls are projected into UMAP from a. f Clonal expansion indexes of T cell subsets from COVID-19 patients and controls are separately displayed. g T cell state transition status among any two clusters is inferred by their shared TCR clonotypes. Each T cell cluster is represented by a unique color. The numbers above the bar indicate the percentages of cells sharing TCRs in those two clusters. severe COVID-19 is likely to lead to increased T cell infiltration. To assess the functional adaptation of migrating T cells, we performed transcriptional analysis between blood vs BALF T cells derived from the same T cell clonotypes (Fig. 5h). We found increased expression of ISGs, CCL4, CXCR4, CXCR6, CD69, GNLY, etc. in BALF cells than those in blood, consistent with an activated and tissue-resident phenotype (Fig. 5g and Supplementary Table S6). Together, these data revealed the active recruitment and peculiar activation of peripheral T cells in human lungs infected by SARS-CoV-2.

Discussion
The scRNA-seq has been recently been applied to study host immune response in COVID-19 by us and others 8,9,13,[17][18][19] . Although these studies helped to reveal several aspects of the COVID-19 pathogenesis, a complete picture has not yet been generated. Here we integrated scRNA-seq analysis of paired blood and BALFs and depicted unprecedented details about the altered immune cell landscape in COVID-19 patients, increasing our mechanistic and systemic understanding of COVID-19 immunopathogenesis, such as cytokine storm and lymphopenia.
IFN production is stunted in SARS-CoV-2-infected cells and COVID-19 patients 7,27 . Consistent with these earlier reports, here we found that, although ISGs in COVID-19 patients were induced in monocyte-macrophages, neither type I nor type III IFNs were produced. However, the higher levels of ISGs in monocyte-macrophages and T cells from mild cases than those in severe COVID-19 still supported a heightened IFN responses linked to resolving the diseases. It is unclear what causes the imbalance of inflammation and IFN production in COVID-19; we assume that the loss of pDCs in severe cases may partly contribute to the diminished IFN production. Notably, in contrast to their blood counterparts, monocyte-macrophages in paired BALFs produced higher levels of cytokines and chemokines, especially from those severe COVID-19 patients. Thus massive tissue-resident production of cytokines/chemokines and lack of IFN induction suggest a crucial role played by local but not by circulating monocyte-macrophages in fueling the cytokine storm during severe COVID-19 and indicate that inhibiting the local cytokine storm might be more effective than those systemic intervention strategies. Moreover, the roles of anti-inflammatory cytokines in COVID-19 pathogenesis, including IL1R2, IL1RN, IL10, and TGFB1, are unclear and should be investigated in the future.
Previously, several studies have reported that MHC II molecules were downregulated in blood monocytes 4,14,18 in patients with severe COVID-19; however, there was still no consensus on the upregulated genes. We found here that genes related to neutrophil activation, including S100A8, S100A9, and S100A12, were expressed at higher levels in monocytes from severe COVID-19 patients than those in mild cases. Interestingly, these upregulated and downregulated genes are makers of monocytic MDSCs whose frequencies are known to increase during various inflammatory conditions 24,25 . Consistent with their potential immune-suppressive status, we also noticed that genes related to cytokine and IFN production were downregulated as well. Thus, contradicting with the inflammatory role by peripheral myeloid cells in severe COVID-19, the loss of mDCs, emergence of MDSC-like monocytes, and reduced cytokine production actually suggested a peripheral immune paralysis. Since the MDSC-like scores associated with lymphopenia and inflammation markers, we speculate a crucial role of MDSC-like monocytes in dampening immune response and amplifying COVID-19 pathogenesis. Indeed, similar results and conclusions were recently published by several independent studies 28-31 when we wrote this report, further corroborating the claim.
The lymphopenia is another prominent feature of immune perturbation of severe COVID-19 2,3,16 . Here we revealed that the COVID-19-associated lymphopenia included not only depletion of CD8 + T cells but also significant loss of innate-like T cells, including MAIT, γδ T, and NKT cells, similarly reported by another study 13,32 . However, the proportion of cycling T cells is increased in COVID-19; whether this reflects SARS-CoV-2-specific T cell response or bystander activation needs further investigation. In addition, we noticed that the frequencies of CD4 + T cells among all CD3 + T lymphocytes were increased, like CD4-GATA3 (Th2), CD4-CCR6 (Th17), and CD4-ICOS (Tfh) cells. Indeed, these data are consistent with earlier reports showing that lymphopenia in severe COVID-19 affects CD4 + T cells less than CD8 + T cells 33 . CD4 + T cell subsets also showed a tendency of higher clonal expansion, suggesting their activation status. We suspected that the disturbed T cell compartments may contribute to COVID-19 immunopathogenesis, e.g., the impaired antiviral responses by innate-like T cells and CD8 + T cells, meanwhile amplifying inflammation and inducing aberrant antibody responses by CD4 + T cells. Moreover, our transcriptional analysis dispute against the evidence of cytokine production, exhaustion, or increased cell death by peripheral T cells in COVID-19 patients noted by several earlier studies 26,34 . Instead, current TCR tracking analysis suggested enhanced recruitment of peripheral CD4 + and CD8 + T cells into lung tissues in COVID-19 patients, where they were induced to made cytokines locally and likely contributed to cytokine storm and peripheral lymphopenia.
In conclusion, we comprehensively delineated the perturbed immune landscapes during SARS-CoV-2 infections from both peripheral blood and infected lungs.
These data reveal potential cellular and molecular mechanisms implicated in COVID-19 immunopathogenesis and identify a peculiar functional dichotomy, peripheral immune paralysis and broncho-alveolar immune hyperactivation, specifically in severe COVID-19.

Ethics statement
This study was conducted according to the ethical principles of the Declaration of Helsinki. Ethical approval was obtained from the Research Ethics Committee of Shenzhen Third People's Hospital (2020-207). All participants provided written informed consent for sample collection and subsequent analyses.

Patients
Thirteen COVID-19 patients were enrolled in this study at the Shenzhen Third People's Hospital. Metadata and patients' samples were collected similarly as previously described 8 . The severity of COVID-19 was categorized to be mild, moderate, severe, and critical according to the "Diagnosis and Treatment Protocol of COVID-19 (the 7th Tentative Version)" by National Health Commission of China (http://www.nhc.gov.cn/yzygj/s7653p/202003/ 46c9294a7dfe4cef80dc7f5912eb1989.shtml).
In this study, we grouped patients with mild and moderate COVID-19 as the mild group and included those with severe and critical diseases as the severe group. Three healthy subjects were enrolled as the control group.

Real-time quantitative PCR (qRT-PCR) for SARS-CoV-2 RNA
In clinical practice, nasal swab, throat swab, sputum, anal swab or BALF could be collected for the SARS-CoV-2 nucleic acid assays. Total nucleic acid was extracted from the samples using the QIAamp RNA Viral kit (Qiagen) and qRT-PCR was performed using a China Food and Drug Administration-approved commercial kit specific for SARS-CoV-2 detection (GeneoDX). Each qRT-PCR assay provided a threshold cycle (Ct) value. The specimens were considered positive if the Ct value was ≤ 37, and otherwise it was negative. Specimens with a Ct value > 37 were repeated. The specimen was considered positive if the repeat results were the same as the initial result or between 37 and 40. If the repeat Ct was undetectable, the specimen was considered to be negative.

Immune cell isolation
For harvesting BALF cells, freshly obtained BALF was placed on ice and processed within 2 h in the Biological Safety Level 3 laboratory. After passing BALF through a 100-µm nylon cell strainer to filter out cell aggregates and debris, the remaining fluid was centrifuged and the cell pellets were re-suspended in the cooled RPMI 1640 complete medium. For PBMC isolation, immune cells from peripheral blood were isolated by Ficoll-Hypaque density gradient centrifugation protocol. For subsequent study, the isolated cells were counted in 0.4% trypan blue, centrifuged, and re-suspended at the concentration of 2 × 10 6 /mL.

Cytokines measurement by cytometric bead array
Twelve cytokines, including IL-1β, IL-6, IL-8, etc., were detected according to the instruction (Uni-medica, Shenzhen, China, Cat. No. 503022). In brief, the supernatant was taken from BALF after 10-min centrifugation at 1000× g. Afterwards, 25 µL Sonicate Beads, 25 µL BALF supernatant or plasma, and 25 µL of Detection Antibodies were mixed and placed on a shaker at 500× rpm for 2 h at room temperature. Then 25 µL of SA-PE was added to each tube directly. The tubes were then placed on a shaker at 500× rpm for 30 min. The data were obtained by flow cytometry (Canto II, BD) and were analyzed use LEGENDplex v8.0 (VigeneTech Inc.).

scRNA-seq library construction
The scRNA-seq libraries were prepared with the Chromium Single Cell VDJ Reagent Kits (10× Genomics; PN-1000006, PN-1000014, PN-1000020, PN-1000005) following the manufacturer's instruction. Briefly, Gel bead in Emulsion (GEM) are generated by combining barcoded Gel Beads, a Master Mix containing 20,000 cells, and Partitioning Oil onto Chromium Chip B. Reverse transcription takes place inside each GEM, after which cDNAs are pooled together for amplification and library construction. The resulting library products consist of Illumina adapters and sample indices, allowing pooling and sequencing of multiple libraries on the next-generation short read sequencer.
scRNA-seq data processing, cell clustering, and dimension reduction We aligned the sequenced reads against GRCh38 human reference genome by Cell Ranger (version 3.1.0, 10× genomics). To remove potential ambient RNAs, we used the remove-background function in CellBender 35 , which removes ambient RNA contamination and random barcode swapping from the raw UMI-based scRNA-seq data. Quality of cells were further assessed by the following criteria: (1) the number of sequenced genes is 200 to 6000; (2) the total number of UMI per cells is > 1000; (3) the percentage of mitochondrial RNA is < 15% per cell.
Data integration, cell clustering, and dimension reduction were performed by Seurat (version 3) 36 . First, we identified 2000 highly variable genes (HVGs), which were used for the following analysis using FindVariableFeatures function. Next, we integrated different samples by Inte-grateData function, which eliminates technical or batch effect by canonical correlation analysis (CCA). Using those HVGs, we calculate a principal component analysis (PCA) matrix with the top 50 components by RunPCA function. The cells were then clustered by FindClusters function after building nearest neighbor graph using FindNeighbors function. The cluster-specific marker genes were identified by FindMarkers function using MAST algorithm. The clustered cells were then projected into a two-dimension space for visualization by a nonlinear dimensional reduction method RunUMAP in Seurat package.

Integrated analysis of peripheral myeloid, NK, and T cells
For cells in PBMCs, we integrated the myeloid compartment including mDC and monocytes or NK and T cells using the similar aforementioned procedure. We re-clustered the myeloid or NK and T cells using the top 20 dimensions of PCA with default parameters. To obtain high-resolution cell clusters for each subset, we set the parameter resolution to 1.2. The cell clusters were annotated by canonical markers.

Integrated analysis of myeloid, NK, and T cells from PBMCs and BALF
Myeloid or NK and T cells from PBMCs and BALF were integrated separately. For myeloid cells, we extracted macrophage and mDC cells in BALF and monocyte and mDC in PBMCs from the corresponding raw count matrix. The extracted cells were integrated using CCA in Seurat (version 3) as mentioned above. For clustering, the resolution parameter was set to 0.6. Similarly, we extracted NK and T cells in BALF and in PBMCs from the corresponding raw count matrix. The extracted cells were integrated using CCA in Seurat (version 3) as mentioned above. For clustering, the resolution parameter was set to 1.5.

Analysis of DEGs
FindMarkers function in Seurat (version 3) with MAST algorithm was used to analyze DEGs. For each pairwise comparison, we run FindMarkers function with parameters of test.use='MAST'. The overlap of differentiated expressed genes among different comparisons was shown by Venn diagram. Genes were defined as significantly upregulated if average natural logarithm fold change (logFC) > 0.25 and adjusted P < 0.01. The genes with logFC < −0.25 and adjusted P < 0.01 were considered significantly down regulated. The selected genes shown on the heatmap have logFC < −0.41 or > 0.41 and adjusted P < 0.01.

GO annotation
ClusterProfiler 37 in R was used to perform GO term enrichment analysis for the significantly upregulated and downregulated genes. Only GO term of Biological Process was displayed.

Pseudo-time trajectory analysis
Trajectories analysis was performed using slingshot 38 for monocyte-macrophages and T cells separately. For T cells in PBMC, naive CD4 and CD8 T cells were set as the start point for CD4 + T cells and CD8 + T differentiation trajectory, respectively. For integrated analysis of monocyte-macrophages in PBMCs and BALF, we deduced the cell trajectory for each individual using the peripheral monocytes as the start point.

Sc-TCR analysis
The amino acid and nucleotide sequence of TCR chains were assembled and annotated by cellranger vdj function in CellRanger (version 3.1.0). Only cells with paired TCRα and TCRβ chains were included in clonotype analysis. Cells sharing the same TCRα-and TCRβ-CDR3 amino acid sequences were assigned to the same TCR clonotype. We assessed the TCR expansion, TCR transition among cell types, and TCR migration between PBMC and BALF using R package STARTRAC (version 0.1.0) 39 . TCR migration between PBMC and BALF were shown using the circos software 40 .

Statistics
The Student's t-test (t.test in R, two-sided, unadjusted for multiple comparisons) was used for pairwise comparisons of the cell proportions between different groups. Statistical difference of TCR expansion index and migration index, between mild and severe disease group, were calculated using t.test in R. Statistical difference of cytokine level between BALF and blood in Fig. 3g were calculated using two-sided wilcox.test in R.