Single-cell analyses and host genetics highlight the role of innate immune cells in COVID-19 severity

Mechanisms underpinning the dysfunctional immune response in severe acute respiratory syndrome coronavirus 2 infection are elusive. We analyzed single-cell transcriptomes and T and B cell receptors (BCR) of >895,000 peripheral blood mononuclear cells from 73 coronavirus disease 2019 (COVID-19) patients and 75 healthy controls of Japanese ancestry with host genetic data. COVID-19 patients showed a low fraction of nonclassical monocytes (ncMono). We report downregulated cell transitions from classical monocytes to ncMono in COVID-19 with reduced CXCL10 expression in ncMono in severe disease. Cell–cell communication analysis inferred decreased cellular interactions involving ncMono in severe COVID-19. Clonal expansions of BCR were evident in the plasmablasts of patients. Putative disease genes identified by COVID-19 genome-wide association study showed cell type-specific expressions in monocytes and dendritic cells. A COVID-19-associated risk variant at the IFNAR2 locus (rs13050728) had context-specific and monocyte-specific expression quantitative trait loci effects. Our study highlights biological and host genetic involvement of innate immune cells in COVID-19 severity.

The dysregulated IFN response has been suggested in COVID- 19 (refs. 8,14,15,29), leading us to evaluate IFN response across PBMCs. We first defined a score of response to Type I IFN and that of response to IFN-γ for each cell based on the expression of a published gene list. We observed high expressions of Type I IFN response genes in NK and CD16 + monocytes and higher response for moderate COVID-19 compared to the other conditions, especially in monocytes and DC (Fig. 1e, Extended Data Fig. 1f and Supplementary Table 3). We also found high expressions of IFN-γ response genes in CD16 + monocytes and conventional dendritic cells (cDC), and higher response for moderate COVID-19, especially in CD16 + monocytes and pDC (Fig. 1f, Extended Data Fig.  1g and Supplementary Table 3). These data are consistent with the previous reports that systemic IFN response is higher for nonsevere disease 14,29 , suggesting the potential importance of innate immune cells in immunopathology of COVID-19.

Monocytes and DC in COVID-19
To characterize transcriptome dynamics in 116,944 cells annotated as monocytes and DC across disease conditions, we subclustered monocytes into five subsets and DC into two subsets according to the expression of canonical gene markers as follows: cMono_S100A, cMono_IL1B, cMono_CCL3, intermediate monocytes (intMono), ncMono, cDC and pDC ( Fig. 2a and Extended Data Fig. 2a).
To reveal the compositional changes between COVID-19 and healthy controls in this subset, we performed differential abundance analysis using Milo 27 , identifying 6,721 neighborhoods, of which 1,265 showed evidence of differential abundance (Fig. 2b). The proportion of ncMono, cDC and pDC declined prominently in COVID-19 patients compared to healthy controls ( Fig. 2c and Extended Data Fig. 2b). Next, we evaluated differential abundance with disease severity, identifying 3,003 neighborhoods, of which 19 showed significant differential abundance (Extended Data Fig. 2c). While few neighborhoods showed significant differential abundance, there was a trend toward a decrease in pDC and ncMono in severe compared to moderate COVID-19 (Extended Data Fig. 2d), implying that ncMono might contribute to immunopathology of COVID-19 severity because decreased cell proportion of ncMono is one of the COVID-19-specific features 13,30,31 .
To gain insight into functions of different cell subsets according to COVID-19 severity, we performed differential expression analyses and a set of Gene Ontology (GO) analyses across the groups (that is, all COVID-19 patients versus healthy, moderate disease versus healthy and severe disease versus healthy) in the five subsets (cMono, intMono, ncMono, cDC and pDC). The top 20 enriched pathways upregulated in COVID-19 versus healthy group were related to innate immunity or antiviral response, and almost all were shared among the five subsets (Extended Data Fig. 2e). We next compared the enriched pathways upregulated in moderate severity versus healthy group and those in severe versus healthy group. The top ten enriched pathways were also almost the same in each subset, while there existed several pathways with different enrichment patterns between the moderate disease and severe group (Fig. 2d). Although 'response to IFN-γ' pathway was enriched in moderate COVID-19 of each subset, it was not evident in three monocyte subsets in severe COVID-19, suggesting that decreased IFN-γ response in monocyte subsets might contribute to the severity of COVID-19. Enrichment of 'response to Type I IFN' pathway was specifically depleted in ncMono of severe COVID- 19.
To further elucidate the mechanism of COVID-19 severity, we conducted differential expression analysis between moderate disease versus severe disease group in ncMono. PLD4, which digests ssRNA and ssDNA 32 , was the most downregulated in severe compared to moderate disease group (Fig. 2e). In addition, the expression of CXCL10, which belongs to IFN-γ-induced gene and critical in response to various infectious pathogens 33 and has been reported to be involved in COVID-19 us important aspects of the immunopathology of COVID-19, the immune response of the host to SARS-CoV-2 still remains unclear.
In addition, genome-wide association studies (GWASs) of COVID-19, highlighted as an achievement by the COVID-19 host genetics initiative (HGI), have demonstrated that the host genetic backgrounds influence susceptibility to and/or severity of . Multiple genetic variants associated with the COVID-19 risk are shared across different populations, while population-specific risk variants also have been reported 20,21 . Considering the different prognoses by ancestry 22 , integrated analysis of transcriptome and genetic data at single-cell resolution by ancestry should provide new insights.
Here we performed a detailed scRNA-seq analysis of over 895,000 peripheral blood mononuclear cells (PBMCs) of 73 patients with COVID-19 as well as 75 healthy controls of Japanese ancestry, and then context-specific and cell type-specific expression quantitative trait loci (eQTL) analysis by integrating scRNA-seq and host genetics data.
The proportion of nonclassical monocytes (ncMono) decreased in COVID-19 patients and RNA velocity analysis revealed the downregulation of the cellular transitions from classical monocytes (cMono) to ncMono in COVID-19 patients. We found that CXCL10 expression was downregulated in ncMono during severe COVID-19, and cell-cell communication analysis inferred that the cellular interactions involving ncMono and plasmacytoid dendritic cells (pDC) were reduced in severe COVID-19. The putative disease genes identified by the GWAS of severe COVID-19 were enriched in monocytes and dendritic cells (DC), and COVID-19-associated variants had context and cell type-specific eQTL effects, with the IFNAR2 variant (rs13050728) in particular having COVID-19-specific and monocytes-specific eQTL effect. In summary, our data linked innate immune cell dysfunction, especially ncMono, with severe COVID-19 and demonstrated the enrichment of host genetic risk in innate immune cells, indicating biological and host genetic critical involvement of innate immune cells in COVID-19 severity.

