key: cord-0724081-kl2w6yiv
authors: Xie, Xiaowei; Cheng, Xuelian; Wang, Gaoxiang; Zhang, Biao; Liu, Mengyao; Chen, Liting; Cheng, Hui; Hao, Sha; Zhou, Jianfeng; Zhu, Ping; Cheng, Tao
title: Single-cell transcriptomes of peripheral blood cells indicate and elucidate severity of COVID-19
date: 2021-02-04
journal: Sci China Life Sci
DOI: 10.1007/s11427-020-1880-y
sha: 131355fffc8f87201a0312e91cd6c84b8fce9d44
doc_id: 724081
cord_uid: kl2w6yiv

The blood and immune system of coronavirus disease 2019 (COVID-19) infected patients are dysfunctional, and numerous studies have been conducted to resolve their characteristics and pathogenic mechanisms. Nevertheless, the variations of immune responses along with disease severity have not been comprehensively documented. Here, we profiled the single-cell transcriptomes of 96,313 peripheral blood mononuclear cells (PBMCs) derived from 12 COVID-19 patients (including four moderate, four severe and four critical cases) and three healthy donors. We showed that proliferative CD8 effector T cells with declined immune functions and cytotoxicity accumulated in the critical stage. By contrast, the quantity of natural killer (NK) cells was significantly reduced, while they exhibited enhanced immune activities. Notably, a gradually attenuated responseto COVID-19 along with disease severity was observed in monocytes, in terms of cellular composition, transcriptional discrepancy and transcription factor regulatory network. Furthermore, we identified immune cell-type dependent cytokine signatures distinguishing the severity of COVID-19 patients. In addition, cell interactions between CD8 effector T/NK cells and monocytes mediated by inflammatory cytokines were enhanced in moderate and severe stages, but weakened in critical cases. Collectively, our work uncovers the cellular and molecular players underlying the disordered and heterogeneous immune responses associated with COVID-19 severity, which could provide valuable insights for the treatment of critical COVID-19 patients. ELECTRONIC SUPPLEMENTARY MATERIAL: Supplementary material is available for this article at 10.1007/s11427-020-1880-y and is accessible for authorized users.

