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

1The coronavirus disease 2019 (COVID-19) pandemic is caused by the betacoronavirus severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2). 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 lesion. 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-19. 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 changes.

2There have been a number of studies that have examined the blood based immune response to SARS-CoV-2 infection–. 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 sites,. 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.

3A total of 20 patients at the Massachusetts General Hospital and 4 patients from Columbia University Irving Medical Center (NYC) who succumbed from SARS-CoV-2 infection underwent autopsy upon consent for clinical care. All patients were confirmed for SARS-CoV-2 infection through qRT-PCR assays performed on nasopharyngeal swab specimens. Clinical summaries of the 24 patients are listed in Supplementary Table 1. The mean age of this cohort was 62.5 years (range 32–89) with 14 males and 10 females. A total of 17 patients had medication records available with 5/17 patients on immunosuppressive medications, including 3 patients on corticosteroids. Most patients received hydroxychloroquine (13/17 = 76%), while none received remdesivir.

4Cases 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-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).

5Based on a quantitative assessment of the RNA-ISH (high viral cases defined as ratio of area of viral infected lung to total lung area ≥ 2%), 11 cases were classified as high viral RNA (Cases 1, 5, 7, 8, 9, 11, 15, 16, 18, C, D) while 12 cases were characterized as low viral RNA (Cases 2, 3, 4, 6, 10, 12, 13, 14, 17, 19, 20, B). High viral RNA correlated with extracellular viral RNA in a hyaline membrane pattern (Chi square p=0.001; Fig. 1b).

6The 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.

7To 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).

8Other 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.

9We 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 10/11, Fisher Exact p = 0.049). Although the limited number of cases precluded correlation with parameters associated with a cytokine syndrome,, 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.

10We next investigated gene expression patterns in separate lobes of the lung (left lower lobe - LLL, left upper lobe - LUL, right lower lobe - RLL, right middle lobe - RML, right upper lobe - RUL) from SARS-CoV-2 cases using total RNA-seq. Unsupervised hierarchical clustering of samples for the 500 most variant genes demonstrated a separation of high viral RNA samples from low viral RNA samples (Fig. 2a). Notably, there was separate clustering of high viral cases (Cases 1, 8, 9, 11, C, D) and low viral cases (Cases 3, 4, 6, 10, B), but interestingly, there was a separate cluster containing lung specimens from both high and low viral cases (Cases 2, 3, 5). Cases 2 and 5 both had intermediate levels of virus when compared to other cases. In addition, the left upper lobe in Case 7 was notable for having very high viral levels unlike the other 4 lobes that had low virus. Gene set enrichment analysis for genes that were different between the three clusters (colored boxed Fig. 2a) was done to understand the host response differences between these lung samples by viral load strata (Supplemental Data 2). The virus high cluster had high expression of interferon stimulated genes (ISGs) including IRF1, IFI44L, IFIT3, and antiviral genes (OAS3, ADAR) with significant enrichment of interferon (IFN) pathway gene sets including HALLMARK_INTERFERON_GAMMA_RESPONSE (FDR = 1.04E-24) and HALLMARK_INTERFERON_ALPHA_RESPONSE (FDR = 1.43 E-21). High virus samples also had enrichment of HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION (FDR = 3.13E-12) consistent with wound healing in these SARS-CoV-2 infected lungs. The mixed virus cluster noted high expression of multiple collagen genes (COL1A1, COL1A2, COL3A1, COL4A1, COL4A2, COL5A1, COL6A1, COL6A2, COL6A3) along with other genes enriched in HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION (FDR = 4.60 E-25) and HALLMARK_MYOGENESIS (FDR = 5.53 E-09), alluding to continued wound healing and early fibrosis. Both high and mixed virus clusters had genes elevated in HALLMARK_COAGULATION (High virus cluster FDR = 2.19 E-10; Mixed virus cluster FDR = 1.46 E-06) with notable increased expression of PECAM1 (CD31) and VWF, two well-known markers for endothelial cells in this gene set, indicating potential implications on pulmonary vasculature and coagulation. Low virus samples had notable elevated surfactant genes (SFTPB, SFTPC, SFTPA1, SFTPA2), mucins (MUC2, MUC3A, MUC4, MUC5AC, MUC16), and keratins (KRT4, KRT5, KRT6A, KRT13). These genes are enriched in epithelial cells and would suggest higher proportion of airway and alveolar lining cells, as opposed to inflammatory cells, in low viral lung specimens. Together, these findings indicate there is a spectrum of viral load not only between patients, but also within the same patient that correlates with different transcriptional profiles of the host response.