Single-cell transcriptional profiling of COVID-19 PBMC
To investigate the immunopathogenesis and host genetics mechanism of SARS-CoV-2 infections in COVID-19 patients, PBMCs were collected from 73 COVID-19 patients and 75 healthy controls of Japanese ancestry at Osaka University. The 73 patients with COVID-19 were classified into two conditions as follows: moderate (n = 9) and severe (n = 64) disease according to the WHO classification 23 (Fig. 1a and Supplementary Table 1). No significant difference was noted in age distribution and sex composition between moderate and severe disease groups. The clinical characteristics are summarized in Supplementary Table 1.
After the unified single-cell analysis pipeline (Methods), we obtained 895,460 high-quality cells from PBMCs of all the samples. Cells were manually annotated based on the RNA expression of known marker genes to discriminate subpopulations 6,7,24,25 . We defined 13 cell subsets ( Fig. 1b and Fig. 1b).
To reveal the compositional changes between COVID-19 and healthy controls, we applied Milo 27 , identifying 39,170 neighborhoods, of which 21,279 showed evidence of differential abundance (FDR < 0.1, Fig. 1c). We found a prominent decrease in T cells, NK cells and DC for COVID-19 patients compared to healthy controls, consistent with the previous reports 13,28 (Fig. 1d and Extended Data Fig. 1c). Next, we evaluated differential abundance with disease severity, identifying 14,350 neighborhoods, of which 76 showed significant differential abundance (Extended Data Fig. 1d). While few neighborhoods showed significant differential abundance, there was a trend toward increased B cells and  Fig. 1a). c, Graph representation of Nhoods identified by Milo. Nodes are Nhoods, colored by their log 2 FC between COVID-19 patients (n = 72) and healthy controls (n = 75) adjusted by age and sex. Nondifferential abundance Nhoods (FDR ≥ 0.1) are colored white, and sizes correspond to the number of cells in a Nhood. Graph edges depict the number of cells shared between adjacent Nhoods. d, Beeswarm plot showing the distribution of adjusted log 2 FC in abundance between COVID-19 patients and healthy controls in Nhoods according to 13 cell types. Colors are represented in the same way as in c. e,f, The module score of Type I IFN response and IFN-γ response in PBMCs. The score was calculated using a gene set termed 'GOBP_ RESPONSE_TO_TYPE_I_INTERFERON' (GO:0034340) and 'GOBP_RESPONSE_TO_ INTERFERON_GAMMA' (GO:0034341), respectively. Upper heatmaps depicting the difference between average scores of 13 cell types and that of all single cells. The module scores of cells in each cell type were compared with the average score of all PBMCs using a two-sided one-sample t-test. Lower heatmaps depicting the difference between average scores of moderate or severe disease group and those of the healthy group in each of 13 cell types (n = 75 healthy, n = 9 moderate, n = 64 severe). The module scores of cells of moderate or severe disease group were compared to those of healthy group in each cell type using a two-sided Welch's t-test. *P uncorrected < 1.0 × 10 −50 , **P uncorrected < 1.0 × 10 −300 . Mono, monocytes; Pro, proliferative; Nhood, neighborhood.  Article https://doi.org/10.1038/s41588-023-01375-1 severity with proteomics analysis 5,13 , was also prominently decreased in severe group.
To analyze how the dynamics of transcriptional activation and cell transition differ in disease status, we performed RNA velocity analysis 34 . The transition potential from intMono to ncMono was observed in healthy controls, while such transition was not observed in COVID-19 (Fig. 2f). To quantitatively compare the differences of estimated cell transition from cMono to ncMono through intMono 35 between COVID-19 and healthy controls, we analyzed the increase of unspliced fractions from cMono to ncMono. We found that condition-specific regression slope was lower in COVID-19 patients with significant interaction between three monocytes and two clinical conditions (P Interaction = 2.0 × 10 −4 ; Fig. 2g and Extended Data Fig. 2f). These data suggest that the decreased proportion of ncMono in COVID-19 patients is a consequence of the downregulation of the cellular transition from cMono to ncMono.

