System-wide transcriptome damage and tissue identity loss in COVID-19 patients.

The molecular mechanisms underlying the clinical manifestations of coronavirus disease 2019 (COVID-19), and what distinguishes them from common seasonal influenza virus and other lung injury states such as acute respiratory distress syndrome, remain poorly understood. To address these challenges, we combine transcriptional profiling of 646 clinical nasopharyngeal swabs and 39 patient autopsy tissues to define body-wide transcriptome changes in response to COVID-19. We then match these data with spatial protein and expression profiling across 357 tissue sections from 16 representative patient lung samples and identify tissue-compartment-specific damage wrought by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) infection, evident as a function of varying viral loads during the clinical course of infection and tissue-type-specific expression states. Overall, our findings reveal a systemic disruption of canonical cellular and transcriptional pathways across all tissues, which can inform subsequent studies to combat the mortality of COVID-19 and to better understand the molecular dynamics of lethal SARS-CoV-2 and other respiratory infections.

Introduction

In March 2020, the World Health Organization (WHO) declared a novel pandemic of coronavirus disease 2019 (COVID-19), an infection caused by the betacoronavirus severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), which is currently attributed to over 320 million cases and over 5.5 million deaths globally (https://coronavirus.jhu.edu). Since the presenting symptoms of COVID-19 resemble those of common viral respiratory infections, a molecular diagnosis is required to distinguish a SARS-CoV-2 infection from influenza and other respiratory illnesses,, and ongoing questions remain about the host responses to SARS-CoV-2 relative to other respiratory pathogens. As severe illness and death continue to impact a segment of COVID-19-positive individuals, urgent questions persist about the molecular drivers of morbidity and mortality associated with SARS-CoV-2 infection. This knowledge could lead to improvements in both the acute treatment and the long-term management of pathological changes in multiple organs.

Prior work has shown that COVID-19 leads to high systemic interferon responses (both alpha and gamma) and that co-infection with other pathogens is relatively rare (3%–15%).4, 5, 6 Yet there are limited data for discriminating the molecular response between different kinds of respiratory infections or pulmonary conditions (e.g., influenza A (IAV) versus SARS-CoV-2 infections) and almost no data on the variegated impact of different pathogens across different human tissues other than clinical observations from intensive care unit (ICU) patients. Delineation of pathogen- and tissue-specific differences is critical for understanding the molecular determinants of mortality associated with COVID-19 and for developing novel diagnostics and therapeutic interventions.

To address this knowledge gap, we used shotgun metatranscriptomics (total RNA sequencing [RNA-seq]) to comprehensively profile human tissues in 39 patients who died from COVID-19 (185 total autopsy samples), including heart, liver, lung, kidney, and lymph nodes, analyzed gene expression, and assessed the system-wide impact of SARS-CoV-2 infection. We also used a spatial protein and transcript mapping platform (GeoMx) to visualize the cartography of the infection in these tissues and to discover disruption of regional and cell-type-specific expression. The spatial transcriptomics data examined 357 total regions of interest (ROIs), which were selected from 13 patients who had SARS-CoV-2, influenza, or bacterial infections and from 3 normal patients as controls, revealing the cellular and regulatory signatures that define these distinct pathological states. Finally, to provide context to earlier stages and sites of infection, we compared these in-depth spatial and tissue-specific transcriptome maps with an independent cohort of nasopharyngeal (NP) swabs from 216 COVID-19-positive patients and 430 COVID-19-negative controls, which revealed a significant and distinct disruption of cellular and transcriptional programs induced by SARS-CoV-2 infection in the patients who unfortunately succumbed to the disease. As a resource for the field, these data have also been placed in an online portal (https://covidgenes.weill.cornell.edu/) for additional data mining and visualization.

Results

System-wide host responses and transcriptome changes by COVID-19

We first used shotgun metatranscriptomics (total RNA-seq) for host and viral profiling on 39 patients who died from COVID-19, including 185 organ-specific tissue samples from the respective autopsies and 22 healthy control samples from organ donor remnant tissues (Figure 1A; Table S1). We examined the COVID-19-specific host responses and transcriptome changes across various organs (heart, kidney, liver, lung, and lymph nodes) to ascertain the differentially expressed genes (DEGs) between COVID-19 and control sample sets (q < 0.01, expression fold change > 1.5 by the differential expression analysis method [DESeq2]; Figure 1B; Table S2). Separately, SARS-CoV-2 RNA reads were aligned to the SARS-CoV-2 genome and the number of reads was discerned across multiple tissues including the lung, lymph node, kidney, liver, and heart (Figure S1) but mostly in the lung. Viral reads were robustly detected in NP swab samples from COVID-19 patients, consistent with previous reports., Normalized coverage values in SARS-CoV-2-positive autopsy tissue samples revealed a detection bias toward the SARS-CoV-2 3′ end sequences, which is consistent with the known viral transcript abundance (Figures S1A and S1B). Reconstruction of the viral genomes revealed known and unknown variants common to many patients and some evidence of intra-host variability (Figure S1C). Coverage limitations notwithstanding, we generally saw evidence of the same viral strain across the organ systems within each patient. We also assessed the variability of the COVID-19 patient samples by clustering based on the SARS-CoV-2 viral loads (and noted as SC2 high and SC2 low), and the viral loads were inversely correlated with the duration of disease and independent of factors such as race, age, or gender (Figures S1D and S1E). In addition to the annotated clinical metadata table (Table S1), we summarized the clinical courses of a few representative COVID-19 patients from hospitalization to intubation to death (Figure S1D).

Figure 1

Tissue- and duration-specific dysregulation of gene expression from SARS-CoV-2 infection

(A) Sample overview by tissue for RNA-seq and GeoMx experiments. Number of regions of interest (ROIs) and patient numbers (n) are summarized, and a representative tissue slide image for GeoMX spatial profiling is presented along with the tissue types.

(B) Volcano plots of the COVID-19-positive versus normal tissues are shown for the five tissues from autopsy: heart (n = 41), kidney (n = 27), liver (n = 40), lung (n = 40), and lymph node (n = 27). Differentially expressed genes (>1.5-fold, q < 0.05, DESeq2) are shown as purple spots (downregulated in COVID-19) and orange (upregulated in COVID-19).

COVID-19 pathway enrichment analysis revealed significant changes (q < 0.01) in pathways for viral infection (regulation of viral genome replication and viral entry into host cell) and immune response (regulation of type 1 interferon response and regulation of tyrosine phosphorylation of Stat protein, Gene Ontology (GO) Regulation of toll-like receptor signaling pathway), and the enrichment varied as a function of the viral load (Figure 2A; Table S2). We observed that each tissue showed its own distinct transcriptional disruption in response to SARS-CoV-2 infection, with the lymph node exhibiting the greatest number of DEGs when compared with controls (in both high- and low-viral-load groups). Of note, both tissue-specific and pan-tissue disruptions of normal-expression programs were observed (Figure 1B), and these were then summarized using gene set enrichment analysis (GSEA; Figures 2A and S2A). Some pathways were consistently dysregulated in all tissues during early infection (SARS-CoV-2 high), such as the G2M checkpoint (q values of 6.4 × 10−19, 3.0 × 10−5, 1.2 × 10−8, 9.7 × 10−8, and 0.002 for lung, liver, kidney, lymph node, and heart, respectively; Table S2), E2F targeting (q values of 2.8 × 10−26, 0.00672, 1.9 × 10−7, 9.9 × 10−9, and 0.047, respectively), and epithelial mesenchymal transition (EMT; q values of 2.1 × 10−22, 3.0 × 10−4, 2.9 × 10−7, 1.4 × 10−7, and 0.0287, respectively), but in the late infection (SARS-CoV-2 low), gene networks showed more inter-tissue heterogeneity in their disrupted pathways, including cytokine activity and inflammatory response. However, in both the SARS-CoV-2 high- and low-viral-load groups, the G2M checkpoint and E2F networks were consistently upregulated, indicating a core, persistent set of dysregulated cell-cycle regulation genes during early and late stages of infection.

Figure 2

Pathways and cell population changes of COVID-19

(A) The pathways that show significant differences in all or the majority of the five tissue types are shown, with statistical significance from GSEA testing across five kinds of tissue (colors in legend), including heart (n = 41), kidney (n = 27), liver (n = 40), lung (n = 40), and lymph node (n = 27). The x axis shows the normalized enrichment score for COVID-19-positive versus control, SARS-CoV-2 high or low versus control, and SARS-CoV-2 high versus Low comparisons.

(B) Cellular deconvolution distribution violin plots for the heart, kidney, liver, and lung. SARS-CoV-2 high (red), low (orange), and normal (blue) tissues are shown using a square root scale. The number of biological replicates is the same as in (A).

The DEGs and GSEA results were then examined for the largest differences between the infection level and stage (SARS-CoV-2 high, early COVID-19 infection versus SARS-CoV-2 low, late COVID-19 infection). Interestingly, the heart tissues showed the largest transcriptional differences, revealing that the later stage of the infection had a much greater impact on cardiac tissues (Figure S2A). To place these results into a larger context and to compare them with the findings of other data sets, we compared DEGs from each tissue with RNA-seq data from NP swabs as well as RNA-seq data from a publicly available data set on monocytes from COVID-19-positive and -negative patients (Figure S2B). While the highest correlations were seen within the same tissue types, most of the tissues with a high viral load showed a statistically significant, positive correlation with the DEGs in the NP swab samples (q < 0.01) when compared with normal/negative patients. In contrast, the later infection (low-viral-load) patients' tissues showed a negative correlation with the NP swab samples, indicating that the systemic impact of SARS-CoV-2 can be missed when not considering the biological impact on different organs. Interestingly, when matched with disease severity (high-, medium-, and low-viral-load groups within NP swab samples), the difference was bigger in the low and medium groups than in the high groups.

To create a more fine-grained analysis of the cellular gene expression states in each tissue, we used the cell deconvolution multiple signal classification algorithm (MuSiC) on each tissue’s RNA-seq data (see STAR Methods). The MuSiC results showed distinct disruptions of the transcriptional programs for each type of tissue in the COVID-19 patients and in the gain or loss of cell types (Figure 2B). Consistent with previous reports, the lung showed a loss of the capillary intermediate cells and alveolar epithelial cell types., Strikingly, we also found decreases in the major cell types in each organ type, suggesting a systemic disruption of the COVID-19 response. The kidney and liver showed a loss of proximal tubule in the kidney and of hepatocyte marker expressions in the liver but an increase in T cells in both organs. Furthermore, the heart showed a near-complete loss of the cell signatures for cardiomyocytes in both SARS-CoV-2 high and low viral loads (Figures 2B and 3A), despite no obvious gross or histologic changes in the heart (Figure 3B). This observation extends from previously reported cardiovascular involvements in COVID-19 and further defines the SARS-CoV-2-specific transcript and cellular changes in lung as well as other organs such as heart, liver, and kidney.,

Figure 3

Cellular disruption and tissue identity loss from SARS-CoV-2 infection

(A) Specific gene expression distributions for cardiomyocyte-related genes (from patient autopsy samples, n = 41).

(B) Sample collection strategy for COVID-19 autopsy samples from complete adult cases.

(A–J) Two cases with representative hematoxylin and eosin (H&E) images are shown: (A–E) Cases 58 and (F–J) 73. (A and F) Lung, (B and G) kidney, (C and H) heart, (D and I) liver, and (E and J) mediastinal lymph node (H&E stain). (B, D, E, I, and J) Original magnification ×50. (A, C, F, G, and H) Original magnification ×100.

(C) Gene expression changes of representative functional markers for each organ: heart (n = 41), kidney (n = 27), liver (n = 40), and lung (n = 40). Five genes related to the organ function were chosen (from the Human Protein Atlas), along with five housekeeping genes on the rightmost side. Darker gray shaded area represents the ranges of the housekeeping gene changes (baseline noise) and error bars show standard error from DESeq2 output, relative to the uninfected control.

We then looked at several markers related to the functions and processes of each organ (Figure 3C). While the lung showed the biggest changes in response to COVID-19, the loss of functional markers was specific to organ type. For example, surfactant proteins (i.e., SFTPA1, SFTPB, SFTPC), which are components of the alveolar lipoprotein complex crucial for gas exchange, were lost in the lung but increased in other organs. Similarly, markers associated with liver function (i.e., liver enzymes and proteins such as ALB, HAO1, and ALDOB) and solute carrier family proteins (i.e., SLC22A13, SLC34A1, uromodulin [UMOD]) were specifically lost in liver and kidney, respectively. In addition to the cardiomyocyte markers (Figure 3A), we also found that markers such as phospholamban (PLN, calcium pump inhibitor), heart and neural crest derivatives expressed 1 (HAND1), and troponin cardiac type (TNNT2) were specifically lost in the heart. This observation suggests that in addition to cell type losses, the disruptions due to COVID-19 affect the biological processes and functions for which each organ system is responsible. Consistent with recent reports related to virus-induced senescence found in COVID-19 lung models, these data also show that each organ system (lung, heart, liver, and kidney) exhibits its own distinct inflammatory response and resultant change of tissue identity.

Spatial and expression profiling of high and low SARS-CoV-2 infection in the lung

For a deeper examination of COVID-19 lung tissue, we then used the GeoMx Digital Spatial Profiling (DSP) platform to perform multiplexed high-resolution spatial transcriptomic profiling of 357 lung tissue ROIs from 16 patients. The ROIs were selected from deceased patients with COVID-19 (n = 8), nonviral acute respiratory distress syndrome (ARDS) (n = 2), influenza-induced ARDS (n = 3), and healthy tissues from individuals without infections as controls (n = 3) using nCounter Multiplex Analysis incorporated with targets for SARS-CoV-2 (i.e., IO360 panel plus COVID-19 Spike-in; Figure 1A). Among COVID-19 patients characterized with total RNA-seq, we identified patients that had high overall SARS-CoV-2 expression (SC2 high) or had low overall expression of SARS-CoV-2 (SC2 low) from their lung-tissue samples, and four representative samples from each group were selected for further downstream analysis. Serial sections were stained with an RNA scope probe against the viral Spike (S) gene, Syto13 (nuclear DNA), macrophages (CD68), immune cells (CD45), and epithelial cells (pan-cytokeratin) along with the GeoMx Cancer Transcriptome Atlas Panel (CTA, 1,811 targets) supplemented with 23 human genes associated with lung biology and two open reading frames (ORFs) from the SARS-CoV-2 genome (Figures 1A and S1). We chose tissue ROIs that captured three structural components of the lung, including vascular, airway, and alveolar regions (Figure S3).

We observed significant differences across tissue types within the lung (DEGs with q < 0.05 and > 1 fold change by DESeq2). For example, vascular ROIs showed an increase in ACTA2 and FLNA, while alveolar regions exhibited a mixture of signals from macrophage and monocyte genes (CD163 and CD68) with collagen (COL1A1, COL1A2, COL6A3). In large airway tissues, genes related to the mucosal layer (MUC4, MUC5AC, MUC5B) as well as cytokine-mediated signaling pathway genes (CCL20, CXCL1, CXCL6, CXCL6, MMP1) and type I interferon genes (IFIT3, ISG15, STAT1, MX1) were upregulated in COVID-19. Regardless of the tissue type, we saw consistent increases in genes such as B2M, CD81, GNAS, HLA.B, HSPA1A, HSPB1, and SERPING1 in COVID-19 lungs (Figure 4A; Table S3). When we summarized these genes using GO terms, tissue-type-specific processes in response to COVID-19 were distinct from overall inflammation (Table S3). While inflammatory responses were present in all tissue compartments, macrophage and monocyte activation were prominent in alveolar regions, while fibrosis occurred near the vascular region. Moreover, type I interferon- and cytokine-related pathways were found significant in large airway tissues, while complement activation was found in vascular tissues.

Figure 4

Spatial transcriptomics identifies tissue- and disease-specific differences

(A) Venn diagram of tissue-specific COVID-19 DEGs (relative to normal, adjusted p values < 0.01, >1-fold) from 16 patients across 357 ROIs in total.

(B) UpSet plot depicting intersections of disease-specific DEGs (p values < 0.05, fold change [|FC|] > 1), also from 16 patients across 357 ROIs in total.

(C) Volcano plot showing differences between normal (n = 3, 64 ROIs) and SARS-CoV-2 high-viral-load samples (n = 4, 86 ROIs) after accounting for compartmental differences. Top genes, in terms of p value or FC are indicated in gray, and COVID-19 Spike-in genes are labeled in black.

(D) Ternary plot of a combined analysis of SARS-CoV-2 high (n = 4, 86 ROIs), low (n = 4, 97 ROIs), and normal (n = 3, 64 ROIs), where genes are projected away from the center based on their marginal means. Genes upregulated in a single group approach that group’s corner. Top genes in terms of p value or FC are labeled. Genes with color were significant with p < 0.05. Genes with outlines are significant after correcting for multiple hypothesis testing.

We also compared differences across injury sources (SARS-CoV-2 infection, influenza, and nonviral ARDS) and found significant differences between normal and SARS-CoV-2-high samples, even after accounting for compartmental variability (Figures 4B and S4). These differences included a decrease in the expression of surfactant genes (SFTPA1, SFTPB, and SFTPC) associated with type 1 and type 2 pneumocytes as well as an enrichment for genes associated with basal cells (TP63) and club cells (SCGB1A1), and several immune markers (e.g., HLA-B, HLA-E, p all < 0.05, mixed effects model) in SARS-CoV-2-high patients compared with both normal and SARS-CoV-2-low patients. Ternary plots of a combined analysis of normal, SARS-CoV-2-high and -low tissues—where transcripts are projected away from the center based on their marginal means—revealed upregulation of several genes enriched in each set of lung tissues. Enrichments included SFTPA1, SFTPB, and SFTPC (alveolar epithelial cell markers) in normal lungs, CLU (lung injury and repair) and S100A9 (enriched in activated macrophages) in SARS-CoV-2-low lungs, and TP63 (basal cell), ID1 (upregulated and a key regulator of lung injury and repair), and interferon-regulated genes including IFI6, IFI27, ISG15, and LY6E in SARS-CoV-2-high lungs (Figures 4C and 4D). Enrichment of interferon-stimulated genes was observed only in SARS-CoV-2-high samples, which corresponds to early stages of the infection, by timing of disease onset and histopathology (Figures S1D and S1E). Moreover, we observed enrichments of CASP3 and ID1, suggesting ongoing cellular injury and repair responses in SARS-CoV-2-high patients. In contrast, we find an enrichment of several markers of pulmonary fibrosis (CLU, COL1A1, COL1A2, and COL3A1) in SARS-CoV-2-low patients (these correspond to the later stages of infection), indicating that there are two distinct stages of infection.

We also examined the spatial transcriptome data for differences between the high- and low-SARS-CoV-2- and the IAV-infected lung samples. This analysis revealed a significant enrichment of THBS1 and NR4A1 in the influenza samples, which are both genes that have been shown to be engaged in response to influenza-induced lung injury (Figure S4B). When lung tissues from SARS-CoV-2-high and -low samples were compared with those from nonviral ARDS, enrichment for S100A8 and S100A9 was observed, which is consistent with the significant enrichment of neutrophils in these samples and was previously suggested to be a driver of COVID-19 pathogenesis (Figure S4A)., Importantly, an enrichment for genes involved in lung injury and repair (ID1), as well as those involved in type I interferon responses including IFI6, IFI27, ISG15, and LY6E in SARS-CoV-2-high lungs, is observed even when compared with IAV infection of the lung and nonviral lung injury. Similarly, we find enrichment of several markers of pulmonary fibrosis (CLU, COL1A1, COL1A2, and COL3A1) in SARS-CoV-2-low samples as compared with nonviral and viral ARDS samples, exemplifying the profound lung injury and fibrosis during later stages of COVID-19.

COVID-19-specific heterogeneity and spatial tropism in the lung

Regardless of specific lung-tissue types, the expression profiles and respective cell-type proportions were enough to distinguish and cluster normal versus COVID-19 (high and low SARS-CoV-2 viral loads) lung ROIs (Figure 5A), which indicates that the SARS-CoV-2 infection is altering the cellular interaction landscape and composition of the lung tissue. To gain additional insights into the heterogeneity and spatial tropism of lung infection independent of ROI origins, we assessed how each ROI compares with the typical healthy tissue. Some ROIs (especially of alveolar regions) in COVID-19 showed significantly diminished similarity to those of healthy lungs (Figures 5B, S5A, and S5B). Interestingly, the loss of this similarity did not result in convergence to other tissue types we observed (i.e., loss of alveolar similarity score did not increase the vascular or large airway scores). To evaluate the cause of this loss of similarities, we discovered genes attributing disease- and tissue-specific clustering and identified 35 genes with the highest contrast (Figure 5C). For example, some of these genes were related to regulatory T cell differentiation in COVID-19 high- and low-viral-load alveolar tissues (PLK1, ATM) and IL6/JAK/STAT3 signaling in SARS-CoV-2-low vascular regions (SOCS1, CXCL11; Table S3). In addition to the known molecular responses to COVID-19, these genes can be used as a set to locate infection response within the lung and disease states.

Figure 5

Evaluation of tissue identity loss and heterogeneity from infection

(A) Principal-component analysis (PCA) to see sample- and disease-clustering. Colors denote disease conditions (64 normal ROIs, 86 and 97 SARS-CoV-2 high- and low-viral-load ROIs), while shapes show the tissue types (alveolar, large airway, and vascular regions).

(B) Similarity score distribution when compared with a tissue-specific (278 alveolar, 45 large airway, and 34 vascular ROIs) healthy reference gene profile (ns, non-significant, ∗p ≤ 0.05, ∗∗p ≤ 0.01, ∗∗∗p ≤ 0.001, and ∗∗∗∗p ≤ 0.0001).

(C) Genes identifying disease- and tissue-specific conditions.

(D) Cell-type proportions in normal versus COVID-19. The median and quartiles are noted by the box plot inside. p value two-tailed t tests were done to compare the means (ns, ∗p ≤ 0.05, ∗∗p ≤ 0.01, ∗∗∗p ≤ 0.001, and ∗∗∗∗p ≤ 0.0001).

(E) Enrichment of cell type- and COVID-19-specific gene signatures from 11 patients across 247 ROIs.

To pinpoint specific changes in different lung tissue types, count estimates of 15 distinct cell types were imputed based on gene expression profiles from the Human Cell Atlas (HCA) adult lung dataset, including neutrophil profiles derived from single-nucleus RNA-seq (snRNA-seq) of lung tumors (see STAR Methods). Consistent with other studies,, we observed that COVID-19 was associated with an increase in tissue-infiltrating immune cells, including T cells, natural killer (NK) cells, monocytes, and macrophages (Figures 5D and S5C). Some immune cell types, such as monocytes, NK cells, and regulatory T cells, showed a statistically significant increase only in the SARS-CoV-2-high conditions. Also, while fibroblasts and endothelial cells increased in both high- and low-viral-load samples (52% and 65% increase, respectively), type 1 and type 2 alveolar epithelial cell proportions decreased (26% and 16% decrease, respectively), reflecting the ongoing tissue remodeling or selective epithelial cell death induced by infection (Figure S5C). We validated these findings by comparing the expression levels with cell-type-specific COVID-19 gene signatures derived from snRNA-seq, which showed enrichments across different cell types (Figure 5E; Table S3). Using an orthogonal approach, we stained the lung tissues with Masson’s trichrome and observed a statistically significant increase (p < 0.01) in cellular collagen-rich areas, confirming the increase in lung fibroblasts (Figure S3A).

Given the observed changes in cell proportions that are induced during SARS-CoV-2 infection, we next examined the impact on the co-occurrence of these changes as a metric of intrapulmonary cellular heterogeneity (Figure 6). Pairwise correlations of all detected cell types under five different conditions (SARS-CoV-2 high and low, IAV, nonviral ARDS, and control) were calculated and visualized (Figure 6A). In normal lungs, we observed three clear cellular correlation clusters: (1) monocytes, fibroblasts, T cells, and NK cells, (2) neutrophils and ciliated cells, and (3) plasmacytoid dendritic cells (pDCs), macrophages, and B cells. While perturbations of these cellular clusters were observed across all injury conditions, the NK cell-T cell correlation was lost only in the SARS-CoV-2 high-viral-load patients and was not present in low-viral-load patients, with the low viral load corresponding to the later stages of infection (Figures 6A and 6B). We then quantified the correlation differences in the lungs' cellular landscapes between SARS-CoV-2-high and -low patients, which indicated that the greatest changes were in the monocyte-T cell correlations and dendritic-neutrophil correlations (Figure S6A), further supporting the view that SARS-CoV-2-specific T cell activity may be disrupted. Of note, it may be possible that the changes in correlation reflect both the generation of long-term memory and T cell-mediated killing of infected epithelial cells.

Figure 6

Cellular tropism and heterogeneity in response to COVID-19

(A) Correlation matrix of pairwise cell-type correlations. Statistically insignificant correlations were not displayed (gray area, p < 0.05); from 16 patients across 357 ROIs in total, as shown in Figure 1A.

(B) Average proportion changes of the cell types relative to normal. The cell types were ordered by average increase in proportions across 357 ROIs. Error bars indicate 0.5∗SD. On the left, stacked bar plot depicts the overall proportions by conditions.

(C and D) Entropy estimates of the (C) tissue types and (D) cell types within the normal (n = 3, 64 ROIs), SARS-CoV-2-high (n = 4, 86 ROIs) and -low (n = 4, 97 ROIs), influenza (n = 2, 46 ROIs), and ARDS (n = 3, 67 ROIs) conditions.

To compare changes across the disease states (COVID-19, influenza, ARDS), we visualized the average proportions of each cell type. Macrophage and neutrophil population levels were found to be much higher in the COVID-19 lungs (both high and low viral loads), while T cell, monocyte, and epithelial cell numbers were much lower than normal. The entropy estimates of the given cell types (Figure 6C), which are a function of variance, showed that macrophages and neutrophils commonly displayed increased heterogeneity (across ROIs and patients) across all injury conditions. Entropy in fibroblast, epithelial, vessel, and T cell populations was greater in SARS-CoV-2-high and -low lung tissue than in IAV infection and ARDS (121.5%, 104.7%, 111.9%, and 124.7% changes, respectively, when compared with the average of SARS-CoV-2 high/low and that of IAV/ARDS; Figure 6D), suggesting an increase in cellular heterogeneity. When measuring and comparing cellular tropism within the tissue, nearly all COVID-19-positive ROIs showed an increase in heterogeneity of the cell populations, with the sole exception being vascular regions with high COVID-19 (Figure 6D). The decrease in heterogeneity in the vascular regions of high-viral-load conditions is mainly from the decrease of fibroblasts and epithelial, T, and NK cells (Figure S6B). The signals from B cells are also specific to large airway tissues, and this observation is consistent with the cell fraction increase in large airway ROIs of COVID-19 samples (Figure S6B). Such changes in heterogeneity are related to specific damages and tissue dysregulations. Single-sample GSEA (ssGSEA) of macrophage, neutrophil, and T cell regulatory pathways showed enrichment in COVID-19, even when compared with influenza and ARDS, including macrophage activation and apoptotic process (1.3- and 2.8-fold increase in averaged ssGSEA scores relative to normal, with p values of 2.91 × 10−8 and 0.001) in patients with high viral loads (Figures S6C and S6E).

Discussion

In this study, we established a clinical analytical pipeline to collect and examine autopsy samples to elucidate and compare the spatial transcriptional landscape induced by SARS-CoV-2, IAV, and nonviral ARDS. By combining transcriptional profiles of 39 patient autopsy tissues from heart, liver, lung, kidney, and lymph nodes, we presented body-wide transcriptome changes in response to COVID-19. Across all tissues, we found system-wide disruptions of cellular and transcriptional pathways and matched the lung data with spatial protein and expression profiling (GeoMx, across 357 tissue sections from 16 representative patient lung samples). We also identified tissue-compartment-specific damage (alveolar, vascular, and large airway compartments within the lung tissue) and the loss of tissue type identity caused by the SARS-CoV-2 infection, which correlated with viral loads (high versus low) and the clinical course of infection.

Patient lung tissue samples containing significant levels of SARS-CoV-2 showed enrichment for genes related to a variety of immune markers specific to certain immune cells and lung injuries as well as for interferon-stimulated genes (e.g., IFI27, IFITM1, and LY6E) and macrophage activation (S100A9, TYMP, and SERPING1). In contrast, patient lung tissue samples containing low levels of SARS-CoV-2 RNA show enrichment for COL1A1 and other markers of pulmonary fibrosis. Compared with other viral-related diseases (influenza), COVID-19 tissue samples still show significant enrichment for genes involved in lung injury and repair, interferon-signaling genes, and pulmonary fibrosis. Of note, COVID-19 (high-viral-load samples), influenza, and ARDS each show differential HLA-B and -C expressions, which are known mediators of NK and T cell activation, and which can mediate host risk of infection; one example is enrichment for HLA-DRB5, whose expression and specific gene polymorphisms are associated with pulmonary fibrosis and severity., Compared across different disease types, all diseases, COVID-19 (low viral load), influenza, and ARDS showed enrichment for DMBT1, a gene known to be upregulated and dysregulated in pulmonary injury and fibrosis., Virus-related diseases (SARS-CoV-2 high viral load and influenza) in particular showed significant changes on the expression of lung epithelial-cell-related transcripts (i.e., ACTB, C1R, and FN1); such changes are known markers of lung-injury gene signatures.,

The spatial analysis platform (Nanostring GeoMx) enabled us to investigate the impact of the disease by incorporating cellular and spatial organization. Consistent with recent reports from bulk cellular profiling, we observed an increase in immune cell types and fibroblasts in COVID-19 but a decrease in alveolar epithelial cells. In SARS-CoV-2 low-viral-load tissues, the proportions of some immune cells (i.e., monocytes, NK cells, or regulatory T cells) were normal, but fibroblasts and vessel cells still exhibited an increase similar to those observed in the high-viral-load samples. Some of these cell types form a “cellular correlation cluster,” and these co-occurrence clusters of cellular changes are uniquely disrupted in COVID-19 (relative to influenza and ARDS), particularly in the COVID-19 high-viral-load sample group. While macrophages and neutrophils showed an increase in entropy across all lung-related injury conditions, NK and T cells showed an increase only in COVID-19 samples. While few studies have interrogated the tissue environments, multiple studies have examined the changes occurring during COVID-19 infection in the peripheral blood and have identified poor T cell responses and T cell dysregulation.,29, 30, 31, 32, 33 Together, these findings highlight the robust and dynamic nature of SARS-CoV-2 engagement with tissue homeostatic processes and that the stage of COVID-19 infection impacts the pathophysiological landscape of the lung.

When these spatial transcriptomics data were compared to the multi-organ bulk RNA-seq data from the autopsy issues, confirmatory as well as additional signatures of COVID-19 disease were found. First, fibroblasts and immune cells (i.e., macrophages) were increased in most tissue types, while tissue-specific cell types, such as alveolar epithelial types 1 and 2 cells in the lung and cardiomyocytes in the heart, showed a decrease in COVID-19 relative to controls. Increases in fibroblasts, endothelial cells, and immune cells may be impacted by a variety of immune activations within each organ, particularly as a long-term response to the infection. This observation may also be related to the decrease in characteristic transcriptomic signatures of the main cell type within each organ, which may contribute to the morbidity and mortality of COVID-19. For example, the reduction in the cardiomyocyte cell fraction correlated with a reduction in several transcripts encoding sarcomeric and contractile proteins (in both high and low viral loads),, representing a persistent transcriptional perturbation and potentially long-term, cardiac-specific impact of COVID-19.

These data support the view of both tissue-specific and time-dependent biological responses to the stages of SARS-CoV-2 infection, and this is buttressed by orthogonal data. For example, DEGs observed from NP swabs showed a high correlation to tissue-specific DEGs in the early stages of infection but very little correlation with later infection. The observed correlations likely reveal monocyte migration and infiltration into tissues during SARS-CoV-2 infection, as seen by others. It has also been reported that monocyte depletion/migration is associated with kidney disease, inducing lupus-like symptoms, which could potentially explain the correlation with kidney tissue. For lymph nodes, there exists evidence in the literature that SARS-CoV-2 will impact the lymph nodes at an early stage of infection, potentially causing T cell lymphopenia and possibly being responsible for focal necrosis seen in the lymph nodes. Nonetheless, the lung, heart, and lymph nodes were the tissues most disrupted by infection.

Overall, these data represent one of the largest autopsy series of COVID-19 disease and synthesize several orthogonal methods, including bulk transcriptomics, digital spatial transcriptomics, multiple imaging technologies, and computational analysis, to build a map of SARS-CoV-2 pathophysiology. Given the ability to combine bulk transcriptomics data from multiple organ types, we find organ-specific changes to immune responses and the loss of tissue functions unique to COVID-19 disease, which can help additional studies and methods for mitigating the systemic damage caused by the SARS-CoV-2 virus across the body.

Limitations of the study

While the size of the cohort we used for spatial profiling is smaller than that used for bulk transcriptomics, we believe an in-depth characterization of patient tissues is crucial in directly validating and supporting the findings from cellular and animal models., Spatial profiling technology captures key aspects of COVID-19 that bulk data cannot, including the locations of transcriptional and cellular changes caused by the disease (especially in late infection) and spatial heterogeneity of cell types. The compartment-specific COVID-19 gene signatures would need further validation (e.g. chromatin states and additional validation), as we rely on computational methodologies and deconvolution techniques, and could also benefit from matched profiling of peripheral blood, examination of strain type, and also differences in vaccination status. Nonetheless, this molecular map of COVID-19 represents a needed cellular and molecular atlas for the community, which can inform future studies into COVID-19 progression and SARS-CoV-2 pathology.

STAR★Methods

Key resources table

Resource availability

Lead contact

Further information and requests for resources and reagents should be directed to the lead contact, Christopher E. Mason (chm2042@med.cornell.edu).

Materials availability

This study does not involve new unique reagents or materials.

Experimental model and subject details

IRB statement

Tissue samples were provided by the Weill Cornell Medicine Department of Pathology. The Tissue Procurement Facility operates under Institutional Review Board (IRB) approved protocol and follows guidelines set by Health Insurance Portability and Accountability Act (HIPAA). Experiments using samples from human subjects were conducted in accordance with local regulations and with the approval of the IRB at the Weill Cornell Medicine. The autopsy samples are considered human tissue research and were collected under IRB protocols 20-04021814 and 19-11021069. All autopsies have consent for research use from next of kin, and these studies were determined as exempt by IRB at Weill Cornell Medicine under those protocol numbers.

Patient sample collection

All autopsies are performed with consent of next of kin and permission for retention and research use of tissue. Autopsies were performed in a negative pressure room with protective equipment including N-95 masks; brain and bone were not obtained for safety reasons. All fresh tissues were procured prior to fixation and directly into Trizol for downstream RNA extraction. Tissues were collected from lung, liver, lymph nodes, kidney, and the heart as consent permitted. For GeoMx, RNAscope, trichrome and histology tissue sections were fixed in 10% neutral buffered formalin for 48 hours before processing and sectioning. These cases had a post-mortem interval of less than 48 hours. For bulk RNA-seq tissues, post-mortem intervals ranged from less than 24 hours to 72 hours (with 2 exceptions - one at 4 and one at 7 days - but passing RNA quality metrics) with an average of 2.5 days. All deceased patient remains were refrigerated at 4°C prior to autopsy performance.

Method details

qRT-PCR

Total RNA was extracted in TRIzol (Invitrogen) according to the manufacturer’s instructions. To quantify viral replication, measured by the expression of sgRNA transcription of the viral N gene, one-step quantitative real-time PCR was performed using SuperScript III Platinum SYBR Green One-Step qRT-PCR Kit (Invitrogen) with primers specific for the TRS-L and TRS-B sites for the N gene as well as ACTB as an internal reference. Quantitative real-time PCR reactions were performed on an Applied Biosystems QuantStudio 6 Flex Real-Time PCR Instrument (ABI). Delta-delta-cycle threshold (ΔΔCT) was determined relative to ACTB levels and normalized to mock infected samples. Error bars indicate the standard deviation of the mean from three biological replicates. The sequences of primers/probes are provided in the key resources table.

RNA-seq analysis

Patient specimens were processed as described in Butler et al. Clinical metadata is summarized in Table S1. Briefly, nasopharyngeal (NP) swabs were collected N using the BD Universal Viral Transport Media system (Becton, Dickinson and Company, Franklin Lakes, NJ) from symptomatic patients. Total Nucleic Acid (TNA) was extracted from using automated nucleic acid extraction on the QIAsymphony and the DSP Virus/Pathogen Mini Kit (Qiagen). Autopsy tissues were collected from lung, liver, lymph nodes, kidney, and the heart and were placed directly into Trizol, homogenized, and then snap frozen in liquid nitrogen. At least after 24 hours these tissue samples were then processed via standard protocols to isolate RNA.

For RNA library preparation, all samples' TNA were treated with DNAse 1 (Zymo Research, Catalog #E1010). Post-DNAse digested samples were then put into the NEBNext rRNA depletion v2 (Human/Mouse/Rat), Ultra II Directional RNA (10 ng), and Unique Dual Index Primer Pairs were used following the vendor protocols from New England Biolabs. Completed libraries were quantified by Qubit and run on a Bioanalyzer for size determination. Libraries were pooled and sent to the WCM Genomics Core or HudsonAlpha for final quantification by Qubit fluorometer (ThermoFisher Scientific), TapeStation 2200 (Agilent), and qRT-PCR using the Kapa Biosystems Illumina library quantification kit.

NYGC RNA sequencing libraries were prepared using the KAPA Hyper Library Preparation Kit + RiboErase, HMR (Roche) in accordance with manufacturer's recommendations. Briefly, 50-200ng of Total RNA were used for ribosomal depletion and fragmentation. Depleted RNA underwent first and second strand cDNA synthesis followed by adenylation, and ligation of unique dual indexed adapters. Libraries were amplified using 12 cycles of PCR and cleaned-up by magnetic bead purification. Final libraries were quantified using fluorescent-based assays including PicoGreen (Life Technologies) or Qubit Fluorometer (invitrogen) and Fragment Analyzer (Advanced Analytics) and sequenced on a NovaSeq 6000 sequencer (v1 chemistry) with 2x150bp targeting 60M reads per sample.

Spatial transcriptomics analysis

Gene Expression profiling of freshly extracted RNA from formalin fixed paraffin-embedded (FFPE) lung samples was performed using the NanoString PanCancer IO360 panel with custom probes for SARS-CoV-2 viral genes. After normalization, high and low COVID-19 clusters were identified by unsupervised analysis, and samples from each cluster were selected for additional profiling. GeoMx Digital Spatial Profiling (DSP) was performed on these samples, and control samples from non-viral ARDS, influenza, and normal lung tissues following standard protocols using the COVID-19 Immune Response Atlas., Samples were stained with immunofluorescent antibodies for CD68, CD45, PanCK, and DNA (Syto-13). Regions profiled included vascular zone, large airway, alveoli zone, and IF-guided segments focused specifically on macrophages. Samples were sequenced on an Illumina NextSeq, processed and filtered for quality as described in supplementary methods. Differential expression was assessed on the resulting normalized data using mixed effect models, accounting for intra-patient heterogeneity to assess differences between SARS-CoV-2 high and low viral load samples, and among distinct tissue structures profiled. Cell deconvolution of the GeoMx data was performed using the SpatialDecon R package. Gene set enrichment analysis (GSEA) was performed to qualify coordinate gene expression changes quantified during differential expression analysis.

Quantification and statistical analysis

Differential gene analysis

RNAseq data was processed through the nf-core/rnaseq pipeline. This workflow involved quality control of the reads with FastQC, adapter trimming using Trim Galore (https://github.com/FelixKrueger/TrimGalore), read alignment with STAR, gene quantification with Salmon, duplicate read marking with Picard MarkDuplicates (https://github.com/broadinstitute/picard), and transcript quantification with StringTie. Other quality control measures included RSeQC, Qualimap, and dupRadar. Alignment was performed using the GRCh38 build native to nf-core and annotation was performed using Gencode Human Release 33 (GRCH38.p13). FeatureCounts reads were normalized using variance-stabilizing transform (vst) in DESeq2 package in R for visualization purposes in log-scale. Cell deconvolution was performed using MuSiC on single cell reference datasets for lung, liver, kidney, and heart.,45, 46, 47, 58 Immune cell deconvolution was performed on lymph node samples using quanTIseq. Differential expression of genes was calculated by DESeq2. Differential expression comparisons were done as either COVID + cases versus COVID- controls for each tissue specifically, correcting for sequencing batches with a covariate where applicable, or pairwise comparison of viral levels from the lung as determined by nCounter data. In the volcano plot protein coding genes were plotted using Gencode classifications using -log10 (adjusted value) and log2 fold-change metrics. Genes with BH-adjusted p value < 0.01 and absolute log2 fold-change greater than 0.58 (at least 50% change in either direction) were taken as significantly differentially regulated. Genes were ranked by their Wald statistic and their log2 fold-change values and used as input for gene set enrichment analysis (GSEA) on the molecular signatures database (MSigDB).,48, 49, 50 Any signature with adjusted p value < 0.01 was taken as significant. List of differentially expressed genes and significantly enriched pathways are reported in Table S2.

Pairwise correlations of cell types by conditions

Correlation matrix visualizes the Pearson correlation coefficient by cell types within each disease condition. Statistically insignificant correlations (p-value > 0.05) are filtered and identified clusters of positive and negative correlation is marked. The correlations from SARS-CoV-2 high and low viral load samples are compared with normal, using R package cocor (v1.1-3). Briefly, the correlation coefficients are tested using Fisher’s r-to-Z transformation to quantify the differences between the two correlations. To quantify correlations, each data point (or correlation coefficient) corresponds to a fisher-tested correlation (z statistics and –log(P-value) for x and y axes, respectively). The entropy calculations were done with the synRNASeqNet R package (v1.0, entropyML function, https://github.com/cran/synRNASeqNet). The deconvoluted cell counts were used as an input to run maximum likelihood entropy calculations.

Similarity analysis

The consensus gene profiles for alveolar, large airway, and vascular healthy samples were built by taking the average gene profiles from healthy ROIs of respective tissue origin. To validate gene profiles, healthy ROIs were randomly sampled (1, 5, and 10 ROIs) and compared with the consensus profiles. With these gene profiles, we assessed the similarity of the profile from each ROI with the reference profile by taking cosine similarity (1 being closer to the reference, 0 being orthogonal, Figure S5A). To identify genes specific to tissue- and disease- states, we performed logistic regression with L1 norm to model the gene expression profiles. The logistic regression was done with glmnet (v4.1-1). The genes with highest coefficients were filtered to identify 35 genes that may distinguish diseased tissue types (Table S3).

GeoMx transcriptomic data normalization and quantification

Spatial transcriptomics analysis

Discerning viral load from bulk nCounter screening—Bulk expression profiling was performed to identify COVID-19 patients with high vs low viral load. To do this, RNA from fresh TRIzol extracted and fixed lung tissue from 29 COVID-19 autopsies, 4 non-COVID-19 lung injury and 3 controls were evaluated by bulk expression analysis using NanoString’s nCounter PanCancer IO360 Panel plus custom probes for eight SARS-CoV-2 viral genes (encoding ORF7a, surface glycoprotein, nucleocapsid phosphoprotein, ORF8, envelope protein, ORF3a, membrane glycoprotein, and ORF1ab) to assess viral content. At least 100 ng of RNA was loaded for hybridization and quantified by the nCounter MAX Analysis System (NanoString Technologies, Seattle, WA, USA).

Transcript counts were normalized to ERCC positive controls and housekeeper reference gene expression prior to analysis. Hierarchical clustering of nCounter results revealed two clusters of high and low severity and four “SARS-CoV-2 high” and four “SARS-CoV-2 low” patients were selected (Figure S1D). These eight samples were analyzed with two influenza-infected patients, three non-viral ARDS patients, and three normal lung control patients using the GeoMx platform.

RNA/NGS slide preparation for GeoMx DSP

For GeoMx DSP slide preparation, we followed the GeoMx DSP slide prep user manual (MAN-10087-04). Briefly, tissue slides were baked in a drying oven at 60 °C for 1 hour and then loaded to Leica Biosystems BOND RX FFPE for deparaffinization and rehydration. After the target retrieval step, tissues were treated with Proteinase K solution to expose RNA targets followed by fixation with 10% NBF. After all tissue pre-treatments were done, tissue slides were unloaded from the Leica Biosystems BOND RX and incubated with RNA probe mix (COVID-19 Immune Response Atlas panel) overnight. The next day, tissues were washed and stained with tissue visualization markers; CD68-647 at 1:400 (Novus Bio, NBP2-34736AF647), CD45-594 at 1:10 (NanoString Technologies), PanCK-532 at 1:20 (NanoString Technologies) and/or SYTO 13 at 1:10 (Thermo Scientific S7575).

GeoMx DSP sample collections

For GeoMx DSP sample collections, we followed the GeoMx DSP instrument user manual (MAN-10088-03). Briefly, tissue slides were loaded on the GeoMx DSP instrument and then scanned to visualize whole tissue images. For each tissue sample, we collected 4 types of functional tissue regions: vascular zone, large airway, alveoli zone, and macrophages. Each tissue region was carefully selected by a board-certified pathologist. Regions of interest (ROIs) were then segmented with corresponding fluorescent tissue markers, when available within the region. Twenty-four to twenty-three GeoMx DSP regions were selected per tissue and collected following UV illumination within the defined segment as described in Merritt et al. Compartments that were segmented within a region of interest were collected separately as unique areas of illumination (AOIs).

GeoMx DSP NGS library preparation and sequencing

Each GeoMx DSP sample was uniquely indexed using Illumina’s i5 x i7 dual-indexing system. 4 μL of a GeoMx DSP sample was used in a PCR reaction with 1 μM of i5 primer, 1 μM i7 primer, and 1× NSTG PCR Master Mix. Thermocycler conditions were 37 °C for 30 min, 50 °C for 10 min, 95 °C for 3 min, 18 cycles of 95 °C for 15 sec, 65 °C for 60 sec, 68 °C for 30 sec, and final extension of 68 °C for 5 min. PCR reactions were purified with two rounds of AMPure XP beads (Beckman Coulter) at 1.2× bead-to-sample ratio. Libraries were paired end sequenced (2 × 75) on a NextSeq550 up to 400 million total aligned reads.

Processing and filtering raw NGS data

Three hundred seventy-nine AOIs plus non-template controls (NTCs) were sequenced, producing about 1.3B reads (with about ∼11% unique). NextSeq-derived FASTQ files for each sample were compiled for each AOI using Illumina’s bcl2fastq program and then demultiplexed and converted to Digital Count Conversion (DCC) files using Nanostring’s GeoMx DnD pipeline (v1). These DCC files were then converted to an expression count matrix using a custom python script. A minimum of 10,000 reads were required for each non-NTC sample (2 AOIs removed). Probes were checked for outlier status by implementing a global Grubb’s outlier test with alpha set to 0.01. The counts for all remaining probes for a given target were then collapsed into a single metric by taking the geometric mean of probe counts. A count of 1 was added to any probe that yielded 0 counts before the geometric mean was taken. For each sample, RNA probe pool specific negative probe normalization factor was generated based on the geometric mean of negative probes in each pool.

Quality control and AOI filtering

Following initial screening above, there were 373 AOIs interrogated using DSP that span 16 patients and three compartments (288 alveolar, 48 large airway, and 37 vascular regions). Of these, 370 AOIs yielded greater than 50 nuclei. The 75th percentile of the gene counts (i.e., geometric mean across all non-outlier probes for a given gene) for each AOI were calculated and normalized to the geometric mean of the 75th percentile across all AOIs to give the upper quartile or Q3 normalization factors for each AOI. The distribution of these Q3 normalization factors were then checked for outliers defined as any AOI greater than two standard deviations from the mean log2 Q3 normalization factor. This criterion removed 15 AOIs that fell below the range and 1 AOI that fell above the range. Following AOI filtering, 358 (of 373, ∼96%) AOIs were used for downstream analyses.

Removal of gene outliers and normalization

Gene outliers were detected by a limit of quantification (LOQ) approach. The LOQ for each AOI was defined as the sample’s negative geometric mean multiplied by its negative standard deviation raised to the power of two. Any target (1,837 total) that was not above LOQ in at least 2% of AOIs were deemed prohibitively low expressors and were removed from analysis. This feature-based filtering approach discarded 171 genes (9.3%) leaving 1,666 genes. Genes were normalized by the Q3 approach as above.

Deconvolution of cell proportions using GeoMx

Cell deconvolution methods followed that of Desai et al. Cell mixing proportions were performed using the R package SpatialDecon using the cell profile matrix based upon the Human Cell Atlas adult lung 10× dataset and appended with a neutrophil profile derived from snRNA-seq of lung tumors. ROIs were selected to be representative across the FFPE slides, based on morphology and immunofluorescence of each tissue.

Differential expression analysis

Two different sets of differential expression (DE) analyses were performed. Common in both sets of models were five groups (COVID-19 SARS-CoV-2 high, SARS-CoV-2 low, influenza, Non-viral ARDS, and Normal) and three compartments (alveolar, large airway, vascular). In DE set one, differences between all 10 pairwise groupings were performed to look for differences between pairwise comparisons and to serve as the basis for downstream gene set enrichment analysis (GSEA). In the second DE set, SARS-CoV-2 high and low viral load groups were compared against one of three non-COVID groupings to identify genes that are up- or downregulated in COVID-19 (sensu lato) relative to non-COVID-19 groups and to identify genes that are consistently differentially expressed in SARS-CoV-2 high and in low separately.

In the first model, each 10 pairwise groupings were considered separately (e.g., SARS-CoV-2 low vs SARS-CoV-2 high). DE analysis was performed by fitting each gene’s normalized log2 expression level using a Linear Mixed Effect model to account for interpatient variation with the R package lmerTest. Patient ID was used as a random effect (random intercept) and grouping, compartment, and grouping-by-compartment interactions were used as fixed effects. Satterthwaite's approximation for degrees of freedom for P-value calculation was used.

In the second set of models, SARS-CoV-2 high and SARS-CoV-2 low AOIs were always included, and DE was used to determine how genes were up- or downregulated compared to each of three non-COVID groups. As such, there were three different “sets” (influenza vs SARS-CoV-2 high vs SARS-CoV-2 low; Non-viral ARDS vs SARS-CoV-2 high vs SARS-CoV-2 low; and Normal vs SARS-CoV-2 high vs SARS-CoV-2 low). For a given gene and a given set, AOIs were first filtered to exclude non-members (i.e., exclude normal and Non-viral ARDS in the set “influenza vs SARS-CoV-2 high vs SARS-CoV-2 low”). The log2 expression of sample for a given gene was fit to a mixed effect model with group (three levels) and compartment (three levels) and their interaction as fixed effects and Patient ID as a random effect. The Least Squares (LS) means or “marginal means” were estimated for each comparison (i.e., log2 means for levels of “group” which are averaged over the levels of other factors in the model). In addition to the LS means, the pairwise P-values for all 3 comparisons were computed.

To visualize the marginal means for each gene relative to SARS-CoV-2 high, SARS-CoV-2 low, and a given normal group, the three-dimensional data were collapsed into two-dimensional ternary plots. Specifically, for a given gene g in set S, the three-element vector of marginal means can be expressed as g. By convention, the order of the elements of g were normal, SARS-CoV-2 high, and SARS-CoV-2 low. Elements of g were then rescaled while preserving their relative relationship by multiplying by a scaling factor and converting from log2 space to linear space. Then g was scaled further by dividing each element by the element with the minimum value. To convert g from a vector of three to a vector of two (representing points along a simplex plane), the new x coordinates were calculated by:where represents the ith element in vector g. Similarly, the y coordinates were calculated by:

A given gene is then assigned a “corner” by how close it is from each of the three simplex’s corners.

Each gene has three p-values associated with it (from the pairwise contrasts above). The contrast with the lowest p-value was selected to represent a given gene (i.e., this corresponded well with the corner that the gene was assigned to). P-values were then adjusted to account for multiple comparisons by using the Benjamini-Hochberg procedure. Differentially expressed genes from this analysis is included in Table S3.

Gene set enrichment

MA plots from the 10 pairwise DE analyses (DE model 1) were used to ensure that low expressors were not accounting for the majority of the large fold changes. Gene Set Enrichment Analysis (GSEA) was conducted using the R package fGSEA with MSigDB Hallmark and Reactome databases. Gene ranks were based on the log2 fold change from the individual DE analyses and gene sets were bound between 15 and 500 genes. P-values for enrichment were estimated by 1,000 permutations of the data. The pathways were then sorted based on adjusted p-value first and then by their Normalized Expression Score (NES).

For the 10 most extreme pathways in each direction as well as the COVID-19 spike-in genes, single sample GSEA (ssGSEA) was performed using the R package gsva with a min and max size of 15 and 500, respectively. Enrichment scores for a given pathway were rescaled by dividing AOI’s enrichment score by the mean across samples and then rescaled between 0 and 1. These rescaled enrichment scores were then visualized for each AOIs' x and y coordinates from their respective FFPE slide.

Histology and imaging analysis

Sections from SARS-CoV-2 high (4), SARS-CoV-2 low (4) and normal lungs (3) cases used for the GeoMx analysis were stained using hematoxylin and eosin and Masson’|'s Trichrome according to standard protocol. Four 20× regions were randomly selected from each slide the Color deconvolution2 algorithm for ImageJ and cellular trichrome rich areas were manually selected and measured for pixel area (cellular areas on red deconvolution, trichome on blue deconvolution). Total image pixel area was used to determine percent fibroblast-rich trichrome positive zones.

Correlation plot comparing different COVID-19 samples and tissues

To generate the correlation plot comparing the global changes for COVID-19 infection between different tissues, we utilized several different RNA-sequence data. The collection of the RNA sequencing data for the NP swab samples were described previously. The NP swab samples were analyzed with different methodologies with one comparing COVID-19 viral infection to the negative patients and the other a regression on continuous variables as a function of SARS-CoV-2 sequence amount. The viral comparison was previously described in Butler et al. and the DESeq2 was utilized to generate the differential expression data.

The monocyte COVID-19 RNA-Seq data, published under the accession GSE159678, was downloaded from SRA and gene expression was quantified using Salmon’s selective alignment approach. The RNA-Seq processing pipeline was implemented using pyrpipe (https://github.com/urmi-21/pyrpipe/tree/master/case_studies/Covid_RNA-Seq). Exploratory data analysis and differential expression analysis were performed using MetaOmGraph. From the differential expression analysis for each group the fold-change values for the genes were filtered with an adjusted p-value < 0.05. A correlation plot between the fold-change values for the significantly regulated genes for each comparison of the COVID-19 samples were plotted using R program corrplot v0.84.

Viral genome analysis

Total RNA-seq reads were classified against a custom, pan-kingdom reference using kraken2. Reads that mapped uniquely to SARS-CoV-2 were aligned to the Wuhan reference genome using bwa-mem. Alignments were deduplicated, assembled, and called for major variants (Variant Allele Frequency > 0.6) using IVAR with a minimum coverage of 5 reads per site.

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Antibodies
Immune Cell Profiling Panel (Core)	Nanostring Technologies, Inc	GMX-PROCONCT-HICP-12, Item 121300101, Lot# 0474026
10 Drug Target Panel	Nanostring Technologies, Inc	GMX-PROMODNCT-HIODT-12, Item 121300102, Lot# 0474029
Immune Activation Status Panel	Nanostring Technologies, Inc.	GMX-PROMODNCT-HIAS-12, Item 121300103, Lot# 0474032
Immune Cell Typing Panel	Nanostring Technologies, Inc	GMX-PROMODNCT-HICT-12, Item 121300104, Lot# 0474035
Cell Death Panel	Nanostring Technologies, Inc	GMX-PROMOD-NCTHCD-12, Lot# 0474050
MAPK Signaling Panel	Nanostring Technologies, Inc	GMX-PROMOD-NCTHMAPK-12, Lot# 0474047
Pl3K/AKT Signaling Panel	Nanostring Technologies, Inc	GMX-PROMOD-NCTHPl3K-12, Lot# 0474053
COVID-19 GeoMx-formatted Antibody Panel including (TMPRSS2, clone EPR3861; ACE2, clone EPR4436; Cathepsin L/V/K/H, clone EPR8011; DDX5, clone EPR7239; and SARS-CoV-2 spike glycoprotein, polyclonal)	Abcam	ab273594, Lot# GR3347471-1
GeoMx Solid Tumor TME Morphology Kit	Nanostring Technologies, Inc	GMX-PRO-MORPH-HST-12; Item 121300310
Alexa Fluor® 647 alpha-Smooth Muscle Actin Antibody, clone 1A4	Novus Bio	IC1420R
Biological samples
Autopsy tissues	Weill Cornell Medicine Department of Pathology	https://pathology.weill.cornell.edu/
Chemicals, peptides, and recombinant proteins
TRIzol	Invitrogen	Cat. #15596026
10% neutral buffered Formalin	Electron Microscopy Sciences	Cat. #15712
DNAse I	Zymo Research	Cat. #E1010
Critical commercial assays
Super-Script III Platinum SYBR Green One-Step qRT-PCR Kit	Invitrogen	Cat. #12594025
BD Univeral Viral Transport Media System	Becton, Dickinson and Company	Cat. #220526
QIAsymphony DSP Virus/Pathogen Mini Kit	Qiagen	Cat. #937036
NEBNext® rRNA Depletion Kit v2 (Human/Mouse/Rat) with RNA Sample Purification Beads	New England BioLabs	Cat. #E7405
NEBNext® Ultra™ II Directional RNA Library Prep Kit for Illumina	New England BioLabs	Cat. #E7760
TapeStation 2200	Agilent Technologies	Cat. #G2964AA
Kapa Biosystems Illumina library quantification kit	Roche	Cat. 07960140001
GeoMx DSP system	Nanostring Technologies, Inc	MAN-10088-03
Deposited data
Raw and analyzed RNA-seq data	This paper	dbGAP: accession #38851 and ID phs002258.v1.p1
Analyzed Nanostring GeoMx data	This paper	GEO: GSE169504
Human reference genome NCBI build 38, Gencode Human Release 33 (GRCH38.p13)	Genome Reference Consortium	http://www.ncbi.nlm.nih.gov/projects/genome/assembly/grc/human/
Raw RNA-seq data	Rother et al.44	GSE159678
Reference scRNA-seq data	Travaglini et al.45; MacParland et al.46; Stewart et al.47; Wang et al.8	https://www.humancellatlas.org/
Molecular Signatures for GSEA (MSigDB)	Liberzon et al.48; Subramanian et al.49; Kuleshov et al.50; Sergushichev et al., 2016	http://www.gsea-msigdb.org/gsea/
Oligonucleotides
Primers for RT-PCR; ACTB-Forward: CGTCACCAACTGGGACGACA	This paper	N/A
Primers for RT-PCR; ACTB-Reverse: CTTCTCGCGGTTGGCCTTGG	This paper	N/A
Primers for RT-PCR; SARS-CoV-2-TRS-L: CTCTTGTAGATCTGTTCTCTAAACGAAC	This paper	N/A
Primers for RT-PCR; SARS-CoV-2-TRS-N: GGTCCACCAAACGTAATGCG	This paper	N/A
Software and algorithms
ImageJ	Schneider et al., 2012	https://imagej.nih.gov/ij/
nf-core/rnaseq pipeline	Ewels et al.51	https://nf-co.re/rnaseq
FastQC	Andrews52	https://www.bioinformatics.babraham.ac.uk/projects/fastqc/
Trim Galore!	N/A	https://github.com/FelixKrueger/TrimGalore
STAR	Dobin et al.53	https://github.com/alexdobin/STAR
Salmon	Patro et al.54	https://salmon.readthedocs.io/en/latest/salmon.html
Picard	N/A	https://github.com/broadinstitute/picard
StringTie	Kovaka et al.55	https://ccb.jhu.edu/software/stringtie/
Samtools	Li and Durbin56	http://samtools.sourceforge.net/
DESeq2 R package	Love et al.57	https://bioconductor.org/packages/release/bioc/html/DESeq2.html
MuSiC R package	Wang et al.58	https://xuranw.github.io/MuSiC/articles/MuSiC.html
quanTIseq R package	Finotello et al.59	https://icbi.i-med.ac.at/software/quantiseq/doc/
Cocor R package	Diedenhofen and Musch60	https://cran.r-project.org/web/packages/cocor/cocor.pdf
synRNASeqNet R package	Luciano Garofano	https://github.com/cran/synRNASeqNet
Other
Resource page to visualize and explore autopsy RNA-seq data	This paper	https://covidgenes.weill.cornell.edu