Main

The host response to severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) infection ranges from asymptomatic infection to severe coronavirus disease 2019 (COVID-19) and death. The leading cause of mortality is acute lung injury and acute respiratory distress syndrome, or direct complications with multiple organ failure1,2,3,4. Clinical deterioration in acute illness leads to ineffective viral clearance and collateral tissue damage1,2,3,4,5. Severe COVID-19 is also accompanied by an inappropriate pro-inflammatory host immune response and a diminished antiviral interferon response6,7,8.

Many molecular and cellular questions related to COVID-19 pathophysiology remain unanswered, including how cell composition and gene programs shift, which cells are infected, and how associated genetic loci drive disease. Autopsies are crucial to understanding severe COVID-19 pathophysiology9,10,11,12, but comprehensive genomic studies are challenged by long post-mortem intervals (PMIs).

Here, we developed a large cross-body COVID-19 autopsy biobank of 420 autopsy specimens, spanning 11 organs, and used it to generate a single-cell atlas of lung, kidney, liver and heart associated with COVID-19 and a lung spatial atlas, in a subset of 14–18 donors per organ. Our atlases provide crucial insights into the pathogenesis of severe COVID-19.

A COVID-19 autopsy cohort and biobank

We assembled an autopsy cohort of 20 male and 12 female donors, of various ages (>30–>89 years), racial/ethnic backgrounds, intermittent mandatory ventilation (IMV) periods (0–30 days) and days from symptom start to death (Fig. 1a, Supplementary Table 1). A biobank was created with a subset of 17 donors. From most donors, we collected at least lung, heart and liver tissue (Fig. 1a, Extended Data Fig. 1a, Supplementary Methods), preserving specimens for single-cell and spatial analysis. We optimized single-cell and single-nucleus RNA sequencing (sc/snRNA-Seq) protocols for Biosafety Level 3 and NanoString GeoMx workflows to spatially profile RNA from different tissue compartments by cell composition or viral RNA (Supplementary Methods).

Fig. 1: Experimental and computational pipeline for a COVID-19 autopsy atlas.
figure 1

a, Sample processing pipeline. Up to 11 tissue types were collected from 32 donors. b, sc/snRNA-Seq analysis pipeline. QC, quality control.

COVID-19 cell atlases

We generated sc/snRNA-Seq atlases of lung (n = 16 donors, k = 106,792 cells/nuclei, m = 24 specimens; donors D1–8, 10–17), heart (n = 18, k = 40,880, m = 19; D1–8, 10–11, 14–17, 27–28, 31–32), liver (n = 15, k = 47,001, m = 16; D1–7, 10–17) and kidney (n = 16, k = 33,872, m = 16; D4–8, 10–12, 14–15, 17, 25–26, 28–30). Although initial tests showed some differences in proportions of cell types between snRNA-Seq and scRNA-Seq, snRNA-Seq performed better overall13 (Extended Data Fig. 1b–d and data not shown) and was used for the remaining samples.

We developed a computational pipeline (Fig. 1b) to tackle unique technical challenges. We used CellBender remove-background14 to remove ambient RNA, which enhanced cell distinction and marker specificity (Extended Data Fig. 1e–h; Supplementary Methods), we rapidly quality-controlled, pre-processed and batch-corrected the data with cloud-based Cumulus15 (Extended Data Fig. 2a–g, Supplementary Methods) and we automatically annotated cells and nuclei by transferring labels from previous atlases (Fig. 2a, Extended Data Fig. 2h, Supplementary Methods). We refined these labels with manual annotation of subclusters in each main lineage (Fig. 2b, Extended Data Fig. 2i–n, Supplementary Methods). The automated annotation approach allowed us to compare against other data resources (without clustering or batch correction), while the manual approach enabled us to refine cell identity assignments with detailed domain knowledge.

Fig. 2: A single-cell and single-nucleus atlas of COVID-19 lung.
figure 2