T cells and NK cells and T cell repertoires in COVID-19
We subclustered 628,715 cells manually annotated as T cells and NK cells from PBMCs and obtained 13 subsets according to RNA expression of canonical markers ( Fig. 3a and Extended Data Fig. 3a).
To reveal the compositional changes between COVID-19 and healthy controls, we performed differential abundance analysis using Milo 27 , identifying 27,182 neighborhoods, of which 10,035 showed evidence of differential abundance (Fig. 3b). The proportion of almost all clusters declined in COVID-19 patients ( Fig. 3c and Extended Data Fig.  3b). We next evaluated differential abundance with disease severity, identifying 7,981 neighborhoods, of which nine showed significant differential abundance (Extended Data Fig. 3c). Although there were few neighborhoods showing significant differential abundance, the proportions of CD4T cells, regulatory T (T reg ) cells and CD56 bright NK (NK_CD56 bright ) cells tended to be lower in severe COVID-19 compared to moderate COVID-19, while those of natural killer T (NKT) cells and NK tended to be higher (Extended Data Fig. 3d).
To gain insight into the clonal relationship among individual T cells across three disease conditions, we performed T cell receptor (TCR) analysis for nine subsets except for γδ T (gdT) cells, NKT, NK and NK_CD56 bright . The detection percentage of TCRs was 73.1% ( Fig. 3d and Extended Data Fig. 3e). Similar to the previous report 6,36 , the large clonal expansions were observed particularly in CD4 + effector T (CD4_Ef) cells, CD8 + effector T (CD8_Ef) cells and mucosal-associated invariant T (MAIT) cells (Fig. 3e,f). The proportion of clonally expanded CD4_Ef increased in moderate disease compared to healthy and severe group ( Fig. 3d and Supplementary Table 4), suggesting that efficient clonal expansion of CD4_Ef might contribute to the prevention of serious COVID-19. We next examined whether expanded clonotypes were shared among each sample, between COVID-19 and healthy controls, and among cell types. The large expanded clonotypes exhibited parallel expansion in several cell types, particularly in CD4_Ef, CD8 + memory T (CD8_memory) cells and CD8_Ef, whereas a vast majority of these expanded clonotypes were unique for individual patients ( Fig. 3f and Extended Data Fig. 3f,g). In contrast, the large expanded clonotypes in MAIT were shared among individuals.
We compared each clonotype from our data against the currently known CDR3 sequences of SARS-CoV-2-specific TCR in VDJdb 37 . A small number of our TCRs (n = 4,143; 0.8%) shared their CDR3 with reported SARS-CoV-2-specific TCRs (Fig. 3g). The distribution of TCR specific to SARS-CoV-2 was relatively uniform across cell types, and the clonal expansions of such TCRs were observed in memory or activated T cells ( Fig. 3g and Extended Data Fig. 3h). Considering the low percentage of SARS-CoV-2-specific TCRs and the low sharing of expanded TCRs among individuals, further accumulation of data on SARS-CoV-2-specific TCRs would be warranted.

B cells and B cell repertoires in COVID-19
We subclustered 123,728 cells manually annotated as B cells and PB from PBMCs and obtained four subsets according to RNA expression of canonical markers ( Fig. 4a and Extended Data Fig. 4a).
To reveal the compositional changes between COVID-19 and healthy controls, we performed differential abundance analysis using Milo 27 , identifying 8,169 neighborhoods, of which 2,120 showed evidence of differential abundance (Fig. 4b). The proportion of B_naive and B_plasma increased in COVID-19 patients, whereas that of B_memory and B_activated decreased ( Fig. 4c and Extended Data Fig. 4b). We also evaluated differential abundance with disease severity, identifying 4,006 neighborhoods, of which none showed significant differential abundance. However, there was a trend toward a higher proportion of B_plasma in severe compared to the moderate disease group, while that of B_activated tended to be lower (Extended Data Fig. 4c,d).
We next performed differential expression and pathway enrichment analysis in four subsets. Pathway enrichment analysis showed that pathways related to antiviral response and immune response were enriched in COVID-19 patients compared to healthy controls (Extended Data Fig. 4e), and that the biological pathways related to IFN response were enriched in three subsets (B_naive, B_memory and B_activated) in the moderate disease group, while those pathways were not in severe group (Extended Data Fig. 4f), consistent with a previous report 5 . Conversely, pathways related to Type I IFN were enriched in B_plasma in severe group, but not in the moderate severity group (Extended Data Fig. 4f). We also conducted DE analysis between moderate versus severe COVID-19. HERC5, which has direct antiviral function by catalyzing ISGylation 38 , was significantly downregulated in B_naive and B_memory of severe compared to moderate COVID-19 (Fig. 4d). IFN-related genes, such as ISG15, IFITM1 and GBP4, were also downregulated in B_naive and B_activated of severe compared to moderate COVID-19 (Fig. 4d).
Finally, we performed B cell receptor (BCR) analysis. The detection percentage of BCRs was more than 80% in B cell subset (Fig. 4e,f). Clonal expansions were most evident in B_plasma, which showed the larger expansion in COVID-19 (Fig. 4e,f and Supplementary Table 5). We next examined whether expanded clonotypes were shared among each sample, between COVID-19 and healthy controls, and among cell types ( Fig. 4g and Extended Data Fig. 4g,h). In contrast to TCR analysis, very few clonotypes were shared between COVID-19 and healthy controls ( Fig. 4g and Extended Data Fig. 4g). In addition, while the expanded clonotypes in B_plasma were not shared with the other B cell subsets, the expanded clonotypes exhibited parallel expansion among B_naive, B_memory and B_activated (Fig. 4g). Together, our results suggest that in COVID-19 disease, a robust antibody response characterized by clonally expanded circulating PB occurs against a background of augmented IFN responses.

