key: cord-0763091-it4pwxvk
authors: Park, Jiwoon; Foox, Jonathan; Hether, Tyler; Danko, David; Warren, Sarah; Kim, Youngmi; Reeves, Jason; Butler, Daniel J.; Mozsary, Christopher; Rosiene, Joel; Shaiber, Alon; Afshinnekoo, Ebrahim; MacKay, Matthew; Bram, Yaron; Chandar, Vasuretha; Geiger, Heather; Craney, Arryn; Velu, Priya; Melnick, Ari M.; Hajirasouliha, Iman; Beheshti, Afshin; Taylor, Deanne; Saravia-Butler, Amanda; Singh, Urminder; Wurtele, Eve Syrkin; Schisler, Jonathan; Fennessey, Samantha; Corvelo, André; Zody, Michael C.; Germer, Soren; Salvatore, Steven; Levy, Shawn; Wu, Shixiu; Tatonetti, Nicholas; Shapira, Sagi; Salvatore, Mirella; Loda, Massimo; Westblade, Lars F.; Cushing, Melissa; Rennert, Hanna; Kriegel, Alison J.; Elemento, Olivier; Imielinski, Marcin; Borczuk, Alain C.; Meydan, Cem; Schwartz, Robert E.; Mason, Christopher E.
title: Systemic Tissue and Cellular Disruption from SARS-CoV-2 Infection revealed in COVID-19 Autopsies and Spatial Omics Tissue Maps
date: 2021-03-09
journal: bioRxiv
DOI: 10.1101/2021.03.08.434433
sha: 6587e220a4dbaaf9940584a1f8b6e0a6ccebbd9b
doc_id: 763091
cord_uid: it4pwxvk

The Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2) virus has infected over 115 million people and caused over 2.5 million deaths worldwide. Yet, the molecular mechanisms underlying the clinical manifestations of COVID-19, as well as what distinguishes them from common seasonal influenza virus and other lung injury states such as Acute Respiratory Distress Syndrome (ARDS), remains poorly understood. To address these challenges, we combined transcriptional profiling of 646 clinical nasopharyngeal swabs and 39 patient autopsy tissues, matched with spatial protein and expression profiling (GeoMx) across 357 tissue sections. These results define both body-wide and tissue-specific (heart, liver, lung, kidney, and lymph nodes) damage wrought by the SARS-CoV-2 infection, evident as a function of varying viral load (high vs. low) during the course of infection and specific, transcriptional dysregulation in splicing isoforms, T cell receptor expression, and cellular expression states. In particular, cardiac and lung tissues revealed the largest degree of splicing isoform switching and cell expression state loss. Overall, these findings reveal a systemic disruption of cellular and transcriptional pathways from COVID-19 across all tissues, which can inform subsequent studies to combat the mortality of COVID-19, as well to better understand the molecular dynamics of lethal SARS-CoV-2 infection and other viruses.

In March 2020, the World Health Organization (WHO) declared a novel pandemic of the coronavirus disease 2019 (COVID- 19) , an infection caused by the betacoronavirus Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2) 1 , currently attributed to over 115 million cases and over 2.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 2,3 , 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 about the molecular drivers of morbidity and mortality associated with SARS-CoV-2 infection persist. This knowledge could lead to improvements in both the acute treatment and long-term management of pathological changes in multiple organs.

(alveolar epithelial cell markers) in normal lungs, CLU (lung injury and repair) and S100A9 (enriched in activated macrophages) in COVID-19 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 COVID-19 high lungs. Enrichment of interferon stimulated genes was only observed in COVID-19 high samples, which correspond to early stages of the infection, by timing of disease onset and histopathology (Extended Data Figure  1 ). Moreover, we observe enrichment of CASP3 and ID1, suggesting ongoing cellular injury and repair responses in COVID-19 high patients. In contrast, we find an enrichment of several markers of pulmonary fibrosis (e.g. CLU, COL1A1, COL1A2, and COL3A1) in COVID-19 low patients (which correspond to the later stages of infection), indicating that these are two distinct stages of infection 7 .

