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  -