The coronavirus disease 2019 pandemic, resulting from the novel severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), has given rise to substantial patients and panics worldwide (Wiersinga et al., 2020; . Up to November 2020, about 47 million patients were diagnosed as pneumonia, including approximately 1.2 million deaths, with a continuous increase (Dong et al., 2020) . Most infection-confirmed patients exhibited mild and moderate symptoms, 10%-20% of which developed into severe pneumonias and about 5% of which eventually deteriorated into critical cases with acute respiratory distress syndrome (ARDS) requiring mechanical ventilation in an intensive care unit (ICU) (Cao, 2020) . More studies are needed to explore the pathogenesis of COVID-19 for better devising effective drugs against SARS-CoV-2 and the sequelae of COVID-19.

The immune system of COVID-19 patients is seriously dysfunctional, characterized by lymphopenia (Guan et al., 2020) and elevation of inflammatory cytokines (Chen et al., 2019; Zhou et al., 2020) . Recently, single-cell transcriptome analysis has contributed much to interpreting the cellular and molecular characteristics of the COVID-19 pandemic. Analysis of the landscape of bronchoalveolar lavage fluid (BALF) has revealed that highly clonally expanded CD8 T cells occurred in moderate patients (Liao et al., 2020) . Moreover, CD8 effector T and natural killer (NK) (Wilk et al., 2020; Zhang et al., 2020) cells displayed higher cytotoxicity and exhaustion in COVID-19 patients than those in healthy individuals.

Regarding myeloid cells, both reduced expression of HLA genes and potential up-regulation of interferon signatures have been demonstrated in the CD14 + monocytes of COVID-19 patients (Wilk et al., 2020) . Severe COVID-19 patients have been characterized by aberrantly abundant FCN1 + inflammatory macrophages (Liao et al., 2020) . A cytokine storm-related monocyte subpopulation has also been identified in severe cases . Additionally, severe COVID-19 patients could be distinguished from mild cases by dysregulated myelopoiesis, especially the occurrence of immature neutrophils (Schulte-Schrepping et al., 2020; Silvin et al., 2020) . Though these observations revealed immuno-associated phenotypes of COVID-19 patients, a systematic interrogation of cellular and molecular variation along with disease severity and their underlying mechanisms is still lacking. In this study, we performed single-cell RNA sequencing of PBMCs derived from 12 COVID-19 patients (grouped by moderate, severe and critical types) and three healthy donors. Firstly, we revealed the disordered landscape of peripheral blood cells from moderate to critical cases. Then, we disclosed the inverse immune responses of proliferative CD8 effector T and NK cells in critical patients. It was notable that monocytes displayed a gradually attenuated response to COVID-19 concomitant with disease severity in multiple aspects. Finally, frequent interaction between CD8 effector T/NK cells and monocytes mediated by inflammatory cytokines was observed in COVID-19 patients.

To interpret the immunological variations by different severities, single-cell transcriptomes of PBMCs were obtained from 12 COVID-19 patients by 10X Genomics. Reference datasets of three healthy donors were obtained from the 10X official website (see Methods). The 12 patients were categorized into three clinical groups: moderate (n=4), severe (n=4) and critical (n=4, and 3 were in ICU) ( Figure 1A ). Detailed clinical characteristics and routine blood tests of 12 enrolled patients are summarized in Figure S1 and Table S1 in Supporting Information. Additionally, routine blood samples of another 45 COVID-19 patients and 25 healthy donors were collected to observe the fluctuation in the number of blood cells. Compared with healthy controls, the number of lymphocytes significantly declined from moderate to critical stage, in accordance with previously reported lymphopenia (Guan et al., 2020) . Of note, the number of monocytes increased in moderate stage but turned to decline in severe and critical stages, while platelets and white blood cells (WBC) fluctuated dramatically in critical COVID-19 patients ( Figure 1B ).

After the low-quality and doublet filter, single-cell transcriptome profiles of 96,313 PBMCs from 15 individuals were integrated and piped into removal of batch effects, followed by dimension reduction and clustering analysis. To accurately annotate the cell types of PBMCs, we aligned each single cell to our Atlas of Blood Cells (ABC) reference (Xie et al., 2020) using "TransferData" from Seurat (Stuart et al., 2019) . The ABC is a comprehensive transcriptome atlas including 32 human blood cell types based on deep sequencing. The majority of PBMCs were T cells and monocytes, and the overall differentiation trajectory coincided with known hematopoietic hierarchy ( Figure S2A in Supporting Information). In addition, expression distributions of canonical signature genes were visualized by uniform manifold approximation and projection (UMAP) to further verify the cell identities of PBMCs. We observed that CCR7, GZMK and GNLY were specifically expressed in T cells, while CD79A and MS4A1 were highly expressed in the B cell population. Monocytes were characterized by the expression of CD14 and FCGR3A ( Figure S2B in Supporting Information). In total, we captured 10 relatively rough transcriptional clusters for comparing the cellular composition between healthy donors and COVID-19 patients ( Figure 1C ). Particularly, we found a cell cluster displaying both monocyte (FCN1, CCL3L) and platelet (PF4 and PPBP) signatures even after strictly doublet filtering ( Figure S2B in Supporting Information), termed as platelet-monocyte hereafter.

Compared with healthy donors, the cellular composition of PBMCs in COVID-19 patients changed dramatically ( Figure  1D ), though the classifications of 10 cell clusters remained consistent from moderate to critical stage ( Figure S2C in Supporting Information). Concretely, the relative abundance of naive T and B cells decreased while that of effector T/NK cells and proliferative T cells increased in COVID-19 patients. Of note, the cell abundance of platelet-monocyte increased significantly in critical stage ( Figure 1D ), consistent with the increase of platelet-monocyte aggregates during acute infections, thus giving rise to the activation of platelet and adverse clinical outcomes in critical stage (Manne et al., 2020) . Thus, based on the roughly clustered 10 cell types, the overall variations of peripheral blood cells coincided with recent reports (Liao et al., 2020; Zhang et al., 2020) .

To further dissect the variations that occurred in NK/T cells across healthy individuals and COVID-19 patients, 55,012 NK/T cells were extracted from PBMCs and re-clustered per the method described above. Combined with ABC reference and expression distributions of canonical markers ( Figure  S3A and B in Supporting Information), we precisely identified eight subtypes of NK/T cells, namely CD8 naive T cells, CD4 naive T cells, CD4 memory T cells, CD4 regulatory T (Treg) cells, CD8 memory T cells, CD8 effector T cells, proliferative CD8 effector T cells and NK cells ( Figure  2A ). We observed that the relative abundance of proliferative CD8 effector T cells gradually increased with disease severity, whereas accumulated NK cells in moderate and severe stages declined in critical stage ( Figure 2B ; Figure S3C in Supporting Information). Next, to elucidate the dysregulated biological processes in each subtype of NK/T cells, we calculated their differentially expressed genes (DEGs) between any two stages, which were then submitted to unsupervised clustering to identify distinct gene modules. Gene ontology enrichment analysis displayed that gradually downregulated genes in the CD8 effector T cells of COVID-19 patients were associated with "IL-17 signaling pathway" and "response to toxic substance". Significantly down-regulated genes in the proliferative CD8 effector T cells of COVID-19 patients were enriched in "lymphocyte/T cell activation" and "positive regulation of immune effector process". We also found a set of genes related to "natural killer cell mediated cytotoxicity" repressed in NK cells from moderate and severe stages, but reactivated in critical stage ( Figure 2C ). Additionally, we observed that CD8 effector T cells exhibited higher cytotoxicity in COVID-19 patients, consistent with a previous report . Of note, proliferative CD8 effector T cells exhibited lower transcription activities of cytotoxicity-related genes in COVID-19 patients, whereas cytotoxic NK cells were reactivated in critical stage, both of which were in accordance with the change tendency acquired from DEGs above ( Figure S3D in Supporting Information).

To resolve the transcription factor regulatory networks underlying NK/T cells of COVID-19 patients, we estimated the regulon activity score (RAS) for enriched transcription factors in each single cell utilizing Single-Cell rEgulatory Network Inference and Clustering (SCENIC) (Aibar et al., 2017) . The regulatory atlas based on RASs basically revealed the differentiation trajectory of NK/T cells ( Figure  S3E in Supporting Information). We observed that RUNX3, EOMES and TBX21 were activated in CD8 memory T cells, CD8 effector T cells and NK cells, while FOXO1, MXI1 and FOXP1 were specifically up-regulated in naive T cells ( Figure S3F in Supporting Information). In consideration of the increased number and dysfunction of CD8 effector T cells in severe and critical COVID-19 patients, we identified differentially activated transcription factors in all CD8 effector T cells between healthy donors and COVID-19 patients. We found that ZBTB7B, which represses the generation of CD8 effector T cells (Wang et al., 2008) , was inactivated in severe and critical stages. While CBFB (Kas- per et al., 2006) , EP300 (Kasper et al., 2006) and RUNX3 (Shan et al., 2017) , which facilitate the generation of CD8 effector T cells, were activated in severe and/or critical stages ( Figure 2D ). In addition, apoptosis-related gene MAF (Peng et al., 2009 ) was activated in the CD8 effector T cells of critical patients, in line with their higher exhaustion scores ( Figure S3D in Supporting Information). Altogether, we observed that the number of proliferative CD8 effector T cells with attenuated function were increased in critical stage, whereas NK cells displayed the inverse variation pattern (decreased number with heightened function), indicating a heterogeneous immune response of NK/T cells to COVID-19. Furthermore, transcription factors regulating the differentiation of CD8 effector T cells were identified, which might explain their expansion in severe and critical stages.

Regarding B cells, they were divided into three subtypes, namely naive B cells, memory B cells and plasma cells. We observed that naive B cells and memory B cells were significantly decreased, while plasma cells were significantly increased in COVID-19 patients ( Figure S3G in Supporting Information). To investigate the transcriptional changes of B cells along with disease severity, we analyzed the DEGs between healthy donors and COVID-19 patients in different B subtypes. We found that the up-regulated genes in COVID-19 patients for both naive B cells and memory B cells were associated with translation (ribosome), RNA splicing and intrinsic apoptotic signaling pathway ( Figure S3H in Supporting Information), basically coincided with the previous report .

Subsequently, 29,694 monocytes isolated from PBMCs were used to analyze the myeloid variations between healthy donors and COVID-19 patients. We identified a sub-cluster of classical monocytes expressing high levels of interferon (IFN)-stimulated genes such as ISG15 and MXI. In addition, HLA-related genes were highly expressed in several subclusters of monocytes (Schulte-Schrepping et al., 2020) ( Figure 3A ; Figure S4A and B in Supporting Information). Next, we found that the cell proportions of HLA + classical monocytes, intermediate monocytes and non-classical monocytes were increased in moderate and severe stages, but were decreased in critical stage ( Figure 3B ), in accordance with the change tendency based on clinical data ( Figure 1B ). Of note, the cell abundance of ISG15 + classical monocytes enhanced significantly in critical patients ( Figure 3B ). To investigate the transcriptional heterogeneity during distinct severities of COVID-19, we normalized the number of DEGs from each stage in contrast to healthy donors, and found that the transcriptional difference in all sub-clusters of monocytes gradually declined from moderate to critical stage ( Figure   3C ; Figure S4C in Supporting Information). Accordingly, gene expression patterns across different severities based on DEGs were elaborated for each sub-cluster of monocytes.

Remarkably, the majority of monocytes (with the exception of non-classical monocytes and dendritic cells) responded to COVID-19 dramatically in moderate stage, but gradually approached healthy states in severe/critical patients. DEGs in gene modules were significantly enriched in interferon-and cytokine-related biological processes ( Figure 3D ; Figure  S4D in Supporting Information). Next, the transcription factor regulatory network underlying monocytes was profiled ( Figure S4E in Supporting Information). CEBPD and RUNX1 were activated in CD14 + monocytes while NFE2L1 and TCF12 functioned in CD16 + monocytes (non-classical monocytes) ( Figure S4F in Supporting Information). By comparing the RASs of transcription factors between healthy individuals and COVID-19 patients, we observed that the activities of NF-κB (NFKB, NFIL3, REL, RELA, RELB) (Ammon et al., 2000; Baek et al., 2009; de Wit et al., 1998) and IFN (IRF1, IRF8, IRF9, PRDM1) (Chistiakov et al., 2018; Doody et al., 2006) gradually declined along with disease progression, for both CD14 + monocytes and CD16 + monocytes ( Figure 3E ; Figure   S4G in Supporting Information). These results suggested that monocytes exhibited a gradually attenuated response to COVID-19 from moderate to critical stage, in the combination of cell number, transcriptional heterogeneity and transcription factor regulatory network.

Considering that many DEGs in NK/T cells and monocytes were enriched in cytokine-mediated signaling pathway, in addition to previously reported cytokine storm, we systematically assessed the transcription levels of inflammatory cytokines across COVID-19 patients with different severities. We found that most inflammatory cytokines such as well-known IL10 (Huang et al., 2020) , CCL3 (Huang et al., 2020) and IL6 displayed higher expression levels in COVID-19 patients, especially in critical stage, than in healthy individuals ( Figure 4A ). The concentrations of cytokines in blood plasma were also measured to further confirm the increased levels of IL10, IL6, CCL3 and CXCL10 in COVID-19 cases ( Figure 4B ). Moreover, there existed favorable cytokines with high expression in B cells (CCL5, IL18, XCL1, XCL2), T cells (CCL4, IL6, IL10, CXCR3, CXCR6) and monocytes (CCL2, IL6, IL10, CXCL8, CXCL10), respectively from moderate, severe and critical stages ( Figure 4C ), suggesting that disease severity of COVID-19 might be distinguished by immune cell-type dependent cytokine signatures.

The immune responses of NK/T cells and monocytes to COVID-19 were distinct, both of which expressed abnormal levels of cytokines. Therefore, we speculated that monocytes might be involved in the regulation of other immune cell types by cytokines. Cell interactions between different cell clusters of PBMCs were evaluated based on the expression levels of cytokine/receptor pairs by way of CellPhoneDB (Efremova et al., 2020) , and the most frequent cell interactions were determined between CD8 effector T/NK cells and monocytes ( Figure 4D ). Next, we found that the interaction intensities between CD8 effector T/NK cells and monocytes (supported by the expression levels and count of cytokine/ receptor pairs) were elevated in moderate and severe cases but were diminished in critical cases ( Figure 4E and F). To trace the pivotal cytokine/receptor pairs, the connection intensities between healthy individuals and COVID-19 cases were compared, and we found that cell interactions involving inflammatory cytokines were significantly increased in COVID-19 patients. Strikingly, CCL3L1 (a chemokine that binds to the HIV receptor CCR5 to inhibit HIV entry (Dolan et al., 2007) ) and DPP4 (a known receptor for Middle East Respiratory Syndrome Coronavirus (MERS-CoV) (Letko et al., 2018)) communicated with each other in critical cases. However, the interaction intensities between CD8 effector T cells and monocytes mediated by CCL3, CCL4, CCL4L2 and IL1B, which were also mentioned recently , strengthened in moderate and severe cases but weakened in critical stage. In addition, two interaction pairs involving CD74 that were reported previously exhibited higher interactions only between CD8 effector T cells and CD14 + monocytes ( Figure 4G ). In conclusion, it was plausible that monocytes interacted with CD8 effector T/NK cells, possibly mediated by inflammatory cytokines, particularly frequently in moderate and severe cases.

In this study, we depicted the disordered single-cell landscape of 96,313 PBMCs derived from 12 COVID-19 patients (from moderate to critical stage) and three healthy donors. We revealed the heterogeneous immune response of CD8 effector T/NK cells to COVID-19. Of note, monocytes displayed a gradually attenuated response to COVID-19 associated with disease severity. Moreover, the frequent communications between CD8 effector T/NK cells and monocytes mediated by inflammatory cytokines in moderate and severe stages became weakened in critical cases. Previous attention has primarily been paid to lymphopenia as the main feature of COVID-19 patients, coincided with the reduced number of naive T and memory T cells in our study. However, the cell proportions of CD8 effector T cells dramatically increased, especially that of proliferative CD8 effector T cells increased along with disease severity. With regard to the reason, we speculated that the inactivation of ZBTB7B (Wang et al., 2008) , repressed by CBFB (Kasper et al., 2006) , might explain the increased number of CD8 effector T cells in severe and critical stages. Activation of EP300 (Kasper et al., 2006) and RUNX3 (Shan et al., 2017) also possibly contributed to the differentiation of CD8 effector T cells in severe COVID-19 patients. The proliferative capacity of CD8 effector T cells was probably associated with their dysfunctional state (Li et al., 2019) , which was consistent with their down-regulated biological function in COVID-19 patients. In addition, the apoptosis-related gene MAF (Peng et al., 2009 ) was activated in the CD8 effector T cells of critical patients. This, combined with their higher exhaustion score, suggested that accumulated CD8 effector T cells in critical stage were more susceptible to apoptosis and exhaustion. In contrast to CD8 effector T cells, NK cells displayed exactly inverse variation pattern in critical stage: the accumulation of NK cells in moderate and severe stages declined, whereas attenuated NK cells were reactivated in critical patients.

Severe COVID-19 patients have been characterized by abnormal myeloid compartments, such as the appearance of cytokine storm-related monocytes and immature neutrophils Schulte-Schrepping et al., 2020; Silvin et al., 2020) . We observed that inflammatory cytokines such as well-known IL10, CCL3 and IL6 definitely displayed higher expression levels in COVID-19 patients, further confirmed by their concentrations in blood plasma. Unexpectedly, monocytes displayed a gradually attenuated response to COVID-19 during disease severity, strongly supported by cell number, biological function and transcription factor regulatory network. Moreover, the frequent interactions between CD8 effector T/NK cells and monocytes in moderate and severe cases also turned to decline in critical stage, possibly affected by the attenuated monocytes. These results indicated that, the heterogeneous immune response of NK/T cells and monocytes to COVID-19 and the cell interaction mediated by inflammatory cytokines between them might play synergistic roles in fighting against COVID-19.