Next, we examined the spatial transcriptome data for differences between the COVID-19 high/low patients and the influenza lung samples. This analysis revealed a significant enrichment of THBS1 and NR4A1 in the influenza samples, both genes that have been shown to be engaged in response to influenza-induced lung injury (Extended Data Figure 2a) . When lung tissues from COVID-19 high and low samples were compared to those from non-viral induced acute respiratory distress syndrome (ARDS), enrichment for S100A8 and S100A9 was observed, which is consistent with the significant enrichment of neutrophils in these samples and previously suggested to be a driver of COVID-19 pathogenesis (Extended Data Figure  2b ) 8, 9 . Importantly, an enrichment for genes involved in lung injury and repair (ID1) as well as those involved in type 1 interferon responses including IFI6, IFI27, ISG15, and LY6E in COVID-19 high lungs is observed even when compared to influenza infection of the lung (Extended Data Figure 2a ) and non-viral lung injury (Extended Data Figure 2b) . Similarly, we find enrichment of several markers of pulmonary fibrosis (e.g. CLU, COL1A1, COL1A2, and COL3A1) in COVID-19 low samples as compared to non-viral and viral ARDS samples, exemplifying the profound lung injury and fibrosis during later stages of COVID-19 infection.

For each ROI, count estimates of 15 distinct cell types were imputed based on gene expression profiles from the Human Cell Atlas (HCA) adult lung dataset, including a neutrophil profile derived from snRNA-seq of lung tumors (see Methods). Consistent with other studies 10, 11 , we observed that COVID-19 was associated with an increase in tissue infiltrating immune cells, including T cells, NK cells, monocytes, and macrophages (Figure 2a) . Some immune cell types, such as monocytes, NK cells, and regulatory T cells, showed a statistically significant increase only in the COVID-19 high condition (Figure 2) . Also, while fibroblasts and endothelial cells increased in both COVID-19 high and low 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 (Extended Data Figure 3 ). Using an orthogonal approach, we stained the lung tissues with Masson's trichrome and observed a statistically significant increase in cellular collagen-rich areas, confirming the increase in lung fibroblasts (Extended Data Figure 4 ). As shown in the cluster analysis these cell type counts and proportions distinguish between the normal vs. COVID-19 (high and low) lungs regardless of lung structure origin (Figure 2b) , indicating that the SARS-CoV-2 infection is altering the cellular landscape and composition of the lung tissue.

Given the observed changes in cell proportions that are induced during SARS-CoV-2 infection, we next examined the impact on the cell-to-cell interaction landscape as a metric of intrapulmonary cellular heterogeneity (Figure 3) . Pairwise correlations of all detected cell types under five different conditions (COVID-19 high and low, flu, ARDS, and normal) were calculated and visualized (Figure 3a) . In normal lungs, we observed three clear cellular correlation clusters (1) monocytes, fibroblasts, T and NK cells, (2) neutrophils and ciliated cells, and (3) pDCs, macrophages, and B cells. While perturbations of these cellular interactions were observed across all injury conditions, the NK-T cell correlation was lost only in the COVID-19 high patients, and not present in COVID-19 low patients, which corresponds to the later stages of infection (Figure 3a-b) . We then quantified the correlation differences in the lungs' cellular landscape between COVID-19 high vs. COVID-19 low patients, which indicated that the greatest changes were in the monocyte-T correlations and dendritic-neutrophil correlations (Extended Data Figure 5a) , further supporting the view that SARS-CoV-2 specific T cell activities may be disrupted. Of note, it may be possible that the changes in correlation reflect both the generation of long-term memory as well as T cell mediated killing of infected epithelial cells 12 .

When next plotted the average proportions of each cell type across the COVID-19 disease states (COVID-19, Flu, ARDS). Macrophages and neutrophils populations were found to be much higher in the COVID-19 lungs (both high and low), while T, monocyte, and epithelial cells were much lower compared to normal. The entropy estimates of the given cell types (Figure 3c) showed that macrophages and neutrophils commonly displayed increased heterogeneity (across ROIs and patients) across all injury conditions. Variances of fibroblast, epithelial, vessel, and T cell populations were bigger in COVID-19 high and low lung tissues compared to flu infection and ARDS, suggesting an increase in the cellular heterogeneity. When comparing the entropies by tissue types, nearly all COVID-19 positive ROIs showed increase in heterogeneity of the cell populations, with the sole exception being vascular regions with high COVID-19 (Figure 3d) . The decrease in entropy in the vascular regions of COVID-19 high condition is mainly from the decrease of fibroblasts, epithelial, T, and NK cells (Extended Data Figure 5b ). 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 (Extended Data Figure 5b ). Single-sample Gene Set Enrichment Analysis (ssGSEA) of macrophage-, neutrophil-, and T cell regulatory pathways showed enrichment in COVID-19, even when compared with flu and ARDS, including macrophage activation and apoptotic process (1.302 and 2.809 fold increase in averaged ssGSEA scores relative to normal, with p-values of 2.91e-08 and 0.001) in COVID-19 high. (Extended Data Figure 6a -c).

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 (Extended Data Figure 7a -b) but mostly in the lung. Viral reads were robustly detected in nasopharyngeal samples from COVID-19 patients, consistent with the published reports 13 . Normalized coverage values in SARS-CoV-2 positive autopsy tissue samples revealed detection bias towards the SARS-CoV-2 3' end consistent with the known life cycle of the virus (Extended Data Figure 7a-b) 13, 14 . Reconstruction of the viral genomes revealed known and unknown variants common to many patients (Extended Data Figure 7c) , and some evidence of intra-host variability.

