Temporal and Spatial Heterogeneity of Host Response to SARS-CoV-2 Pulmonary Infection

The relationship of SARS-CoV-2 lung infection and severity of pulmonary disease is not fully understood. We analyzed autopsy specimens from 24 patients who succumbed to SARS-CoV-2 infection using a combination of different RNA and protein analytical platforms to characterize inter- and intra- patient heterogeneity of pulmonary virus infection. There was a spectrum of high and low virus cases that was associated with duration of disease and activation of interferon pathway genes. Using a digital spatial profiling platform, the virus corresponded to distinct spatial expression of interferon response genes and immune checkpoint genes demonstrating the intra-pulmonary heterogeneity of SARS-CoV-2 infection.


INTRODUCTION
The coronavirus disease 2019  pandemic is caused by the betacoronavirus severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2)1. Although there has been significant progress in understanding the factors involved with SARS-CoV-2 cellular infectivity, the relationship of SARS-CoV-2 lung infection and severity of pulmonary disease manifestations is not fully understood. Immune responses to viral infection have evolved to clear the pathogen, and differences in these responses amongst patients is likely to affect clinical outcomes.
Autopsy series have revealed that the predominant pattern of lung injury in COVID-19 patients is diffuse alveolar damage, characterized by hyaline membrane formation and in most cases, a presumed healing phase of this lesion2. However, these studies are limited in their ability to elucidate the complex immune response in SARS-CoV-2 pulmonary infection. An initial brief report of single cell RNA-seq analysis of bronchoalveolar lavage fluid from 9 patients noted an abundance of inflammatory monocyte-derived macrophages, lower CD8+ T cell infiltration, and high cytokine levels (IL-8, IL-6 and IL-1β) in patients with severe COVID-193. This suggested that macrophage driven responses and a "cytokine storm" were potentially preventing adequate T cell response to SARS-CoV-2 in patients with severe disease. Another study of 3 patients that focused on the peripheral blood response to SARS-CoV-2 found elevated IL-1 pathway cytokines and subsequent decrease in peripheral T cells, potentially linking intrapulmonary immune response with systemic changes4.
There have been a number of studies that have examined the blood based immune response to SARS-CoV-2 infection5-7. Tissue based examination has the potential to provide a more accurate assessment of SARS-CoV-2 related immune signatures, particularly if the immune cells are restricted to the affected organs. The ability to visualize SARS-CoV-2 at a tissue level provides unique information on the cell types infected by the virus and the spatial relationship of infected cells with immune and non-immune cells in the microenvironment. This provides a strategy to elucidate the roles of direct viral cytopathic effect and cellular injury from aberrant immune reaction, both within the lung and at extra-pulmonary sites8,9. Herein, we examined autopsy material from 24 COVID-19 patients collected at two institutions. The results demonstrate heterogeneous levels of SARS-CoV-2 RNA which correlate with duration of disease and show a range of host immune response patterns as well as significant spatial heterogeneity of both viral load and immune response. Cases were evaluated with RNA in situ hybridization (RNA-ISH) using a SARS-CoV-2 RNA specific probe targeting the S gene applied to multiple at least 2 different lobes of the lung and selected extrapulmonary organs. RNA-ISH positive cases noted intracellular staining detectable with a predominance in pneumocytes (Fig. 1a). Robust extracellular staining in hyaline membranes was detectable in 11 of 23 cases (Fig. 1a). One sample (Case A) failed by RNA-. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted August 2, 2020. . https://doi.org/10.1101/2020.07. 30.20165241 doi: medRxiv preprint ISH due to sample quality. Preservation of RNA quality was confirmed by GAPDH RNA-ISH. A similar pattern of reactivity was noted on an immunohistochemical stain for the SARS-CoV nucleocapsid protein ( Supplementary Fig. 1).
The histological analysis of the tissue sections used for viral RNA-ISH revealed a mixed picture with 5 cases showing only acute diffuse alveolar damage pattern of injury; in 9 patients the acute injury was accompanied by interstitial and/or airspace organization. Two cases lacked acute injury, instead showing only organizing-pattern injury. Low viral RNA correlated with evidence of interstitial/airway organization (p=0.049). Of note, the two patients with the purely organizing pattern of injury showed extremely low viral RNA (0% and 0.01%). Cases with low viral RNA also showed higher numbers of Napsin A positive cells (p=0.02) and trend towards higher keratin positive cells (p=0.08), supporting the repopulation of pneumocytes/bronchial epithelial cells. Collectively, early pattern of lung injury (exudative diffuse alveolar damage) correlated with high viral RNA while the presence of organization and intact pneumocytes correlated with low viral RNA.
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted August 2, 2020. . https://doi.org/10.1101/2020.07.30.20165241 doi: medRxiv preprint . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted August 2, 2020. . https://doi.org/10.1101/2020.07.30.20165241 doi: medRxiv preprint Fig. 1 Detection of SARS-CoV-2 in Human Autopsy Samples. a, Paraffin embedded sections from the lung of Case1 show abundant SARS-CoV-2 extracellular RNA-ISH signal (red) predominantly localization to the hyaline membranes (arrow). The inset shows the corresponding hematoxylin and eosin stained section shows histologic features of exudative diffuse alveolar damage with prominent hyaline membranes. b, Percentage of viral load in the lung as determined by a quantitative analysis of SARS-CoV-2 RNA-ISH. c, Expression heatmap of RNA-seq aligned counts of genes in the SARS-CoV-2 genome Log2(RPM) from autopsy cases. Consistent with the quantitative analysis on the RNA-ISH platform, cases 9, 1, C, 11, D and case 8 showed the highest viral load. The non-pulmonary organs were virtually negative for virus, except two bowel tissues. d, Swimmers plot highlighting difference in duration of illness between viral high and viral low cases. e, Quantitative protein expression of keratin and Napsin A by immunohistochemistry performed on lung sections between viral high and low cases. Viral low cases showed higher number of keratin and Napsin A positive cells, both surrogate markers of pulmonary pneumocytes. Boxand-whisker plot, center line, median; box limits, upper and lower quartiles; whiskers, range. Pvalue two tailed t-test.
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted August 2, 2020. . https://doi.org/10.1101/2020.07. 30.20165241 doi: medRxiv preprint To validate the RNA-ISH data, we performed molecular confirmation through quantitative RT-PCR (qRT-PCR)( Supplementary Fig. 2) and Total RNA-sequencing (RNA-seq) (Fig. 1c) for Cases 1-11, Cases A-D, and 5 non-COVID-19 autopsy lung specimens (Negative Control).
There was concordance of orthogonal techniques of viral detection, although as expected, qRT-PCR had the highest sensitivity of detecting SARS-CoV-2 in tissues (Supplementary Data1, Supplementary Table 3). RNA-seq also identified negative sense viral reads in most specimens with detectable virus indicating active viral replication in the lung (Supplementary Table 6).
Interestingly, antisense viral transcripts averaged 12.4% of total viral transcripts across all samples with detectable viral transcripts. The sense percentage was higher in the viral high cases, although this difference was not statically significant (mean 91.1% vs 81.6%; p=0.24).
Other organs examined in these cases including the heart, liver, jejunum, bowel, bone marrow, adipose tissue, skin and kidney were negative for detectable SARS-CoV-2 by RNA-ISH (Supplementary Data 1). Total RNA-seq detected viral RNA in two bowel samples (Cases 8 and 11). All other extrapulmonary samples were negative for SARS-CoV-2 by total RNA-seq. The absence of detectable virus by RNA-ISH suggests that viral RNA detected in non-pulmonary organs may represent viral RNA in blood, or that viral load is beneath the limit of detection by ISH.
We then evaluated if RNA-ISH viral load were associated with clinical parameters in our cases (Supplementary Table 1). Patients classified as high viral RNA showed a shorter duration of disease (mean 7.2 days, SD 3.02) than patients with low viral load (mean 19 days, SD 4.9) ( Fig.   1d; t-test p=0.0001). Additionally, patients classified as high viral RNA load had significantly fewer days in the hospital (mean 3.6 days, SD 2.2) than patients with low viral load (mean 14 days, SD 7.9) (t-test p=0.003). These differences could be because of a higher proportion of patients not receiving intubation and mechanical ventilation in the high viral RNA cases (4/9 vs . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted August 2, 2020. . https://doi.org/10.1101/2020.07.30.20165241 doi: medRxiv preprint 10/11, Fisher Exact p = 0.049). Although the limited number of cases precluded correlation with parameters associated with a cytokine syndrome,10,11 all parameters, except for absolute lymphocyte counts, showed higher numbers in low viral cases, although none of these were statistically significant (Supplementary Table 2). There were no correlations between viral load and age, gender or immunosuppression medication. Collectively, the data suggest a predictable pattern of viral infection manifested histologically by prominent diffuse alveolar damage with robust reactivity by RNA-ISH in the early phase and a shift towards organizing fibrosis with more sparse reactivity by RNA-ISH in the later phase. is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

Significant Innate Immune Response in High SARS-CoV-2 Infected Lungs
The copyright holder for this preprint this version posted August 2, 2020. We then analyzed RNA-seq data based on patient lung viral load and not individual samples as done with clustering. Using a cutoff of an average of 50 reads per million (RPM) of total viral gene expression across lung samples, we separated patients into high and low virus that was concordant with RNA-ISH results with the exception of Case 7, which again had one lung lobe (LUL) having high expression levels of virus (pink star) and 3 additional lobes (RLL, RML, LLL) that had detectable virus but below 50 RPM. Unbiased differential expression analysis (FDR < 0.01) identified 338 host genes that were higher in the high viral cases and 5710 genes higher in the low viral cases (Supplemental Data 3). The genes expressed higher in high viral cases . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted August 2, 2020.  Supplementary Fig. 3).
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted August 2, 2020. . https://doi.org/10.1101/2020.07.30.20165241 doi: medRxiv preprint . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted August 2, 2020. . https://doi.org/10.1101/2020.07.30.20165241 doi: medRxiv preprint Figure 2: Interferon gamma response genes define high SARS-CoV-2 cases a, Unsupervised hierarchical clustering of 500 most variant genes across lung specimens from SARS-CoV-2 infected patients separating into high, mixed, and low viral RNA cases. b, Gene expression heatmap of interferon gamma response genes differentially expressed between high and low viral RNA cases. The high viral cases were enriched with higher interferon response genes.
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted August 2, 2020. . https://doi.org/10.1101/2020.07.30.20165241 doi: medRxiv preprint We next looked at factors shown to be potentially involved in SARS-CoV-2 infectivity including TMPRSS2, ACE2, and CD147 (BSG)12,13. CD147/BSG had the highest mean expression (120.5 RPM) and ACE2 the lowest (4.2 RPM) with frequent undetectable levels in COVID-19 lung samples (Fig. 3a). We validated our immune cell RNA-seq analysis with IHC and quantification of immune cell types, which again showed very high numbers of CD163 positive cells in all samples ( Fig. 3b and c). There was a non-significant trend toward higher CD163 (p=0.07), CD3 (p=0.06), and CD4 (p=0.14) in low viral cases. To determine if there were differences in the type of CD163 cells, we performed deconvolution of RNA-seq data ( Fig. 3d and Supplementary Fig. 5), which showed lower cellular fraction of uncommitted (M0) macrophages (p = 0.01) in low viral cases and higher M1 polarized macrophages in high viral cases (p = 0.015). There was also a trend towards M2 polarization seen in low viral cases (p = 0.07).
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted August 2, 2020. . https://doi.org/10.1101/2020.07.30.20165241 doi: medRxiv preprint . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