a, Automatic prediction identifies 28 cell subsets across compartments. UMAP embedding of 106,792 harmonized sc/snRNA-Seq profiles (dots) from 24 tissue samples of 16 lung donors with COVID-19, coloured by automatic annotations (legend). b, Epithelial cell subsets. UMAP embedding of 21,661 epithelial cell or nucleus profiles, coloured by manual annotations, with highly expressed marker genes (boxes). c, d, Cell composition and expression differences between COVID-19 and healthy lung. c, Cell proportions (x axis: mean, bar; 95% confidence intervals, line) in each automatically annotated subset (y axis) in COVID-19 snRNA-Seq (red, n = 16), healthy snRNA-Seq (grey, n = 3) and healthy scRNA-Seq (n = 8, blue). Cell types shown have a COVID-19 versus healthy snRNA-Seq false discovery rate (FDR) of <0.05 (Dirichlet multinomial regression). d, Significance (−log10(P), y axis) versus magnitude (log2(fold change), x axis) of differential expression of each gene (dots; horizontal dashed line, FDR < 0.05) between COVID-19 and healthy lung from a total of 2,000 AT2 cells and 14 studies (two-sided test; Supplementary Methods). e, f, An increased pre-alveolar type 1 transitional cell state (PATS)20,21,22 program in pneumocytes in COVID-19 versus healthy lung. e, Distribution of PATS signature scores (y axis) for 17,655 cells from COVID-19 and 24,000 cells from healthy lung pneumocytes (x axis). P < 2.2 × 10−16 (one-sided Mann–Whitney U test). f, UMAP embedding of 21,661 epithelial cell profiles (dots) coloured by signature level (colour legend, lower right) for the PATS (top) or intrapulmonary basal-like progenitor (IPBLP) cell (bottom) programs. g, Model of epithelial cell regeneration in healthy and COVID-19 lung. In healthy alveoli (top), AT2 cells self-renew (1) and differentiate into AT1 cells (2). In COVID-19 alveoli (bottom), AT2 cell self-renewal (1) and AT1 differentiation (2) are inhibited, resulting in PATS accumulation (3) and recruitment of airway-derived IPBLP cells to alveoli (4).

A cell census of the COVID-19 lung

Automatic annotation defined 28 subsets of parenchymal, endothelial and immune cells (Fig. 2a, Supplementary Table 2, Supplementary Methods), with further manual annotation within subgroupings (Fig. 2b, Extended Data Figs. 2, 4, Supplementary Methods). Deconvolution of bulk RNA-Seq from the same samples largely agreed (Extended Data Fig. 3a, b, Supplementary Methods), and our two annotation strategies had 94% agreement (Extended Data Fig. 3c–e).

Among immune cells we distinguished six cell myeloid subsets: CD14highCD16high inflammatory monocytes with antimicrobial properties and five macrophage subsets (Extended Data Figs. 2j, 4b) that were enriched for scavenger receptors, toll-like receptor ligands, inflammatory transcriptional regulators or metabolism genes; four B cell and plasma cell subsets: BLIMP1high plasma cells16,17, BLIMP1intermediate plasma cells, B cells and JCHAIN-expressing plasmablasts (Extended Data Figs. 2k, 4b); five T and natural killer (NK) cell subsets: two CD4+ subsets: regulatory T cells (Treg cells) and a metabolically active subset; one CD8+ subset; and two T or NK cell subsets (Extended Data Figs. 2l, 4b), including one with cytotoxic effector genes. The dearth of neutrophils (Fig. 2a, 419 cells) is likely due to freezing or limitations of droplet-based sc/snRNA-Seq13.

We identified seven endothelial cell (EC) subsets18,19 (Extended Data Figs. 2m, 4b): arterial, venous and lymphatic, capillary aerocytes, capillary EC-1 and capillary EC-2 and a mixed subset (Supplementary Methods), and three stromal subsets: fibroblasts, proliferative fibroblasts and myofibroblasts19 (Extended Data Fig. 2n, Supplementary Table 3).

There were eight epithelial subsets, including club/secretory cells, AT1 cells, AT2 cells, and proliferative AT2 cells (Fig. 2b). One subset corresponded to a previously described AT2 to AT1 transitional cell state (KRT8+ pre-alveolar type 1 transitional cell state (PATS); PATS/ADI/DATP)20,21,22 (Fig. 2b).

Changes in lung cell composition