We next used shotgun metatranscriptomics (total RNA-seq) for host and viral profiling on 39 patients that died from COVID-19, including 190 organ-specific tissue samples from the respective autopsies and healthy controls from organ donor remnant tissues (n=3). 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 high, COVID-19 low, and control sample sets (q-value <0.01 expression fold-change >1.5-fold, DESeq2, Figure 4a-d) . Pathway enrichment analysis The DEGs and GSEA results were then examined for the largest differences between infection level and stage (COVID-19 high, early infection vs. COVID-19 low, late 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 4) . To place these results into greater context, we compared DEGs from each tissue to RNA-seq data from nasopharyngeal (NP) swabs, previously described in Butler et al 4 , as well as RNA-seq data from a publicly available dataset on monocytes from COVID-19 positive and negative patients (Extended Data Figure 10 ). The majority of the tissues with high viral load were significantly positively correlated with the DEGs in the NP swab samples (q-value <0.01), compared to normal/negative patients. However, in contrast, the later infection (low viral load) patients' tissues showed 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.

Despite significant transcriptional differences in the lung, there were no clear gross pathologic differences in bulk tissue or cell surface proteins (Extended Data Figure 8a, 8b) . Thus, to create a more fine-grained analysis of the cellular gene expression state in each tissue, we used the cell deconvolution algorithm MuSiC on each tissue's RNA-seq data. The MuSiC results showed distinct disruptions of the transcriptional programs for each tissue in the COVID-19 patients and gain or loss of cell types (Figure 4b) . Specifically, the kidney and liver showed a loss of proximal tubule in the kidney and hepatocyte marker expressions in the liver (Extended Data Figure 4a, 10) , but an increase in T cells found in the kidney and liver. The lung showed a loss of the capillary intermediate cells and alveolar epithelial cells. Finally, the heart showed a striking, near-complete loss of the cell signatures for cardiomyocytes (Figure 4c) , in both COVID-19 high and low samples (Extended Data Figure 4a, 9) , despite no obvious gross or histologic changes in the heart (Extended Data Figure 8b ).

Each tissue was examined for transcript switches across nine isoform disruption categories, such as loss of function, disrupted open reading frames, and nonsense mediated decay (FDR <0.05, Figure 5a ). Splicing alterations were most pronounced in the heart and lung tissues (7/9 isoform types significantly altered), but less significantly in the lymph nodes (2/9) and liver (0/9) again highlighting a particularly disruptive set of transcriptional changes in the heart and lung. Next, alternative splicing was analyzed within each tissue (Figure 5a ) from the ISAR data, showing significant exon skipping in the lung and increased exon skipping events in the heart, and more alternate start and stop sites within the heart, lung, and lymph nodes (FDR<0.05). No significant isoform switching or alternative splicing was detected in the kidney, despite significant disruption in expression profiles.

To put these differences into broader context, we then examined the intersection of specific splicing-disrupted genes across tissues and pathways (Figure 5c,d) . The greatest number of genes with consequential isoform switches (significant genes from Figure 5a ,b) was found in the lymph node (1353) and heart (1029), followed by lung (231) and liver (27) . Most splicedisrupted genes were unique to a given tissue, but the greatest overlap was observed between the heart and the lymph node (162 genes). Next, pathway enrichment was calculated with PathwaySplice 19 , using canonical pathways from the Reactome database (Figure 5d ). In the tissues with significant isoform switches, the greatest number of significantly disrupted pathways (adjusted p-value < 0.05) occurred in the heart (97), followed by the lymph node (72) and lung (65) , indicating a stronger pathway disruption effect in the heart than the lymph node and lung. Of note, the top dysregulated pathways in these tissues were related to viral responses, including the export of viral ribonucleoproteins and antiviral mechanism by interferon-stimulated genes, reflecting a broad, but tissue-specific, transcriptional response to viral infection (Extended Data Figure 10 ).