In summary, our current work uncovered a disordered immune profile along with disease severity of COVID-19, the characteristics of which were especially manifested in critical patients. This work could provide valuable insights for the treatment of COVID-19 patients.

A total of 12 COVID-19 patients diagnosed as moderate (n=4), severe (n=4) or critical (n=4) stage from Tongji Hospital (Wuhan, China) were involved in this study. The categorization of clinical groups (moderate, severe and critical) was defined based on Guidelines for Diagnosis and Treatment of Coronavirus Disease 2019 issued by the National Health Commission of China (5th edition). Moderate, fever, respiratory tract symptoms and pneumonia manifestations on imaging; severe, any who meet any of the following criteria: respiratory frequency ≥30 min -1 , oxygen saturation ≤93% at a resting state, arterial partial pressure of oxygen/oxygen concentration ≤300 mmHg, lung in-filtrates>50% within 24 to 48 h on imaging; critical, meeting any of the following criteria: respiratory failure requiring mechanical ventilation, septic shock, organ failure that requires monitoring and treatment in the ICU. The blood samples were collected from patients after hospitalization. Sample collection, pre-processing and laboratory operations were approved by the Tongji Hospital Research Ethics Committee. Written informed consents were obtained from each enrolled patient during the COVID-19 pandemic.

Clinical information including medical history, symptoms and routine blood samples of 12 COVID-19 patients were collected from medical records complying with WHO guidance. Additionally, routine blood samples of another 45 COVID-19 patients were obtained from Tongji Hospital to compare the cell number of each cell type between moderate (n=15), severe (n=15) and critical (n=15) patients. The routine blood samples of 25 healthy individuals were included as control groups, which were recruited from the Institute of Hematology and Blood Disease Hospital, Chinese Academy of Medical Sciences (Tianjin, China).