Fig. 3 Monocytes and Macrophages Dominate SARS-CoV-2 Lung
Immune Response a, Expression heatmap of genes involved with SARS-CoV-2 viral entry, macrophage/monocytes, and MHC Class I. b, Immunohistochemical quantitation of immune cell subsets in lung samples for CD163, CD3, CD4, CD8, CD56, CD20, and CD138. Box-and-whisker plot, center line, median; box limits, upper and lower quartiles; whiskers, range. P-value two tailed t-test. c, Representative images of CD163 IHC staining in virus high and virus low case. d, Cell fraction (%) of immune cell estimated in lung tissue obtained by deconvolution of bulk RNA-seq data using CIBERSORTx. Bar = mean. P-value two tailed t-test.
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted August 2, 2020. is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted August 2, 2020. . https://doi.org/10.1101/2020.07.30.20165241 doi: medRxiv preprint unsupervised hierarchical clustering of ROIs using ISGs showed that SARS-CoV-2 positive ROIs were highly enriched for IFN response (Fig 4e). This demonstrates that the IFN response genes are preferentially expressed in regions of SARS-CoV-2 virus and not as a generalized inflammatory response in lung. Moreover, we hypothesize that SARS-CoV-2 negative ROIs in Cluster 1 represent resolving viral clearance and that Cluster 2 potentially represents never infected regions of lung or unresponsive lung regions with SARS-CoV-2 virus.
To determine if these changes are driven by specific samples, we performed differential expression of ROIs within each lung sample and identified consistent commonly expressed genes in each lung sample (Fig. 5a). Not surprisingly, the number of differentially expressed genes within lung samples was driven by the amount of SARS-CoV-2 heterogeneity within samples (i.e. Case 1 and 7 had more genes differentially expressed between SARS-CoV-2 positive vs. negative ROIs). Of the genes that were differentially expressed in multiple patients, protein was also elevated in these cases, which was also seen at the RNA level in both the GeoMx DSP analysis (Fig. 5b) and in bulk RNA-seq as an IFN gamma response gene (Fig. 2b).
Similarly, PD-L1 protein was elevated in Cases 1 and 7 (Fig. 5c) with PD-L1 RNA expression (CD274) elevated in Cluster 1 (Fig. 4e). IHC staining and whole slide quantification of PD-L1 and IDO1 in samples did not identify clear differences ( Supplementary Fig. 6), which illustrates the importance of analyzing ROIs to uncover intrapulmonary relationships between SARS-CoV-2 and host response. Notably, IDO1 staining was predominantly found in endothelial cells, which suggests a relationship between viral infection and changes in pulmonary vascular response.
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted August 2, 2020. . https://doi.org/10.1101/2020.07.30.20165241 doi: medRxiv preprint Finally, we attempted to dissect the intrapulmonary heterogeneity of the IFN response by looking at different classes of IFN genes that were co-expressed with each other (Fig. 5d). . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted August 2, 2020. . https://doi.org/10.1101/2020.07.30.20165241 doi: medRxiv preprint . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted August 2, 2020.  . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted August 2, 2020. . https://doi.org/10.1101/2020.07.30.20165241 doi: medRxiv preprint . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

