Tocilizumab treatment in severe COVID-19 patients attenuates the inflammatory storm incited by monocyte centric immune interactions revealed by single-cell analysis

40 41 Coronavirus disease 2019 (COVID-19) has caused more than 40,000 deaths 42 worldwide. Approximately 14% of patients with COVID-19 experienced severe 43 disease and 5% were critically ill. Studies have shown that dysregulation of the 44 COVID-19 patients’ immune system may lead to inflammatory storm and cause severe 45 illness and even death. Tocilizumab treatment targeting interleukin 6 receptor has 46 shown inspiring clinical results of severe COVID-19 patients. However, the immune 47 network with Tocilizumab treatment at single cell resolution has not been uncovered. 48 Here, we profiled the single-cell transcriptomes of 13,289 peripheral blood 49 mononuclear cells isolated at three longitudinal stages from two severe COVID-19 50 patients treated with Tocilizumab. We identified a severe stage-specific monocyte 51 subpopulation and these cells centric immune cell interaction network connected by the 52 inflammatory cytokines and their receptors. The over-activated inflammatory immune 53 response was attenuated after Tocilizumab treatment, yet immune cells including 54 plasma B cells and CD8 T cells still exhibited an intense humoral and cell-mediated 55 anti-virus immune response in recovered COVID-19 patients. These results provided 56 critical insights into the immunopathogenesis of severe COVID-19 and revealed 57 fundamentals of effectiveness in Tocilizumab treatment. 58

network with Tocilizumab treatment at single cell resolution has not been uncovered. Next, we explored the transcription factors (TFs) that may be involved in the 152 regulation of inflammatory storm in monocytes. We used SCENIC 15 and predicted 15 153 TFs that may regulate genes which were enriched in severe stage-specific monocyte 154 (Fig 2g). By constructing a gene regulatory network, we found 3 of them, namely ETS2, 155 NFIL3 and PHLDA2 were regulating the cytokine storm relevant genes (Extended Data 156 Fig.2d). In addition, we found the expressions of ETS2, NFIL3 and PHLDA2 and their 157 target genes were enhanced in severe-specific monocyte subpopulation (Fig 2h), 158 suggesting these 3 TFs may regulate the inflammatory storm in monocytes. 159 Given that monocytes in the severe stage may be involved in the regulation of a 160 variety of immune cell types, we used the accumulated ligand/receptor interaction 161 database 16 CellPhoneDB (www.cellphonedb.org) to identify alterations of molecular 162 interactions between monocytes and all the immune cell subsets (Supplementary Table  163 6). We found 15 pairs of cytokines and their receptors whose interaction were 164 significantly altered in severe versus recovery and healthy stages (Fig 3a). Consistent 165 with previous study 4 , we found monocytes interacted with CD4 + T cells and plasma B 166 cells in patients at severe stage through the ligand/receptor pairs of IL-6/IL-6R. This 167 interaction, together with many other cytokine storm relevant cell-cell communications 168 were then extensively attenuated after the treatment of Tocilizumab (Fig 3b), suggesting 169 that this drug may functioning by effectively blocking the inflammatory storm in severe 170