The blood samples were collected on the admission of COVID-19 patients within 24 h. Peripheral blood mononuclear cells (PBMCs) were isolated using a Ficoll-Paque (GE Healthcare) solution according to standard density gradient centrifugation methods. Cells were harvested and resuspended in freezing media (90% FBS, 10% DMSO) and frozen using a Freezing Container (Nalgene® Mr. Frosty® Cryo 1°C, Thermo Fisher Scientific) in a -80°C freezer for less than 20 days until use. The number and viability of cells were measured using a TC20 automated cell counter (Biorad). Dead cells (cell viability less than 80%) were removed by magnetic bead purification (Miltenyi Biotech) according to the manufacturer's protocol. The single-cell transcriptome libraries of 12 patients were constructed following the instructions of single-cell 3′ solution v3 reagent kit (10X Genomics). Briefly, cells were resuspended with PBS containing 0.04% BSA and loaded onto a chromium single cell chip B along with partitioning oil, reverse transcription reagents and gel beads. After generating single-cell gel bead-in-emulsions (GEMs), reverse transcription was performed using a Verity Thermal Cycler (Life Technologies). Barcodes were added to released RNAs from lysed cells through reverse transcription in individual GEM, followed by fragmentation, end repair, polyA tailing and adaptor ligation based on the standard protocol. The cDNA purification and size selection were achieved by SPRI select beads (Beckman Coulter), whose quality was then assessed using an Agilent 2100 Bioanalyzer. Finally, the libraries were sequenced on a MGISEQ-2000 sequencer as paired-end 150 bp reads by Beijing Genomics Institute (BGI, Shenzhen, China).