Finally, to examine the impact that COVID-19 has on T cell receptors in the patients that could be found in each tissue, we assembled CDR3 sequences from the patient-specific RNA-Seq data and derived a set of T Cell Expressed Motifs (TCEMs); TCEMs are amino acid sequences that can be recognized as an antigen by a T Cell when presented by an MHC. We mapped these TCEMs to the SARS-CoV-2 reference genome, to quantify the likely response of the T cells to the virus (Figure 6) . When compared the abundance of TCEMs that mapped to the SARS-CoV-2 genome, the abundance of TCEMs that mapped to SARS-CoV-2 was significantly higher in cases than in controls (3.7x increase using one-sided T-test, p-value of 1.24e-12, Figure 6a and Supplementary Table 3) .

We next considered variable (non-uniform) TCEM abundance across different tissue types. By taking the average abundance of TCEMs that mapped to each gene in the SARS-CoV-2 genome across all patients for each tissue type (lung, liver, heart, kidney, and lymph nodes, Figure 6b ), we found that TCEMs matching SARS-CoV-2 are differentially abundant by tissue types. The majority (2592/3174, 81.7%) of unique TCEMs were found in samples from lymph nodes, consistent with their importance in immune function. Despite the largest number of unique TCEMs, those from the lymph nodes were typically low abundance compared to other tissues, which may be evidence of T Cell maturation before effective clones are found and distributed across the body. By contrast TCEMs from heart samples were highly abundant and less diverse, indicating a possible targeted response, which also provides further evidence of the distinct cardiovascular COVID-19 phenotype of heart damage.

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, influenza, and ARDS. First, we found that the lung and other tissue samples demarcated COVID-19 patients into two groups, based on the amount of SARS-CoV-2 RNA present: COVID-19 high versus COVID-19 low. Lung samples examined using two techniques (Nanostring nCounter or GeoMx DSP) both support the correlation that high SARS-CoV-2 viral loads correlate with early COVID-19 disease while low SARS-CoV-2 viral loads correlate with late COVID-19 disease. The nature of the COVID-19 disease and the impact on the lung (and other tissues) appears to be defined principally by the presence of SARS-CoV-2 viral RNA.

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 enriched 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 to other viral related diseases (influenza), COVID-19 tissue samples still show enrichment for genes involved in lung injury and repair and interferon signaling genes, but the COVID-19 tissue samples also show enrichment of several markers of pulmonary fibrosis. Of note, COVID-19 (high), flu, and ARDS each show differential HLA-B and -C expression, which are known mediators of natural killer and T cell activation 20, 21 and which can mediate host risk of infection 4 . Enrichment for HLA-DRB5, whose expression and specific gene polymorphisms are associated with pulmonary fibrosis and severity 22, 23 . Comparing across different disease types, all diseases-COVID-19 (low), influenza, and ARDS showed enrichment for DMBT1, a gene known to be upregulated and dysregulated in pulmonary injury and fibrosis 24, 25 . Virus-related diseases (COVID-19 high and flu) particularly showed significant change of expression in lung epithelial cell-related transcripts (i.e. ACTB, C1R, and FN1), and such changes are known markers of the lung injury gene signature 26 .