In comparison with normal lung from a matching region (Fig. 2c, Supplementary Methods), numbers of AT2 cells were significantly decreased (false discovery rate (FDR) = 2.8 × 10−15, Dirichlet multinomial regression; Supplementary Methods), possibly reflecting virally induced cell death23,24,25. Numbers of dendritic cells (FDR = 0.004), macrophages (FDR = 3.6 × 10−10), NK cells (FDR = 0.018), fibroblasts (FDR = 0.013), lymphatic endothelial cells (FDR = 0.00058) and vascular endothelial cells (FDR = 0.00011) all increased.

Cell proportions varied between donors (Extended Data Fig. 5a, b). Whereas variation was not significantly correlated with PMI, age or sex, IMV was positively correlated with epithelial cell fraction (FDR = 0.007; Spearman ρ = 0.765) and negatively correlated with T and NK cell fraction (FDR = 0.041; ρ = −0.62). Fewer days on a ventilator may indicate a rapidly deteriorating condition. This is corroborated by the nominally significant positive correlation between epithelial cell fraction and days from symptom start to death (ρ = 0.671, P = 0.004, but FDR = 0.053).

Induced programs in epithelial cells

There were widespread, cell-type-specific transcriptional changes in lung cell types associated with COVID-19 (Extended Data Fig. 5c, Supplementary Methods), most notably in CD16+ monocytes (1,580 upregulated genes), lymphatic endothelial (578), vascular endothelial (317), AT2 (309) and AT1 (307) cells. Within AT2 cells, there was higher expression (P < 0.0004) of genes associated with host viral response (Fig. 2d), including those for programmed cell death (STAT1), inflammation and adaptive immune response (Supplementary Table 4). Lung surfactant genes were downregulated, consistent with reports from in vitro studies21.

Failed paths for AT1 cell regeneration

The PATS program signature was increased in COVID-19 pneumocytes (P < 2.2 × 10−16, one-sided Mann–Whitney U test) (Fig. 2e, Extended Data Fig. 5d). This progenitor program is induced during lung injury20,21,22 (for example, idiopathic pulmonary fibrosis), consistent with fibrosis in severe COVID-1926,27. These studies also highlight fibroblast expansion, which we also observed (Fig. 2c).

A subset of PATS program cells, distinct from KRT5+TP63+ airway basal cells, expressed canonical (KRT8/CLDN4/CDKN1A) and non-canonical (KRT5/TP63/KRT17) PATS markers (Fig. 2f, Extended Data Fig. 5d, Supplementary Table 3). These may be TP63+ intrapulmonary basal-like progenitor (IPBLP) cells, which were identified in H1N1 influenza mouse models28 and act as an emergency cellular reserve for severely damaged alveoli29. The putative IPBLP cells express interferon virus defence and progenitor cell differentiation genes (Supplementary Table 3). Thus, multiple emergency pathways for alveolar cell regeneration are activated in lung (Fig. 2g, Discussion).

Changed cell composition with viral load

To determine viral load and associated host responses, we analysed donor- and cell-type-specific distribution of SARS-CoV-2 reads (Fig. 3a, b, Extended Data Fig. 6a–d, Supplementary Methods). Reads spanned the entire SARS-CoV-2 genome, with bias towards positive-sense alignments. A few cells had reads aligning to all viral segments, including the negative strand (Extended Data Fig. 6e), potentially indicating productive infection. Virus detection was not technically driven (Extended Data Fig. 6f–i), and inter-donor variation was consistent with SARS-CoV-2 qRT–PCR on bulk RNA (Extended Data Fig. 6j–l, Supplementary Methods). Viral load was negatively correlated with days from symptom start to death (Fig. 3c), as previously reported30,31. Bulk RNA-Seq yielded nine unique complete viral genomes from nine donors with high viral loads (Extended Data Fig. 6m, Supplementary Methods); all genomes carried the D614G allele. We identified no other common respiratory viral co-infections (Extended Data Fig. 6n). Total viral burden per sample (including ambient RNA; Supplementary Methods) positively correlated with proportions of mast cells, specific macrophage subsets, venular endothelial cells and capillary aerocyte endothelial cells (Extended Data Fig. 6o–u).

