key: cord-0786621-e2qmd3kj authors: Ma, Y.; Qiu, F.; Deng, C.; Li, J.; Huang, Y.; Wu, Z.; Zhou, Y.; Zhang, Y.; Xiong, Y.; Yao, Y.; Zhong, Y.; Qu, J.; Su, J. title: Single cell sequencing analysis uncovers genetics-influenced CD16+monocytes and memory CD8+T cells involved in severe COVID-19 date: 2022-02-07 journal: nan DOI: 10.1101/2022.02.06.21266924 sha: 764d0bbd6de0d8f9d4023e9ddd836525a757b726 doc_id: 786621 cord_uid: e2qmd3kj Background: Understanding the host genetic architecture and viral immunity contributes to the development of effective vaccines and therapeutics for controlling the COVID-19 pandemic. Alterations of immune responses in peripheral blood mononuclear cells play a crucial role in the detrimental progression of COVID-19. However, the effects of host genetic factors on immune responses for severe COVID-19 remain largely unknown. Methods: We constructed a powerful computational framework to characterize the host genetics-influenced immune cell subpopulations for severe COVID-19 by integrating GWAS summary statistics (N = 969,689 samples) with four independent scRNA-seq datasets (N = 606,534 cells). Results: We found that 34 risk genes were significantly associated with severe COVID-19, and the number of highly-expressed genetics-risk genes increased with the severity of COVID-19. Three cell-subtypes that are CD16+monocytes, megakaryocytes, and memory CD8+T cells were significantly enriched by COVID-19-related genetic association signals. Notably, three causal risk genes of CCR1, CXCR6, and ABO were specifically expressed in these three cell types, respectively. CCR1+CD16+monocytes and ABO+ megakaryocytes with significant up-regulated genes including S100A12, S100A8, S100A9, and IFITM1 confer higher risk to the cytokine storms among severe patients. CXCR6+ memory CD8+ T cells exhibit a notable polyfunctionality of multiple immunologic features, including elevation of proliferation, migration, and chemotaxis. Moreover, we observed a prominent increase in cell-cell interactions of both CCR1+ CD16+monocytes and CXCR6+ memory CD8+T cells in severe patients compared to normal controls among both PBMCs and lung tissues, and elevated interactions with epithelial cells could contribute to enhance the resident to lung airway for against COVID-19 infection. Conclusions: We uncover a major genetics-modulated immunological shift between mild and severe infection, including an increase in up-regulated genetic-risk genes, excessive secreted inflammatory cytokines, and functional immune cell subsets contributing high risk to severity, which provides novel insights in parsing the host genetics-influenced immune cells for severe COVID-19. cells that had UMIs more than 10,000 and expressed reads of mitochondrial genes greater than 10% 135 were removed). The CellCycleSoring function was applied to remove the effects of confounding 136 factors. Principal component analysis (PCA) was carried out to extract principal components (PCs) 137 that could explain most of datasets via using high variable genes. Top 20 PCs were utilized to 138 conduct uniform manifold approximation and projection (UMAP) to embed the dataset into two 139 dimensions. Subsequently, we constructed a shared nearest-neighbor graph (SNN) using the 140 FindNeighbors function based on the top 20 PCs, and applied a graph-based modularity-141 optimization algorithm from the Louvain method [31] on this SNN for clustering the dataset with 142 the cluster resolution set to 0.5. We used the RunHarmony function with PCA reduction method 143 from harmony R package [32] to integrate samples to correct batch effects. The 144 FindConservedMarkers function in Seurat was implemented to find differential expressed genes 145 for determining cellular identity. Well-defined markers were used to annotate clusters, and 146 uncharacterized clusters in the first round of clustering were extracted to run the second round of 147 clustering (Supplemental Table S2 ). A total of 606,534 cells with 563,856 PBMC cells and 42,678 148 BALF cells were yielded from 300 samples based on the four independent scRNA-seq datasets 149 (Supplemental Table S1 and Figure 1A ). To allow comparison across samples and datasets, we 150 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted February 7, 2022. ; https://doi.org/10.1101 https://doi.org/10. /2022 Page 7 of 39 used a common dictionary of gene symbols to annotate genes and these unrecognized symbols were Hierarchical clustering analysis 172 To examine the similarity of the transcriptome profiles between cell types across different 173 COVID-19 severities (Supplemental Figure S4) , we merged the counts of UMI for each cell type 174 according to normal, mild, moderate, and severe COVID-19. In order to normalize gene expression, 175 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted February 7, 2022. ; https://doi.org/10.1101 https://doi.org/10. /2022 Page 8 of 39 we divided the counts of UMI for each gene by the counts of total UMI for all genes in each cell 176 type and then multiplied by 100,000, as refer to the method in a previous study [18] . A median 177 expression value of greater than 0.5 was used to calculate the relative change in each gene 178 expression by dividing it by the median value for each gene, and the Pearson correlation coefficient 179 (PPC) of the relative change in gene expression was used for current hierarchical clustering analysis. the LD information among SNPs extracted from GWAS summary data on ) , MVN( D . Namely, T follows a mixture distribution of independent 2 1  198 random variables. A total of 19,138 genes were included in the current analysis. We used the 199 Benjamini-Hochberg false discovery rate (FDR) method, in which a gene with a FDR ≤ 0.05 (P ≤ 200 6.8×10 -5 ) was interpreted as significant, to adjust for multiple testing. Pathway enrichment analysis 203 We applied the built-in functions of MAGMA [35] , using the results from GWAS summary 204 statistics as its input, to examine genome-wide enriched biological pathways for severe COVID-19. 205 We calculated competitive P values by examining the results that the combined effect of genes 206 within a pathway is significantly greater than the combined effect of all other genes, and 10,000 207 permutations was used to adjust competitive P values. Additionally, we leveraged the over- and a gene with a value of FDR ≤ 0.05 (P ≤ 3.8×10 -5 ) is considered to be significant. . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted February 7, 2022. ; https://doi.org/10.1101 https://doi.org/10. /2022 In silico permutation analysis 245 To explore the concordance of results from both MAGMA analysis (Gene set #1: N = 944, P 246 ≤ 0.05) and S-MultiXcan analysis (Gene set #2: N =1,274, P ≤ 0.05), we performed an in silico 247 permutation analysis which consisted 100,000 times (N Total) random selections [41, 42] . We first 248 calculated the number of overlapped genes between Gene Set #1 and #2 (N Observation = 302), then 249 employed the total number of genes in S-MultiXcan analysis as background genes (N Background = 250 22,326). By randomly selecting the same number of genes as Gene set #2 (N = 1,274) from the 251 background genes, and after repeating it 100,000 times, we calculated the number of overlapped 252 genes between Gene Set #1 and the sample we selected each time (N Random). Finally, we calculated 253 the empirically permuted P value using the following formula: P = ≥ , and 254 empirical P value ≤ 0.05 is considered to be significant. 257 We conducted a drug-gene interaction analysis for identified genetics-risk genes by using 258 protein-chemical interactions in the context of STRING-based PPI networks [43] and STITCH-259 based drug annotation information (v5.0, http://stitch.embl.de/) [44] . Only experimentally-260 validated gene-drug interactions with ranked confidence score were selected for constructing a 261 drug-gene interaction network. To examine the potential therapeutic effects of highly-expressed 262 genes in each immune cell, we conducted an enrichment analysis of 43 druggable categories based 263 on the DGIdb database (https://www.dgidb.org/druggable_gene_categories). Additionally, we 264 collected 1,263 human druggable proteins, which are therapeutic targets of clinical stage or 265 approved drugs, from a previous study [26] . Among them, 704 proteins are targets for potential 266 COVID-19-relevant drugs based on registers of clinical trials for COVID-19, approved 267 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted February 7, 2022. ; https://doi.org/10.1101 https://doi.org/10. /2022 immunomodulatory/anticoagulant drugs, or have biological functions associated with SARS-CoV-268 2 infection (Supplemental Table S11 SNPs within the 1 Mb window based on the 1,000 Genome Project European Phase 3 panel [33] . 287 We restricted the analysis to SNPs in the autosomes, and any SNPs with MAF ≤ 5% were excluded. . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted February 7, 2022. ; https://doi.org/10. 1101 /2022 Page 13 of 39 The major histocompatibility complex region (Chr6: 25-35 Mbp) was also excluded due to the 289 extensive LD in this region. Defining cell state scores 292 We leveraged cell state scores to assess the immunological degree of each immune cell type 293 expressed a pre-curated expression gene set [11, 14, 29] 311 To identify potential cellular interactions of CCR1 + CD16+monocytes and CXCR6 + memory 312 CD8+T cells with other immune cells, we utilized the CellChat R package [47] for inferring the 313 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted February 7, 2022. As shown in Figure 1 , we devised a computational framework to parse the host genetics-336 modulated immune cell subpopulations implicated in severe COVID-19. It included three main 337 parts: 1) integrative analysis that combined GWAS summary statistics with scRNA-seq data to 338 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted February 7, 2022. ; https://doi.org/10. 1101 /2022 Page 15 of 39 genetically map single-cell landscape for severe COVID-19 ( Figure 1A) ; 2) identifying genetics- Figure S1 ). We leveraged well-known marker genes to assign these clusters to 13 distinct cell types, CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted February 7, 2022. CD4+T cells, and NK cells were significantly decreased with the increased severities 363 (Supplemental Figure S5 ). 366 Through performing a meta-analysis of 21 independent GWAS studies from the COVID-19 367 Host Genetic Consortium, eight genomic loci were identified to be associated with hospitalized 368 COVID-19 at a genome-wide significant level, including 1p22.2 (rs2166172, P = 2.74×10 -8 ), 369 3p21.31 (rs35081325, P = 3.32×10 -58 , and rs33998492, P = 3.59×10 -14 ), 6p21.33 (rs143334143, P 370 = 1.28×10 -10 ), 7p11.2 (rs622568, P = 2.57×10 -8 ), 9q34.2 (rs505922, P = 2.24×10 -9 ), 12q24.13 371 (rs2269899, P = 3.24×10 -8 ), 19p13.3 (rs2109069, P = 6.4×10 -13 ), and 21q22.11 (rs13050728, P 372 =1.91×10 -11 ) ( Figure 2A , Supplemental Table S3, Figure S6 , and Materials S2). Among these eight 373 loci, three loci, 1p22.2, 6p21.33 and 7p11.2, were newly identified. It should be noted that there 374 were two independent genetic association signals (Index SNPs: rs35081325 and rs33998492) in the 375 3p21.31 locus for severe COVID-19 ( Figure 2B and Supplemental Figure S7A -C). Using the 376 Variant2Gene (V2G) algorithm [48], we prioritized CXCR6 as a candidate causal gene for 377 rs35081325 and causal gene CCR1 for rs33998492 (Supplemental Method S1). 378 Furthermore, the index SNP of rs505922 (P = 2.24×10 -9 ) in the 9q34.2 locus is highly LD with 379 the reported SNP of rs657152 (R 2 = 0.874) [27] and rs8176719 (R 2 = 0.876) [25] . Based on the top-380 ranked V2G score for rs505922, we prioritized ABO as a potential causal gene contributing 381 susceptibility to severe COVID-19. By performing a MAGMA gene-level association analysis, we 382 observed that 25 genes including CXCR6, CCR1, IFNAR2, IL10RB, and OAS1 were significantly 383 associated with severe COVID-19 (FDR < 0.05, Supplemental Figure S8 and Table S4 ). GWAS-384 based pathway enrichment analysis revealed that 19 biological pathways, including cytokine-385 cytokine receptor interaction, influenza A, and TNF signaling, were significantly associated with 386 hospitalized COVID-19 (Supplemental Figure S9 and Table S5 ). . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted February 7, 2022. ; https://doi.org/10.1101/2022.02.06.21266924 doi: medRxiv preprint 389 To obtain combined signals from multiple tissues [49], we leveraged S-MultiXcan to meta-390 analyze the tissue-specific associations from 49 tissues in GTEx (see Methods), which showed that 391 the genetically predicted expressions of 16 genes were significantly associated with severe COVID-392 19 (FDR < 0.05, Figure 2C and Supplemental Table S6 ). Of note, 14 of 16 genes (87.5%) were 393 identified to be significant in MAGMA analysis (Supplemental Figure S10A -B). Through 394 conducting S-PrediXcan analysis of blood and lung tissues that were linked with SARS-CoV-2 395 infection, we found eight genes whose genetically-regulated expression were significantly 396 associated with severe COVID-19 (FDR < 0.05, Supplemental Table S7 ). Using in silico 397 permutation analysis, we further observed that there existed a high consistence among results from 398 MAGMA, S-PrediXcan, and S-MultiXcan analyses (P < 1.0×10 -5 , Supplemental Figure S11A -C). 399 The aforementioned multiple genomic analyses identified 34 risk genes that showed supportive 400 evidence of involvement in the etiology of COVID-19 (Supplemental Figure S12A -B). 403 The result of a network-based enrichment analysis suggested that 22 of 34 risk genes were 404 significantly enriched in a PPI subnetwork (P = 2.85×10 -13 , Figure 2D ), which is consistent with 405 the consensus that disease-related genes are more densely connected [50, 51] . To functionally 406 characterize the drug targets of these genes, we conducted a drug-gene interaction analysis and 407 identified 11 genes including CCR1, IFNAR2, IL10RB, and OAS1 were targeted by at least one 408 known drug ( Figure 2D and Supplemental Figure S14 CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted February 7, 2022. [25, 52, 53] . 415 Based on the expression profile of dataset #1, we conducted a hierarchical clustering analysis 416 of these identified risk genes on COVID-19 severity, and found that these risk genes predisposed 417 be highly-expressed in severe patients compared to normal group (Permuted P = 0.023, Figure 2F 425 To identify genome-wide genetics-influenced immune cells for severe COVID-19, we first 426 leveraged a regression-based polygenic model to integrate GWAS summary data on severe 427 COVID-19 with single-cell transcriptomic profiles (dataset #1) according to different COVID-19 428 severities (See methods). We found that CD16+monocytes were significantly associated with three 429 phases of COVID-19, mature B cells showed remarkable associations with mild COVID-19, 430 megakaryocytes were significantly associated with moderate and severe COVID-19, and memory 431 CD8+T cells showed significant associations with severe COVID-19 (simulated P < 0.05, Figure 432 3A). Further, we used a generalized linear regression model to validate these severe COVID-19-433 associated cell types by conditioning on the 10% most specific genes for each type, and consistently 434 found that CD16+monocytes and megakaryocytes showed notable associations with severe 435 COVID-19 (P < 0.05, Supplemental Method S2). These results indicated that CD16+monocytes, 436 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted February 7, 2022. ; https://doi.org/10. 1101 /2022 Page 19 of 39 megakaryocytes, and memory CD8+T cells were more vulnerable to the influence of genetic 437 components on severe-stage patients. 438 Based on the specificity algorithm used in MAGMA, we noticed that the top specific cell type 439 of CCR1 was CD16+monocytes, CXCR6 was most specifically expressed in memory CD8+T cells, 440 and ABO was specific to megakaryocytes (Supplemental Figure S15A) , recalling that CXCR6, 441 CCR1 and ABO were prioritized to be candidate causal genes for severe COVID-19 based on the 442 V2G score in above genetics-based analysis. Compared with other cell types, CCR1 was primarily 443 expressed in CD16+monocytes (24.77%), CXCR6 was mainly expressed in memory CD8+T cells 444 (40.29%), and the ABO-expressed cells were highly specific to megakaryocytes (54.63%) 445 (Supplemental Figure S15B and Table S9 ). To gather additional empirical support, we analyzed the 446 combined dataset of both datasets #2 and #3 as a validation and found CCR1, CXCR6, and ABO 447 showed a consistent specificity in the three cell types (Supplemental Figure S16 ). 448 Given that the primary goal of current study was to characterize genetics-influenced peripheral 449 immune cell types for severe COVID-19, the majority of our subsequent detailed analyses would 450 be concentrated on three immune cell subpopulations: CCR1 + CD16+monocytes, ABO + 451 megakaryocytes, and CXCR6 + memory CD8+T cells ( Figure 3B ). 452 453 CCR1 + CD16+monocytes and ABO + megakaryocytes contributing higher risk to cytokine storm 454 The accumulating lines of evidence [29, 54] have suggested that subsets of monocytes and 455 megakaryocytes might be the major resources of inflammatory storm. We sought to examine 456 whether CCR1 + CD16+monocytes and ABO + megakaryocytes play more important roles in 457 cytokine storm among severe patients. As for CCR1 + CD16+monocytes, we found that the 458 inflammatory cytokine score was significantly higher than that of CCR1 -CD16+monocytes (P = 459 2.5×10 -7 , Figure 4A ). Consistently, the combined score of both cytokine-cytokine receptor 460 interaction and chemokine signaling pathway was prominently higher in CCR1 + CD16+monocytes 461 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted February 7, 2022. ; https://doi.org/10. 1101 /2022 Page 20 of 39 (P < 2.2×10 -16 , Supplemental Figure 17A ). Compared with CCR1 -CD16+monocytes, there were 462 351 significantly highly-expressed genes in CCR1 + CD16+monocytes, such as inflammatory and 463 cytokine genes of IL1B, IL27, CXCL10, CXCL8, CD14, and OSM (FDR < 0.05, Figure 4B and 464 Supplemental Table S11), which have been documented to be associated with the inflammatory 465 response and chemotaxis of immune cells among 15, 55, 56] . Functionally, 466 19 KEGG pathways were significantly overrepresented by the 351 highly-expressed genes (FDR < 467 0.05, Figure 4C and Supplemental Table S12), including cytokine-cytokine receptor interaction and 468 chemokine signaling pathway, reminiscing that most of them were identified in above genetics-469 based pathway analysis. Additionally, these highly-expressed genes among CCR1 + CD16+monocytes have a remarkably higher proportion of druggable genes and COVID-19-471 associated druggable genes (P ≤ 0.01, Supplemental Figure S17 and Table S13 ). 472 The cell percentage of CCR1 + CD16+monocytes showed a notable elevation among moderate 473 and severe patients compared with normal controls (P < 0.001), with no significant difference 474 between mild patients and normal controls (P = 0.1, Figure 4D ). Furthermore, the inflammatory 475 cytokine scores among CCR1 + CD16+monocytes were significantly elevated with increased 476 severities (Trend P = 0.0013, Figure 4E ). In comparison with normal controls, mild, moderate, and 477 severe patients displayed significantly up-regulated expressions (up-DEGs) with 14, 169, and 190 478 genes respectively (FDR < 0.05, Figure 4F and Supplemental Figure S17D ). Notably, there existed 479 a high correlation between up-DEGs of moderate and severe patients (r = 0.937, P < 2.2×10 -16 ; 480 Figure 4G ), such as S100A8, S100A9, and IFITM1 ( Figure 4H-4J) , indicating a similar expression 481 pattern between moderate and severe patients. Accumulating release of massive amounts of 482 calprotectin (S100A8/S100A9) in monocytes contributes to inflammatory response among severe 483 16, 29] . 484 Furthermore, these 190 up-DEGs were significantly enriched in disease-terms associated with 485 viral infection and inflammation and 17 functional GO-terms (FDR < 0.05, Figure 4K , 486 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted February 7, 2022. Figure S17E and Tables S14-S15), including interferon alpha/beta signaling and 487 interferon gamma signaling. These interferon-related genes including IRF3, IRF2, IFI6, IFITM1, 488 ISG15, and ICAM1 may induce autoinflammatory and autoimmune conditions contributing to the 489 innate immune cells against SARS-CoV-2 infection [57, 58] . Of note, a high proportion of 63.68% 490 among 190 up-DEGs such as CXCL8, IFITM1, S100A8, and S100A9 were annotated into 15 491 potential druggable gene categories (Supplemental Figure S17F -L and Table S16 ). These results 492 indicated that interferon-related genes among CCR1 + CD16+monocytes have instrumental effects 493 in exacerbating inflammation among severe patients. 494 In addition, we found that ABO + megakaryocytes had a significantly higher inflammatory 495 cytokine score than that in ABOcells (P < 0.001, Supplemental Figure S18A -B). Compared with 496 ABOmegakaryocytes, 424 genes were significantly highly-expressed in ABO + megakaryocytes 497 (FDR < 0.05, Supplemental Figure S18C and Table S17 ). These 424 highly-expressed genes were 498 significantly enriched in systemic lupus erythematosus, alcoholism, and platelet activation (FDR < 499 0.05, Supplemental Figure S18D and Table S18 ). Similar to CCR1 + CD16+monocytes, the cell 500 percentage of ABO + megakaryocytes was significantly elevated among moderate and severe 501 patients (P < 0.01, Supplemental Figure S18E ). Among ABO + megakaryocytes, 20 and 35 up-DEGs 502 were notably associated with moderate and severe patients, respectively (FDR < 0.05, Supplemental 503 Figure S18F-G). There was a highly overlapped rate of these up-DEGs between moderate and 504 severe COVID-19 groups, including ACP1, S100A8, and A100A9 (18/20 = 90%, Supplemental 505 Figure S18F-N) . These 35 up-DEGs were annotated to 12 druggable gene categories and 506 significantly enriched in several disease terms (Supplemental Figure S18H and Tables S19-S20), 507 such as shock and thrombocytopenia, which were reported to be associated with . 508 Overall, these results suggest that both CCR1 + CD16+monocytes and ABO + megakaryocytes 509 contribute higher risk to inflammatory storm among severe patients. 510 511 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted February 7, 2022. Compared with CXCR6memory CD8+T cells, we found that scores of cytokine, chemokine, IFN-519 ɑ/β response, T cell activation, proliferation, and migration were significantly higher among 520 CXCR6 + memory CD8+T cells (P < 0.05, Figure 5A -D and Supplemental Figure S19A -C). There 521 were 158 highly-expressed genes among CXCR6 + memory CD8+T cells in comparison with 522 CXCR6cells (FDR < 0.05, Figure 5E ). These highly-expressed genes were significantly enriched 523 in two biological pathways of cytokine-cytokine receptor interaction and inflammatory bowel (FDR 524 < 0.05, Supplemental Figure S19D and Table S21 ). The chemokine signaling pathway showed a 525 suggestive enrichment (P < 0.05). These highly-expressed genes contained numerous pro-526 inflammatory cytokine and chemokine genes, such as CCR1, CCR2, CCR5, CCR6, CCL3L1, 527 IFNGR1, IL18R1, IL23R, MYC, and TNFSF14, which may be associated with the activation of 528 memory CD8+T cells. 529 Furthermore, the cell proportion of CXCR6 + memory CD8+T cells was significantly higher 530 among both mild and moderate COVID-19 than that among normal group (P < 0.05), whereas the 531 cell proportion of CXCR6 + memory CD8+T cells among severe COVID-19 was remarkably lower 532 than that among normal group (P = 0.012, Figure 5F ). Consistently, we found that the scores of 533 chemokine, T cell activation, and migration were increased with the increasing patient severities 534 among CXCR6 + memory CD8+T cells (Trend P < 0.05, Figure 5G -I), and that lower cytotoxicity 535 score and exhaustion score were observed among moderate-to-severe patients (Trend P < 0.05, 536 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted February 7, 2022. 548 To gain refined insights into CCR1 + CD16+monocytes and CXCR6 + memory CD8+T cells, we 549 examined the cellular interactions among cell populations in PBMCs and BALFs according to the 550 COVID-19 disease status using the CellChat algorithm [47] . For CCR1 + CD16+monocytes in 551 PBMCs, we found a notable increase in cell-to-cell interactions with other immune cells among 552 severe patients than that in normal controls (P < 0.05, Figure 6A and Supplemental Figure S20 ). 553 There was no statistical difference in cellular communications of CCR1 -CD16+monocytes with 554 other cells between normal and COVID-19 patients (P > 0.05, Figure 6B ). Compared with normal 555 controls, CCR1 + CD16+monocytes showed elevated interactions with megakaryocytes, memory 556 CD8+T cells, NK, effector CD8+T cells, and CD14+monocytes among severe patients 557 (Supplemental Figure S20 ). There were 14 ligand-receptor interactions observed to be remarkably 558 dominated among severe patients ( Figure 6C ), including ANXA1-FPR1, ITGB2-ICAM2/CD226, LGALS9-CD44, SELPLG-SELL/SELP, APP-CD74, and THBS1-CD36/CD47. . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted February 7, 2022. With regard to CXCR6 + memory CD8+T cells in PBMCs, the predicted cell-to-cell interactions 561 showed a prominent elevation with increased severities of COVID-19 (P < 0.05, Figure 6D ). Figure 6J ). Notably, there was a 60% increase in cellular interactions 581 between CCR1 + CD16+monocytes and epithelial cells compared with that of CCR1 -582 CD16+monocytes (Supplemental Figure S21B) . We also found a 33.33% increase in the 583 interactions between CXCR6 + memory CD8+T cells and epithelial cells compared with that of 584 CXCR6memory CD8+T cells (Supplemental Figure S21C) , such as enhanced ligand-receptor 585 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted February 7, 2022. CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted February 7, 2022. [11, 29, 54] showed that the influence caused by monocytes and megakaryocytes in 615 inflammatory storms is noteworthy among severe COVID-19 patients. We found that CCR1 + 616 CD16+monocytes and ABO + megakaryocytes showed a significantly increased propensity to cause 617 inflammatory storms among severe patients. The observations suggest highly-expressed interferon-618 related genes, including S100A8, S100A9, S100A12, CD14, CXCL8, IGSF6, IRF3, IFI6, IFITM1, 619 and IFITM3 among the two cell subsets contribute to exacerbate inflammation among severe 620 patients. The inflammatory mediator of EN-RAGE encoded by S100A12 was significantly 621 correlated with , and S100A8, S100A9, IRF3, IFI6, IFITM1, and IFITM3 have been CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted February 7, 2022. ; https://doi.org/10.1101 https://doi.org/10. /2022 Page 28 of 39 cases [15] . We demonstrated that CXCR6 + memory CD8+T cells mounted highly effective immune 660 responses to against COVID-19, highlighting the remarkable biological plasticity in subsets of 661 memory CD8+T cells differentiation. 662 The power of this study is limited by the lack of matched genetic data and scRNA-seq data in CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted February 7, 2022. . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted February 7, 2022. Competing interests: The authors declare no competing interests. . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted February 7, 2022. ; CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted February 7, 2022. ; https://doi.org/10.1101 https://doi.org/10. /2022 Page 33 CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted February 7, 2022. ; https://doi.org/10.1101 https://doi.org/10. /2022 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted February 7, 2022. ; . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted February 7, 2022. Genetics-related risk pathways S100A8, S100A9, S100A12, IFITM1 • Genetics-influenced immune cell subpopulations for severe COVID-19 CCR1 CCL3L1 CD40LG IFNGR1 IL18R1 IL18RAP IL23R MYC SLAMF1 TNFSF14 CCR7 CD14 CD40 IL2R RGS1 TCF7, HLA-DPA1, STOML2, VMP1, ITM2A, TARS, WNK1, KRAS, TMEM205, IRF7, LAMTOR5, TMPO, PUF60, MYDGF, TRGC2, WAS, TUBA4A, TROVE2 An interactive web-based dashboard to track COVID-19 in 828 real time Characteristics of and Important Lessons From the Coronavirus 830 Disease 2019 (COVID-19) Outbreak in China: Summary of a Report of 72 314 Cases 862 et al: Single-cell landscape of immunological responses in patients with COVID-19 COVID-19 severity correlates with airway epithelium-866 immune cell interactions identified by single-cell analysis Subsets Discriminate Severe from Mild COVID-19 IRF3 and type I interferons fuel a fatal response to 988 myocardial infarction Impaired type I interferon activity and inflammatory 991 responses in severe COVID-19 patients Thrombocytopenia is associated with severe 993 coronavirus disease 2019 (COVID-19) infections: A meta-analysis Multifunctional T-cell analyses to study response and progression 997 in adoptive cell transfer immunotherapy Origin and differentiation of human 1000 memory CD8 T cells after vaccination Granzyme H destroys the function 1002 of critical adenoviral proteins required for viral DNA replication and granzyme B 1003 inhibition Multi-1005 platform omics analysis reveals molecular signature for COVID-19 pathogenesis, 1006 prognosis and drug target discovery MHC class II transactivator CIITA induces cell 1009 resistance to Ebola virus and SARS-like coronaviruses CXCR6 regulates localization of tissue-resident 1012 memory CD8 T cells to the airways Interstitial-resident memory CD8(+) T cells sustain 1015 frontline epithelial memory in the lung Relationship between the ABO Blood Group and the COVID-19 Susceptibility. Clin 1018 Infect Dis Incidence of thrombotic 1021 complications in critically ill ICU patients with COVID-19 Associated with COVID-19 Pneumonia Detected with Pulmonary CT Angiography Sequence-based prediction of SARS-1028 CoV-2 vaccine targets using a mass spectrometry-based bioinformatics predictor 1029 identifies immunogenic T cell epitopes Inflammatory monocytes activate memory inflammation CD14+,CD16+ blood monocytes and joint inflammation 1037 in rheumatoid arthritis STAT2 deficiency and susceptibility 1040 to viral illness in humans Understanding memory CD8(+) T cells c-Myc is a universal amplifier of expressed genes in lymphocytes and 1045 embryonic stem cells Phenotypical and functional alteration of unconventional 1048 T cells in severe COVID-19 patients Heritability enrichment of specifically expressed genes 1051 identifies disease-relevant tissues and cell types Ziegler-Heitbrock L: The CD14+ CD16+ blood monocytes: their role in infection and