Profiling of lung SARS-CoV-2 and influenza virus infection dissects virus-specific host responses and gene signatures

Patient characteristics

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 non-viral and non-respiratory diseases (Extended Data table 3).

Histopathological characterisation

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.

RNAscope® and ROI selection for digital spatial transcriptional profiling

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 (figure 2a) and these cores were comprehensively evaluated using the GeoMx DSP platform (shaded regions, figure 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).

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 1

Schematic of the study.

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 2

RNA-FISH and digital spatial profiling morphology marker visualisation of nCoV2019 mRNA positive cores. ROI selection was guided by RNAscope-FISH staining for nCoV2019 spike mRNA. (a) RNA-FISH staining on COVID TMA showing the cores LN1 and LN3 which were highly positive for nCoV2019 spike mRNA (red). Nuclei are shown in blue (DAPI staining). Scale bar, 1000 µm. (b) Representative ROIs selected across LN1 and LN3 SARS-CoV-2 virus positive cores on immunohistological paraffin embedded sections stained for monocytes (CD68, yellow), T cells (CD3ε, red), lung epithelial cells (PanCK, green) and nuclei (DAPI, blue). Scale bar, 500 µm. c. High resolution image of high viral core LN1 where each region of interest has been determined for the majority cell type (Hyaline membrane/T2 pneumocyte). d. High resolution image of high viral core LN3 where each region of interest has been determined for the majority cell type (Hyaline membrane, T2 pneumocyte, macrophages, bronchiolar epithelial cells).

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®.

Histopathology and dominant features of ROIs

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.

