Main

The underlying immune responses to severe acute respiratory syndrome (SARS)-CoV2 remain unclear. Abnormal immune responses have been reported in coronavirus animal infection models1,2 and lymphopenia is a prominent feature of severe COVID-19 (ref. 3). Bronchoalveolar lavage fluid (BALF) contains microenvironment information on bronchioles and lung alveoli. Here, we performed single-cell RNA sequencing (scRNA-seq) on BALF cells from three patients with moderate COVID-19 (M1–M3), six patients with severe/critical infection (S1–S6), three healthy controls (HC1–HC3) and a publicly available BALF (HC4)4 sample (Supplementary Tables 1 and 2). Clustering analysis showed 31 distinct clusters composed of macrophages (CD68), neutrophils (FCGR3B), myeloid dendritic cells (mDCs) (CD1C, CLEC9A), plasmacytoid dendritic cells (pDCs) (LILRA4), natural killer (NK) cells (KLRD1), T cells (CD3D), B cells (MS4A1), plasma cells (IGHG4) and epithelial cells (TPPP3, KRT18) (Fig. 1a), identified by signature genes (Extended Data Fig. 1a,b). Major cell types, including mDCs, pDCs, mast cells, NK cells, T cells and B cells, contained cells from most samples, whereas macrophages showed specific enrichment in different groups (Extended Data Fig. 1c–e). BALFs of patients with severe/critical COVID-19 infection contained higher proportions of macrophages and neutrophils and lower proportions of mDCs, pDCs and T cells than those with moderate infection (Fig. 1b and Extended Data Fig. 1f).

Fig. 1: Dysregulated bronchoalveolar immune landscapes in patients with severe COVID-19 infection.
figure 1

a, UMAP presentation of major cell types and associated clusters in BALFs (n = 13). b, The proportions of the major BALF immune cell types among healthy controls (HC) and patients with moderate (M) and severe (S) COVID-19 infection (using two-sided Student’s t-test for pairwise comparisons). c, UMAP projection of four macrophage groups among controls (n = 4) and patients (moderate, n = 3; severe, n = 6). d, Expression of FCN1, SPP1 and FABP4 by BALF macrophages from each sample.

Full size image

We re-clustered 20 macrophage subclusters to further dissect their heterogeneity (Extended Data Fig. 2a). According to recent classification criteria4 and typical markers (Extended Data Fig. 2b), the macrophages were grouped by FCN1, SPP1 and FABP4 expression patterns including the FCN1hi group 1, FCN1loSPP1+ group 2, SPP1+ group 3 and FABP4+ group 4 macrophages (Extended Data Fig. 2c). The macrophage compartments differed largely in both composition and expression of FCN1, SPP1 and FABP4 in different cell groups (Fig. 1c and Extended Data Fig. 2d). FABP4 was preferentially expressed by controls and patients with moderate infection, while FCN1 and SPP1 were highly expressed by patients with severe/critical infection (Fig. 1d). We conducted differentially expressed gene (DEG) analysis (Extended Data Fig. 2e), gene ontology (GO) analysis and gene set enrichment analysis (GSEA) (Extended Data Fig. 2f) between cell groups. Group 1 expressed the peripheral monocyte-like markers S100A8, FCN1 and CD14, and group 2 expressed high levels of the chemokines CCL2, CCL3, CXCL10 and other genes. Both groups had gene expression patterns suggestive of classic M1-like macrophages. Group 3 expressed the immunoregulatory genes A2M, GPR183 and CCL13 and the profibrotic genes TREM2, TGFB1 and SPP1, suggestive of alternative M2-like macrophages, a reparative but profibrotic subset5. Group 4 expressed the alveolar macrophage (AM) genes FABP4, APOC1, MARCO and genes associated with lipid metabolic functions6 and was enriched in BALFs of controls and moderate cases but not in severe cases (Extended Data Fig. 2f,g). Using SCENIC, we found enhanced expression of the transcription factor genes NFκB, STAT1, STAT2 and interferon (IFN) regulatory factors in group 1 and 2, M2-promoting TFEB, NR1H3, PPARA and CREB1 activities in group 3 and AM-specific PPARG and CEBPB in group 4 macrophages (Extended Data Fig. 2h and Methods). These findings suggest that a highly proinflammatory macrophage microenvironment is present in the lungs of patients with severe COVID-19, which is consistent with previous knowledge of macrophage populations during steady-state, inflammation and recovery7.