The spatial analysis platform (GeoMx) enabled a novel analysis of the impact of the disease across entire tissues. Consistent with recent reports from bulk cellular profiling, we observed an increase in the immune cell types and fibroblasts in COVID-19, high but a decrease in alveolar epithelial cells 27 . In COVID-19-low condition, the proportions of some immune cells (i.e. monocytes, NK cells, or regulatory T cells) were similar to normal, but the fibroblasts and vessel cells still exhibited an increase in COVID-19. Some of these cell types form a "cellular correlation cluster," and these clusters of cell-cell interactions are uniquely disrupted in COVID-19 (relative to flu and ARDS), particularly in COVID-19 high. While macrophages and neutrophils showed an increase in entropy across all lung-related injury conditions, NK and T cells only showed an increase in COVID-19 conditions. 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 [28] [29] [30] [31] . 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 omics 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, the fibroblasts and immune cell types such as macrophages increased in most tissues, while alveolar epithelial types 1 and 2 cells in the lung showed a decrease in COVID-19 (relative to controls). Interestingly, other organ types such as heart or kidney also showed similar trends in COVID-19. Increase in fibroblasts, endothelial, and immune cells may be impacted by a variety of immune activations within each organ, particularly as s long-term response to the infection. This observation may also be related to the decrease in the 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 cardiomyocyte cell fraction and disrupted splicing isoforms within heart tissues was concomitant with a reduction in several transcripts encoding sarcomeric and contractile proteins, regardless of viral level 32,33 , representing a persistent transcriptional perturbation and cardiac-specific impact of COVID-19.

These data support the view of distinct 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 high correlation to tissue specific DEGs in the early stages of infection, but very little correlation later infection. To explain this phenotype, recent evidence has suggested that monocytes can migrate into tissues as the SARS-CoV-2 infection progresses 15 , and this response can potentially explain the different correlations. It has also been reported that monocyte depletion/migration is associated with kidney disease, inducing lupus like symptoms 16 , 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 responsible for focal necrosis seen in the lymph nodes 17 . Nonetheless, the lung, heart, and lymph nodes were the tissues most disrupted by infection.

In addition to our observations, there are several papers indicating the connections of cytokine storms with macrophage activation syndromes and/or other imbalance of the immune system environment 34, 35 , which also can be variable across tissues. To explore these connections, we examined the relative abundance in TCEMs mapped to the SARS-CoV-2 genome in our RNAseq results. Organ-wide TCEM mapping also revealed a wide variety of TCEM sites in heart, liver, or lung and showed high similarity within the body sites or patients. Indeed, the ssGSEA from the lung specific dataset of COVID-19 on the T cell mediated pathways revealed distinct from that of ARDS, flu, and the normal samples (Extended Data Figure 6c ). The expression signature was particularly prominent for T cell proliferation, chemotaxis, and activation, and this impact was most evident in the lymph node and the heart.

Overall, these data represent one of the largest autopsy series of COVID-19 and synthesizes several orthogonal methods, including: bulk transcriptomics, digital spatial transcriptomics, multiple imaging technologies, and novel computational analysis (TCEMs, ISAR) to build a map of SARS-CoV-2 pathophysiology. The vast majority of COVID-19 data previously published has been from NP swab or peripheral blood, and here we provide evidence that these data cannot capture the key, tissue-specific impacts of the disease, especially in late infection. Moreover, given the ability to also use these data to map the distinct viral variants as the pathogen moves through the body, our data also provide additional support of intra-host viral diversity for SARS-CoV-2, including variations specific to distinct tissues 36 , which can help guide future studies and treatments. This molecular map of COVID-19 serves as an atlas for the community and can inform future studies into COVID-19 progression and SARS-CoV-2 pathology. also further funding was provided by KBR, Inc. Sequencing of some samples was performed at the New York Genome Center (NYGC) as part of the COVID-19 Genomic Research Network (CGRN) with funds generously provided by NYGC donors. We would like to thank the Epigenomics Core Facility at Weill Cornell Medicine, the Scientific Computing Unit (SCU), as well as the Starr Cancer Consortium: I9-A9-071, I13-0052, the NIH: R21AI129851,  R01MH117406, R01CA249054 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 Supplementary Table  2 .

Patient specimens were processed as described in Butler et al., 2020 4 . 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 (10ng), 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.

