key: cord-268112-zk8z8z4i
authors: Zhu, Linnan; Yang, Penghui; Zhao, Yingze; Zhuang, Zhenkun; Wang, Zhifeng; Song, Rui; Zhang, Jie; Liu, Chuanyu; Gao, Qianqian; Xu, Qumiao; Wei, Xiaoyu; Sun, Hai-Xi; Ye, Beiwei; Wu, Yanan; Zhang, Ning; Lei, Guanglin; Yu, Linxiang; Yan, Jin; Diao, Guanghao; Meng, Fanping; Bai, Changqing; Mao, Panyong; Yu, Yeya; Wang, Mingyue; Yuan, Yue; Deng, Qiuting; Li, Ziyi; Huang, Yunting; Hu, Guohai; Liu, Yang; Wang, Xiaoqian; Xu, Ziqian; Liu, Peipei; Bi, Yuhai; Shi, Yi; Zhang, Shaogeng; Chen, Zhihai; Wang, Jian; Xu, Xun; Wu, Guizhen; Wang, Fusheng; Gao, George F.; Liu, Longqi; Liu, William J.
title: Single-cell sequencing of peripheral blood mononuclear cells reveals distinct immune response landscapes of COVID-19 and influenza patients
date: 2020-07-19
journal: Immunity
DOI: 10.1016/j.immuni.2020.07.009
sha: 
doc_id: 268112
cord_uid: zk8z8z4i

SUMMARY COVID-19 is a severe infectious disease that is a current global health threat. However, little is known about its hallmarks compared to other infectious diseases. Here, we report the single-cell transcriptional landscape of longitudinally collected peripheral blood mononuclear cells (PBMCs) in both COVID-19 and influenza A virus (IAV)-infected patients. We observed increase of plasma cells in both COVID-19 and IAV patients, and XAF1-, TNF- and FAS-induced T cell apoptosis in COVID-19 patients. Further analyses revealed distinct signaling pathways activated in COVID-19 (STAT1 and IRF3) vs. IAV (STAT3 and NFκB) patients and substantial differences in the expression of key factors. These factors include relatively increase of IL6R and IL6ST expression in COVID-19 patients, but similarly increased IL-6 concentrations compared to IAV patients, supporting the clinical observations of increased pro-inflammatory cytokines in COVID-19 patients. Thus, we provide the landscape of PBMCs and unveil distinct immune response pathways in COVID-19 and IAV patients.

(MS4A1 + IGHG1 + ); two were annotated as NK cells (KLRF1 + ) and one population was 167 labeled DCs (CD1C + LYZ + ) and monocytes (LYZ + CD68 + ) (Figures 1D and S1A) . Most 168 of the clusters consisted of cells from multiple patients indicating common immune 169 traits among patients. In addition, PBMC samples from patients did not express ACE2 170 and TMPRSS2 receptors and did not exhibit viral reads, indicating that SARS-CoV-2 171 may not infect PBMCs ( Figure S1B) . 172 173