11We 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 were enriched for IFN response (HALLMARK_INTERFERON_GAMMA_RESPONSE FDR = 2.66 E-81; HALLMARK_INTERFERON_ALPHA_RESPONSE (FDR = 1.40 E-75) with 60 IFN gamma response genes (Fig. 2b). HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION (FDR = 3.90 E-14) and HALLMARK_COAGULATION (FDR = 2.54 E-08) were also enriched in high versus low cases (Supplementary Fig. 3).

12We next looked at factors shown to be potentially involved in SARS-CoV-2 infectivity including TMPRSS2, ACE2, and CD147 (BSG),. 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). TMPRSS2 had moderate expression in all lung samples (mean 27.2 RPM), and interestingly, was significantly lower in virus high compared to low samples (−1.72 Fold; FDR 0.0036). Recent single cell RNA-seq work in normal lungs has shown low expression of ACE2 and higher expression of TMPRSS2 in a subpopulation of pneumocytes, consistent with our findings. This difference may reflect selective destruction of TMPRSS2 positive pneumocytes in high relative to low viral infections and repopulation by pneumocytes. Analysis of immune cell and inflammatory markers by RNA-seq showed a predominance of macrophage markers in all cases compared to other immune cells with markedly high expression of CD163 (Fig. 3a and Supplementary Fig. 4). Analysis also showed markedly increased expression of multiple MHC Class I related genes (HLA-A, HLA-B, HLA-C, B2M, TAP1) in viral high versus low cases, which is consistent with high interferon activity.

13We 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).

14Given our observation that individual patients demonstrated anatomical heterogeneity with respect to histopathology and viral load, an in-depth assessment of the immune response within infected airspaces was undertaken. This was accomplished by the viral RNA directed sampling of paraffin embedded tissue sections using the GeoMxTM Digital Spatial Profiler (DSP) in 6 lobes from 5 patients. Quantitative analysis of RNA and protein in discreet regions of interest (ROIs) was performed separately on SARS-CoV-2 viral positive and negative areas as identified by RNA-ISH (Fig. 4a). Deconvolution of immune subsets using RNA expression profiles of each ROI noted significant spatial heterogeneity of the immune response across cases (Fig. 4b). Macrophages were more abundant in areas with SARS-CoV-2 virus present (red circles) compared to non-viral areas (blue circles) in 4 of the 6 tissues with detectable virus. However, Case 2-lung LUL with very low virus had regions of high macrophage infiltrate without detectable virus. Interestingly, Case 1-lung LLL was very different from the other specimens with high dendritic cells (DC), T-cells, NK cells, and to a lesser extent plasmacytoid DCs (pDCs) in SARS-CoV-2 virus areas without monocyte or macrophage infiltrates. This contrasts the other cases as well as Case 1-lung LUL, which highlights the disparate immune responses to SARS-CoV-2 even within the same patient.

15We then used an unbiased tSNE clustering of all ROIs and identified two distinct clusters (Fig. 4c). Cluster 1 included the majority of the RNA-ISH positive ROIs and a smaller proportion of ISH negative ROI, while Cluster 2 included the majority of RNA-ISH negative ROI, control SARS-CoV-2 negative lung samples, and a minority of RNA-ISH positive ROI. Differential expression analysis between clusters noted high enrichment of ISGs (ISG15, IRF7, IRF9) in Cluster 1. Other genes overexpressed in Cluster 1 were chemokine and cytokine related genes (CCL15, CCL18, TNF, IL11RA, XCR1), the immune checkpoint protein PD-1 (PDCD1), and cancer-related genes (SSX1 and MAGEA4). Consistent with the bulk RNA sequencing data, 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.

16To 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, (Fig. 5b) the majority were again related to IFN response including gene expression of ISGs (GBP1, IFITM1, IFI6) and MHC Class I (HLA-A, HLA-B, HLA-C, HLA-F, B2M, TAP1). Analysis of protein expression in these same samples (Fig. 5c) revealed elevated immune checkpoint proteins in virus high ROIs including CTLA4 and PD-L1 in Case 1 and 7 lung samples. IDO1 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.

17Finally, 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). Genes that are secreted or associated with extracellular signaling molecules including CXCL9, CXCL10, and IDO1 were expressed at much higher levels in SARS-CoV-2 positive ROIs compared to other genes related to MHC Class I (HLA-A, HLA-B, HLA-C, HLA-F, B2M) or genes related to antiviral cellular response (IFI6, IFIT3, IFITM1). The high levels of IDO1 identified in endothelial cells by IHC (Supplementary Fig. 6) and the known expression of CXCL9 and CXCL10 in endothelial cells suggests these differences may be related to local viral effects on vasculature, but other cell types including monocytes/macrophages can be the source of the expression of these genes. The spatial patterns in these IFN response genes indicate that there are differences in local host response to the virus compared to the more diffuse impact of IFN pathway activation in the lung. Altogether, there is clear intrapulmonary heterogeneity of SARS-CoV-2 infection with a non-uniform IFN response by the immune system with activation of immunoregulatory pathways (PD1/PD-L1, CTLA4, IDO1).

