key: cord-0858376-aynkjngy authors: Kulasinghe, Arutha; Tan, Chin Wee; dos Santos Miggiolaro, Anna Flavia Ribeiro; Monkman, James; SadeghiRad, Habib; Bhuva, Dharmesh D.; da Silva Motta Junior, Jarbas; Vaz de Paula, Caroline Busatta; Nagashima, Seigo; Baena, Cristina Pellegrino; Souza-Fonseca-Guimaraes, Paulo; de Noronha, Lucia; McCulloch, Timothy; Rodrigues Rossi, Gustavo; Cooper, Caroline; Tang, Benjamin; Short, Kirsty R.; Davis, Melissa J.; Souza-Fonseca-Guimaraes, Fernando; Belz, Gabrielle T.; O'Byrne, Ken title: Profiling of lung SARS-CoV-2 and influenza virus infection dissects virus-specific host responses and gene signatures date: 2021-10-21 journal: Eur Respir J DOI: 10.1183/13993003.01881-2021 sha: 90b530f38e9b072a3e8d9a692bc28f4bf2eddb8e doc_id: 858376 cord_uid: aynkjngy The severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) that emerged in late 2019 has spread globally, causing a pandemic of respiratory illness designated coronavirus disease 2019 (COVID-19). A better definition of the pulmonary host response to SARS-CoV-2 infection is required to understand viral pathogenesis and to validate putative COVID-19 biomarkers that have been proposed in clinical studies. Here, we use targeted transcriptomics of FFPE tissue using the Nanostring GeoMX™ platform to generate an in-depth picture of the pulmonary transcriptional landscape of COVID-19, pandemic H1N1 influenza and uninfected control patients. Host transcriptomics showed a significant upregulation of genes associated with inflammation, type I interferon production, coagulation and angiogenesis in the lungs of COVID-19 patients compared to non-infected controls. SARS-CoV-2 was non-uniformly distributed in lungs (emphasising the advantages of spatial transcriptomics) with the areas of high viral load associated with an increased type I interferon response. Once the dominant cell type present in the sample, within patient correlations and patient-patient variation had been controlled for, only a very limited number of genes were differentially expressed between the lungs of fatal influenza and COVID-19 patients. Strikingly, the interferon-associated gene IFI27, previously identified as a useful blood biomarker to differentiate bacterial and viral lung infections, was significantly upregulated in the lungs of COVID-19 patients compared to patients with influenza. Collectively, these data demonstrate that spatial transcriptomics is a powerful tool to identify novel gene signatures within tissues, offering new insights into the pathogenesis of SARS-COV-2 to aid in patient triage and treatment. Since its emergence in 2019, severe acute respiratory syndrome-associated coronavirus-2 (SARS-CoV-2) has caused a broad spectrum of disease, ranging from asymptomatic infections to acute respiratory distress syndrome (ARDS). Despite a lower fatality rate than other related coronaviruses, such as SARS-CoV-1 and middle east respiratory syndrome coronavirus (MERS-CoV), SARS-CoV-2 is highly transmissible and to date has resulted in over >175M infections and >3.5M deaths worldwide and rising [1] . COVID-19 is the clinical manifestation of SARS-CoV-2 infection. ARDS develops in 42% of COVID-19 patients presenting with pneumonia and accounts for a significant number of deaths associated with COVID-19 [2, 3] . ARDS is a form of hypoxemic respiratory failure defined by the presence of diffuse alveolar damage (DAD) commonly associated with bacterial pneumonia, sepsis, pancreatitis or trauma. The pathogenesis of ARDS is typically attributed to inflammatory injury to the alveolar-capillary membrane which leads to pulmonary oedema and respiratory insufficiency. In response to injury, alveolar macrophages orchestrate inflammation of the lung by recruiting neutrophils and circulating macrophages to the site of injury leading to the death of alveolar epithelial cells and further impairing lung function [4] . Severe COVID-19 is often also characterised by concurrent vascular disease and coagulopathy. This includes breakdown of the vascular barrier, oedema, endothelialitis and thrombosis. Thrombotic and microvascular complications are frequently recorded in deceased patients, suggesting that vascular pathology is a major driver of severe disease [5] . While these mechanisms contribute to disease progression, it is unclear which mechanisms are the predominant factors in causing fatality. Transcriptomic analysis of patient tissue offers a unique opportunity to not only understand the pathogenesis of SARS-CoV-2 but also to assist in biomarker validation. To date, the major search for clinical biomarkers, including transcriptomic analyses on COVID-19 patient samples, have been performed on whole blood. Blood samples from infected patients are highly accessible but it is not clear whether blood parameters accurately reflect either the viral load or the level of tissue damage that clinical treatments aim to control [6] . The lung is the primary site of viral replication, yet only a limited number of transcriptomic studies have been performed on the lungs of COVID-19 patients. Such studies have used transcriptomic analyses on only a small patient cohort (1-2 patients) [7, 8] or relied on 'bulk' sequencing of whole lung tissue, where any insight into tissue heterogeneity in this complex organ is lost [5, 9, 10] . In contrast, Nanostring GeoMX TM Digital Spatial Profiling (DSP) gene expression panels applied directly on tissue sections take into account the spatial location of transcriptomic features and are able to directly sample intra-organ heterogeneity [11, 12] . This provides a much deeper picture of the cellular changes driven by viral infection, and is a powerful way to define the characteristics of the host response to the virus and the spatial relationships between them. Indeed, such preliminary transcriptomic analysis of COVID-19 patient lungs appears promising [13] . However, such analyses require the application of advanced bioinformatics techniques to control for the numerous potential confounders present in gene expression data sampled in such a complex experimental design. In addition, there is a clear need to distinguish the host SARS-CoV-2 transcriptomic response (and associated spatial heterogeneity) from other infectious triggers of ARDS such as influenza virus, which may be circulating concurrently with SARS-CoV-2. To address these questions, we investigated the pathogenesis of SARS-CoV-2 in lung tissue using spatial transcriptomics approaches from rapid autopsy samples taken from 10 COVID-19 patients, five 2009 pandemic H1N1 (pH1N1) influenza virus patients and four control patient samples collected at a single institution. The Nanostring GeoMX TM DSP platform combined with bioinformatic modelling allowed deconvolution of individual variability from viral and host factors. These data showed that viral load has a heterogenous impact within a patient on the pulmonary transcriptomic response, and there is significant overlap in the gene expression profiles between influenza virus and SARS-CoV-2 infected samples. However, several unique transcriptomic signatures were detected in the lungs of COVID-19 patients. This study demonstrates the strengths of spatial transcriptomics coupled with powerful linear modelling bioinformatic techniques to broadly identify key pathways involved in viral infections and to clearly differentiate the involvement of distinct cellular and molecular pathways in different infections. Tissue microarray (TMA) cores were prepared from autopsied pulmonary tissue from 10 SARS-CoV-2 and 5 pH1N1 patients who died from respiratory failure (ARDS)(Supplementary Table 1 -3). Control material was obtained from 4 uninfected patients (Supplementary Table 3 For both groups, initially, the imaging exams such as X-rays and CTs were analyzed to identify the pulmonary segments with more severe lung injury. Once we confirmed radiographically evident and representative lesions (especially the left lobes pulmonary segments, since they can be removed in a more agile and practical mini-thoracotomy technique), we collected the sample through an anterior mini-thoracotomy on the fourth or fifth intercostal space. The time elapsed between the patient's death, obtaining the consent form by the relatives, and removing the fragments through the mini-thoracotomy did not exceed 4 hours (COVID-19 and H1N1 groups). In all these cases (COVID-19 and H1N1 group), the samples were similarly sized (3x3 cm) and were delicately handled and resected using surgical scissors. The samples were then fixed into a 10% concentration formalin solution and kept in it for at least 24 hours until blocking and slicing for microscopic analysis. The lung formalin-fixed paraffin-embedded (FFPE) samples were stained with hematoxylin and eosin (H&E) to observe the histopathological aspects and find the appropriate and representative areas for punch and construct the TMA. Clinical data were obtained from medical records during hospitalization in the ICU (at Hospital Marcelino Champagnat in Curitiba, Brazil, for the COVID-19 group and at Hospital de Clínicas in Curitiba, Brazil for the H1N1 group). Also included in this study was a Control group of lung samples from patients who died from neoplastic or cardiovascular diseases (not involving lung lesions). The samples of the Control group were obtained by the conventional autopsy technique. The time elapsed between the patient's death, obtaining the consent form by the relatives, and removing the fragments through the conventional autopsy did not exceed 6 hours. These Control group samples were then fixed into a 10% concentration formalin solution and kept in it for at least 24 hours until blocking and slicing for microscopic analysis. The lung formalin-fixed paraffin-embedded (FFPE) samples were stained with hematoxylin and eosin (H&E) to observe the histopathological aspects and find the proper punch areas to construct the TMA. Clinical data were obtained from medical records during hospitalization at Hospital de Clínicas in Curitiba, Brazil. Tissue Preparation and Histopathology. Thirty 5µm thick serial sections were cut from the TMA blocks onto positively charged slides (Bond Apex) and sections were stained with hematoxylin and eosin (H&E) and Masson's trichrome stain. Brightfield images were obtained using an Aperio (Leica Biosystems, US) slide scanner for assessment of histopathology by a pathologist. TMA cores were assessed for overall histological pattern. Regions of interest (ROIs) were semi-quantitatively analysed for alveolar haemorrhage and oedema, hyaline membrane formation, acute inflammation, type 2 pneumocyte hyperplasia, squamous metaplasia, capillary congestion, fibroblastic foci, and interstitial inflammation. Bond-RX autostainer (Leica Biosystems, US) with antibody targeting SARS-CoV-2 spike protein (Abcam, ab272504) at 2g/ml. Heat induced epitope retrieval was performed in buffer ER1 at 100°C for 20 minutes, and signal visualized with 3,3'-Diaminobenzidine (DAB) substrate. Slides were imaged using a Zeiss Axioscanner (Carl Zeiss, Germany) and IHC was scored by a pathologist for bronchiolar epithelium, type 2 pneumocytes, interstitial lymphocytes and alveolar macrophage compartments. RNAscope ® probes (ACDbio, US) targeting SARS-CoV-2 spike mRNA (nCoV2019, #848561-C3), ACE2 host receptor mRNA (#848151-C2), and host serine protease TMPRSS2 bead-to-sample ratio for PCR product purification. Paired-end sequencing (2×75) was performed using NextSeq550 up to 400M total aligned reads. Fastq files were processed by DND system and uploaded to GeoMX DSP system where raw and Q3 normalized counts of all targets were aligned with regions of interest (ROIs). Transcriptomic Data. Data used in this study result from an mRNA assay conducted with the NanoString's GeoMx Immune atlas panel using the GeoMX Digital Spatial Profiler (DSP). The data were measurements of RNA abundance of over 1800 genes, including 22 addin COVID-19 related genes, 4 SARS-CoV-2 specific genes and 2 negative control (SARS-CoV-2 Neg, NegProbe) genes and 32 internal reference genes. Samples were acquired from 3 TMA slides each containing 10 tissue cores from 10 patients. The three slides correspond to either COVID-19, H1N1 or Control Tissue biopsies. Transcriptomic measurements were made on regions of interests within each core. Control samples are COVID-19-free and H1N1-free, but they are not necessarily disease-free samples; they originate from patients with other nonviral diseases. In total, 60 ROIs were analysed (47 COVID-19, 8 H1N1 and 5 "Control"). Factors considered in this dataset include disease type (COVID-19, H1N1 and Control), patient of origin, dominant tissue type (Type 2 pneumocytes, bronchiolar epithelium, hyaline membranes, macrophages) and viral load (based on RNAscope scores). Data exploration and quality checks were conducted on the Q3 normalised count data generated from the CTA DSP mRNA assay. Relative log expression (RLE) plots were used assess the presence of unwanted variation in the form of batch effects [14] . Data were normalised using the trimmed mean of M-values (TMM) method [15] using all genes in the panel. Specifically, log2-transformed transcript abundance data were mediancentred for each gene, and then within each sample the difference between the observed and population median of each gene was calculated. Principal components analysis (PCA) of the samples was conducted to identify variability related to specific factors in the dataset and experimental design. Differential expression (DE) analysis, Gene Ontology and KEGG pathway enrichment analysis were carried out using R/Bioconductor packages edgeR (v3.30.3) [16, 17] and limma (v3.44.3) [18] . Briefly, differential expression was modelled using linear models with various experimental factors as predictors. Variation in gene expression was modelled as the combination of a common dispersion that applies to all genes and a gene-specific dispersion. Limma was used to model and estimate the variation of each gene by borrowing information from all other genes using an empirical Bayes approach, thereby allowing estimation of the common and gene-wise variation with as few as 2 replicates per condition. Once the linear model was fit to a given experimental design, various contrasts of interest were used to query the data for differential expression. The resulting statistic was an empirical Bayes moderated t-statistic which was more robust than a t-statistic from a classic t-test. Based on the observed differences between cores and tissue structures from the PCA analyses, considerations were made to allow for similarity that exists for regions originating from the same core using duplicationCorrelations in limma [19] . The flexible modelling framework afforded by linear models was used to take into account differences between heterogenous tissue structures (i.e. bronchiolar epithelial samples vs the rest of the samples) by including them as covariates in the models. Two main factors were investigated in this analysis, namely disease type (COVID-19, H1N1 and Control) and viral load. For disease type, two comparisons were investigated: COVID-19 against Control with cores and dominant tissue type as covariates, and COVID-19 against H1N1 samples, with cores and dominant tissue type as covariates. For these two contrasts, the voom-limma with duplicationCorrelations pipeline [20] was used to fit linear models. The TREAT criteria was then applied [21] (p-value <0.05) to perform statistical tests and subsequently calculate the t-statistics, log-fold change (logFC), and adjusted p-values for all genes. For viral load, the voom-limma with duplicationCorrelations pipeline [20] (p-value <0.05) was used to fit linear models for the comparison between high and low viral load regions originating from COVID-19 patient cores (disease type as a covariate). The contrast was applied on two resolutions, firstly by grouping ROIs from the same cores with the same degree of viral load (cores/patients-based approach); secondly by treating each regions/ROIs sample as an independent observation (ROI-based approach). Gene set enrichment analysis (GSEA) were performed using the fry approach from the limma package. Gene sets from the Molecular Signatures Database (MsigDB) Hallmarks [22, 23] , C2 (curated gene sets) and C7 (immunologic signature gene sets) categories were tested using fry. Pathways from the KEGG pathway database were tested using the kegga function in the limma package and gene ontology enrichment was assessed using goana from the limma package [19] . Lung tissues were obtained by rapid autopsy from 10 SARS-CoV-2 infected patients, five pH1N1 influenza patients and four non-virally infected control patients. The SARS-CoV-2 cohort was made up of four females and six males with mean age of 76 years (ranging from 46-93 years). These patients exhibited multiple comorbidities including obesity, diabetes and heart disease and had received varying treatment following hospital admission. Patient survival after admission ranged from 8 to 38 days, with time on mechanical ventilation ranging from 3 to 21 days (Supplementary Table 1 ). Blood profiling performed on admission and again immediately prior to death showed elevated D-dimer levels for 7 of 10 patients consistent with vascular coagulopathy (Supplementary Table 2 ). The pH1N1 influenza virus patient cohort comprised five males with a mean age of 45.8 years (ranging from 31-53 years). All patients had received mechanical ventilation from the time of admission until death. The uninfected control cohort consisted of four males with a mean age of 36.25 (1 ranging from 18-60 years) whose tissues were donated following death from nonviral and non-respiratory diseases (Extended Data Table 3 ). Tissue microarrays (TMAs) were prepared from pulmonary formalin-fixed paraffin-embedded cores of each of the three cohorts (one core per patient) and blinded histological examination was performed by a pathologist (Supplementary Table 4 ). Acute phase diffuse alveolar damage (DAD) characterised by hyaline membranes, type 2 pneumocyte hyperplasia and alveolar septal thickening was a prominent feature of SARS-CoV-2 cores. pH1N1 tissues also showed features of DAD, with alveolar haemorrhage evident in two of the cores. These findings were not seen in control cores. Fibrosis was not a prominent feature, apart from one SARS-CoV-2 core that exhibited moderate interstitial and alveolar fibrosis with concordant histological organising pneumonia. Two pH1N1 cores showed mild interstitial fibrosis that were characterised by acute phase or late organising phase DAD. Transcriptional profiling of target cores was guided by the detection of SARS-CoV-2 spike mRNA. Two cores ( Figure 2 ; LN1 and LN3) exhibited strong signals for viral mRNA (Fig. 2a) and these cores were comprehensively evaluated using the GeoMx DSP platform (shaded regions, Fig. 2b ). Twenty-five regions of interest (ROIs) were selected from these two RNAscope® positive cores and 21 ROIs were selected from remaining eight cores of SARS-CoV-2 infected patients (minimum of one ROI per patient) for which viral mRNA was below detection by RNAscope. Eight ROIs were selected from the lungs of five pH1N1 patients (minimum of one ROI per patient), and five ROIs were selected from the lungs of control patients (minimum of one ROI per patient). SARS-CoV-2 RNAscope ® abundance was semi-quantitatively assessed for each ROI. Representative scoring criteria from 0 to 3 are shown in Supplementary Fig. 2 where positive scores were allocated only to regions within LN1 and LN3 (Supplementary Table 5 ). In addition, concordant spike protein immunohistochemistry (IHC) appeared to be specific for type 2 pneumocytes and bronchiolar epithelial regions, however, some staining was also observed within interstitial lymphocytes and alveolar macrophages ( Supplementary Fig. 1 ). Interestingly, spike protein IHC indicated virus in cores other than LN1 and LN3 (Supplementary Table 5 ) that did not produce a detectable signal using RNAscope ® . Low level non-specific staining was observed in some regions of pH1N1 and control tissues Supplementary Table 5 . We thus focussed our ROI selection on cores that indicated high mRNA integrity and viral detection by RNAscope ® . Assessment of histology by a pathologist allowed ROIs to be grouped by tissue architecture prior to further analysis. The predominant cell types found within each ROI were annotated in addition to scoring of their histopathology and spike protein IHC (Supplementary Table 5 ). ROIs could be broadly categorised into regions dominated by bronchiolar epithelial structures, type 2 pneumocytes, hyaline membranes and macrophages ( Supplementary Fig. 2 ). Alveolar oedema, acute inflammation and squamous metaplasia were not observed in diseased or control cores, however, type 2 pneumocyte hyperplasia and capillary congestion were present in SARS-CoV-2 pulmonary tissues. Q3 housekeeping normalised data was corrected for systematic bias by TMM to allow for comparisons between SARS-CoV-2, pH1N1 influenza virus and control samples ( Supplementary Fig. 3 ). Principal components analysis (PCA) plots were then used to explore the variability in SARS-CoV-2, pH1N1 influenza virus and control lung samples and to determine whether they could be clearly separated simply based on the type of infecting pathogen and to identify factors that might confound differential expression analysis. Principal components (PCs) capture orthogonal dimensions of variability in the data in descending order of contribution. Non-infected samples clearly separate from other samples on PCs 1 and 3, while SARS-CoV-2 and pH1N1 influenza samples remained associated on these axes (Fig. 3a) . Indeed, across PCs 1-4, SARS-CoV-2 and pH1N1 influenza samples substantially overlapped. This suggests that the transcriptomic profile of SARS-CoV-2 and pH1N1 influenza virus tissues were not highly variable at low viral load, and that a common transcriptome was induced by these two viruses in fatal cases. To further examine the major sources of variation in the dataset, we labelled PCA plots by multiple different parameters including sample core (Fig. 3b) , dominant cell type (Fig. 3c ) and viral load (Fig. 3d) . When samples were labelled by core, PC 3 and PC 4, were able to distinguish patient specific effects, as shown by a clear separation of cores LN1 and LN3, with clustering of other cores with multiple ROIs (LN10) (Fig. 3b ). This suggests that within patient correlations and patient-patient variation has an important underlying effect on the observed transcriptome and that this needs to be taken into account in any subsequent analysis Additionally, samples separating on PC2, and clearly clustered from other samples when plotting PC1 and 2, were grouped based on the dominant cell type found within the sample core ( Fig. 3c) . Specifically, ROIs where bronchiolar epithelial cells were the dominant cell type clustered together independent of whether the patient had pH1N1 influenza virus or SARS-CoV-2 ( Fig. 3c) , or belonged to the non-infected controls. This implied that the tissue heterogeneity captured by dominant cell type was a substantial contributor to the variability in our experiment. While epithelial cells are the primary cellular target of both SARS-CoV-2 and influenza virus [24] , our data did not identify viral epithelial tropism as a key regulator of signal as the cores dominated by bronchial epithelial cells did not exhibit the highest viral load (as defined by RNA Scope for SARS-CoV-2) (Fig. 3d) . However, the high SARS-CoV-2 viral load cores separate from other virally-infected samples and control samples on PC3 (Fig. 3d) . Together, these analyses suggest that unbiased and accurate comparative transcriptomics analyses between uninfected, SARS-CoV-2 and pH1N1 influenza lung samples must take account of the dominant cell type present in samples, inter-patient variations and within-patient correlations. The ability to include structural and spatial covariates demonstrates the key advantage in using a spatial transcriptomics as compared with bulk RNA sequencing, but also necessitates careful construction of the model used to compare expression across the ROIs sampled. together with CSF3, colony stimulating factor, normally associated with mobilization of cells from the bone marrow (Supplementary Table 6 ). A number of tumor-associated genes including GAGE1 and MAGEC2 were downregulated, with these genes being associated with cell survival. Interestingly, PCK1, which is required for tissue gluconeogenesis, and MPL, that regulates the generation of megakaryocytes and thus platelet formation were also highly reduced [25] . A more diverse array of upregulated genes was found to be differentially expressed in the lungs of SARS-CoV-2 patients compared with those from patients who died of non-viral causes (Supplementary Table 7 ). These included genes associated with antigen presentation (e.g. B2M, HLA-DRA, HLA-C and HLA-E), the Type I interferon (IFN) response (e.g. LY6E and IFI27), fibrosis and epithelial cell growth (e.g. COL3A1, FN1, KRT7, COL1A1, KRT18 and KRT19) and the complement cascade (C1QA and C1R). Consistent with these observations, gene set enrichment analysis showed that hallmark genes of biological processes such as reactive oxygen species, complement, IL6/JAK/STAT3 signalling, apoptosis, p53 signalling IFN-γ response, hypoxia, IFN-α response, coagulation, TGF-β signalling, IL-2/STAT5 signalling, angiogenesis, TNF-α signalling via NF-B, and inflammatory response were all significantly upregulated in SARS-CoV-2 infected lungs (Table 1 ). Gene set enrichment using blood coagulation, hypoxia responses and angiogenesis related genes from nanoString's nCounter® PanCancer Progression Panel further confirmed the upregulation of pathways associated with blood coagulation and angiogenesis responses ( Table 2 ). These analyses show that SARS-CoV-2 samples significantly upregulate genes associated with blood coagulation, angiogenesis, the IFN-α and IFN-γ response and positively regulate cytokine production and cytokine stimulus during the immune response (Fig. 5 ). These data are consistent with previous descriptions of a 'cytokine storm' response in patients with severe SARS-CoV-2, often accompanied by various coagulopathies [26] . They are also in agreement with the notion that a late stage or delayed type I IFN response may be observed in concert with a pro-inflammatory response [27] . The pathogenesis of SARS-CoV-2 is a complex interplay of host immunopathology and viral load. To determine how viral load influenced the host transcriptomic signal, we performed a differential expression analysis of SARS-CoV-2 samples stratifying based on viral load. Patient sample cores were classified as either 'high' or 'low' based on the presence of viral mRNA by RNAscope ® . Thus, two cores which were scored positively by RNAscope ® , LN1 and LN3, were designated as 'high' while the remaining 8 cores from SARS-CoV-2 patients were scored as 'low'. This analysis showed that only 17 up-regulated and 7 down-regulated genes were differentially expressed between the high viral load SARS-CoV-2 infected cores and the other SARS-CoV-2 cores (Fig. 6 , Table 3 ). Consistent with our RNAScope analyses, SAR-CoV-2 specific genes S and ORF1ab were the most strongly upregulated genes in the high viral load group (1.916 and 1.861 log2 fold change, respectively). This is also consistent with previous observations that in at least some COVID-19 patients, these are the most highly expressed viral genes in the lung [28] . All other upregulated genes with >1 log2 fold change in expression were associated with the type I IFN response (CXCL10, IFIT3, ISG15, MX1, GBP1 and IFI6) (Fig. 6a , Table 3 ). The above analysis assumed that every ROI derived from the same core would have the same degree of viral infection (i.e. high or low). However, viral infection is unlikely to be uniform across tissues. Indeed, none of the bronchiolar epithelial ROIs in the two viral high cores were classified as high (Fig. 3d , Supplementary Table 5 ). To examine the impact of this, we considered each ROI as an independent sample which was classified as either 'high' (score > 0) or 'low' (score = 0 ) by RNAscope positivity. By stratifying cores with this criterion and accounting for the dominant cell type as well as correlations between ROIs from the same core, we identified 33 differentially upregulated genes and 132 downregulated genes in high viral load ROIs compared to low viral load ROIs (Supplementary Table 8 ). Once again, S and ORF1ab were the most strongly upregulated genes in the high viral load group (2.415 and 2.365 log2 fold change respectively). Other upregulated genes of note were those associated with the type I IFN response (e.g. RSAD2, OASL, TLR8, IL1RN and IKBKE), chemokines associated with IFN (e.g. CXCL11, CCL8 and CCL7) and the SARS-CoV-2 receptor ACE2. A more diverse array of genes were downregulated in high viral load samples including those associated with the complement cascade (e.g. C3), ribosomal proteins (e.g. RPS6 and RPL23) and antioxidants (e.g. PRDX5) ( Figure 6B ). Taken together, these data suggest that a pronounced type I interferon gene signature is associated with SARS-CoV-2 RNA expression and the spatial context provided by this technology enables finer details of viral load associations to be determined. Table 7 ). Two heat shock protein family members (HSPA1A and HSPA) were also identified as significantly down-regulated in COVID-19 tissues. These two genes were also significantly upregulated in pH1N1 tissues compared to uninfected tissues (data not shown), suggesting a potential influenza specific signature. To understand the relevance of these genes to COVID-19 and pH1N1, we have compared these differentially expressed genes to other published spatial transcriptomics approaches both experimentally and analytically and that is shown by the varying sizes of the gene sets ( Figure 8 ). It is clear that our approach is able to identify a focussed set of gene signatures specific to COVID-19 infection as compared to other studies (IFI27 and LY6E common across all studies, and IFI6 common in 4 out of the 5 studies). Together, these data indicate that fatal pH1N1 influenza virus and SARS-CoV-2 lung samples have only subtly different host transcriptomes but may be potentially distinguished by specific gene signatures. This is further clarified in the pH1N1 vs Control (Supplementary Fig 6) and GSEA analysis ( Supplementary Fig 7) which shows many common activation pathways for both COVID-19 and pH1N1 infected tissue when compared to control. These data suggest that many DE pathways identified in the COVID-19 sample may be indicative of a more generic viral-induced ARDS, whilst only a limited number of DE genes may be unique to COVID-19 itself. Understanding the biological functions, networks, and host-pathogen interactions that impact organ and disease development requires both cellular information and a spatial context. Analyses of bulk RNA sequencing or scRNA sequencing provides a global overview of an organ's response to a pathogen, typically identifying broad inflammatory pathways. However, such approaches fail to identify subtle individual cellular changes that are spatially distinct. This has particular impact when considering innate and adaptive immune cells that may be crucial for understanding pathogen-specific responses, and to distinguish the profile of one pathogen from another. In contrast, spatially resolved transcriptomes of virally-infected tissues offers the possibility of disentangling the individual infected cells, contributions of viral load, cellular responses and hence patient-to-patient variability. The clinical spectrum of COVID-19 is highly variable. Association of viral load with disease severity would provide a mechanism for early stratification of patients for treatment options. Indeed, analyses of nasopharyngeal swabs suggests that nasopharyngeal SARS-CoV-2 RNA is independently correlated with disease severity [29] , similar to earlier observations with SARS-CoV-1 [30] . However, the association between viral replication in the lung and severe disease remains more complicated. Here, we observed that 8 out of 10 patients who died because of COVID-19 had no detectable viral RNA in lung biopsies as determined by RNAscope. Whilst this rate of viral RNA detection may be lower than expected, it is important to recognise that viral antigen deposition in the lung is heterogeneous, and therefore these data may simply reflect random sample selection. There was no notable difference in time from admission to death in patients who were positive vs. those who were negative by viral RNA-Scope. However, we do not have information for these patients on time between symptom onset and hospital admission which may have influenced viral clearance. Nevertheless, it is important to note that all patients in this study were PCR positive for SARS-CoV-2 RNA via nasopharyngeal swab. Therefore, if the absence of viral RNA in the lungs of 8 patients reflects viral clearance this was restricted to the lower respiratory tract. Transcriptomic analysis showed that areas of high viral load were associated with a pronounced type I IFN response, consistent with other preliminary spatial analysis of the lungs of COVID-19 patients [13] and assumedly due to increased viral pathogen-associated molecular patterns (PAMPs) available for type I IFN stimulation. Nevertheless, even in samples where no viral RNA was detected, a pronounced type I IFN gene signature could still be observed in the pulmonary tissue. These observations add weight to the growing understanding of the role type I IFNs in SARS-CoV-2 infections [31] . Recent findings in a COVID-19 tissue atlas study found upregulation of IFN-α and IFN-γ response genes and oxidative phosphorylation pathways in the lungs of COVID-19 patients vs control [32] . Current evidence suggests that an early and robust induction of type I interferons can help control viral replication and help ensure a mild infection. In contrast, a delayed induction of type I IFN (i.e. after the peak of viral replication) is largely irrelevant for viral control as viral titres have already declined at this point in the infection but may help perpetuate the detrimental pro-inflammatory response and lung damage [31, 33] . The transcriptomic data presented here thus speak to the potential dual-role of type I IFNs in SARS-CoV-2 infection. Pulmonary transcriptomic analysis is a powerful tool to delineate the pathogenesis of SARS-CoV-2 and identify how it differs compared to other respiratory pathogens such as influenza virus. Previous studies have suggested that the lungs of COVID-19 patients display an increased incidence of alveolar capillary microthrombi and thrombosis with microangiopathy compared to those of influenza patients [5] . Consistent with these observations, an earlier report using bulk lung RNA analyses to identify 69 differentially-expressed angiogenesis-related genes in COVID-19 patients, but not in influenza patients [13] . In the present study, we observed that both COVID-19 and influenza patients had an upregulation of genes associated with coagulation and angiogenesis, these genes were not differentially expressed between the two virus-infected patient groups. These contrasting data most likely reflects the fact that all samples in the present study were collected at the end stage of disease and accordingly, earlier differences in influenza and COVID-19 disease pathogenesis and severity were unable to be examined. However, these data do highlight both the difficulties and importance of finding specific biomarkers that are sensitive enough to differentiate COVID-19 from other respiratory viral infections even at end-stage disease. Several potential clinical and transcriptional biomarkers for triaging patients with COVID-19 have been explored. Clinical biomarkers include C-reactive protein, serum amyloid A, interleukin-6, lactate dehydrogenase, neutrophil-to-lymphocyte ratio, D-dimer, lymphocytes and platelet count [34, 35] . However, the majority of these studies have focussed on biomarkers in patient blood. While peripheral blood samples are a convenient and highly accessible site for clinical sampling, a biomarker that could be rapidly identified in the respiratory tissue (e.g. via a nasal swab) may be more accessible in low resource settings. IFI27 levels in the blood have been previously identified to have a 88% diagnostic accuracy (AUC) and 90% specificity in discriminating between influenza and bacterial infections in patients [36] . Interestingly, previous analyses of IFI27 in the blood of COVID-19 patients revealed that IFI27 was upregulated in SARS-CoV-2 infection [37] [38] [39] . Here, we identified IFI27 as differentiated upregulated in the lungs of both COVID-19 patients vs control patients and in the highly restricted set of genes differentially expressed between COVID-19 and influenza patients. These data raise the exciting possibility that IFI27 may not only represent a biomarker for severe COVID-19 but that it may also help differentiate this disease from other clinically similar viral infections. This becomes particularly important as the 'second wave' of COVID-19 in the US and Europe is set to overlap with the winter influenza season. Validation of this gene in nasal specimens, as well associations with disease severity will be required to confirm IFI27 as a gene signature that is useful in stratifying COVID-19 patients. This study was subject to several important limitations. Firstly, all data was derived from a small sample cohort derived from a single study site, including unbalanced patient vs control groups, and it remains to be determined how much these data can be extrapolated to other patient populations. Importantly, future studies with a larger number of patients should include an equal number of males and females as male sex is independently associated with increased COVID-19 severity [40] . Furthermore, additional studies are required across a broader range of patients (i.e. those with mild and moderate disease) to determine the therapeutic value of measuring IFI27 levels during infection. Moreover, there have been rapid advances using antibody guided ROI capture, whole transcriptome assays, and cellular deconvolution to infer cell types. However, despite these limitations these data reveal the unprecedented power of spatial profiling combined with detailed multiparameter bioinformatic analyses to dissect the key variables that contribute to differential gene expression across highly variable patient cohorts and the heterogeneous distribution of virus and immune responsiveness within tissues. Nanostring data that support the findings of this study will be deposited with the Gene Expression Omnibus (GEO) repository and made available on publication. All the data to evaluate the conclusions in this paper are present either in the main text or in supplementary materials. Table 2 Subepithelial lymphocytes 0 0 0 3 0 1 0 0 2 0 0 18 COVID LN3 Type 2 pneumocytes 0 0 2 3 0 1 0 2 1 3 2 41 COVID LN3 Type 2 pneumocytes 0 3 2 3 0 1 0 2 2 3 FN1 TYMP GPX1 HIF1A VEGFA GPI LAMA5 ITGA5 JUN PTEN TNFRSF12A VEGFB MCAM FGFR1 NRP1 EPHA2 CASP8 SETD2 ROBO4 SYK FAP TIE1 PDGFA THY1 JAM3 PIK3CA IL18 RORA MAP3K7 ANGPT1 ANGPTL4 PIK3CG MMRN2 ZC3H12A NOS3 CX3CL1 DLL4 KDR TNFSF12 JAG1 ITGB3 FGFR2 CXCR3 CEACAM1 PTGS2 VEGFC FGF18 EGF ANGPT2 SOX17 S100A7 FGF9 TGFB2 TAL1 CXCL8 ID1 H 1 N 1 v s Control ) World Health Organization Coronavirus Disease (COVID-19) Dashboard Pathophysiology of COVID-19-associated acute respiratory distress syndrome: a multicentre prospective observational study COVID-19 acute respiratory distress syndrome (ARDS): clinical features and differences from typical pre-COVID-19 ARDS Pathogenesis of influenzainduced acute respiratory distress syndrome. The Lancet Infectious Diseases Pulmonary Vascular Endothelialitis, Thrombosis, and Angiogenesis in Covid-19 Transcriptomic characteristics of bronchoalveolar lavage fluid and peripheral blood mononuclear cells in COVID-19 patients Imbalanced Host Response to SARS-CoV-2 Drives Development of COVID-19 Spatial mapping of SARS-CoV-2 and H1N1 Influenza Lung Injury Identifies Differential Transcriptional Signatures Dynamics of IgG seroconversion and pathophysiology of COVID-19 infections Multiplex digital spatial profiling of proteins and RNA in fixed tissue High-Plex and High-Throughput Digital Spatial Profiling of Non-Small-Cell Lung Cancer (NSCLC) Temporal and Spatial Heterogeneity of Host Response to SARS-CoV-2 Pulmonary Infection RLE plots: Visualizing unwanted variation in high dimensional data A scaling normalization method for differential expression analysis of RNA-seq data Differential expression analysis of multifactor RNA-Seq experiments with respect to biological variation edgeR: a Bioconductor package for differential expression analysis of digital gene expression data limma powers differential expression analyses for RNA-sequencing and microarray studies Use of within-array replicate spots for assessing differential expression in microarray experiments voom: precision weights unlock linear model analysis tools for RNA-seq read counts Testing significance relative to a fold-change threshold is a TREAT Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profiles The Molecular Signatures Database (MSigDB) hallmark gene set collection Distribution of ACE2, CD147, CD26, and other SARS-CoV-2 associated molecules in tissues and immune cells in health and in asthma, COPD, obesity, hypertension, and COVID-19 risk factors Coagulopathy in COVID-19: Focus on vascular thrombotic events Mouse model of SARS-CoV-2 reveals inflammatory role of type I interferon signaling Differential Expression of Viral Transcripts From Single-Cell RNA Sequencing of Moderate and Severe COVID-19 Patients and Its Implications for Case Severity SARS-CoV-2 viral load predicts COVID-19 mortality Coagulation disorders in coronavirus infected patients: COVID-19, SARS-CoV-1, MERS-CoV and lessons from the past Type I and Type III Interferons -Induction, Signaling, Evasion, and Application to Combat COVID-19 The type I interferon response in COVID-19: implications for treatment The role of biomarkers in diagnosis of COVID-19 -A systematic review Biomarkers and outcomes of COVID-19 hospitalisations: systematic review and metaanalysis A novel immune biomarker IFI27 discriminates between influenza and bacteria in patients with suspected respiratory infection The differential immune responses to COVID-19 in peripheral and lung revealed by single-cell RNA sequencing A single-cell atlas of the peripheral immune response in patients with severe COVID-19 Systems biological assessment of immunity to mild versus severe COVID-19 infection in humans Diabetes and Overweight/Obesity Are Independent, Nonadditive Risk Factors for In-Hospital Severity of COVID-19: An International, Multicenter Retrospective Meta-analysis Spatial mapping of SARS-CoV-2 and H1N1 lung injury identifies differential transcriptional signatures The NUSSI. Circuits between infected macrophages and T cells in SARS-CoV-2 pneumonia Temporal and spatial heterogeneity of host response to SARS-CoV-2 pulmonary infection Covid-1 Covid-10 Covid-14 Covid-18 Covid-19 Covid-21 Covid-22 Covid-26 Covid-28 Covid-29 Covid-30 Covid-1 Covid-10 Covid-14 Covid-18 Covid-19 Covid-21 Covid-22 Covid-26 Covid-28 Covid-29 Covid-30 STAT1 HIF1A IFNGR1 TNFRSF1A IL1R1 JAK1 PARP9 PTPRC CXCR4 CSF1 IFNGR2 TNFAIP3 IRF7 PTPN11 KLF4 IFIH1 YTHDF2 RUNX1 IRF3 TLR2 MAVS IFNAR2 CASP1 IRAK1 CCL5 NLRC5 SOCS3 IKBKG CASP8 RIPK1 PYCARD DDX58 BIRC3 SOCS1 JAK2 HSPA1A SYK FADD IKBKB IRAK3 CYLD SIGIRR EDN1 PIAS4 IL1R2 CHUK RIPK2 AXL TICAM2 ANGPT1 TRAF1 TNF IKBKE TLR4 IRAK2 TBK1 CD24 TRAF2 IL1RN TREM2 TXK PPARG IL6 APOA1 F2RL1 IRGM TSLP WNT5A LIFR IFNG ARG1 IL3 IL7 -1 0 1 2 3 4 CD74 B2M HLA-F HLA-E CD81 MIF IL1R1 FCER1G HK1 MAPK3 MAPK14 TNFRSF14 MAVS LILRB1 DDX58 CD36 RSAD2 IL18 TRAF6 MAP3K7 FZD5 TLR4 TRAF2 IL18R1 TICAM1 NOD2 IL1B BCL10 HLA-G IL6 GATA3 CD160 NR4A3 F2RL1 TNFSF4 WNT5A XCL1 NLRP3 -1 0 1 2 3 4 CD74 B2M LY6E HLA-C IFI27 C1S TXNIP IFITM1 STAT2 IFITM2 PSMB9 MX1 CD47 TAP1 ISG15 IFIT3 CXCL10 LAMP3 PARP9 UBA7 PSMB8 CSF1 BST2 IFI35 IRF9 IRF7 OAS1 GBP4 IFIH1 EIF2AK2 IL4R PARP12 CASP1 TRIM21 CASP8 IRF2 IFIT2 RSAD2 DHX58 TMEM140 SAMD9 HERC6 RIPK2 SELL CCRL2 GBP2 IRF1 CXCL11 OASL ISG20 IL15 IL7 -1 0 1 2 3 4 CD74 B2M LY6E IFI27 C1R HLA-A HLA-B SERPING1 C1S TXNIP STAT1 STAT3 STAT2 APOL6 HIF1A HLA-DMA IFITM2 PSMB9 CDKN1A CFB MX1 ICAM1 TNFSF10 TAP1 ISG15 IFIT3 CXCL10 PSMB10 PML PSMB8 SOD2 BST2 IFI35 HLA-DQA1 IRF9 CCL2 TNFAIP3 IRF7 CXCL9 FCGR1A IL10RA OAS3 GBP4 CD40 FAS IFIH1 EIF2AK2 OAS2 TAPBP IL4R PARP12 CSF2RB NFKBIA IFNAR2 CASP1 CCL5 IRF8 NLRC5 CMKLR1 SOCS3 PIM1 TRIM21 IDO1 CASP8 RIPK1 PSMB2 CD38 DDX58 IL15RA CASP7 SOCS1 JAK2 IRF2 NFKB1 IFIT2 RSAD2 FPR1 CASP3 DHX58 SLAMF7 VCAM1 MYD88 HERC6 RIPK2 IRF4 CD86 ITGB7 GZMA CD274 IL2RB CD69 SELP IRF1 PTGS2 CXCL11 PLA2G4A OASL ISG20 TNFAIP6 HLA-G NOD1 IL6 STAT4 KLRK1 IFIT1 CCL7 XCL1 IRF5 IL15 IL7 -