The general patterns of PBMC cell populations were comparable across patients 175 (Figure 2A and Table S2 promoting Ig synthesis. XBP1 is a positively-acting transcription factor in the 189 CREB-ATF family that is expressed at a high amounts in plasma cells, and is crucial 190 for increasing protein synthesis in plasma cells (Shaffer et al., 2004) . IRF4 regulates 191 immunoglobulin class switch recombination and sustained and higher concentrations 192 of IRF4 are known to promote the generation of plasma cells (Ochiai et al., 2013) . 193

found that the expression of CD2AP on activated CD4 + T cells was elevated 195 compared with healthy controls (Figure 2E) . The adaptor molecule CD2AP in CD4 + T 196 cells modulates the differentiation of follicular helper T cells, and improves protective 197 antibody responses in viral infection (Raju et al., 2018) . In addition to supporting 198 plasma cell function, TNFSF14 expression was also increased in activated CD4 + T 199 cells and cytotoxic CD8 + T cells. These factors promote T cell activation as well as T 200 cell recruitment to tissues from peripheral blood. KDM5A, which encodes an 201

H3K4me3 demethylase that is required for NK cell and T cell activation, is 

Gene Ontology (GO) analyses were performed to gain insight into functions of 209 different cell subsets between COVID-19 patients and healthy controls. Genes in the 210 group "Response to type I interferon signaling" were enriched in T, B and NK cell 211 subsets of D1 and D4, but not D16 samples (Figures 3A and S3A 3A and S3A) , which may be due to a higher proportion of plasma cells in B cell 219 clusters because high demand of protein synthesis was required during antibody 220 production. Other signaling pathways, such as "Regulation of chromosome 221 organization", "DNA conformation change", etc. in SARS-CoV-2 infection were also 222 upregulated. The roles of these genes will need further investigation. 223 224 Differential expressed genes (DEGs) in these transcriptomic profiles were then 225 compared between COVID-19 patients and healthy controls. ISGs, which are vital to 226 early viral control (Schoggins and Rice, 2011), were identified in patients infected with 227 SARS-CoV-2 on D1 and D4 (Figures 3B and S3B ), in line with the enrichment for 228 "Response to type I interferon signaling" pathways in our GO analysis (Figures 3A 229 and S3A). Among these ISGs, ISG15, IFI44L, MX1, and X-linked inhibitor of 230 apoptosis (XIAP)-associated factor 1 (XAF1) were up-regulated in T, B, NK, and DC 231 cell subsets (Figures 3B and S3B) . The expression of these four genes was also 232 much higher in COVID-19 patients compared to healthy controls at the bulk level 233 ( Figure 3C) . We then examined transcription dynamics of these genes during the 234 disease process. To achieve this, we divided the disease processes of COVID-19 235 patients from symptom onset to discharge into 4 stages ( Table S3) . We identified 6 236 time-dependent expression patterns ( Figure 3D ) and investigated their biological 237 significance ( Figure 3E ). Cluster 3 contain 158 genes that had decreased expression 238 levels over time. The functions of these genes were significantly enriched in biological 239 processes associated with interferon responses, indicating that the transcriptional 240 regulation of these genes is dynamic, and that they were activated at early time points 241 and deactivated at late time points (Figure 3E ). Cluster 1 contains 38 genes whose 242 expression levels were elevated from stage 2. GO enrichment analysis showed that 243 their functions were significantly enriched in translation and protein synthesis-related 244 pathways. This is consistent with the timing of antibody production (Thevarajan et al., drives apoptosis under stress (Jeong et al., 2018) . TP53-mediated apoptosis was also 258 enhanced by XAF1 via post-translational modification (Zou et al., 2012) .Therefore, 259 expression of genes linked to XAF1-mediated apoptosis, including IRF1, TP53, 260 BCL2L11 and CASP3, was analyzed ( Figure S3D) . Expression of IRF1, TP53, and 261 CASP3 was increased in T, B and NK cell subsets in COVID-19 patients compared to 262 controls, while BCL2L11 exhibited different patterns in different cell subsets. In 263 addition to the XAF1 related apoptosis pathway, expression of genes in other 264 apoptosis-linked pathways, including TNF-and Fas pathways (Elmore, 2007) , was 265 examined in both COVID-19 patients and healthy controls (Figures 3F and S3E) . The 266 expression of TNFSF10 (TRAIL) and its receptor TNFRSF10A were increased in T 267 cells from COVID-19 patients relative to healthy controls. Other TNF path member, 268

including TNFRSF1B, were also relatively up-regulated in COVID-19 patients. As for 269 the FAS path, the expression of FAS, FASLG, and FADD were up-regulated in T cells 270 of COVID-19 patients, though not significantly ( Figure 3F ). In B and NK cell subsets, 271 TNFSF10 and FADD were notably increased in COVID-19 patients with other genes 272 increased mildly, except that Fas in B cells and TRADD in NK cells were subtly 273 decreased ( Figure S3E) . Taken together, we find that up-regulated genes relevant to 274 the XAF1-, TNF-, and Fas pathways may lead to increased T cell apoptosis in 275 COVID-19 patients. (Ludwig and Planz, 2008) . Activation of NF-κB plays a key role in regulating the 292 proinflammatory innate immune response and adaptive immune response. Previous 293 work showed that IL-6 and IL-10 are highly increased in severe IAV patients (Yu et al., 294 2011), known activators of STAT3 signaling. Expression of STAT3 was elevated in IAV 295 patients compared to COVID-19 patients and healthy controls, and seemed to 296 correlate with time post-admission. In addition, RUNX3 expression was up-regulated 297 in activated CD4 + T cells of IAV patients compared to the other groups ( Figure 4B) . It 298 has been suggested that RUNX3 induction is a key step for CD4 + T cells to acquire 299 cytotoxic activity, whereas another study showed that RUNX3 induced by IAV 300 infection through the NF-κB pathway promoted apoptosis in airway epithelial cells 301 (Gan et al., 2015) . The role of RUNX3 in the T cell-mediated immune response 302 remains to be determined. Several pro-viral host factors such as CHD1, BCLAF1, and 303 PHF3 that contribute to viral infection, replication and immune evasion, were also 304 up-regulated in IAV patients (Figures 4A and S4A ) (Shapira et al., 2009) . To see how 305 response of these genes changes over time, we divided the disease processes of IAV 306 patients into 2 stages (Table S3) . We did not find any gene exhibiting statistically 307 significant time-dependent regulation, although some them had trends toward 308 increased or decreased expression (Figures S4B and S4C) . For COVID-19 patients, STAT1, a major transcription factor activated in response to 311 interferon, was up-regulated in activated CD4 + T cells, cytotoxic CD8 + T cells, naïve T 312 cells and DCs (Figures 4A and S4A) . Several proinflammatory factors, including TNF 313 and TNFSF14, were elevated in activated CD4 + T cells, cytotoxic CD8 + T cells, MAIT 314 cells, and NK cells (Figures S4A and S4E) , suggesting enhanced effector function 315 and memory cell development (Desai et al., 2018) . Expression of IL6 was not 316 detectable in PBMCs of all patients or healthy controls. We measured plasma 317 concentrations of IL-6 during hospitalization of these COVID-19 patients and after 318 leaving the hospital ( Figure S4D ). Plasma IL-6 were above the normal range (0.0-7.0 319 pg/mL) (Chen et al., 2020) in 2 out of 5 patients when admitted to the hospital, 320 dropped back to normal during hospitalization, and remained stable after recovery. 321

There was a rising phase of IL-6 in most COVID-19 patients indicating an active 322 inflammatory response, which is also observed in IAV patients (Yu et al., 2011) . The 323 patient with severe symptoms had a much higher level of IL-6 when hospitalized, and 324 it took longer to return to normal plasma IL-6 levels, which corresponded to disease 325 severity. Finally, we examined temporal changes in their expression (Figures S4B and S4C) . 353 Although some interferon-responsive transcription factor genes, such as STAT1 and 354 IRF3, tend to have decreased expression levels over time, none of these patterns 355 were statistically significant due to the limited number of patients. Future studies with 356 more patients and more time points will address this question. (Xiong et 399 al., 2020). In addition to XAF1-induced apoptosis, the extrinsic pathway of apoptosis, 400 including TNF-α/TNFR1 and Fas/FasL path (Elmore, 2007) were also found to be 401 See also Figure S1 and Table S1 . See also Figure S2 and Table S2 . See also Figure S3 and Table S3 . All clinical information including demographic data, medical history, symptoms, signs, 595 and laboratory data were collected from patient medical records. Laboratory data, 596 included blood routine, lymphocyte subsets, infection-related biomarkers, 597 inflammatory cytokines. The total number of leukocytes, percentage of neutrophils 598 and percentage of lymphocytes in peripheral blood were counted by hemocytometer. 599

The number and percentage of lymphocyte subsets was analyzed using the 600 FACSCanto flow cytometer for COVID-19 patients on admission. C-reactive protein 601 and Lactic acid was detected by the Beckman automatic biochemical analyzer. 602

Interleukin 6 (IL-6) was detected using the ROCHE Elecsys IL-6 assay. 

Cell clustering was performed by Seurat (v3.1) R toolkit (Butler et al., 2018) 659 (https://github.com/satijalab/seurat). Genes expressed in less than 3 cells were 660 filtered out and cells with less than 300 or more than 6000 genes detected were 661 excluded. In order to deal with the batch effect, the "NormalizeData" and 662 "FindVariableGene" function were performed respectively for each sample (n=23). 663

After that, these 23 batches were integrated together using "FindIntergrationAnchors" 664 and "IntegrateData" function with dims parameter set to 20. Then, the integrated 665 dataset was scaled and PCA was calculated. The first 20 PCs were used to construct 666 a SNN network and a graph-based clustering approach, louvain algorithm, was 667 applied to identify cell clusters with the resolution set to 1. Finally, UMAP was applied 668 to visualize the clustering result in 2D space. To further regress out the influence of 669 blood cells, we excluded this cluster of cells as well as 10 hemoglobin genes in the 670 expression matrix (HBA1,HBA2, HBB, HBG1, HBG2, HBQ1, HBD, HBM, HBE1 and 671 HBZ), and performed clustering again using the same method described above. 672 673

To annotate each cluster to a specific immune cell type, we selected some classic 675 markers for immune cells and used violin plots ( Figure 1D ) and UMAP feature plots 676 ( Figure S1 ) to annotate cell types. The following genes were used for cell type 

We performed DEG analysis using the "FindMarkers" function in Seurat R package. 684

The cell populations we wanted to compare were input as ident.1 and ident. (Table S3) . First, the average expression 695 of each gene was calculated for each stage. Next, we add 0.000001 for every gene to 696 avoid 0 in the expression matrix which was not acceptable for Mfuzz. Then, the 697 "filter.std(min.std=0)", "standardize()" and "mestimate()" functions were performed for 698 preprocessing according to the tutorial. After that, we clustered the genes into 6 699 different expression programs (Figure 3D) . For the GO analysis of the 6 programs, 700 we first excluded genes with the maximum expression less than 1. Only GO terms 701 with adjusted P value <= 0.05 were shown in the figure (Figure 3E) . ribonucleoprotein complex assembly proteasome-mediated ubiquitin-dependent protein catabolic process regulation of symbiosis, encompassing mutualism through parasitism proteasomal protein catabolic process ERAD pathway response to interferon-gamma response to interferon-alpha interferon-gamma-mediated signaling pathway regulation of multi-organism process response to type I interferon cellular response to type I interferon type I interferon signaling pathway regulation of mononuclear cell proliferation regulation of lymphocyte proliferation leukocyte migration neutrophil mediated immunity neutrophil activation involved in immune response neutrophil degranulation lymphocyte differentiation interleukin-6 production regulation of interleukin-6 production neutrophil activation protein targeting establishment of protein localization to membrane nuclear-transcribed mRNA catabolic process protein targeting to membrane protein localization to endoplasmic reticulum establishment of protein localization to endoplasmic reticulum nuclear-transcribed mRNA catabolic process, nonsense-mediated decay protein targeting to ER cotranslational protein targeting to membrane translational initiation SRP-dependent cotranslational protein targeting to membrane regulation of chromosome organization covalent chromatin modification histone modification histone methylation mRNA splicing, via spliceosome RNA splicing, via transesterification reactions with bulged adenosine as nucleophile RNA splicing, via transesterification reactions regulation of mRNA metabolic process regulation of mRNA processing RNA splicing antigen receptor-mediated signaling pathway endoplasmic reticulum to cytosol transport retrograde protein transport, ER to cytosol proteasomal protein catabolic process IRE1−mediated unfolded protein response cellular response to topologically incorrect protein cellular response to unfolded protein endoplasmic reticulum unfolded protein response protein folding response to endoplasmic reticulum stress response to topologically incorrect protein response to unfolded protein regulation of intracellular protein transport regulation of viral process regulation of viral genome replication regulation of histone modification RNA localization viral life cycle cellular response to type I interferon type I interferon signaling pathway positive regulation of chromosome organization chromatin remodeling viral genome replication regulation of chromosome segregation negative regulation of chromosome organization response to type I interferon regulation of mRNA processing regulation of chromatin organization regulation of RNA splicing double−strand break repair regulation of mRNA metabolic process sister chromatid segregation histone modification covalent chromatin modification defense response to virus regulation of DNA metabolic process response to virus RNA splicing, via transesterification reactions mRNA splicing, via spliceosome RNA splicing, via transesterification reactions with bulged DNA conformation change regulation of chromosome organization RNA splicing BTF3  SAP18  ZFP36  ZFP36L2  IL10RB  IL6R  IL6ST  IL27RA  IL17RA  IL2RA  IL10RA  IL7R  IFNAR1  TNFSF14  STAT1  IRF3  TNFSF10  TNF  IL23A  IFNG  STAT3  CHD1  REL  RUNX3  BCLAF1  JUND  NFKB1  PHF3  CD69  IL4R Expression level Expression level

Integrating single-cell 723 transcriptomic data across different conditions, technologies, and species

Epidemiological and clinical characteristics of 99 cases of 2019 novel coronavirus pneumonia in 727

China: a descriptive study

The TNF superfamily molecule LIGHT 730 promotes the generation of circulating and lung-resident memory CD8 T cells following an acute 731 respiratory virus infection

STAR: ultrafast universal RNA-seq aligner

Apoptosis: a review of programmed cell death

Transcription factor Runx3 is induced by influenza A 736 virus and double-strand RNA and mediates airway epithelial cell apoptosis

From "A"IV to "Z"IKV: attacks from emerging and re-emerging pathogens

Clinical 740 features of patients infected with 2019 novel coronavirus in Wuhan

Angiotensin II plasma levels are linked to disease severity and predict fatal outcomes in 743 H7N9-infected patients

XAF1 forms a positive 745 feedback loop with IRF-1 to drive apoptotic stress response and suppress tumorigenesis

Angiotensin-converting enzyme 2 is a functional 749 receptor for the SARS coronavirus

T-cell immunity of SARS-CoV: 751 Implications for vaccine development against MERS-CoV

Influenza viruses and the NF-kappaB signaling pathway -towards a 753 novel concept of antiviral therapy

IL-6/IL-6 receptor system 755 and its role in physiological and pathological conditions

Transcriptional regulation of germinal center B 758 and plasma cell fates by dynamical control of IRF4

ISG15 in antiviral immunity and beyond

The adaptor molecule

CD4 T cells modulates differentiation of follicular helper T cells during chronic LCMV infection

TNF-induced signaling in apoptosis

Annual estimates of the burden of seasonal influenza 768 in the United States: A tool for strengthening influenza surveillance and preparedness. Influenza 769 Other Respir

Interferon-stimulated genes and their antiviral effector 771 functions

XBP1, downstream of Blimp-1, expands the secretory apparatus and 774 other organelles, and increases protein synthesis in plasma cell differentiation

A physical and regulatory map of host-influenza interactions reveals pathways in H1N1 infection

XAF1 mediates tumor necrosis factor-alpha-induced apoptosis and X-linked inhibitor of 781 apoptosis cleavage by acting through the mitochondrial pathway

Breadth of concomitant immune responses prior to patient 785 recovery: a case report of non-severe COVID-19

A novel coronavirus outbreak of global 787 health concern

Structural and functional basis of SARS-CoV-2 entry by using human ACE2

Clinical features of 69 cases with Coronavirus 792 Disease

Interleukin-6 and its receptors: a highly regulated and 794 dynamic system

Nowcasting and forecasting the potential domestic and 796 international spread of the 2019-nCoV outbreak originating in Wuhan, China: a modelling study

Transcriptomic characteristics of bronchoalveolar lavage fluid and peripheral blood mononuclear 800 cells in COVID-19 patients

High expression of 802 ACE2 receptor of 2019-nCoV on the epithelial cells of oral mucosa

Effective treatment of severe COVID-19 patients with tocilizumab

ChIPseeker: an R/Bioconductor package for ChIP peak 807 annotation, comparison and visualization

Intensive cytokine induction in pandemic H1N1 influenza virus infection accompanied by 39 robust production of IL-10 and IL-6

Combined effects of XAF1 and TRAIL on the 812 apoptosis of lung adenocarcinoma cells

Human T-cell immunity against the 814 emerging and re-emerging viruses

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

Discovery of a novel coronavirus associated with the recent pneumonia outbreak in 820 humans and its potential bat origin. bioRxiv

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

XIAP-associated factor 1 (XAF1), a novel target of p53, enhances 826 p53-mediated apoptosis via post-translational modification