Changes of intercellular interactions in PBMC across clinical status
To map the cellular interaction differences between COVID-19 and healthy controls, we inferred all possible intercellular communications by the expression of ligand-receptor pairs in both cell populations using CellPhoneDB 39 and NATMI 40 . CellPhoneDB and NATMI revealed strong interactions particularly among monocytes and DC in both COVID-19 and healthy groups ( Fig. 5a and Extended Data Fig. 5a).
In addition to simple edge count analysis, we examined the differences in the cell-connectivity-summary networks based on mean expression weight between COVID-19 and healthy controls using NATMI 40 . The cellular interactions involving pDC as the sender in COVID-19 patients were lower than in healthy controls, and those involving B_plasma as a receiver in COVID-19 patients were lower than in healthy controls (Fig. 5b).    We next investigated the differences in the cell-connectivitysummary networks based on mean expression weight between moderate and severe COVID-19. Almost all intercellular interactions were reduced in severe compared to moderate disease group (Fig. 5c and Extended Data Fig. 5b). Notably, the intercellular interactions signaling from pDC and those signaling to ncMono and pDC were reduced in severe COVID-19, implying that the dysfunction of the intercellular interactions involving these two subsets might be related to the severity of COVID-19.
Finally, we explored the intercellular interactions centered on CXCL10 in ncMono, which was significantly downregulated in severe COVID-19 (Fig. 2e). CXCL10 is an IFN-γ-induced gene and exerts its biological effects by binding to CXCR3 (ref. 33). Therefore, we investigated IFNG/IFNGR interactions of ncMono as receiver and CXCL10/CXCR3 interaction of ncMono as sender using NATMI (Fig. 5d). Activated T cells and NK cells showed the strong interactions of IFNG/IFNGR with ncMono, which were enhanced as the severity of COVID-19 increased. On the other hand, pDC and to a lesser extent activated T cells showed the CXCL10/CXCR3 interaction with ncMono in moderate COVID-19, which was not observed in healthy controls and severe COVID-19 (Fig. 5d). These intercellular interaction analyses computationally infer the possibility that dysfunction of ncMono and the consequently reduced interaction of CXCL10/CXCR3 might be one of the factors responsible for driving COVID-19 severity.

Host genetics risks of the severity of COVID-19 in PBMC
Elucidating interaction between host genetics and transcriptional dynamics resolves causal biological mechanism of infection. To evaluate genome-wide host genetics risk of COVID-19 and identify subpopulations of disease-associated cells in PBMCs, we integrated information from our scRNA-seq data with polygenic signals from COVID-19 GWAS using scDRS 41 .
First, we computed a disease score for each cell observed at our COVID-19 scRNA-seq datasets according to the COVID-19 case-control GWAS summary statistics of COVID-19 HGI (round 6; ref. 18), and projected the scores on Uniform manifold approximation and projection (UMAP; Extended Data Fig. 6a). The disease scores at the COVID-19 scRNA-seq dataset were similar across any cell types when using the GWAS summary statistics of self-reported COVID-19 (C2, n Case = 112,612), while the cells annotated as monocytes and DC showed the higher disease scores than the other cell types when GWAS cases were restricted to severe ones (that is, hospitalized COVID-19 (B2, n Case = 24,274) and very severe COVID-19 (A2, n Case = 8,779); Fig. 6a). When comparing disease score for each phenotype, disease scores from severe GWASs were higher than those from self-reported GWAS, particularly at cells annotated as monocytes and DC, which was more prominent in very severe GWAS (Fig. 6b, Extended Data Fig. 6b and Supplementary Table 6). Next, we assessed associations between the six major cell types (Methods) and three COVID-19 GWAS phenotypes, and also within-cell type association heterogeneity using scDRS 41 . No cell type was enriched in the self-reported infection GWAS, whereas monocytes were associated with very severe GWAS, and DC was associated with hospitalization and very severe GWASs (Fig. 6c and Supplementary Table 7), demonstrating that polygenic risks involved in the severity of COVID-19 were enriched in the cells responsible for innate immunity. In addition, three cell types showed heterogeneity in association with hospitalization and critical illness (Fig. 6c). To investigate the host genetics association of monocytes and DC with severe COVID-19 in more detail, we examined cell type-disease association and its heterogeneity of the five innate immune subsets (cMono, intMono, ncMono, cDC and pDC). All of these subsets were associated with GWAS for severe disease, with stronger associations in very severe cases ( Fig. 6c and Supplementary Table 7). Some of them showed significant heterogeneity within the subset (Fig. 6c and Supplementary Table 7). The subset analysis of manually annotated 25 clusters revealed the significant association of MAIT with GWAS for severe disease in addition to monocytes and DC subsets (Extended Data Fig. 6c and Supplementary Table 7). MAIT functions as innate sensors of viral infection 42 , again implicating the involvement of host genetic risk of severe COVID-19 with innate immunity.