Fig. 3: SARS-CoV-2 RNA+ single cells are enriched for phagocytic and endothelial cells.
figure 3

a, b, Many SARS-CoV-2 RNA+ single cells do not express known SARS-CoV-2 entry factors. UMAP embedding of all 106,792 lung cells or nuclei (as in Fig. 2a), showing either a, only the 40,581 cells from seven donors containing any SARS-CoV-2 RNA+ cell, coloured by viral enrichment score (Supplementary Methods, red: stronger enrichment) and by SARS-CoV-2 RNA+ cells (black points), and marked by annotation and FDR of enrichment (legend) or b, all 106,792 cells/nuclei, coloured by expression of SARS-CoV-2 entry factors (co-expression combinations with at least 10 cells are shown). Dashed lines, major cell types. c, Reduction in SARS-CoV-2 RNA with prolonged symptom onset to death interval (Spearman ρ = −0.68, P < 0.005, two-sided test). Symptom onset to death (x axis, days) and lung SARS-CoV-2 copies per nanogram input RNA (y axis) for each donor (n = 16). d, Expression changes in SARS-CoV-2 RNA+ myeloid cells. Significantly differentially expressed (DE) host genes (log-normalized and scaled digital gene expression, rows; cutoff: FDR < 0.05 and log2(fold change) > 0.5) across SARS-CoV-2 RNA+ (n = 158) and SARS-CoV-2 RNA myeloid cells (n = 790) (columns).

Genes upregulated in biopsy samples with high versus low or no viral load (Supplementary Methods) included viral response and innate immune processes (log2(fold change) > 1.4, Wald test, FDR-corrected P < 0.05; Extended Data Fig. 6v, Supplementary Table 4) and significantly overlapped with those in bulk RNA-Seq of post-mortem COVID-19 lungs in another study32 (FDR = 3.12 × 10−6, Kolmogorov–Smirnov test). Downregulated genes (log2(fold change) < 1.4, Wald test, FDR-corrected P < 0.05) were involved in surfactant metabolism dysfunction and lamellar bodies (secretory vesicles in AT2 cells33).

Lung cells enriched for SARS-CoV-2 RNA

Myeloid cells were the cell category most enriched for SARS-CoV-2 RNA (158 cells after correction for ambient RNA, FDR < 0.013; Fig. 3a, Extended Data Fig. 6w–y, Supplementary Methods), with particular enrichment in CD14highCD16high inflammatory monocytes (FDR < 0.005) and LDB2highOSMRhighYAP1high macrophages (FDR < 0.02; Extended Data Figs. 6x, 7a, b), although enrichment scores in individual donors varied. There was elevated, but non-significantly enriched, viral RNA in endothelial cells, with the capillary EC-2 (cluster 3, FDR < 0.017) and lymphatic endothelial cells (cluster 7, FDR < 0.006) enriched compared with other endothelial subsets (Fig. 3a, Extended Data Figs. 6w, y, 7c, d). There were also SARS-CoV-2 RNA+ cells among mast cells, and B and plasma cells, and viral RNA reads in multiple other cell types (Fig. 3a, Extended Data Fig. 6w). Notably, SARS-CoV-2 RNA+ cells did not co-express the entry factors ACE2 and TMPRSS2, or other hypothesized entry cofactors (Fig. 3b, Extended Data Fig. 7e–h).

Immune programs in SARS-CoV-2 RNA+ cells

SARS-CoV-2 RNA+ cells had distinct transcriptional programs compared with SARS-CoV-2 RNA counterparts, with differentially expressed genes (FDR < 0.05; Supplementary Methods) in epithelial and myeloid cells, including PPARGhighCD15Lhigh macrophages and CD14highCD16high inflammatory monocytes (Supplementary Table 5). Genes upregulated in epithelial SARS-CoV-2 RNA+ cells were enriched for TNF, AP1 and chemokine and cytokine signalling, SARS-CoV-2-driven cell responses in vitro32, and keratinization pathways, which may reflect an injury response (Extended Data Fig. 7i). Genes upregulated in myeloid SARS-CoV-2 RNA+ cells were those associated with chemokine and cytokine signalling, and responses to interferon, TNF, intracellular pathogens and viruses (Fig. 3d, Extended Data Fig. 7j–m, Supplementary Table 5), as previously described34,35. Cytokines and viral host response genes were upregulated in both CD14highCD16high inflammatory monocytes and PPARGhighCD15Lhigh macrophages (Extended Data Fig. 7m, Supplementary Table 5), including CXCL10 and CXCL11, which were upregulated in nasopharyngeal swabs35 and bronchoalveolar lavages34.

