Molecular mapping of interstitial lung disease reveals a phenotypically distinct senescent basal epithelial cell population

Abstract

Compromised regenerative capacity of lung epithelial cells can lead to cellular senescence, which may precipitate fibrosis. While increased markers of senescence have been reported in idiopathic pulmonary fibrosis (IPF), the origin and identity of these senescent cells remain unclear, and tools to characterize context-specific cellular senescence in human lung are lacking. We observed that the senescent marker p16 is predominantly localized to bronchiolized epithelial structures in scarred regions of IPF and systemic sclerosis-associated interstitial lung disease (SSc-ILD) lung tissue, overlapping with the basal epithelial markers Keratin 5 and Keratin 17. Using in vitro models, we derived transcriptional signatures of senescence programming specific to different types of lung epithelial cells and interrogated these signatures in a single-cell RNA-Seq data set derived from control, IPF, and SSc-ILD lung tissue. We identified a population of basal epithelial cells defined by, and enriched for, markers of cellular senescence and identified candidate markers specific to senescent basal epithelial cells in ILD that can enable future functional studies. Notably, gene expression of these cells significantly overlaps with terminally differentiating cells in stratified epithelia, where it is driven by p53 activation as part of the senescence program.

Keywords: Aging; Cellular senescence; Fibrosis; Pulmonology.

Figure 1. Expression of the senescent marker p16 (CDKN2A) is induced in IPF and SSc-ILD, where it localizes to basal epithelial cells in bronchiolized epithelium.

(A) Microarray analysis of CDKN2A gene expression (mean ± SD, control n = 8, IPF n = 40) in control and IPF lung tissue. ***P < 0.005 (unpaired 2-tailed Student’s t tests). (B) RNA-Seq data for CDKN2A gene expression (mean ± SD, control n = 4, IPF n = 10, SSc-ILD n = 3) in control, IPF, and SSc-ILD lung tissue explants. ***P < 0.005 (Tukey’s multiple comparisons test). nRPKM, normalized reads per kilobase gene model per million total reads. (C) Top 50 correlates with CDKN2A gene expression in control and IPF lung tissue, sample set from A. (D) H&E (top row) and immunohistochemical staining of control, IPF, and SSc-ILD lung tissue sections for p16 protein. Highlighted regions from H&E shown at higher magnification from immunostainings of serial sections. (E) Serial sections of lung tissue stained for p16 and KRT17 proteins. (F) Costaining of IPF lung tissue sections for p16/SP-C and p16/KRT5 proteins. Scale bars: 200 μm (H&E), 50 μm (enlarged) (D); 50 μm (E); 100 μm (F).

Figure 2. Establishing in vitro senescent lung epithelial cell models.

(A) SAECs stained for senescence-associated β-galactosidase activity (SAβ-Gal) at 5 days after treatment (control: DMSO/H₂O₂ at 50 nM, senescence-inducing: doxorubicin/H₂O₂ at 100 nM). (B) Immunofluorescence (IF) staining for Ki-67 and Actin in SAECs at 5 days after treatment. (C) Quantification of SAEC area from representative fields (mean ± SD, n = 10) after induction of senescence by doxorubicin (day 5). ***P < 0.005 (unpaired 2-tailed Student’s t tests). (D) CDKN2A (mean ± SD, n = 4 biological replicates) gene expression in nonsenescent (DMSO/H₂O₂ 50 nM) and senescent (doxorubicin/H₂O₂ 100 nM) cultures at day 5. ***P < 0.005 (unpaired 2-tailed Student’s t tests). (E) IF staining for p16 and actin in SAECs at day 5 after treatment. (F) SASP gene expression, CSF2, IL6, and CXCL8, (mean ± SD, n = 4 biological replicates) in nonsenescent (DMSO/H₂O₂ 50 nM) and senescent (doxorubicin/H₂O₂ 100 nM) cultures at day 5. ***P < 0.005 (unpaired 2-tailed Student’s t tests). (G) Quantification of secreted SASP proteins, GM-CSF, IL-6, IL-8, (mean ± SD, n = 3 biological replicates) in supernatants of treated SAEC cultures over a 48-hour period (day 5 to day 7 after initiation of treatment). ***P < 0.005 (unpaired 2-tailed Student’s t tests). (H) Relative cell number in SAEC cultures at day 4 cultured in SAGM basal media, conditioned media from control cultures (CM-Cont), or conditioned media from senescent SAEC cultures (CM-Sen) (mean ± SD, n = 3). **P < 0.05, ***P < 0.005 (Tukey’s multiple comparisons test). (I) SAβ-Gal stainings of SAEC cultures at day 4. (J) Quantification of percentage of senescent cells in SAEC cultures treated with conditioned media at day 4 (mean ± SD, n = 3). ***P < 0.005 (Tukey’s multiple comparisons test). Scale bars: 100 μm (A and I), 40 μm (B), 40 μm (E).

Figure 3. Cell type dictates senescent transcriptional phenotype.