Context and cell-type-specific eQTLs of COVID-19-related variants
To gain a better understanding of transcriptional variability and dynamics regulated by the GWAS-identified COVID-19-associated variants, we examined eQTL effects of the replicated variants at the GWAS of COVID-19 in the Japanese population 20 , separately for COVID-19 (n = 67) and healthy controls (n = 75; Extended Data Fig. 7a).
First, we performed eQTL analysis for the six major cell types (Methods). COVID-19-associated variants showed different cell type distributions with significant eQTL effects between the COVID-19 patients and healthy controls, demonstrating context-specific eQTL effects ( Fig. 7a and Supplementary Table 8). Among them, monocytes of the COVID-19 patients had significant eQTL effects in multiple variants (FDR < 0.02). Given the previous analyses demonstrating the involvement of monocytes with COVID-19 severity in this dataset, we examined eQTL effects in the three subsets of monocytes. The two variants with eQTL effects in monocytes of the COVID-19 patients had significant eQTL effects specifically in cMono (FDR = 2.6 × 10 −6 for ABO and FDR = 0.017 for IFNAR2), and no such eQTL effect of the IFNAR2 variant was observed in the healthy controls (FDR = 0.66; Fig. 7b and Supplementary Table 9).

Discussion
Here we reported comprehensive analyses of single-cell transcriptome and TCR and BCR of PBMCs from the COVID-19 patients and healthy controls in Japanese, integrated with host genetics data. Our data presented the dysfunction of monocytes or DC, particularly ncMono, in severe COVID-19, the enrichment of host genetics COVID-19 risks in monocytes or DC, as well as COVID-19 context-specific eQTL effect of the IFNAR2 variant (rs13050728) in monocytes. These highlighted the biological and host genetic critical involvement of innate immune cells in COVID-19 severity.
While previous studies on IFN and SARS-CoV-2 have focused on Type I IFN due to their robust capacity to interfere with viral replication 43 , several studies have indicated that IFN-γ is also an essential component in the severity of COVID-19 (refs. 29,44,45). We found the depleted enrichment of IFN-γ response pathways in monocyte subsets Article https://doi.org/10.1038/s41588-023-01375-1 from patients with severe COVID-19, again demonstrating the importance of IFN-γ response for COVID-19. One of the COVID-19-specific features is decreased cell fraction of ncMono in severe disease 13,31 , whereas its fraction generally increases in patients with sepsis and inflammatory disease 46 . This phenomenon was also observed in this study, and we revealed the downregulation of the cellular transitions from cMono to ncMono in COVID-19 patients using RNA velocity analysis. In addition to decreased cell fraction, differential expression analysis revealed the severely downregulated expression of CXCL10, which is IFN-γ-induced gene and reported to be involved in COVID-19 severity with proteomics analysis 5,13 , in ncMono for severe COVID-19, and cell-cell communication analysis inferred the possibility that CXCL10/CXCR3 interaction between ncMono and pDC was depleted whereas ncMono firmly received IFN-γ signal from activated T cells in severe COVID-19. Thus, these findings indicate that ncMono might contribute to immunopathology of COVID-19 severity via decreased cell fraction as well as biological dysfunction. However, the mechanisms of the reduced differentiation to ncMono and the dysfunction  Human genetic background has been demonstrated to influence the pathogenesis of COVID- 19 (refs. 17-20). Functional analysis has mostly focused on the variant at LZTFL1 on 3p21, which showed the strongest severity association in Europeans but conferred a rare frequency of the risk allele in East Asian 47,48 , suggesting the importance of additional studies in non-Europeans. We found that the putative disease genes identified by the GWAS of severe phenotypes showed cell type-specific expressions in monocytes and DC. We also showed that eQTL effects of the COVID-19-associated variants replicated in the Japanese GWAS were context and cell type-specific, with the IFNAR2 variant, in particular, having COVID-19-specific and monocytes-specific eQTL effect, indicating the host genetic involvement of innate immune cells in COVID-19 severity. Because the context-specific eQTL effects were confirmed, collecting more COVID-19 cases and comparing them with healthy controls will help to clarify the pathogenesis of COVID-19 from the perspective of the host genome.
Collectively, our results motivate us for a detailed examination of ncMono function in the context of COVID-19, as well as to increase sample size to perform integrated analysis with genetic data on a larger scale.

Online content
Any methods, additional references, Nature Portfolio reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41588-023-01375-1.

scRNA-seq computational pipelines and analysis
The R package Seurat (v4.1.0) was used for data scaling, transformation, clustering, dimensionality reduction and most visualization 26 . Data were scaled and transformed using the SCTransform() function, and linear regression was performed to remove unwanted variation due to cell quality (% mitochondrial reads). For integration, we identified 3,000 shared highly variable genes (HVGs) using SelectIntegra-tionFeatures() function. Principal component analysis (PCA) was run on gene expression, followed by batch correction using harmony (v0.1) 53 . UMAP dimension reduction was generated based on the first 30 harmony-adjusted principal components 54 . A nearest-neighbor graph using the first 30 harmony-adjusted principal components was calculated using FindNeighbors() function, followed by clustering using FindClusters() function. Cellular identity was determined by finding differentially expressed genes (DEGs) for each cluster using FindMarkers() function with parameter 'test.use=wilcox', and comparing those markers to known cell type-specific genes (Extended Data Fig. 1a). At the first round of clustering, we identified 13 cell subsets (Fig. 1b). To identify clusters within each major cell type, we performed a second round of clustering on monocytes/DC (CD14 + monocytes, CD16 + monocytes, cDC and pDC), T/NK cells (CD4T, T reg , CD8T, MAIT, proliferative T cells (Pro_T) and NK) and B cells (B and PB), separately, and then we obtained 25 cell subsets (Extended Data Figs. 1a, 2a, 3a and 4a). The final annotation of PBMCs was compared to the PBMC annotation of Azimuth using Seurat 26 . For each of the 25 clusters in our data, the percentage of cells mapped to each Azimuth annotation was calculated (Extended Data Fig. 1b).