A spatial atlas of COVID-19 lung

To provide tissue context, we used Nanostring GeoMx Digital Spatial Profiling (DSP) for transcriptomic profiling from regions of interest (Supplementary Methods) in 14 donors, including three deceased healthy donors (‘healthy’) (Extended Data Fig. 1a). Regions of interest spanned a range of anatomical structures and viral abundance on the basis of SARS-CoV-2 RNA hybridization signals; when possible, we segmented them to PanCK+ and PanCK, and inflamed and normal-appearing alveoli areas of illumination (AOIs) to capture RNA (Fig. 4a, Extended Data Figs. 8a, 9a, Supplementary Methods). We acquired high-quality profiles (Extended Data Fig. 8b) from matched AOIs on the basis of distance to morphological landmarks (Supplementary Methods). SARS-CoV-2 RNA expression varied by donor, with elevated levels in four donors (Extended Data Fig. 8c, d, Supplementary Methods), consistent with viral qRT–PCR and sc/snRNA-Seq. Given the good agreement between a targeted 1,811-gene panel and a whole transcriptome (WTA) panel (18,335 genes) (Extended Data Fig. 8e–g, Supplementary Table 6), we focused our analyses on WTA data. For D8–12, 18–24, we contrasted donors with COVID-19 and healthy donors and COVID-19 epithelial and non-epithelial AOIs; for D13–17, we focused on distinct anatomical regions and inflamed versus normal-appearing regions within donors.

Fig. 4: Composition and expression differences between COVID-19 and healthy lungs and between infected and uninfected regions within COVID-19 lungs.
figure 4

a, Example of analysed regions. Top: RNAscope (left) and immunofluorescent staining (right) of donor D20 with collection regions of interest (ROIs) and matched areas in white rectangles. Bottom: one ROI (yellow rectangle) from each scan (left and middle) and the segmented collection areas of illumination (AOIs) (right). b, Cell composition differences between PanCK+ and PanCK alveolar AOIs and between AOIs from COVID-19 (n = 9, 190 AOIs) and healthy (D22–24, 38 AOIs) lungs. Expression scores (colour bar) of cell type signatures (rows) in PanCK+ (left) and PanCK (right) alveolar AOIs (columns) in whole transcriptome (WTA) data from different donors (top colour bar). c, Differential gene expression in COVID-19 versus healthy lung. Left: significance (−log10(P), y axis) and magnitude (log2(fold change), x axis) of differential expression of each gene (dots) in WTA data between PanCK+ alveoli AOIs from COVID-19 (n = 78) and healthy (n = 18) lung. Right: significance (−log10(q)) of enrichment (permutation test) of different pathways (rows). d, e, Changes in gene expression in SARS-CoV-2 high versus low AOIs within COVID-19 lungs in WTA data. d, SARS-CoV-2 high and low alveolar AOIs. PanCK+ alveolar AOIs (dots) rank ordered by their SARS-CoV-2 signature score (y axis) in WTA data, and partitioned to high (red), medium (grey) and low (blue) SARS-CoV-2 AOIs. e, Significance (−log10(P), y axis) and magnitude (log2(fold change), x axis) of differential expression of each gene (dots) in WTA data between SARS-CoV-2 high and low AOIs for PanCK+ alveoli (AOIs: 17 high, 3 medium, 58 low). Horizontal dashed line, FDR = 0.05; vertical dashed lines, |log2(fold change)| = 2. FC, fold change. The top 10 differentially expressed genes by fold change are marked.

Inflammatory activation in alveoli

Deconvolution of major cell type composition (Fig. 4b, Extended Data Fig. 8h, Supplementary Table 7, 8, Supplementary Methods) showed inferred AT1 and AT2 cells dominating the PanCK+ compartments and greater cellular diversity in the PanCK compartment. COVID-19 PanCK AOIs had increased fibroblast and myofibroblast scores compared with controls, in line with parallel spatial studies36,37.