Memory T cell responses have been observed in patients with SARS8 and are necessary for resolving SARS-CoV infection in mice9,10. We identified six major clusters of T and NK lymphocytes (Extended Data Fig. 3a) based on expression of canonical genes (Extended Data Fig. 3b,c). Lower CD8+ T cell and higher proliferating T cell proportions were observed in patients with severe/critical infection than in patients with moderate infection (Extended Data Fig. 3d,e). These CD8+ T cells expressed higher levels of GZMA, GZMK and FASLG and the tissue-residence markers ITGA1, CXCR6 and JAML (Extended Data Fig. 3f), with upregulation of genes involved in activation, migration and cytokine-related pathways in moderate cases and translation initiation, cell homeostasis and nucleoside metabolic pathways in severe cases (Extended Data Fig. 3g).

We further assessed T cell clonal expansion using single-cell T cell receptor sequencing (scTCR-seq) of seven T cell clusters from BALF from patients with COVID-19 (Fig. 2a and Extended Data Fig. 4a). TCR clonotype was obtainable from 66.9% (26.6–87.3%) of T cells in each patient and from 60.9% (45.2–84.2%) of cells in each T cell cluster (Supplementary Table 2 and Extended Data Fig. 4b). On the uniform manifold approximation and projection (UMAP), projections of clonally expanded T cells showed an aggregative distribution, indicating transcriptional homogeneity in patients with moderate infection versus a diffusive distribution indicating transcriptional heterogeneity in patients with severe/critical infection (Fig. 2b). Cell numbers and clonal expansion differed between CD8+ T cell clusters in patients with moderate versus severe/critical infection (Fig. 2c). More than 50% of ZNF683+ CD8+ T cells from patients with moderate infection belonged to expanded clones, showing a higher amplification index (> 5 cells) than those from patients with severe/critical infection (Fig. 2c). Among the top 20 individual clones (Supplementary Table 3), the majority of cells from patients with moderate infection belonged to the ZNF683+ cluster, but exhibited various expression patterns in patients with severe/critical infection (Extended Data Fig. 4c). ZNF683+-expanded clones that were enriched in patients with moderate infection could represent SARS-CoV-2-specific CD8+ T cells. These expanded clones preferentially expressed tissue-residence genes, including XCL1, CXCR6 and ITGAE (Extended Data Fig. 4d). A higher score of core tissue-residence genes11 was seen in expanded CD8+ T cells versus nonexpanded cells in patients with moderate infection and in total CD8+ T cells in moderately versus severely/critically infected patients (Fig. 2d). Collectively, our data suggest that CD8+ T cells in BALFs from patients with severe/critical infection were less expanded, more proliferative and more phenotypically heterogeneous, whereas a larger proportion of CD8+ T cell effectors with tissue-resident and highly expanded features were present in BALFs from patients with moderate infection.

Fig. 2: T cell clonal expansion in BALF from patients with COVID-19.
figure 2

a, The clustering of BALF-derived T cells in patients with COVID-19 (n = 8, S3 excluded). The dotted lines partition the UMAP into CD4+ T cells (upper), CD8+ T cells (lower), innate T cells (left) and proliferating T cells (right). b, UMAP projection of clonally expanded T cells from patients with COVID-19. c, Clonal expansion status of T cell subsets in patients with moderate and severe/critical infection. Clonotype counts (left) and cell counts (right) are shown. d, Tissue-resident signature score by the expanded versus nonexpanded BALF CD8+ T cells in patients wuth moderate infection and by CD8+ T cells in patients with moderate versus severe/critical infection. In b–d, moderate was n = 3 and severe was n = 5.

Full size image

Finally, we measured several cytokines and chemokines in BALF of patients with COVID-19 (Supplementary Tables 1 and 4). Compared to patients with moderate COVID-19 infection, patients with severe/critical infection had much higher levels of inflammatory cytokines, particularly interleukin (IL)-8, IL-6 and IL-1β, in their BALFs (Extended Data Fig. 5a). Three patients (S1, S8 and S9) with critical infection displayed persistently high levels of cytokines. In scRNA-seq data, IL1B, IL6, TNF and various chemokines (CCL2, CCL3, CCL4 and CCL7) were expressed at higher levels in lung macrophages from patients with severe COVID-19 infection. CXCL9, CXCL10 and CXCL11 levels were much higher in both COVID-19 groups than in healthy people, but CXCL16, whose product binds CXCR6, was more highly expressed in patients with moderate than severe infection (Extended Data Fig. 5b,c). We also examined the corresponding receptors to other upregulated chemokines in patients with COVID-19 (Extended Data Fig. 5d). These data suggest that lung macrophages in patients with severe COVID-19 infection may contribute to local inflammation by recruiting inflammatory monocytic cells and neutrophils though CCR1 and CXCR2, whereas macrophages in patients with moderate COVID-19 infection produce more T cell-attracting chemokines through engaging CXCR3 and CXCR6.