RNAseq data was processed through the nf-core/rnaseq pipeline 41 . This workflow involved quality control of the reads with FastQC 42 , adapter trimming using Trim Galore! (https://github.com/FelixKrueger/TrimGalore), read alignment with STAR 43 , gene quantification with Salmon 44 , duplicate read marking with Picard MarkDuplicates (https://github.com/broadinstitute/picard), and transcript quantification with StringTie 45 . 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 46 . Cell deconvolution was performed using MuSiC on single cell reference datasets for lung, liver, kidney, and heart [47] [48] [49] [50] [51] . Immune cell deconvolution was performed on lymph node samples using quanTIseq 52 . 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 p-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 53 . 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) [54] [55] [56] [57] . Any signature with adjusted p-value < 0.01 was taken as significant. List of differentially expressed genes and significantly enriched pathways are reported in Supplementary Table 1 .

Differential isoform usage was estimated with IsoformSwitchAnalyzeR 58 . Briefly, Salmon isoform count matrices from every sample were imported using importRdata 59 , using Gencode v33 exon annotations and nucleotide sequences. Single-isoform genes were filtered, as well as any genes below an expression cutoff of 1. Isoform switch testing was performed using isoformSwitchTestDEXSeq with default parameters 60, 61 . All tissue types were tested separately, with COVID versus Control as the experimental condition. Differentially spliced genes were annotated as such: coding potential was assigned using CPAT 62 , protein domains were predicted with Pfam 63 , signal peptides were predicted with SignalP 64 , and intrinsically disordered regions were predicted using IUPred2a 65 , all using default parameters. Final isoform switch lists were generated using a differential isoform fraction cutoff of 0.1, a coding cutoff of 0.725, and non-coding ORFs removed. Gains and losses of predicted splicing events were predicted with extractSwitchingSummary and the functional significance of events was predicted with extractConsequenceSummary. Finally, differentially spliced pathways were estimated using PathwaySplice 19 , a gene feature-aware pathway analysis tool, using no GO size limitation, the Wallenius method of hypergeometric bias adjustment, and an exon-counting bin size of 20. Pathways were annotated against the Reactome Pathway Database 66 .

We used the method described by Danko et al 67 to identify potential T Cell Expressed Motifs (TCEMs) in transcriptomic data from this study. Briefly, that method includes assembling IgH sequences using MiXCR 68 , then selecting 5 amino acid k-mers from the CDR3 regions of these sequences using three specific patterns specified by Bremel et al 69 . The whole set of 5aa sequences, with abundances, is the TCEM repertoire for a sample. We normalized the abundances of our TCEM repertoires by the total number of TCEMs of each pattern in a sample to yield relative abundances.

We mapped the TCEM repertoire from each sample to canonical SARS-CoV-2 protein sequences. We took every 5 amino acid subsequence from each protein and matched these sequences to the 5aa sequences in our TCEM repertoire. Any 5aa sequence which did not uniquely identify a single SARS-CoV-2 protein was discarded. We mapped each remaining 5aa sequence to a position on the SARS-CoV-2 genome by offsetting the position within the protein to the position by the position of the protein on the genome. Some 5aa sequences occurred multiple times in the same protein. For these sequences we arbitrarily used the highest (most 3') coordinate. normal tissues are shown for the five tissues from autopsy: heart, kidney, liver, lung, lymph node. Differentially expressed genes (>1.5-fold, q-values <0.05, DESeq2) are shown as purple spots (down-regulation) and orange (up-regulation). e. Select pathways that show significant differences in different tissues are shown, with the statistical significance shown as a color range (legend, GSEA permutation test and adaptive multilevel splitting Monte Carlo method), and non-significant differences shown in grey. f. The pathways that show significant differences in all or majority of the five tissues are shown, with statistical significance from GSEA testing across five tissues (color in legend) and x-axis for the normalized enrichment score for COVID-19 positive vs. Control, COVID-19 high or low vs. Control, and COVID-19 high vs. COVID-19 low comparisons. Extended Data Figure 9 . All dysregulated pathways across tissue types. Pathways found to be dysregulated based on isoform switch analysis are shown for each tissue type examined. Only pathways in which at least one tissue had log fold change > 1 and adjusted p-value < 0.05 were considered. Odds ratio is shown per pathway, and pathways are colored by adjusted pvalues (with no shading for non-significant effects). Control COVID Low COVID High 

• • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• •• • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •

Temporal dynamics in viral shedding and transmissibility of COVID-19

Clinical Characteristics of Coronavirus Disease 2019 in China

Clinical course and risk factors for mortality of adult inpatients with COVID-19 in Wuhan, China: a retrospective cohort study

Shotgun Transcriptome and Isothermal Profiling of SARS-CoV-2 Infection Reveals Unique Host Responses, Viral Diversification, and Drug Interactions

Co-infections in people with COVID-19: a systematic review and meta-analysis

Activation and evasion of type I interferon responses by SARS-CoV-2

Temporal and spatial heterogeneity of host response to SARS-CoV-2 pulmonary infection

Induction of alarmin S100A8/A9 mediates activation of aberrant neutrophils in the pathogenesis of COVID-19

Elevated Calprotectin and Abnormal Myeloid Cell Subsets Discriminate Severe from Mild COVID-19

A distinct innate immune signature marks progression from mild to severe COVID-19

Single-cell landscape of bronchoalveolar immune cells in patients with COVID-19

T cell-mediated immune response to respiratory coronaviruses

Detection of SARS-CoV-2 in Different Types of Clinical Specimens

The Architecture of SARS-CoV-2 Transcriptome

Monocytes and macrophages, targets of SARS-CoV-2: the clue for Covid-19 immunoparalysis

Patrolling monocytes promote the pathogenesis of early lupus-like glomerulonephritis

Pathological inflammation in patients with COVID-19: a key role for monocytes and macrophages

The Landscape of Isoform Switches in Human Cancers

PathwaySplice: an R package for unbiased pathway analysis of alternative splicing in RNA-Seq data

HLA-C as a mediator of natural killer and T-cell activation: spectator or key player?

Differential expression of human major histocompatibility class I loci

Wholeexome sequencing identifies susceptibility genes and pathways for idiopathic pulmonary fibrosis in the Chinese population

Genome-wide association meta-analysis identifies five modifier loci of lung disease severity in cystic fibrosis

Single-Cell Transcriptomic Analysis of Human Lung Provides Insights into the Pathobiology of Pulmonary Fibrosis

Quantitative proteomic characterization of lung tissue in idiopathic pulmonary fibrosis

Mechanical forces induce an asthma gene signature in healthy airway epithelial cells. Sci Rep

The spatio-temporal landscape of lung pathology in SARS-CoV-2 infection

Sex differences in immune responses that underlie COVID-19 disease outcomes

Longitudinal analyses reveal immunological misfiring in severe COVID-19

Antigen-Specific Adaptive Immunity to SARS-CoV-2 in Acute COVID-19 and Associations with Age and Disease Severity

Profiling of immune dysfunction in COVID-19 patients allows early prediction of disease progression

Pathological features of COVID-19-associated myocardial injury: a multicentre cardiovascular pathology study

The cytokine storm and COVID-19

COVID-19 cytokine storm: the interplay between inflammation and coagulation

Hidden genomic diversity of SARS-CoV-2: implications for qRT-PCR diagnostics and transmission

Multiplex digital spatial profiling of proteins and RNA in fixed tissue

SpatialDecon: Deconvolution of mixed cells from spatial and/or bulk gene expression data

Fast gene set enrichment analysis

cocor: A Comprehensive Solution for the Statistical Comparison of Correlations

The nf-core framework for community-curated bioinformatics pipelines

FastQC: A Quality Control Tool for High Throughput Sequence Data

STAR: ultrafast universal RNA-seq aligner

Salmon provides fast and bias-aware quantification of transcript expression

Pertea M Transcriptome assembly from long-read RNA-seq alignments with StringTie2

Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2

Bulk tissue cell type deconvolution with multi-subject single-cell expression reference

A molecular cell atlas of the human lung from single-cell RNA sequencing

Single cell RNA sequencing of human liver reveals distinct intrahepatic macrophage populations

Spatiotemporal immune zonation of the human kidney

Single-cell reconstruction of the adult human heart during heart failure and recovery reveals the cellular landscape underlying cardiac function

Molecular and pharmacological modulators of the tumor immune contexture revealed by deconvolution of RNA-seq data

Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing

Molecular signatures database (MSigDB) 3.0. Bioinformatics

Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles

Enrichr: a comprehensive gene set enrichment analysis web server 2016 update

An algorithm for fast preranked gene set enrichment analysis using cumulative statistic calculation

IsoformSwitchAnalyzeR: analysis of changes in genomewide patterns of alternative splicing and its functional consequences

Differential analyses for RNA-seq: transcript-level estimates improve gene-level inferences

Detecting differential usage of exons from RNA-seq data

limma powers differential expression analyses for RNA-sequencing and microarray studies

Coding-Potential Assessment Tool using an alignment-free logistic regression model

The Pfam protein families database

SignalP 5.0 improves signal peptide predictions using deep neural networks

IUPred2A: context-dependent prediction of protein disorder as a function of redox state and protein binding

We thank the patients, their families, and healthcare workers fighting the COVID-19 pandemic.