Gene-barcode matrices for 12 COVID-19 patients were generated by Cell Ranger (v.3.1.0) referring to the human genome (GRCh38). Filtered gene-barcode matrices of three healthy donors (connect_5k_pbmc_NGSC3_ch1, 5k_pbmc_v3 and pbmc_10k_v3) were downloaded from the 10X Genomics official website. The Seurat (Stuart et al., 2019) package (version 3.0.2) implemented in R (version 3.6.0) was applied to achieve quality-filter, batch effects removal, dimension reduction and clustering analysis based on gene expression matrices. First, only single cells with detected genes more than 300 and total counts in the range from 2,000 to 20,000 were submitted to subsequent analysis. Doublets identified by DoubletFinder (version 2.0.3) (McGinnis et al., 2019) with default parameters were discarded. Then, "IntegrateData" was used to integrate PBMCs from different donors with the default 2000 features. Based on the union of highly variable genes calculated by "Find-VariableFeatures", "RunPCA" and "RunUMAP" (Becht et al., 2019; McInnes et al., 2018) were performed to project single cells into a two-dimensional space. Finally, cell clusters were identified by "FindNeighbors" and "FindClusters" at the resolution of 0.42 and visualized by UMAP. To elaborately dissect the heterogeneity, NK/T cells and monocytes were extracted from PBMCs and re-clustered following the procedures as described above. The annotations of cell types referring to ABC (Atlas of blood cells) were achieved by "TransferData" from Seurat.