Comparing COVID-19 alveolar AOIs with control lungs from deceased healthy donors, there was upregulation of IFNα and IFNγ response genes and oxidative phosphorylation pathways (Fig. 4c, Extended Data Fig. 8i–k, Supplementary Table 6), similar to bulk RNA-Seq of highly infected tissue (IFIT1, IFIT3, IDO1, GZMB, LAG3, NKG7 and PRF1) and to SARS-CoV-2+ myeloid cells (TNFAIP6, CXCL11, CCL8, ISG1 and GBP5) and consistent with PANoptosis in a COVID-19 model38. Conversely, TNF, IL2–STAT5 and TGFβ signalling as well as apical junction and hypoxia were downregulated. Decreased TNF signalling expression in PanCK+ alveoli contrasts with its increase in SARS-CoV-2+ epithelial cells in snRNA-Seq and with reported38 synergy between TNF and IFNγ in mouse models of COVID-19.

Comparison of inflamed and normal-appearing AOIs within the same alveolar biopsy samples of COVID-19 lungs (Extended Data Fig. 9, Supplementary Table 9, D13–D17), showed that upregulated genes were enriched for innate immune and inflammatory pathways39,40, including neutrophil degranulation (FDR = 5.2 × 10−17) and IFNγ (FDR = 3.4 × 10−15) and interleukin (FDR = 1.4 × 10−13) signalling. TNF pathway expression was elevated in inflamed tissue, albeit not significantly (FDR = 0.097). Claudins and tight junction pathways were downregulated, corroborating a disrupted alveolar barrier, as in influenza41,42. Cilium assembly genes were enriched when comparing bronchial epithelial AOIs and matched normal-appearing alveoli (Extended Data Fig. 9d, Supplementary Table 9).

Comparison of SARS-CoV-2 high and low AOIs (Fig. 4d, e, Extended Data Fig. 8l, m, Supplementary Methods) revealed induction of the viral ORF1ab and S genes and upregulation of chemokine genes (CXCL2 and CXCL3) and immediate early genes in the PanCK+ compartment, consistent with snRNA-Seq (Supplementary Table 9, Extended Data Fig. 7i). NT5C, which encodes a nucleotidase with a preference for 5′-dNTPs, is consistently upregulated in SARS-CoV-2-high AOIs (Fig. 4e, Extended Data Fig. 8m, Supplementary Table 9). This gene is not known to have a role in lung injury and should be further studied.

COVID-19 effect on heart, kidney and liver

We next profiled liver, heart and kidney by snRNA-Seq with automated and manual annotation of parenchymal, endothelial and immune cells (Supplementary Methods, Extended Data Figs. 10, 11). Although other studies have reported viral reads in COVID-19 non-lung tissues43, we detected very few viral RNA reads in all three tissues, most of which could not be assigned to nuclei (Extended Data Fig. 11l); this absence was confirmed by NanoString DSP and RNAscope (data not shown).

Focusing on heart, both cell composition and gene programs changed between COVID-19 and healthy heart. There was a significant reduction in the proportion of cardiomyocytes and pericytes, and an increase in vascular endothelial cells (Extended Data Fig. 11e). Genes upregulated (FDR < 0.01) in cardiomyocytes, pericytes or fibroblasts (Extended Data Fig. 11g–i, Supplementary Table 10) included PLCG2, the cardiac role of which is unknown but which was induced in all major heart cell subtypes (Extended Data Fig. 11j), and AFDN, which is upregulated in endothelial cells (Extended Data Fig. 11k), and which encodes a junction adherens complex component44 that is necessary for endothelial barrier function. Upregulated pathways include oxidative-stress-induced apoptosis in pericytes, cell adhesion and immune pathways in cardiomyocytes, and cell differentiation processes in fibroblasts (Supplementary Table 10).

COVID-19 cell types related through GWAS

