Chronic lung diseases are associated with gene expression programs favoring SARS-CoV-2 entry and severity

Abstract

Patients with chronic lung disease (CLD) have an increased risk for severe coronavirus disease-19 (COVID-19) and poor outcomes. Here, we analyze the transcriptomes of 611,398 single cells isolated from healthy and CLD lungs to identify molecular characteristics of lung cells that may account for worse COVID-19 outcomes in patients with chronic lung diseases. We observe a similar cellular distribution and relative expression of SARS-CoV-2 entry factors in control and CLD lungs. CLD AT2 cells express higher levels of genes linked directly to the efficiency of viral replication and the innate immune response. Additionally, we identify basal differences in inflammatory gene expression programs that highlight how CLD alters the inflammatory microenvironment encountered upon viral exposure to the peripheral lung. Our study indicates that CLD is accompanied by changes in cell-type-specific gene expression programs that prime the lung epithelium for and influence the innate and adaptive immune responses to SARS-CoV-2 infection.

Fig. 1. Percentage of cells expressing SARS-CoV-2 receptor genes in lung cell types in different diagnosis subgroups.

a Percentage of cells expressing ACE2 and TMPRSS2 in all cell types. Numbers are the total number of ACE2 + or TMPRSS2 + cells in each cell type in the dataset. b, c Percentage of cells expressing ACE2 (b) and TMPRSS2 (c) in each diagnosis group in the epithelial cell types. d Venn diagram shows overlapping of cells co-expressing the proposed receptors (ACE2, BSG and NRP1) and the protease TMPRSS2. e, f Percentage of cells co-expressing receptors and TMPRSS2 split by cell type and diagnosis group. Plots were generated with mean values of percentage of cells per individual samples, and data are presented as mean values ± SEM. Significant differences between diagnosis groups were calculated using Tukey_HSD test, p value < 0.05: *p value < 0.01: **p value < 0.001: ***p value < 0.0001: ****.

Fig. 2. Expression profile of SARS-CoV-2 mediators and response genes in the epithelial cell population.

a Binary heatmap representing a manually curated list of genes associated with SARS-CoV-2. Orange elements indicate genes with increased expression and white elements indicate genes with decreased expression in CLD samples; Not detected: gene expression was not detected in either of the two tested populations (CLD vs. Control). Differential expressed genes (DEGs) between CLD and control samples (FDR ≤ 0.1) are outlined in black. b Violin plot depicts gene expression level in CLD and control of the two SARS-CoV-2 proteases TMPRSS2 and CTSL. c SARS-CoV-2 entry module score in different cell types, SARS-CoV-2 mediators included ACE2, BSG (CD147), NPR1, HSPA5 (GRP78), TMPRSS2, CTSL, ADAM17, FURIN. The outliers were removed in this plot, please see Supplementary Fig. 8a with outliers included. Boxes: interquartile range, lower and upper hinges correspond to the first and third quantiles, upper and lower whisker extends from the hinge to the largest values or smallest values of 1.5 x interquartile range; *p value < 0.05, **p value < 0.01, ***p value < 0.001, ****p-value < 0.0001, Tukey_HSD post-hoc test. ACE2 and ITGB6 protein expression in IPF lung sections. IPF lung sections stained for ACE2: d small airway, e large airway and f lung parenchyma. IPF lung sections stained for αvβ6: g small airway, h large airway and i lung parenchyma. j Semi-quantitative evaluation of ACE2 scoring among control (n = 12 for each tissue) and IPF (n = 62 for each tissue) sections (both the percentage of staining and staining intensity of ACE2 expression; 0-Negative; 1–0–⩽10%; 2-11–⩽25%; 3-⩽26%), data are presented as mean values ± SEM. Significant differences between IPF and control were calculated using Tukey HSD test, p value < 0.05 *. Scale bar = 100 µm. For d–j: a total of 12 normal lung samples and 62 IPF samples were used.

Fig. 3. CLD AT2 cells exhibit baseline differences in gene expression profile coping with viral infection.

a Significant gene expression correlation in AT2 cells between TMPRSS2 and ACE2, BSG (CD147) and NPR1 in COPD and IPF samples, each dot represents the average expression level of the genes of interest per sample, pairwise gene correlation analysis was done using a fitting linear model and p value was calculated using Anova. b Boxplot shows differences in gene expression of selected SARS-CoV-2 response genes in the AT2 cell types among different diagnosis groups, Boxes: interquartile range, lower and upper hinges correspond to the first and third quantiles, upper and lower whisker extends from the hinge to the largest values or smallest values of 1.5 × interquartile range; **p value-adj ≤ 0.05 (negative binomial test, corrected for Age, Ethnicity, Smoking_status and Dataset). c Upset plot shows shared differential expression genes (DEGs) between different comparisons: ACE2− CLD vs. Control, ACE2 + CLD vs. Control, CLD ACE2 + vs. ACE2-, Control ACE2 + vs. ACE2- and ACE2 correlated genes in the AT2 cells. d Upregulation of two genes uniquely differentially expressed in the CLD ACE2 + vs. ACE2−. e Spearman gene correlation analysis identified genes correlated with ACE2 expression in AT2 ACE2 + cells in different diagnosis groups, p-value was adjusted using Benjamini-Hochberg corrections, dashed lines indicate the 99th percentile of Spearman rho values.

Fig. 4. Analysis of SARS-CoV-2 candidate immune response genes in immune cells.

a Quantification of cell types as a percent of all immune cells in control and diseased lungs, numbers represent the total numbers of the cell type per individual samples, data are presented as mean values ± SEM. b Binary heatmap representing a manually curated list of genes associated with SARS-CoV-2. Orange elements indicate genes with increased expression and white elements indicate genes with decreased expression; Not detected: gene expression was not detected in either of the two tested populations (CLD vs. Control); DEGs with FDR ≤ 0.1 are outlined in black. c Differential expression analysis for SARS-CoV-2 immune candidate genes in cDCs, Macrophages and Monocytes. *p-adjusted value < 0.1, **p-adjusted value < 0.05, p-adjusted value was Bonferroni adjusted from Seurat FindMarkers differential expression analysis using a negative binomial test and corrected for Age, Ethnicity, Smoking_status and Dataset. d Compared to the healthy control samples, HLA type II score is higher in all disease groups (especially Other-ILD). In the T cell population, cytotoxicity scores (e) and exhaustion scores (f) are higher in the disease samples than in control samples. In a, d, e, and f: Boxes: interquartile range, lower and upper hinges correspond to the first and third quantiles, upper and lower whisker extends from the hinge to the largest values or smallest values of 1.5 x interquartile range; Tukey_HSD post-hoc test: *p value < 0.05, **p value < 0.01, ***p value < 0.001, ****p value < 0.0001. See Supplementary Fig. 10 for plots with outliers included for d–f.

Fig. 5. Model of alterations in the diseased lung related to SARS-CoV2 pathogenesis.

(1) In the IPF lung, there is a proximalization of the distal airway. ACE2 + epithelial cells cluster in the small airways though total ACE2 + cell numbers are similar to control. (2) The viral entry score (accounting for all described putative receptors and proteases) is increased in diseased lungs. (3) Diseased epithelial cells have alterations in key SARS-CoV-2 response genes/pathways. (4) In the CLD lung, there is increased expression of cytotoxicity and exhaustion genes in immune cell populations and alterations in viral response pathways (interferon, antigen presentation). Figure created in Biorender.com.