Fig. 5 Spatial Distribution of Innate Immune Response Linked with Presence of SARS-CoV-2 Virus.
a, Differential expression of all genes between SARS-CoV-2 positive/negative regions within each tissue. Horizontal position shows genes' log2 fold-change, with points farther right having higher expression in SARS-CoV-2 positive regions. Vertical position shows -log10(p-value), which increases with statistical significance. Red points show genes that were consistently upregulated, blue points show genes that were consistently down-regulated in SARS-CoV-2 positive regions across the 6 patient tissues. Each gene has the same point color in all 6 panels. b, Genes with consistent differential expression between SARS-CoV-2 positive/negative regions across all tissues. Only consistently up/down-regulated genes (red/blue in panel a) are shown. Grid color shows log2 fold-change, with red indicating higher expression in virus-positive regions. Results with p < 0.01 are given color. Columns are ordered by hierarchical clustering. c, Protein differential expression between SARS-CoV-2 positive/negative regions. Only proteins with FDR < 0.05 in at least one tissue are shown. Grid color shows log2 fold-change, with red indicating higher expression in virus-positive regions. Results with p < 0.01 are given color. Columns are ordered by hierarchical clustering. d, Spatially-resolved expression of viral and interferon signaling genes. Rows show distinct gene sets; columns show distinct tissues. Pie position shows the physical location of regions within each tissue. Wedge volume shows a gene's background-subtracted expression; within each row, all genes are scaled to have the same maximum. Wedge color denotes whether a region was classified SARS-CoV-2 positive or negative by ISH.
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted August 2, 2020. . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted August 2, 2020. . https://doi.org/10.1101/2020.07. 30.20165241 doi: medRxiv preprint Our results are consistent with recent data that suggests SARS-CoV-2 appears to be responsive to IFN-I treatment in vitro and more sensitive than SARS-CoV18. Initial clinical data with IFN beta-1b, lopinavir-ritonavir, and ribavirin has shown some promise with decreased time to clearance of SARS-CoV-2 from nasopharyngeal swab testing19. However, the results contrast with recent efforts that show minimal amounts of IFNs in the peripheral blood of patients with severe COVID-195,6. These two studies have suggested that type-I IFN deficiency in the blood could be a hallmark of severe COVID-19, although neither study evaluated pulmonary tissue.
The presence of robust expression of numerous ISGs in bronchoalveolar lavage fluid from SARS-CoV-2 supports this hypothesis20. While the precise role of the IFN response and its ability to clear virus remains an open question, our work underscores the possibility of a robust IFN response in the majority of patients with severe COVID-19 infection and that these IFN producing cells may remain localized to the lung. The study also highlights the challenges associated with using blood as a surrogate to assess pulmonary disease in viral pneumonia, and this is further compounded by the temporal and spatial heterogeneity of the virus and corresponding immune response.
The use of the GeoMXTM digital spatial profiler has provided unprecedented spatial transcriptomic and proteomic analysis of the intra-pulmonary heterogeneity of SARS-CoV-2 infection. Within individual lobes, infected airspaces with significant interferon response pathway activation were juxtaposed to uninvolved lung tissue providing an opportunity to understand regional variation in lung tissue response to the virus. This intra-patient heterogeneity also extended to immune subsets, best exemplified in Case 1 in which a prominent population of dendritic cells, plasmacytoid dendritic cells, NK, and T-cells in one lobe was distinct from the more common monocyte/macrophage heavy infiltrate seen in another lobe of the same patient. Analysis of virus high ROIs across all cases showed significant interferon . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted August 2, 2020. . https://doi.org/10.1101/2020.07.30.20165241 doi: medRxiv preprint response gene expression compared to virus low ROIs, but we note there were two outlier genes SSX1 and CCL15 that were very high in virus high ROIs. SSX1 is known mostly as a partner of the SS18-SSX1 translocation commonly found in synovial sarcoma, and interestingly, has been found to be upstream of SHCBP1 expression21, a gene involved with paramyxovirus viral replication and suppression of the interferon response22, 23. This suggests that SSX1 may be induced by SARS-CoV-2 as an adaptive response to interferon, but further mechanistic studies will be needed to evaluate this hypothesis. CCL15, also known as leukotactin-1, is highly expressed in M1 compared to M2 macrophages24 and is a chemokine involved with demonstrated an improvement of 28-day mortality among those who were receiving either invasive mechanical ventilation or oxygen alone but not among those receiving no respiratory support28. Notably, dexamethasone was associated with a reduction in 28-day mortality among those with symptoms for more than 7 days but not among those with a more recent symptom onset, findings consistent with high viral load in the early phase of the disease. Our data also supports the likely greater benefit of antiviral therapy such as remdesivir in the early phase of . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted August 2, 2020. . https://doi.org/10.1101/2020.07.30.20165241 doi: medRxiv preprint the disease given high viral levels associated with diffuse alveolar damage. The data also raise the provocative hypothesis that many patients clear SARS-CoV-2 within two weeks of the onset of the disease, largely driven by a robust IFN response. However, the small, but detectable presence of virus RNA in low viral cases would indicate potential continued benefit of antivirals in later stages of the disease. The presence of T-cell inhibitory molecules (CTLA4, PD-L1, IDO1) in high virus cases also suggests that immune checkpoint inhibitors could have a role in patients with ineffective immune response, though this is highly speculative. The RNA-ISH and qRT-PCR assay could potentially be performed as complementary diagnostics on bronchoalveolar lavage samples to distinguish the two phases of the disease as well as discriminate active infection from asymptomatic carrier with another pulmonary process.
However, one notable caveat is the heterogeneity in viral load, necessitating that multiple lobes must be sampled.
In summary, detailed molecular analysis of multiple lung lobes from patients with severe COVID-19 infection highlights two phases in patients who succumb to the disease. While antiviral agents are likely to be most beneficial earlier in the disease, the timing and type of immune modulation, activating or inhibiting, must be carefully considered given the heterogeneous immune responses observed in these patients. Our findings highlight the importance in serially assessing tissue pulmonary SARS-CoV-2 RNA levels and the immune response as well as attempt to identify corresponding surrogate markers in the blood, although the heterogeneity in the viral load and immune response across the lobes of the lung may prove a challenge.
Additional mechanistic and biomarker driven studies of SARS-CoV-2 are needed to optimize patient selection and the timing of treatment administration to address this inherent spectrum of COVID-19 presentation. The current work provides the foundation to evaluate a larger series of autopsies characterizing the spatiotemporal relationship of viral load and host microenvironment . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted August 2, 2020. . https://doi.org/10.1101/2020.07.30.20165241 doi: medRxiv preprint response, and these findings will help inform the design of current and future interventional trials.