"FindAllMarkers" was applied to identify cell cluster-specific marker genes, while "FindMarkers" was used to detect the DEGs between any two given stages. The filtered criteria for marker genes and DEGs were fold-change ≥1.5 or ≤0.67 and adjusted P-value ≤0.05. Gene ontology enrichment analysis on DEGs was performed by Metascape (Zhou et al., 2019) .

To compare the cytotoxicity and exhaustion state of CD8 effector T/NK cells between healthy donor and COVID-19 patients, we calculated the UMIs fraction of cytotoxicityand exhaustion-associated gene sets, as the total UMIs of those genes divided by the total UMIs of each single cell (Li et al., 2019) . Twelve well-defined cytotoxicity markers (PRF1, IFNG, GNLY, NKG7, GZMB, GZMA, GZMH, KLRK1, KLRB1, KLRD1, CTSW and CST7) and six welldefined exhaustion-related genes (LAG3, TIGIT, PDCD1, CTLA4, HAVCR2 and TOX) were used to evaluate the cytotoxicity and exhaustion state of NK/T cells derived from (Guo et al., 2018; Zhang et al., 2020) .

The identification of gene patterns among COVID-19 patients with different severities was performed by way of SOTA (Herrero et al., 2001) function from the clValid package. The number of gene modules was adjusted in terms of their variation tendencies. Then, genes belonging to each pattern were used to carry out the gene ontology enrichment analysis by Metascape (Zhou et al., 2019) .