18This study supports two phases of disease evolution in patients with severe COVID-19 pneumonia: 1) high viral RNA with abundant extracellular virus, and a histologic picture of exudative diffuse alveolar damage, and 2) low (or undetectable) viral RNA, a mixed histologic picture, but dominated by an organizing form of diffuse alveolar damage. The high viral cases were associated with a shorter disease duration, high expression of IFN pathway genes, and activation of wound healing and endothelial genes indicative of tissue damage and initiation of fibrosis. Although these could represent two distinct disease phenotypes, we favor a model of an early phase of the disease characterized by acute cellular injury, high viral load, and a robust IFN response followed by a late form of disease associated with tissue organization on histology and expression profiling, viral clearance, and a waning IFN response.

19The IFN response is the primary early defense against viruses including SARS family of viruses. Viral detection by various pathogen recognition mechanisms stimulates the production of type I and type III IFN resulting in the expression of ISGs via the JAK-STAT signaling pathway. In preclinical models, SARS-CoV-2 viral infection of human intestinal epithelial cells elicits a robust IFN response that is efficient at controlling viral replication and de novo virus production. The current study demonstrates a more nuanced view of the expression of ISGs and ties their expression with the presence of virus and duration of disease; notably, we were able to document the association at the levels of individual air spaces. Collectively, the data suggest the IFN response in SARS-CoV-2 correlated with clearance of the virus from the lung. While the IFN gene signature characterize the patients with high viral load, there was no clear immune or cytokine signature in patients with a low virus and in a later phase of the disease. Instead, in low viral cases we find indication of resolution of wound healing, presence of M0 and M2 macrophages, and return of pulmonary epithelial cells.

20Our 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-CoV. 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 testing. However, the results contrast with recent efforts that show minimal amounts of IFNs in the peripheral blood of patients with severe COVID-19,. 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 hypothesis. 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.

21The 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 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 expression, a gene involved with paramyxovirus viral replication and suppression of the interferon response,. 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 macrophages and is a chemokine involved with chemotaxis of neutrophils, monocytes, and lymphocytes. The presence of higher M1 macrophages in SARS-CoV-2 high cases would indicate that CCL15 is important in the immune response to the virus. Our analysis of proteins enriched in high virus ROIs revealed enriched expression of multiple immune checkpoint molecules including CTLA4, PD-L1, and IDO1, which supports an immune microenvironment that is inhibitory to T-cell activation. Prior studies have reported T cell exhaustion or increased expression of inhibitory receptors on peripheral blood T cells of COVID-19 patients supporting these findings of T-cell suppression,,.

22Determining the phase of the disease in the lungs, specifically the peak and nadir, may have therapeutic implications. Given the clear interferon response to the virus in early infection, immunosuppressive therapy would potentially be most effective in the decelerating viral load phase of the disease. A recent preliminary report of dexamethasone use in COVID-19 patients 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 support. 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 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.

23In 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 response, and these findings will help inform the design of current and future interventional trials.

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

25There was no custom code or mathematical algorithm used for this work. All software for RNA-seq and NanoString GeoMX data analysis is described in the methods and all software will be provided if requested.

26Analysis 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 calculate the number of positive cells per square millimeter of tissue. Areas of pulmonary parenchyma with positive SARS-CoV-2 RNA ISH signal were annotated manually. The percentage of ISH signal was calculated as follows: area of lung parenchyma with positive signal/area of total lung parenchyma.

27All the IHC staining was performed on Leica Biosystem Bond III with Leica Bond Polymer Refine Detection (Catalog number DS9800)

28We performed total RNA-seq and quantitative reverse transcriptase polymerase chain reaction (qRT-PCR). RNA extraction from FFPE slides was done using the FormaPure Total nucleic acid extraction kit (C16675, Beckman Coulter) according to manufacturer instructions. Three 5 um thin tissue sections from areas devoid of acute inflammation were used per sample. SARS-CoV-2 RNA was detected in extracts from FFPE samples with qRT-PCR by following the guidelines of the Center for Disease Control and Prevention for the qRT-PCR diagnostic panel. A total of 1 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. 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 package was used for differential expression analysis between samples. Cell type deconvolution from gene expression was performed using CIBERSORTx 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.

29Autopsy tissues from COVID-19 infected patents were processed following the GeoMx DSP slide prep user manual (MAN-10087–04). Autopsy FFPE slides were baked in oven at 60°C for at least 1 hour, and then rehydrated and blocked by Nanostring block buffer for 1 hour. CD68–594 (Novus Bio, NBP2–34736AF647), CD45–647 (Novus Bio, NBP2–34527AF647), and PanCK-488 (eBioscience, 53-9003-82) were added on the sections along with the Nanostring protein 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.

30Autopsy 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 × 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 bead-to-sample ratio. Libraries were paired-end sequenced (2×75) 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.

31Per 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 above-background 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 log2-transformed normalized data.

32Probes 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.

33For each data point (gene × ROI), the expected background was estimated as the geometric mean of the negative control probes from the appropriate probe pool (CTA or COVID-19 spike-in) 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.

34Cell mixing proportions were estimated using the SpatialDecon R library, which performs mixture deconvolution using constrained log-normal regression. The algorithm was run using a 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 tumors. 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.

35Only 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.

36For 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.