Data normalisation strategy and exploratory analysis of spatial transcriptomic data

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 (figure 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.

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 3

Principal components identify variability and factors in the transcriptomic data. Principal components (PCs) capturing orthogonal dimensions of variability in the transcriptomic data in descending order of contribution (i.e. PC1-PC2, PC3-PC2 and PC3-PC4) were plotted stratifying based on the following factors in the experimental data: (a) Disease types, (b) cores/patients with disease types, (c) dominant tissue types with disease types and (d) Viral Load in cores/patients. Disease type in panels a, b and c labelled as shapes with Control as circles, COVID-19 as triangles and H1N1 as squares. Viral load labelled in panel b as shapes with high viral load as solid circle and low viral low as unfilled triangle. Cores and patients labelled in panel b and c using differing colors. Dominant tissue types in panel c labelled as differing fill colors (BRO_EPI, bronchiolar epithelium; HYALINE, hyaline membranes; MACRO, macrophages; T2_PNEU, Type2 pneumocytes; OTHERS, other classifications).

To further examine the major sources of variation in the dataset, we labelled PCA plots by multiple different parameters including sample core (figure 3b), dominant cell type (figure 3c) and viral load (figure 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) (figure 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 (figure 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 (figure 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) (figure 3d). However, the high SARS-CoV-2 viral load cores separate from other virally-infected samples and control samples on PC3 (figure 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.

Severe COVID-19 infection results in a global pro-inflammatory response and downregulation of immune cell effector and regeneration pathways

Differential gene expression analysis was subsequently performed on all sample cores, taking into account the effect of the dominant cell type and the appropriate within and inter-patient variations. 286 genes were identified as significantly upregulated and 20 genes were significantly downregulated in the lungs of SARS-CoV-2 patients compared with non-virally infected controls (figure 4, Supplementary Tables 6, 7).

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 4

COVID-19 infection drives pro-inflammatory response and suppresses immune cell effector and regeneration. Distribution of differentially expressed genes as a function of the average transcript expression and fold change (log2) identified in COVID-19 samples versus uninfected (control) samples. Down-regulated genes included immune related, cytokine, tumor/cell-survival associated genes and cell regeneration genes, inferring a suppression of immune cell effect and regeneration. Up-regulated genes are enriched with pro-inflammatory genes including Type I interferon (IFN) response and fibrosis genes. Representative genes from these processes are highlighted. Differential expression genes generated using the voom-limma pipeline with limma:duplicationCorrelations and applying TREAT criteria with absolute log fold change >1.2 with p-value <0.05.

Down regulated genes included immune-related genes such as B lymphocyte antigen CD19, cytokine genes such as IL7, IL9 and IL13 and components of the complement cascade (e.g. C9) together with CSF3, colony stimulating factor, normally associated with mobilisation 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 (figure 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].

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 5

Gene set enrichment analyses (GSEA) reveals upregulated cytokine responses with accompanying coagulopathies in COVID-19 infection. GSEA for differentially expressed genes comparing COVID-19 and uninfected samples were conducted using the estimated log2 fold change and p-values with the distribution of significantly differentiated genesets visualised as a function of the log2 fold change (logFC) and genes sorted by the logFC (waterfall plot). The genesets visualised are custom genesets for “angiogenesis response”, “blood coagulation”, “hypoxia” responses based on nanoString's nCounter® PanCancer Progression Panel were identified by GSEA to be differentially upregulated in COVID-19 samples. (see table 2), MSigDB Hallmark gene set for IFN-α and IFN-γ and MSigDB Gene Ontology (GO) gene sets for “Regulation of response to cytokine stimulus’ and “Positive regulation of cytokine production involved in immune response’. Genesets ordered by gene counts. Refer to comprehensive list in table 1. GSEA conducted using limma-fry with FDR <0.05. The direction and relative change in expression are shown. Details of each geneset are provided in Supplementary Fig 4 and Supplementary Fig 5.

Limited differential gene expression associated with SARS-CoV-2 viral load

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 (figure 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) (figure 6a, table 3).

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 6

Resolution provided by spatial data reveals type I interferon gene signature associates with COVID-19 viral load. Distribution of differentially expressed genes as a function of the average transcript expression and fold change (log2) identified in high viral load versus low viral load COVID-19 samples. Analyses were conducted at two resolutions, firstly by (top) grouping ROIs from the same cores with the same degree of viral load (core patients-based approach); secondly by (bottom) treating each regions/ROIs sample as an independent observation (ROI-based approach). In both approaches, SAR-CoV-2 specific genes S and ORF1ab were strongly upregulated in the high viral load group, consistent with results of the RNAscope. Consistent with previous reports, type I IFN response genes were associated with SAR-CoV-2 RNA expression. Notably the core-based approach reveals limited differential gene expression with the finer resolution ROI-based approach revealing additional differential changes in complement cascade, ribosomal protein and antioxidants genes. Representative genes are highlighted. Differential expression genes derived using voom-limma pipeline with limma:duplicationCorrelations and applying absolute log2 fold change >1.0 with p-value <0.05.

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 (figure 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.

Severe COVID-19 is associated with a limited, differential transcriptome when compared with severe influenza

Only a limited number of genes were identified as differentially expressed between COVID-19 and pH1N1 influenza samples (2 downregulated and 4 upregulated genes, figure 7). Of the six genes significantly differentially regulated in COVID-19 samples, three were associated with the type I IFN response (LY6E, IFI27 & IFI6), two were heat shock protein family members (HSPA6 & HSPA1A) and one, CT45A1, when combined with various growth factors, is associated with cell survival and tumorigenesis. Interestingly, all three of the interferon-inducible genes were also significantly upregulated in COVID-19 samples compared to those without a viral infection (Supplementary 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 versus 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.

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 7

Limited transcriptomic differences associated with COVID-19 infection compared with pH1N1. Distribution of differentially expressed genes as a function of the average transcript expression and fold change (log2) identified in COVID-19 samples versus pH1N1 samples were visualised with the 6 differentially expressed genes labelled and classified based on associated biological processes. The genes associated with the type I Interferon response, heat shock protein family members and (associated) with cell survival and tumorigenesis are shown. This exclusive list of genes reveals the subtle differences between COVID-19 and pH1N1 infected transcriptome and presents a potential disease specific transcriptomic signature for distinguishing the two disease types. Differential expression genes derived using voom-limma pipeline with limma:duplicationCorrelations and applying TREAT criteria with absolute log fold change>1.2 with p-value <0.05.

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 8

Comparison of COVID specific gene sets identified across multiple studies. Upset plot describing the overlaps between gene sets across 5 different studies. The size of the gene sets varies across the studies with COVID specific gene set identified in this study (COVID19 versus pH1N1) having the smallest number of genes (6). Notably IFI27 and LY6E are common across all the studies while IFI6 is common amongst 4 out of the 5 studies. Note: red font: upregulation, blue font: downregulation in this study. Studies compared include Margaroli et al. 2021 [41], Grant et al. 2021 [42], Desai et al. 2020 [43] and Butler et al. 2021 [44].