First, regulon activity scores (RASs) of enriched transcription factors in each single cell were acquired utilizing SCENIC (Aibar et al., 2017) software based on UMI matrix. The RAS matrices were then submitted to Seurat to investigate the regulatory networks underlying NK/T cells and monocytes. Meanwhile, the transcription factors with high RAS were visualized by UMAP. Finally, differentially activated transcription factors in each cell cluster between healthy donors and COVID-19 patients with different severities were identified by "FindMarkers".

To resolve the cellular communications between different cell clusters of PBMCs, CellphoneDB (Efremova et al., 2020) was implemented to detect the cell interactions on the basis of normalized UMI (UMI*10000/sum) from healthy, moderate, severe and critical stage, respectively. Significant cell interactions (P value<0.05 and mean>1.0) were retained to compare the communications between monocytes and CD8 effector T/NK cells in healthy donor and COVID-19 patients with different severities.

Both Student's t-test and nonparametric Wilcoxon test were conducted by R language to compare the difference between two groups. Reported P values were from two-sided tests and P values <0.05 were considered to be significant. R language was applied to plot heat maps (by pheatmap package), box plots, bar plots, line plots, point plots, etc. (by ggplot2 package).

Blood plasma was collected as early as possible during hospitalization and preserved in a -80°C freezer until use. Eight cases for each clinical group (moderate, severe and critical) and five healthy donors were involved in this measurement. The concentrations of 48 cytokines and chemokines were measured using Bio-Plex ProHuman Cytokine Screening Panel (48-Plex#12007283, Bio-Rad) on a Luminex FLEXMAP 3D system in accordance with the manufacturer's instructions.

The author(s) declare that they have no conflict of interest. For studies involving human, the authors also state that they conformed with the Helsinki Declaration of 1975 Declaration of (as revised in 2008 concerning Human and Animal Rights, and that they followed out the policy concerning Informed Consent as shown on Springer.com.

SCENIC: single-cell regulatory network inference and clustering

Differential expression of the transcription factor NF-κB during human mononuclear phagocyte differentiation to macrophages and dendritic cells

Identification of novel transcriptional regulators involved in macrophage differentiation and activation in U937 cells

Dimensionality reduction for visualizing single-cell data using UMAP

COVID-19: immunopathology and its implications for therapy

Clinical and immunologic features in severe and moderate forms of coronavirus disease

The impact of interferon-regulatory factors to macrophage differentiation and polarization into M1 and M2

Regulation of p100 (NFKB2) expression in human monocytes in response to inflammatory mediators and lymphokines

CCL3L1 and CCR5 influence cell-mediated immunity and affect HIV-AIDS pathogenesis via viral entry-independent mechanisms

An interactive web-based dashboard to track COVID-19 in real time

BLIMP-1 is a target of cellular stress and downstream of the unfolded protein response

CellPhoneDB: inferring cell-cell communication from combined expression of multi-subunit ligand-receptor complexes

Clinical characteristics of coronavirus disease 2019 in China

Single-cell analysis of two severe COVID-19 patients reveals a monocyte-associated and tocilizumab-responding cytokine storm

Global characterization of T cells in non-small-cell lung cancer by single-cell sequencing

A hierarchical unsupervised growing neural network for clustering gene expression patterns

Clinical features of patients infected with 2019 novel coronavirus in Wuhan

Conditional knockout mice reveal distinct functions for the global transcriptional coactivators CBP and p300 in T-cell development

Adaptive evolution of MERS-CoV to species variation in DPP4

Dysfunctional CD8 T cells form a proliferative, dynamically regulated compartment within human melanoma

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

Platelet gene expression and function in patients with COVID-19

DoubletFinder: doublet detection in single-cell RNA sequencing data using artificial nearest neighbors

Umap: Uniform manifold approximation and projection for dimension reduction

c-Maf increases apoptosis in peripheral CD8 cells by transactivating Caspase 6

Single cell RNA sequencing of 13 human tissues identify cell types and receptors of human coronaviruses

Severe COVID-19 is marked by a dysregulated myeloid cell compartment

The transcription factor Runx3 guards cytotoxic CD8 + effector T cells against deviation towards follicular helper T cell lineage

Elevated calprotectin and abnormal myeloid cell subsets discriminate severe from mild COVID-19

Comprehensive integration of single-cell data

The zinc finger transcription factor Zbtb7b represses CD8-lineage gene expression in peripheral CD4

Pathophysiology, transmission, diagnosis, and treatment of coronavirus disease 2019 (COVID-19)

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

Single-cell transcriptomic landscape of human blood cells

Evolution of the novel coronavirus from the ongoing Wuhan outbreak and modeling of its spike protein for risk of human transmission

Single-cell landscape of immunological responses in patients with COVID-19

Pathogenic T-cells and inflammatory monocytes incite inflammatory storms in severe COVID-19 patients

Metascape provides a biologist-oriented resource for the analysis of systems-level datasets

 (81421002, 81730006, 81890990, 81870086, 81670106, 81922002), and CAMS Initiative for Innovative Medicine (2017Medicine ( -I2M-1-015, 2019Medicine ( -I2M-1-006, 2017.

The supporting information is available online at https://doi.org/10.1007/s11427-020-1880-y. The supporting materials are published as submitted, without typesetting or editing. The responsibility for scientific accuracy and content remains entirely with the authors.