Differential abundance analysis
We used Milo 27 (v1.2.0) to test for the differential abundance of cells within defined neighborhoods, between two conditions (that is, COVID-19 versus healthy controls or moderate COVID-19 versus severe . We first used the buildGraph function to construct a KNN graph with k = 30, using 30 principal components (d = 30). Next, we used the make neighborhoods function to assign cells to neighborhoods based on their connectivity over the KNN graph. For computational efficiency, we subsampled 5% for all PBMCs and T cells and 10% for mono/DC and B cells. To test for differential abundance, Milo fit an NB GLM to the counts for each neighborhood, accounting for different numbers of cells across samples using TMM normalization, and use the QL F-test with a specified contrast to compute a P value for each neighborhood. We included age, sex, days since symptom onset and duration of systemic steroids treatment at the time of specimen collection (the last two variables included only in moderate COVID-19 versus severe COVID-19) as covariates in testNhoods function. To control for multiple testing, we adapted the spatial FDR implemented in Milo and used 10% spatial FDR as a threshold for significance. The spatial FDR and log 2 foldchange (FC) of number of cells between two conditions in each neighborhood were used for visualization.

IFN response scoring
To evaluate IFN response across PBMCs, we downloaded a gene set termed 'GOBP_RESPONSE_TO_TYPE_I_INTERFERON (GO:0034340)' and 'GOBP_RESPONSE_TO_INTERFERON_GAMMA (GO:0034341)' from MSigDB. IFN response scores were evaluated using AddModuleScore() function implemented in Seurat with default parameter. To identify cell

Ethics and specimen collection of PBMC for scRNA-seq
Peripheral blood samples were obtained from patients with COVID-19 (n = 73) and healthy controls (n = 75) at Osaka University Hospital. Patients with COVID-19 were further categorized into groups of moderate (n = 9) and severe (n = 64) according to disease severity based on the highest score on the World Health Organization (WHO) Ordinal Scale for Clinical Improvement ever-present (WHO, R&D Blueprint-new coronavirus-COVID-19 Therapeutic Trial Synopsis, 2020; ref. 23). Almost all cases were patients who were transferred from nearby general hospitals because of severe or potentially severe illness during treatment and already initiated with systemic corticosteroids therapy at others hospitals according to RECOVERY study 49 . The detailed clinical data are summarized in Supplementary Table 1. Part of the subjects (n COVID-19 = 30, n Control = 31) are described elsewhere 20 . One patient with COVID-19 had the karyotype abnormality and was excluded from the analyses including sex as a covariate. This study strictly follows the principles according to the Declaration of Helsinki, with written informed consent obtained from all participants before sample collection according to regular principles. Ethical approvals were gained from the Ethics Committees of Osaka University. There was no compensation for participants.

Preparation of single-cell suspensions
For both patients with COVID-19 and healthy controls, blood was collected into heparin tubes and PBMCs were isolated using Leucosep (Greiner Bio-One) density gradient centrifugation according to the manufacturer's instructions. Blood was processed within 3 h of collection for all samples, and stored at −80 °C until use.

Droplet-based single-cell sequencing
Single-cell suspensions were processed through the 10X Genomics Chromium Controller following the protocol outlined in the Chromium Single Cell V(D)J Reagent Kits (v1.1 Chemistry) User Guide. Chromium Next GEM Single Cell 5′ Library & Gel Bead Kit v1.1 (PN-1000167), Chromium Next GEM Chip G Single Cell Kit (PN-1000127) and Single Index Kit T Set A (PN-1000213) were applied during the process. Oil droplets of encapsulated single cells and barcoded beads (GEMs) were subsequently reverse-transcribed in a Veriti Thermal Cycler (Thermo Fisher Scientific), resulting in cDNA tagged with a cell barcode and unique molecular index (UMI). Next, cDNA was then amplified to generate single-cell libraries according to the manufacturer's protocol. Quantification was made with an Agilent Bioanalyzer High-Sensitivity DNA assay (Agilent, High-Sensitivity DNA Kit, 5067-4626). Subsequently amplified cDNA was enzymatically fragmented, end-repaired and polyA tagged. Cleanup/size selection was performed on amplified cDNA using SPRIselect magnetic beads (Beckman-Coulter; SPRIselect, B23317). Next, Illumina sequencing adapters were ligated to the size-selected fragments and cleaned up using SPRIselect magnetic beads. Finally, sample indices were selected and amplified, followed by a double-sided size selection using SPRIselect magnetic beads. Final library quality was assessed using an Agilent Bioanalyzer High-Sensitivity DNA assay. Samples were then sequenced on an Illumina NovaSeq 6000 as paired-end mode.

Alignment, quantification and quality control of scRNA-seq
Droplet libraries were processed using Cell Ranger 5.0.0 (10X Genomics). Sequencing reads were aligned with STAR (v2.7.2a) 50 using the GRCh38 human reference genome. Filtered expression matrices generated using Cell Ranger count were used to perform the analysis. We excluded cells that had fewer than the first percentile of UMIs or greater than 99th percentile of UMIs in each sample. We also excluded cells with <200 genes expressed or >10% of reads from mitochondrial genes or hemoglobin genes. Additionally, putative doublets were removed using Scrublet (v0.2.1) 51 and scds (v1.10.0) 52 for each sample.
Article https://doi.org/10.1038/s41588-023-01375-1 types with high IFN response, we calculated the average module scores across each of 13 cell types and compared them with that of all PMBCs. The module scores of cells in each cell type were also compared with the average score of all PBMCs using two-sided one-sample t-test. To evaluate IFN response of COVID-19 patients in each cell type, we calculated the average module scores across each of 13 cell types by three clinical statuses and compared the scores of moderate or severe group with those of healthy group. The module scores of cells of moderate or severe group were also compared to those of healthy group in each cell type using two-sided Welch's t-test.

Differential gene expression analysis and GO enrichment
Differential gene expression analysis was performed among (1) all COVID-19 patients versus healthy controls, (2) moderate patients versus healthy controls, (3) severe patients versus healthy controls and (4) moderate versus severe patients. Donor pseudobulk matrices were first created by aggregating gene counts for each cell type, within each sample. Genes were considered for the analysis if they were expressed (UMI count > 0) in more than 10% of cells per cell type. Samples with more than five cells in a cell type were considered in the analysis of the corresponding cell type. Differential gene expression testing was performed using an NB GLM implemented in the Bioconductor package edgeR (v3.32.0) 55 . We included age, sex, days since symptom onset and duration of systemic steroids treatment at the time of specimen collection (the last two variables included only in moderate COVID-19 versus severe COVID-19) in the model as covariates. Statistically significant DEGs were defined with FDR < 0.05 and FC > 2. To find the function of upregulated DEGs, we used the function compareCluster (fun = "enrichGO," pvalueCutoff = 0.05, pAdjustMethod = "BH," OrgDb = "org.Hs.eg.db", ont = "BP") of Clusterprofiler (v3.14.3) 56 .

Estimation of RNA velocity
Spliced and unspliced transcripts were quantified using dropEst (v0.8.6) 57 . Monocytes and DC clusters were evaluated by RNA velocity analysis using scVelo (v0.2.3) 34 , separately for each of the two conditions. dropEst-derived counts were processed, filtered and normalized before velocity estimation on the basis of the top 2,000 HVGs with at least 20 UMI for both spliced and unspliced transcripts across all cells. The moments facilitated the RNA velocity estimation implemented in function scv.tl.velocity with mode set to 'dynamical'. The estimated velocities were used to construct a velocity graph representing the transition probabilities among cells by function scv.tl.velocity_graph. Finally, the velocity graph was used to embed the RNA velocities into the UMAP by the function scv.pl.velocity_embedding_stream. The fractions of unspliced counts were adopted to quantitatively compare the differences in estimated cell transition in monocytes between COVID-19 and healthy groups. The average unspliced ratio of each sample for each of the three monocyte clusters was calculated. The difference of increase in unspliced ratio across three monocyte clusters (ordered from cMono < intMono < ncMono) was evaluated using a linear regression model in each of the two conditions. The interaction between three monocyte clusters and two conditions (ordered from healthy < COVID-19) was also evaluated using a linear regression model with cluster, condition and cluster × condition as covariates.

TCR and BCR analysis
Droplet-based sequencing data for TCR sequences and BCR sequences were aligned and quantified using 5.0.0 (10X Genomics) against the GRCh38 human VDJ reference genome. Filtered annotated contigs for TCR sequences and BCR sequences were analyzed using Scirpy (v0.10.0) 58 . For TCR analysis, we selected T cells that were annotated as following nine cell types via single-cell RNA-seq analysis: CD4_naive, CD4_memory, CD4_Ef, T reg , CD8_naive, CD8_memory, CD8_Ef, MAIT and Pro_T (Extended Data Fig. 3a). Only cells with both TCR α-chain (TRA) and TCR β-chain (TRB) remained for the downstream analysis.
Each unique TRA(s)-TRB(s) pair was defined as a clonotype. Similarly, for BCR analysis, we selected B cells which were annotated as following four cell types via scRNA-seq analysis: B_naive, B_memory, B_activated and B_plasma (Extended Data Fig. 4a). Only cells with both heavy chain (IGH) and light chain (IGK or IGL) were kept for further analysis. Each unique IGH(s)-IGK/IGL(s) pair was defined as a clonotype.
For TCR and BCR data, clonotypes were defined based on CDR3 amino acid sequences with receptor_arms = 'any', metric = 'alignment' and default cutoff of ten. If one clonotype was present in at least two cells, cells harboring this clonotype were considered to be clonal and the number of cells with such pairs indicated the degree of clonality of the clonotype. Using barcode information, T cells with prevalent TCR clonotypes and B cells with prevalent BCR clonotypes were projected on UMAP embedding. To evaluate differences of clonal state between disease status in each cluster, the average clonal expansion rate of each sample was evaluated using two-sided Welch's t-test, where only cells mapped with TCRs/BCRs were included in the analysis.
We downloaded VDJdb 37 , a curated database of TCR sequences with known antigen specificities, and then investigated the TCR that was specific to SARS-CoV-2, based on CDR3 amino acid sequences with the same parameter as above.

Cell-cell interaction analysis in PBMC
At first, to reduce the influence of individual samples contributing a larger number of cells and to speed up computation, we capped the number of cells per sample at 2,500 randomly sampled cells. This was done using the SubsetData() function in Seurat. Genes were adopted if they were expressed in more than 1% of all PBMCs. Putative cell-cell interactions of COVID-19 and healthy were quantified using CellPhoneDB (v2.0.0) and NATMI with default settings 39,40 . To investigate the differences in cell-cell interactions between COVID-19 versus healthy and between moderate versus severe COVID-19, we evaluated FC of mean expression weight using DiffEdges.py implemented in NATMI with default settings. As CXCL10/CXCR3 interaction is not listed on connectomeDB2020, we manually added CXCL10 as ligand and CXCR3 as receptor. Then, cell-cell interactions of CXCL10/CXCR3 and IFNG/IFNGR were evaluated for each of the three conditions by mean expression weight using NATMI.

Polygenic GWAS signals on PBMC
We used MAGMA (v1.07) 59 to compute gene-level association P values and z-score from GWAS summary statistics of COVID-HGI (round 6) 18 . We used a reference panel based on individuals of European ancestry in the 1000 Genomes Project and used a 10-kb window around the gene body to map SNPs to genes. We selected the top 100 genes based on MAGMA P values as putative disease genes.
We used scDRS (v1.0.1) 41 to quantify the aggregate expression of putative disease genes derived from GWAS summary statistics using MAGMA (each putative disease gene is weighted by its GWAS MAGMA z-score and inversely weighted by its gene-specific technical noise level in the single-cell data) in each cell of COVID-19 scRNA-seq data to generate cell-specific raw disease scores. A 1,000 sets of cell-specific raw control scores were calculated from matched control gene sets (matching the gene set size, mean expression and expression variance of the putative disease genes). Then, we normalized the raw disease score and raw control scores for each cell, producing the normalized disease score and normalized control scores. To compute the scores described above, we used the function scdrs compute-score (--n_ctrl = 1000, -cov-file = age, sex, number of genes per cell and disease severity).
For downstream analysis, we performed cell type-level analyses to associate predefined cell types to disease and assess heterogeneity in association to disease across cells within a predefined cell type using the function scdrs perform-downstream with default settings. To correct multiple testing, FDR was calculated via the Benjamini-Hochberg method across all pairs of cell types and three GWAS phenotypes. To compare scDRS disease scores between COVID-19 phenotypes in six Article https://doi.org/10.1038/s41588-023-01375-1 major cell types, the differences in average disease scores of each sample were evaluated using two-sided paired t-test (adjusted for multiple comparisons using Bonferroni's correction).

Genotype data, quality control and genotype imputation
We performed GWAS genotyping of COVID-19 cases and healthy controls using Infinium Asian Screening Array (Illumina) through collaboration with Japan COVID-19 Task Force (https://www.covid19-taskforce. jp/en/home/). We applied stringent quality control filters to the samples (sample call rate < 0.98, related samples with PI_HAT > 0.175 or outlier samples from East Asian clusters in PCA with HapMap project samples), and variants (variant call rate < 0.99, deviation from Hardy-Weinberg equilibrium with P < 1.0 × 10 −6 , or minor allele count < 5). We used SHAPEIT4 software (version 4.2.1) 60 for haplotype phasing of autosomal genotype data. After phasing, we used Minimac4 software (version 1.0.1) 61 for genome-wide genotype imputation. We used the population-specific imputation reference panel of Japanese (n = 1,037) combined with 1,000 Genomes Project Phase3v5 samples (n = 2,504) 62 .

Single-cell eQTL analysis
We applied pseudobulk approach for single-cell eQTL analysis. First, we performed single-cell-level normalization using scran (v1.18.5) 63 . Gene expression per cell type per sample was then calculated as the mean of log2-transformed counts-per-cell-normalized expression across cells. Samples with more than five cells in a cell type were considered in the analysis of the corresponding cell type. We examined eQTL effects of the replicated variants in Japanese COVID-19 GWAS on V2G, the highest gene prioritized by the V2G score of Open Target Genetic, separately for COVID-19 and healthy controls 20,64 . rs35081325 and rs77534576 were excluded from eQTL analysis due to low allele frequency and low expression of V2G (FLJ45513), respectively. For PCA, genes were adopted if they were expressed in more than 1% of all PBMCs.
In the eQTL analysis of COVID-19-associated variants, dosage effects of the variants on the gene expression mean were evaluated using linear regression models with the top two PCs of the genotype data, the top two PCs of the gene expression, age, sex, days since symptom onset, duration of systemic steroids treatment at the time of specimen collection and disease severity (the last three variables included only in COVID-19 analysis) as covariates. In the interaction eQTL analysis of the IFNAR2 variant (rs13050728), the top two PCs of the genotype data, the top two PCs of the gene expression, age and sex were included as covariates. R statistical software (version 4.0.2) was used for the analysis. To correct multiple testing, FDR was calculated via the Benjamini-Hochberg method across all pairs of cell types and five variants, separately for COVID-19 and healthy controls.

Statistics and reproducibility
No statistical method was used to predetermine sample size. No data were excluded from the analyses. We did not use any study design that required randomization or blinding.

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