(A) Four-way comparison of differentially expressed genes in senescent SAEC cultures, compared with DMSO control cultures, induced by doxorubicin versus H₂O₂. (B) Four-way comparison of differentially expressed genes, compared with DMSO control cultures, in doxorubicin-treated SAEC versus doxorubicin-treated BEC cultures. (C) Four-way comparison of differentially expressed genes, compared with DMSO control cultures, in doxorubicin-treated SAEC versus doxorubicin-treated NHBE cultures. (D) Four-way comparison of differentially expressed genes, compared with control cultures, in senescent IMR90 fibroblasts induced via irradiation versus senescent WI-38 fibroblasts after replicative senescence. (E) Four-way comparison of differentially expressed genes, compared with control cultures, in doxorubicin-treated SAEC versus senescent IMR90 fibroblasts induced via irradiation. (F) Overlap in upregulated gene expression among 4 senescent lung epithelial culture models. (G) Overlap in upregulated gene expression among 3 senescent fibroblast models. (H) Overlap between the consensus senescent lung epithelial signature (F) and consensus senescent fibroblast epithelial signature (G). IR, irradiation-induced senescence; RS, replicative senescence. (A–E) P = Spearman’s correlation, 95% confidence intervals in parentheses.

Figure 4. Consensus senescence signature identifies replicative senescence in primary human lung epithelial cell cultures.

(A) SAβ-Gal staining of BEC cultures at passages 2 and 9. (B) Percentage of SAβ-Gal⁺ cells (mean ± SD, n = 5 biological replicates) in BEC cultures at passages 2 and 9. ***P < 0.05 (unpaired 2-tailed Student’s t tests). (C) CDKN2A gene expression (mean ± SD, n = 3 biological replicates) in BEC cultures at P2 and P9 and P2 doxorubicin-treated cultures. ***P < 0.05 (Tukey’s multiple comparisons test). (D–F) Gene expression (mean ± SD, n = 3 biological replicates) of SASP, proliferation markers, and basal cell markers in BEC cultures at P2 and P9 and P2 doxorubicin-treated cultures. **P < 0.05, ***P < 0.005 (Tukey’s multiple comparisons test). (G) Uniform manifold approximation and projection (UMAP) plot showing scRNA-Seq cluster assignments in BEC culture at P4, with the overlaid black lines showing inferred trajectory between clusters as a minimum spanning tree with cluster E2 as the starting point. (H) Distribution of cells in each BEC culture cluster that were annotated as being in the G1, S, or G2/M phase of the cell cycle. (I) Violin plots showing signature scores of BEC culture clusters for both the epithelial-derived (EPI) and fibroblast-derived (FB) senescence-associated gene sets; the black horizontal bar denotes the mean expression value within the given cluster. (J) Violin plots showing the normalized expression of select senescence- and proliferation-associated genes in each BEC culture cluster; each point represents an individual cell, and the black horizontal bar denotes the mean expression value within the given cluster. Scale bars: 100 μm (A).

Figure 5. Characterization of the epithelial populations in IPF and SSc-ILD.

(A) UMAP plot of epithelial cell populations from human lung explants after clustering. PNEC, pulmonary neuroendocrine cells. (B) Percentage of cells in each cluster shown as individual boxplots for the control, IPF, and SSc-ILD samples; y axis is on a square root scale. The box plots depict the minimum and maximum values (whiskers), the upper and lower quartiles, and the median. The length of the box represents the interquartile range. (C) Dot plot of selected canonical cell type markers and marker genes identified by across-cluster differential expression; color (yellow to red) denotes mean normalized expression for the cluster; dot size denotes the percentage of cells within the cluster for which any expression was detected, with no dot shown for percentage detected values under 15%.

Figure 6. Identification of a senescent BEC population enriched in the fibrotic lung.

(A) Boxplots showing the distribution and mean signature scores for both epithelial cell–derived and fibroblast-derived consensus senescence gene sets; each point reflects the mean signature score for an individual subject within the given cluster. The box plots depict the minimum and maximum values (whiskers), the upper and lower quartiles, and the median. The length of the box represents the interquartile range. (B) Heatmap showing the log fold change in normalized expression between the Basal-1 and Basal-2 populations for genes in the union set of control, IPF, and SSc-ILD differentially expressed genes across these 2 clusters; genes upregulated in Basal-2 are shown in red, and genes upregulated in Basal-1 are shown in blue. (C) IHC staining of control, IPF, and SSc-ILD lung tissue sections for LY6D. (D) IHC staining of IPF lung tissue serial sections for LY6D and p16. (E) Scoring of clusters E1–6, from Figure 4, against Basal1/2 gene expression signatures. Scale bars: 100 μm (C), 40 μm (D).

Figure 7. Basal cell senescence is transcriptionally related to squamous terminal differentiation.

(A) Top correlates with LY6D gene expression in lung tissue from control and IPF patient explants. (B) Top correlates with LY6D across various tissues. (C) Phase images of normal human epidermal keratinocyte (NHEK) cultures under basal, differentiation-inducing, and senescence-inducing culture conditions (on left) and IF staining for Ki-67 in NHEK cultures in proliferating, differentiated, and senescent cultures (on right). (D) Comparison of differential gene expression between terminally differentiated keratinocytes and senescent keratinocytes compared with undifferentiated cultures (mean ± SD, n = 3). *P < 0.05 (ratio of means t test). (E) GSEA Hallmark pathway enrichment in the Basal-2 population based on differential gene expression with associated FDR q values and overlap enrichment k:K (# of overlapping genes: # of genes in pathway signature). (F) Gene expression of squamous and SASP markers in gene-edited NHBE cultures at 24 hours (mean ± SD, n = 3). ***P < 0.005 (unpaired 2-tailed Student’s t test). Scale bars: 100 μm (phase in C), 40 μm (IF in C).