Data Availability Statement
All RNA-seq data has been uploaded to NCBI GEO GSE150316. All NanoString GeoMX data will be uploaded to NCBI GEO and released before publication.

Code Availability Statement
There was no custom code or mathematical algorithm used for this work. All software for RNAseq and NanoString GeoMX data analysis is described in the methods and all software will be provided if requested.

Analysis of patient autopsy material was reviewed and approved by the Partners Human
Research IRB (Protocol #: 2020P001001). Autopsy consent was per clinical care as directed by the patient or health care proxy. We evaluated hematoxylin and eosin (H&E) stained sections from FFPE tissue from the lungs, heart, liver, intestine, bone marrow, adipose tissue and kidney. We performed immunohistochemistry (IHC) for CD3, CD8, CD20, CD163, CD123, IDO1, PD-L1, Napsin A, keratin and SARS-CoV N protein on an immunohistochemical platform.
RNA-ISH was performed on FFPE sections using a SARS-CoV-2 RNA specific probe on an automated Leica BOND RX (RNAscope 2.5 LS Probe V-nCoV2019-S, #848568, and RNAscope 2.5 LS Reagent Kit-RED, #322150; Advanced Cell Diagnostics). Slides were imaged using a Leica Aperio CS-O slide scanning microscope at 40x magnification. Image quantification was performed using Halo software (Indica Labs). Tissue regions of interest were annotated by hand, excluding any folds or debris. The Multiplex IHC module was used to . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted August 2, 2020. is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted August 2, 2020. . https://doi.org/10.1101/2020.07.30.20165241 doi: medRxiv preprint ng of RNA input was used per reaction. 1-step qRT PCR was performed with the GoTaq Probe 1-Step RT-qPCR kit (A6120, Promega) using the CDC-approved 2019-nCoV RUO primer-probe kit (10006713, IDT). For Total RNA-sequencing, The Smarter Stranded Total RNA-Seq kit v2 (634413, Illumina) was used with 10 ng RNA input, according to the manufacturer instructions to generate libraries. Dual-indexed pooled libraries were sequenced on the Illumina NextSeq 500 platform using a 150 cycles kit with paired end read mode. For RNA expression analysis, the initial quality control of sequencing data was carried out using the tool FASTQC and alignment of sequencing reads to the reference genome was carried out using STAR aligner, version 2.7 29. We used the genome annotation and GTF for SARS-CoV-2 available on NCBI. A joint annotation was created by adding the COVID19 genome to the HG38 genome sequence and the GTF sequence. A new index for STAR aligner was created using this new annotation. Post alignment using this new annotation, the duplicate reads were marked using PICARD and removed using SAMtools. The resulting BAM files were used to quantify the read counts per gene using HTSeq-count program. The downstream analysis was carried out in the R statistical programming language including hierarchical clustering. The DESeq2 package30 was used for differential expression analysis between samples. Cell type deconvolution from gene expression was performed using CIBERSORTx31 using the LM22 signature matrix and batch correction for bulk sorted reference profiles. Plots were made using the heatmap.2 function in the gplots package in R. All RNA-seq data has been uploaded to NCBI GEO GSE150316.
CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted August 2, 2020. . https://doi.org/10.1101/2020.07.30.20165241 doi: medRxiv preprint cocktail for overnight incubation. The slides were washed and stained with Syto83 (ThermoFisher, S11364) on the next day. 20X fluorescent images were scanned after loading the slides to GeoMx machine. Regions of interest (ROIs) in Alveoli were selected in both COVID19 high and COVID19 low regions based on the COVID-19 ISH staining in the serial section. Oligos from antibodies were cleaved and collected into 96-well plates. Then these oligos were hybridized with NanoString barcodes overnight and read with an nCounter machine.
Digital accounts of each antibody in each ROI were generated for data analysis.

GeoMx DSP for CTA profiling
Autopsy FFPE tissues from COVID-19 infected patents were processed following the GeoMx DSP slide prep user manual (MAN-10087-04). Autopsy slides were baked in oven at 60°C for at least 1 hour, and then deparaffinized and hydrated by Leica Biosystems BOND RX. Proteinase K was added prior to the incubation of incubated with RNA probe mix (CTA and COVID-19 spike-in panel). After overnight incubation, slides were washed with buffer and stained with CD68-594 (Novus Bio, NBP2-34736AF647), CD45-647 (Novus Bio, NBP2-34527AF647), and PanCK-488 (eBioscience, 53-9003-82) and Syto83 (ThermoFisher, S11364) for 1 hour, and loaded to the GeoMx DSP machine to scan 20X fluorescent images. Regions of interest (ROIs) were placed by aligning to the ROIs placed during protein profiling. Oligos were cleaved and collected into 96-well plates. Oligos from each AOI was uniquely indexed using Illumina's i5 x i7 dual-indexing system. 4 uL of a GeoMx DSP sample was used in the PCR reaction. PCR reactions were purified with two rounds of AMPure XP beads (Beckman Coulter) at 1.2x beadto-sample ratio. Libraries were paired-end sequenced (2x75) on a NextSeq550 up to 400M total aligned reads. Fastq files were further process by DND system and raw and Q3 normalized counts of all CTA targets in each AOI were obtained through Nanostring Azorius program.
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted August 2, 2020. . https://doi.org/10.1101/2020.07.30.20165241 doi: medRxiv preprint

Analysis of GeoMx protein data
Per manufacturer's recommendations, the data was normalized by scaling to the negative control IgG probes, which reflect the rate at which antibodies bind to each region. The Ms IgG2a probe was excluded from this calculation due to poor concordance with the other IgGs. 6 ROIs were removed due to low signal. 10 proteins were excluded from analysis due to lack of abovebackground signal. For each tissue and each protein, differential expression vs. SARS-CoV-2 presence/absence was evaluated with an unpaired, heteroscedastic t-test of the protein's log2transformed normalized data.

GeoMx RNA data normalization
Probes for each gene were screened for outliers and their raw counts collapsed to a single estimate of counts using their geometric mean. Each ROI's data was scaled to have the same 75th percentile.

Estimation of background in GeoMx RNA data
For each data point (gene x ROI), the expected background was estimated as the geometric mean of the negative control probes from the appropriate probe pool (CTA or COVID-19 spikein) within the ROI in question. These expected background values were used as an input in mixed cell deconvolution and for background-subtraction in the plot of spatially-resolved expression.

Deconvolution of cell proportions from GeoMx RNA data
Cell mixing proportions were estimated using the SpatialDecon R library, which performs mixture deconvolution using constrained log-normal regression. The algorithm was run using a . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted August 2, 2020. . https://doi.org/10.1101/2020.07.30.20165241 doi: medRxiv preprint cell profile matrix derived from the Human Cell Atlas adult lung 10X dataset and appended with a neutrophil profile derived from snRNA-seq of lung tumors32. The neutrophil profile was scaled to have the same 75th percentile expression value as the average cell type's profile in the Human Cell Atlas lung dataset.

Differential Expression Analysis of GeoMx RNA data
Only genes that rose 2-fold above background in at least one ROI were considered. For each tissue and each gene, differential expression vs. SARS-CoV-2 presence/absence was evaluated with an unpaired, heteroscedastic t-test of the gene's log2-transformed normalized data. For analysis of expression within tissues, genes were defined as consistently up-regulated if they 1. had a log2 fold-change > 0.2 and a Benjamini-Hochberg FDR < 0.1 in at least 2 tissues, and 2. never had a log2 fold-change < 0 and a Benjamini-Hochberg FDR < 0.1 in any tissue. Genes were defined as consistently down-regulated by an equivalent rule.
For analysis across patient samples, expression was modeled using linear mixed effect models allowing for random slope and intercept terms per patient sample. P-values were estimated using Satterthwaite's method for approximation, and adjusting using the Benjamini-Hochberg FDR.

COMPETING INTERESTS
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted August 2, 2020. . https://doi.org/10.1101/2020.07.30.20165241 doi: medRxiv preprint