Finally, we aimed to identify genes and cell types associated with COVID-19 risk by integrating our atlas data with genome-wide association studies (GWAS)45 for common46 variants associated with COVID-19 (Supplementary Methods). Among 26 genes proximal to six COVID-19 GWAS regions (Supplementary Table 11, Supplementary Methods), 14 genes had higher average expression in the lung (P < 0.05, t-test; Extended Data Fig. 12a–d), 21 had significant (FDR < 0.05) expression specificity in at least one lung cell type, including FOXP4 (chromosome (chr.) 6, AT1 and AT2 cells), and CCR1 and CCRL2 (chr. 3, macrophages) (Extended Data Fig. 12e, Supplementary Table 11), and 18 were differentially expressed (FDR < 0.05) in COVID-19 compared with healthy lung (for example, SLC6A20 in goblet cells, CCR5 in CD8+ T cells and Treg cells, and CCR1 in macrophage and CD16+ monocytes (Extended Data Fig. 12f, Supplementary Table 11).

We related heritability from GWAS of COVID-19 severity traits to either cell type programs (genes enriched in a cell type in each tissue) or disease progression programs (genes differentially expressed between COVID-19 and controls in a cell type) in each tissue using sc-linker47 (Supplementary Methods). AT2 (4.8× heritability enrichment, P = 0.04), CD8+ T cells (4.4×, P = 0.009) and ciliated cell programs in the lung, proximal convoluted tubule and connecting tubule programs in kidney, and cholangiocyte programs in liver attained nominal (but not Bonferroni-corrected) significance (Extended Data Fig. 12g, h, Supplementary Table 11). Of all disease progression programs, only the club cell program (single-cell level model) had nominally significant heritability enrichment (10.5×, P = 0.04 for severe COVID-19) (Extended Fig. 12g, Supplementary Table 11).

The highest number of driving genes was observed for lung AT2 cells and spanned several loci, hinting at a polygenic architecture linking AT2 cells with severe COVID-19 (Supplementary Methods, Supplementary Table 11). Implicated GWAS proximity genes include OAS3 in lung AT2 and club cells, and SLC4A7 in lung CD8+ T cells (Supplementary Table 11), as well as genes at unresolved significantly associated GWAS loci (Extended Data Fig. 12i), such as FYCO1 (AT2, ciliated, club; chr. 3p), NFKBIZ (AT2; chr. 3q) and DPP9 (AT2; chr. 19) (Supplementary Table 11).

Discussion

We built a biobank of severe COVID-19 autopsy tissue and atlases of COVID-19 lung, heart, liver and kidney (Extended Data Fig. 12j), complementing a sister lung atlas48.

Among the changes in lung cell composition in COVID-19, is a reduction in AT2 cells and the presence of PATS and IPBLP-like cells, indicating that multiple regenerative strategies are invoked to re-establish alveolar epithelial cells lost to infection. A serial failure of epithelial progenitors to regenerate at a sufficient rate, first by secretory progenitor cells in the nasal passages and large and small airways, followed by alveolar AT2 cells, PATS and IPBLP cells, may eventually lead to lung failure.

Viral RNA in the lung varied, was negatively correlated with time from symptom start to death, and was primarily detected in myeloid and endothelial cells (as in nonhuman primates49); spatial analysis supports high virus levels at the earlier stages of infection36,37,50. Epithelial cells were not enriched in high viral RNA samples or in SARS-CoV-2+ cells, consistent with their excessive death. Cell-associated SARS-CoV-2 unique molecular identifiers may represent a mixture of replicating virus, immune cell engulfment and virions or virally infected cells attached to the cell surface. We did not detect viral RNA in the heart, liver or kidney, but observed other changes, including broad upregulation of PLCG2, a target of Bruton’s tyrosine kinase (BTK), in the heart51.

Combining our profiles with GWAS of COVID-19, we related specific cell types to heritable risk, especially AT2, ciliated and CD8+ T cells and macrophages, as well as genes in multi-gene regions underlying the association. This analysis can improve as GWAS grows and atlases expand.

Our study was limited by a modest number of donors without pre-selection of features, the terminal time point, limited distinction between viral RNA and true infection, and technical confounders such as PMIs. Nevertheless, our methods would enable studies in diverse diseased or damaged tissues. Future meta-analyses will further enhance its power and provide crucial resources for the community studying host–SARS-CoV-2 biology.

Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this paper.