Together with previous studies3,12,13,14, our data suggest that cytokine storm is associated with disease severity in COVID-19. Our study has several shortcomings, including limited and heterogeneous patients and a lack of longitudinal samples collected before and after infection. In addition, the B cell response could not be analyzed due to low cell numbers. The effects of age, pre-existing conditions and immunoregulatory therapies could not be fully assessed. Nevertheless, this landscape of bronchoalveolar immune cells reveals aberrant macrophage and T cell responses underlying immunopathogenesis in COVID-19.

Methods

Patients

Ethics statement

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

Patients and clinical sample collection

All 13 patients with COVID-19 in this study were enrolled from the Shenzhen Third People’s Hospital from January to February, 2020. Disease severity was defined as moderate, severe and critical, according to the ‘Diagnosis and Treatment Protocol of COVID-19 (the 7th Tentative Version)’ by the National Health Commission of China issued on 3 March 2020 (http://www.nhc.gov.cn/yzygj/s7653p/202003/46c9294a7dfe4cef80dc7f5912eb1989.shtml).

Briefly, moderate cases have fever, respiratory symptoms and pneumonia evidenced by computed tomography (CT) imaging. This study enrolled three patients with moderate infection with bilateral pneumonia, labeled M1 to M3. Patients with severe infection were diagnosed on the basis of one of the following criteria: (1) respiratory distress with respiratory rate ≥ 30 times min−1; (2) fingertip oxygen saturation ≤ 93% at resting state; (3) ratio of partial pressure of arterial oxygen to fraction of inspired oxygen (PaO2/FiO2) ≤ 300 mm Hg (1 mm Hg = 0.133 kPa); and (4) obvious progression of lesions in 24–48 h shown by pulmonary imaging > 50%. Two patients with severe infection were enrolled, labeled S1 and S7. These two patients showed PaO2/FiO2 ≤ 300 mm Hg. The patients with critical infection (S2–S6 and S8–S10) were diagnosed on the basis of one of the following criteria: (1) respiratory failure and an artificial airway required for invasive mechanical ventilation; (2) shock; and (3) combined failure of other organs that required intensive care unit monitoring. All of the eight patients with critical infection received invasive mechanical ventilation.

The demographic characteristics of the nine patients studied by scRNA-seq are listed in Supplementary Table 1, which includes three moderate cases (M1–M3), one severe case (S1) and five critical cases (S2–S6). The median age was 57 years, and the participants included six male and three female patients. All nine patients had Wuhan exposure history and had cough and/or fever as the first symptom. Diagnosis of SARS-CoV-2 was based on clinical symptoms, exposure history, chest radiography and were SARS-CoV-2 RNA-positive using commercial quantitative PCR with reverse transcription (qRT–PCR) assays in the sputum, nasal swab and/or BALF. Influenza A and B virus infection were excluded at enrollment. The patient M2 had clinical symptoms and CT evidence showing bilateral pneumonia, although the SARS-CoV-2 RNA test was negative. Patient M2 had close contact with patient S2 and was confirmed to have anti-SARS-CoV-2 antibodies using a SARS-CoV-2 total-IgG antibody detection kit (chemiluminescent immunoassay, nCoV-04100C, innoDx). The cutoff value was a cutoff index of 1 for the kit. Nineteen days after the onset of illness, her plasma IgG cutoff index value was increased to 25.39, showing a strong positive antibody signal and confirming the diagnosis of COVID-19. For assessment of cytokine and chemokine levels, BALF samples from additional patients with COVID-19 were examined. Clinical information of the additional patients is shown in Supplementary Table 4. In addition, three healthy controls, previously characterized and labeled as HC1 to HC3, were included in the study. These donors were confirmed to be free of tuberculosis, tumor and other lung diseases through CT imaging and other laboratory tests. The demographic characteristics of the three healthy controls are shown in Supplementary Table 4.

qRT–PCR assay for SARS-CoV-2 RNA

A throat swab, sputum, nasal swab, anal swab and BALF were collected from the patients at various time points after hospitalization and were sent to the diagnosis laboratory. 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, indicating the number of cycles surpassing the threshold for a positive test. 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.

Isolation of BALF cells

Approximately 20 ml of BALF was obtained and placed on ice. BALF was processed within 2 h and all operations were performed in a BSL-3 laboratory. After passage of BALF through a 100-µm nylon cell strainer to remove clumps and debris, the supernatant was centrifuged and the cells were re-suspended in cooled RPMI 1640 complete medium. BALFs of patients with severe/critical COVID-19 infection contained a larger number of cells (1.25 × 105 to 2.25 × 106 cells ml−1) than those from moderate cases (6.63 × 103 to 5.6 × 104 cells ml−1) and healthy controls (1.17 × 104 to 2.1 × 104 cells ml−1). Then cells were counted in 0.4% Trypan blue, centrifuged and re-suspended at a concentration of 2 × 106 ml−1 for further use.

Cytokine measurement by cytometric bead array

The cytokines including IL-1β, IL-2, IL-4, IL-5, IL-6, IL-8, IL-10, IL-12p70, IL-17, IFN-α, IFN-γ, tumor necrosis factor-α were detected according to the instruction (Uni-medica, cat. no. 503022). In brief, BALF was centrifuged at 1,000g for 10 min and supernatant was taken. Then 25 µl of sonicated beads, 25 µl of BALF supernatant and 25 µl of detection antibodies were mixed and placed on a shaker at 500 r.p.m. for 2 h at room temperature. Then 25 µl of SA-PE was added to each tube directly. The tubes were placed on a shaker at 500 r.p.m. for 30 min. The data were obtained by flow cytometry (Canto II, BD) and were analyzed with LEGENDplex v.8.0 (VigeneTech).

ScRNA-seq and TCR-seq capturing, library construction and sequencing

Single-cell capturing and downstream library constructions were performed using Chromium Single Cell 3g Reagent kits v.2 (10x Genomics; PN-120237, PN-120236 and PN-120262) for HC1 and HC2, kit v.3 (10x Genomics; PN-1000075, PN-1000073 and PN-120262) for HC3 and Chromium Single Cell V(D)J Reagent kits (10x Genomics; PN-1000006, PN-1000014, PN-1000020, PN-1000005) for all patients with COVID-19, according to the manufacturer’s protocol. Briefly, a total 11 µl of single cell suspension (containing 20,000 cells), 40 µl of barcoded gel beads and partitioning oil were loaded to Chromium Chip A/B to generate single-cell gel bead-in-emulsion. The polyadenylated transcripts were reverse-transcribed inside each gel bead-in-emulsion afterwards. Full-length cDNA along with cell barcode identifiers were PCR-amplified and sequencing libraries were prepared and normalized. The constructed library was sequenced on a BGI MGISEQ-2000 or Illumina platform.

Publicly available healthy control data

In addition to data for the three healthy control BALFs, scRNA-seq data from BALF from one additional healthy control were acquired from the Gene Expression Omnibus (GEO) database under accession code GSE128033 (ref. 4), which contains data of one fresh BALF (GSM3660650) from a lung transplant donor generated using a 3′ V2 chemistry kit on Chromium Single Cell Controller (10x Genomics). For the purpose of cross-reference, control data (labeled HC4) were analyzed together with our own dataset.

ScRNA-seq data alignment and sample aggregating

The Cell Ranger Software Suite (v.3.1.0) was used to perform sample de-multiplexing, barcode processing and single-cell 5′ unique molecular identifier (UMI) counting. To detect SARS-CoV-2 reads, a customized reference was built by integrating human GRCh38 and SARS-CoV-2 genome (severe acute respiratory syndrome coronavirus 2 isolate Wuhan-Hu-1, complete genome, GenBank MN908947.3). Specifically, splicing-aware aligner STAR15 was used in FASTQs alignment. Cell barcodes were then determined based on the distribution of UMI counts automatically. The following criteria were then applied to each cell of all nine patients and four healthy controls: gene number between 200 and 6,000, UMI count > 1,000 and mitochondrial gene percentage < 0.1. After filtering, a total of 66,452 cells were left for the following analysis. Finally, a filtered gene-barcode matrix of all samples was integrated with Seurat v.3 (ref. 16) to remove batch effects across different donors. In parameter settings, the first 50 dimensions of canonical correlation analysis (CCA) and principal-component analysis (PCA) were used.

Dimensionality reduction and clustering

The filtered gene-barcode matrix was first normalized using ‘LogNormalize’ methods in Seurat v.3 with default parameters. The top 2,000 variable genes were then identified using the ‘vst’ method in Seurat FindVariableFeatures function. Variables ‘nCount_RNA’ and ‘percent.mito’ were regressed out in the scaling step and PCA was performed using the top 2,000 variable genes. Then UMAP was performed on the top 50 principal components for visualizing the cells. Meanwhile, graph-based clustering was performed on the PCA-reduced data for clustering analysis with Seurat v.3. The resolution was set to 1.2 to obtain a finer result.

Macrophage, NK and T cell re-integration

Macrophages and T cells were re-integrated and re-clustered, respectively, with Seurat v.3. Specifically, macrophages of all samples were integrated using the first 50 dimensions of canonical correlation analysis and PCA. The parameter k.filter was set to 115 to account for the small cell number in sample M3. In the clustering step, parameter resolution was set to 0.8. For all NK and T cells, samples HC2 and S3 were excluded in the integration (Extended Data Fig. 3) due to small T cell numbers in both samples. The parameter k.filter was set to 175 in the integration step. Moreover, the T cells in patients with COVID-19 were re-integrated by excluding sample S3, leaving eight samples in the following TCR analysis (Fig. 2 and Extended Data Fig. 4). The parameter k.filter was set to 140 in this integration step.

Differential gene expression analysis

MAST17 in Seurat v.3 (FindAllMarkers function) was used to perform differential gene expression analysis. For each cluster of macrophage and T cells, DEGs were generated relative to all of the other cells. DEGs of expanded versus nonexpanded CD8+ T cells in moderate cases (n = 3) were calculated using MAST with default parameters. A gene was considered significant with adjusted P < 0.05 (P values were adjusted by false discovery rate in MAST). To compare the differential expression in CD8+ T cells between moderate and severe cases, we first selected upregulated genes in T cells compared to macrophages derived in patients with severe and moderate COVID-19 infection, respectively, using the MAST method with default parameters, to remove higher background noises of macrophage-specific transcripts in T cells. Only upregulated genes in T cells were considered for further DEG analysis between moderate and severe group comparisons. The DEG results are summarized in Supplementary Table 5.

Regulatory network inference

A single-cell regulatory network for four macrophage groups was constructed with SCENIC18. Specifically, GRNBoost2 (https://github.com/tmoerman/arboreto) in pySCENIC v.0.9.19 was applied to infer gene regulatory networks from raw count data. Then potential direct-binding targets (regulons) were selected based on DNA-motif analysis. Finally, gene regulatory network activity for individual cells was identified. To find the regulators for macrophage groups, the regulon activity was averaged. A regulon-group heat map was generated with the R pheatmap package.

Gene functional annotation

For DEGs, GO and GSEA19 were performed with clusterProfiler20, which supports statistical analysis and visualization of functional profiles for genes and gene clusters. In GSEA, 50 hallmark gene sets in MSigDB21 were used for annotation.

TCR V(D)J sequencing and analysis

Full-length TCR V(D)J segments were enriched from amplified cDNA from 5′ libraries via PCR amplification using a Chromium Single-Cell V(D)J Enrichment kit according to the manufacturer’s protocol (10x Genomics). TCR sequences for each single T cell were assembled by Cell Ranger vdj pipeline (v.3.1.0), leading to the identification of CDR3 sequence and the rearranged TCR gene. Cells with both TCR α- and β-chains were kept and cells with only one TCR chain were discarded. For cells with more than one α- or β-chain, the chain with higher UMIs or reads was kept. Then the clonality of the different T cells defined by scRNA-seq data was explored.

Calculations of the composite signature scores by single cells

Tissue-resident core signature scores was calculated for CD8+ T cells by a modified version of Single-Cell Signature Scorer22, which considers a core list of 31 upregulated and downregulated genes by tissue-resident T cells11.

Statistics

The Student’s t-test (t.test in R 3.6.2, two-sided, unadjusted for multiple comparisons) was used to compare tissue-resident core signature scores of the CD8+ T cells belonging to expanded versus nonexpanded clonotypes and patients with moderate versus severe COVID-19 infection. Effect size (Cohen’s d) of the comparisons was calculated with R lsr v.0.5. Differences of median percentage between healthy controls, moderate and severe groups of all cell types were compared using a Student’s t-test (two-sided, unadjusted for multiple comparisons) with R ggpubr v.0.2.5.

Reporting Summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.