COVID-19 patients. 171
In addition, we also discovered many other ligand/receptor pairs involved in a 172 broader spectrum of immune cell communications enriched at the severe stage. For 173 example, TNF-α and its receptors, which connected monocytes with many other 174 immune cells. Others like IL-1β and its receptor, which connected monocytes with 175 CD8 + T cells. Chemokines, such as CCL4L2, CCL3 and CCL4 and their receptors were 176 also found enriched at severe stage. These cytokines and their receptors may serve as 177 potential drug targets to treat COVID-19 patients at severe stage, and some of their inhibitors are undergoing clinical trials in multiple places around the world 179 (Supplementary Table 7). Collectively, these findings illustrated the molecular basis of 180 cell-cell interactions at the peripheral blood of COVID-19 patients, leading to a better 181 understanding of the mechanisms of inflammatory storm of the disease. 182 Robust multi-factorial immune responses can be elicited against viral infection, 183 such as avian H7N9 disease 17,18 . A recent report has found effective immune responses 184 from a non-severe COVID-19 patient 19 . However, it is not clear whether the anti-virus 185 immune response would be affected after Tocilizumab treatment. Therefore, the anti-186 virus immune responses from the humoral and cell-mediated immunity in COVID-19 187 patients at severe stage were compared with recovery stage and healthy controls. As 188 expected, we found plasma B cells were barely existing in healthy controls (Fig. 4a). 189 By contrast, significantly higher proportion of plasma B cells was exclusively increased 190 in both severe and recovery stages (Fig. 4a,  and a subset of CD8 + T cells with proliferation characteristics (cluster 12) (Fig. 4c, d). 198 Among them, the cells from the severe patients were mainly distributed in the effector 199 CD8 + T cell cluster (Fig. 4c, d). We then conducted pairwise comparisons to identify 200 DEGs of CD8 + T cells between the severe, recovery and healthy stages (  Table 9; P < 10 -10 ). 203 Meanwhile, genes involved in "leukocyte mediated cytotoxicity" and "positive 204 regulation of cytokine production" were highly enriched in CD8 + T cells from COVID-205 19 patients at both severe and recovery stage (Fig. 4g, Supplementary Table 9; P < 10 -206 5 ). We further detected elevated expression of the 108 and 449 genes involved in these 207 GO terms (Fig. 4h, i, Supplementary Table 10). Together, these results demonstrated 208 the critical evidence that robust adaptive immune responses against SARS-CoV-2 209 infection exist in severe stage and remain after Tocilizumab treatment. 210 The immune system is crucial to fight off viral infection 20,21 . Recent studies have 211 illustrated that monocytes may be the main cause of exacerbation and even death of 212 COVID-19 patients through inflammatory storms 4 . In this study, we discovered a 213 specific monocyte subpopulation that may lead to the inflammatory storm in patients We generated single-cell transcriptome library following the instructions of single-cell 260 3' solution v2 reagent kit (10x Genomics). Briefly, after thawing, washing and counting 261 cells, we loaded the cell suspensions onto a chromium single-cell chip along with 262 partitioning oil, reverse transcription (RT) reagents, and a collection of gel beads that contain 3,500,000 unique 10X Barcodes. After generation of single-cell gel bead-in-264 emulsions (GEMs), RT was performed using a C1000 Touch TM Thermal Cycler (Bio-265 Rad). The amplified cDNA was purified with SPRIselect beads (Beckman Coulter). 266 Single-cell libraries were then constructed following fragmentation, end repair, polyA-267 tailing, adaptor ligation, and size selection based on the manufacturer's standard 268 parameters. Each sequencing library was generated with unique sample index. Libraries 269 were sequenced on the Illumina NovaSeq 6000 system. numbers between 500 and 6,000 and mitochondrial UMIs less than 10%. For cells from 282 healthy donors, we retained cells with detected gene numbers between 300 and 5,000 283 and mitochondrial UMIs less than 10%. Subsequently we normalized gene counts for 284 each cell using the "NormalizeData" function of Seurat with default parameters. 285 In downstream data processing, we used canonical correlation analysis (CCA) and 286 the top 40 canonical components to find the "anchor" cells between patients and healthy 287 controls. We then used the "IntegrateData" function in Seurat to integrate all the cells. To search for the differentially expressed genes (DEGs), we first assign the negative 296 elements in the integrated expression matrix to zero. We used Wilcoxon rank-sum test 297 to search for the DEGs between each pair of the 3 stages of cells (i.e. severe stage, 298 recovery stage and healthy control). We applied multiple thresholds to screen for DEGs, 299 including mean fold-change >2, P value <0.001, and were detected in >10% of cells in 300 at least one stage. 301 We defined stage A specific-DEGs as the intersections between the DEGs in stage 302 A versus stage B and the DEGs in stage A versus stage C. We defined stage A and B 303

common-DEGs as the intersections of the DEGs in stage A versus stage C and the 304
DEGs in the stage B versus stage C, minus the DEGs between stage A and B. In this 305 way, we obtained the specific-DEGs for each stage, and the common-DEGs for each 306 pair of the 3 stages. We then uploaded these DEG groups to the Metascape 24 website 307 (https://metascape.org/gp/index.html#/main/step1), and used the default parameters to 308 perform Gene Ontology (GO) analysis for each stage. 309 310

Motif enrichment and regulatory network 311
We adopted SCENIC 15 (version 1.1.2) and RcisTarget database to build the gene 312 regulatory network of CD14 + monocytes. Since the number of CD14 + monocytes from 313 healthy control (N = 9,618) was more than those from the severe and recovery stages 314 (N = 1,607), to balance their contributions in the motif analysis, we randomly sampled 315 2,000 CD14 + monocytes from the healthy control for calculation. We selected 13,344 316 genes that were detected in at least 100 monocytes or included in the DEGs of the 3 317 stages as the input features for SCENIC. With default parameters, SCENIC generated 318 the enrichment scores of 427 motifs. We used the student's t-test to calculate the P values of these motifs between severe stage and healthy control, and selected severe-320 specific enriched motifs with fold change >1.5 and P value < 10 -100 . 321 We then applied the enrichment scores of the severe-specific enriched motifs and 322 the expressions of their targeted genes to Cytoscape 25 to construct a connection map for

Data Availability 339
We are uploading the scRNA-seq data of PBMCs from the 2 severe COVID-19 patients 340 to the Genome Sequence Archive at BIG Data Center and the accession number will be 341 available upon request. We also used the scRNA-seq data of PBMCs from 2 healthy 342 donors, which can be downloaded from the 10X genomics official website.   NK I L 6 _ I L 6 R T N F _ V S I R T N F _ T N F R S F 1 B T N F _ T N F R S F 1 A I L 1 B _ A D R B 2 C C L 2 0 _ C X C R 3 F L T 3 L G _ F L T 3 C C L 4 L 2 _ F F A R 2 C C L 4 L 2 _ V S I R C C L 4 L 2 _ P G R M C 2 C C L 3 L 1 _ C C R 1 C C L 4 _ C N R 2 C C L 4 _ S L C 7 A 1 C C L 3 _ C C R   CD83  IGHG1  IGHM  IGKC  IL1B  IL18   TNFRSF18  MAP3K8  CLEC7A  PIK3R6   FCGR2B  FCGR3A  KIR2DL1  KIR2DL4  KIR3DL1  KIR3DL2  KLRB1  KLRC1  KLRD1  KLRF1   NCR1  NCR3  FASLG  GZMB  PRF1  TNF