Interpretable machine learning uncovers epithelial transcriptional rewiring and a role for Gelsolin in COPD

Abstract

Transcriptomic analyses have advanced the understanding of complex disease pathophysiology including chronic obstructive pulmonary disease (COPD). However, identifying relevant biologic causative factors has been limited by the integration of high dimensionality data. COPD is characterized by lung destruction and inflammation, with smoke exposure being a major risk factor. To define previously unknown biological mechanisms in COPD, we utilized unsupervised and supervised interpretable machine learning analyses of single-cell RNA-Seq data from the mouse smoke-exposure model to identify significant latent factors (context-specific coexpression modules) impacting pathophysiology. The machine learning transcriptomic signatures coupled to protein networks uncovered a reduction in network complexity and new biological alterations in actin-associated gelsolin (GSN), which was transcriptionally linked to disease state. GSN was altered in airway epithelial cells in the mouse model and in human COPD. GSN was increased in plasma from patients with COPD, and smoke exposure resulted in enhanced GSN release from airway cells from patients with COPD. This method provides insights into rewiring of transcriptional networks that are associated with COPD pathogenesis and provides a translational analytical platform for other diseases.

Keywords: COPD; Cell biology; Cytoskeleton; Pulmonology.

Figure 1. scRNA-Seq defines cell type profiles in a mouse model of COPD.

C57BL/6J mice (n = 3 per group; 5,000 cells per mouse) were exposed to 6 months of air or CS. Harvested lung tissue was dissociated into single cell suspensions and processed individually for scRNA-Seq as described. (A) Experimental workflow and data analysis using standard Seurat cell clustering versus LOVE functional gene clustering. (B) Uniform Manifold Approximation and Projection (UMAP) representation of Seurat cell clustering identifying 34 cell clusters among the air and smoke exposed groups, notated with color and number labels (n = 4 individually sequenced mice per group). (C) Canonical cell marker genes identify distinct cell populations within myeloid, lymphoid, stromal, epithelial, and endothelial cell clusters. Groups are separated by exposure groups: air (red) or smoke (teal), n = 3 mice per group. Gene expression is reported as z-score by color histogram.

Figure 2. CS results in transcriptional rewiring within epithelial cells in the mouse lung.

(A) Heat maps of z-scores showing the top 5 unique downregulated and upregulated genes between air and smoke exposure groups for each epithelial cell cluster. Groups are separated by exposure groups: air (red) or smoke (teal). Cell-type text color corresponds to the top 5 genes in the same color on the Y-axis label. Gene expression is reported as z-score by color histogram. (B and C) Epithelial cells were analyzed using a functional gene clustering model (LOVE) with functional gene clusters created for each exposure group. (B) Epithelial cells with air exposure, (C) Epithelial cells with 6 months smoke exposure. Functional groups were identified within air and smoke exposure groups. Gene nodes are shown as triangles, and predicted, experimentally validated interactor genes are shown as circles. Red-yellow shading represents the average gene expression intensity across all of the cells in that group. n = 3 mice per group.

Figure 3. CS exposure resulted in transcriptional rewiring in AT1 and AT2 cell populations.

scRNA-Seq data from mice exposed to 6 months of air or CS (n = 3 per group; 5,000 cells per mouse) were analyzed using SLIDE for dataset A. (A) Standalone significant latent (marginal) factors for AT1 cells are in teal and interacting latent factors are in red. Genes comprising each latent factor by cell type are reported in the table and network connectivity maps with genes that characterized CS exposure in red/squares and air in blue/circles. (B) Latent factor genes and network map for AT1 cells, (C and D) SLIDE latent, and interacting factors identified for AT2 cells with latent factor genes shown that characterize CS (red/squares) versus air (blue/circles).

Figure 4. CS exposure resulted in transcriptional rewiring in ciliated cells.

scRNA-Seq data from mice exposed to 6 months of air or CS (n = 3 per group; 5,000 cells per mouse) were analyzed using SLIDE for dataset A. (A) Standalone significant latent (marginal) factors for ciliated cells (CCs) are in teal and interacting latent factors are in red. (B) Genes comprising each latent factor for CCs are reported in the table and network connectivity maps with genes that characterized CS exposure in red/squares and air in blue/circles.

Figure 5. A Gsn-centric latent factor gene set can predict the smoke-exposed group within other epithelial cell types.

scRNA-Seq data from mice exposed to 6 months of air or CS (n = 3 per group; 5,000 cells per mouse) were analyzed using SLIDE for dataset B for AT2 cells. (A) Standalone significant latent (marginal) factors for AT2 cells are in teal and interacting latent factors are in red. Genes comprising each latent factor by for the AT2 cells are reported in the table and network connectivity maps with genes that characterized CS exposure in red/squares and air in blue/circles. Of note, no air-associated latent factors (blue/circles) were present in this analysis. Cross prediction analysis was completed for between dataset A and B to determine if the CS treatment group could be identified. Area under the curve (AUC) is reported. Statistical comparison by 2-tailed Student’s t test with Mann-Whitney test for data in B–D. P values are noted. (B) AT2 dataset B predicting the CS group from AT2 cells in dataset A, (C) AT2 dataset B predicting the CS group from ciliated cells in dataset B, (D) AT2 dataset A predicting the CS group from ciliated cells in dataset A.

Figure 6. Gsn is enriched and increases in ciliated cells in mouse lung after CS and human COPD lung.

Expression of Gsn was plotted from scRNA-Seq data from (A) Dataset A, epithelial cells isolated from mouse lung, 6 month smoke exposure model (n = 3 mice per group; 5,000 cells per mouse) and (B) human control versus COPD lung epithelial cells. Data are shown as mean gene expression in each group and fraction of cells with expression. Cells were isolated from lung samples from patients with COPD (GOLD stage IV, n = 6, ages 58–68, 5 males and 1 female) and normal nonsmoker donor controls (n = 4, ages 56–68, 3 males, 1 female). (C and D) Human lungs from people in the control group or patients with COPD was stained by IF and imaged on a confocal microscope n = 3–4 participants per group (8 images per participant). Data represent normalized mean grey value ± SEM. Statistically significant P values are noted. Statistics by 2-tailed Student’s t test with Mann-Whitney post test. Representative images are shown for (C) airway epithelium stained for GSN (magenta) and EPCAM (green). Scale bars: 50 μm (left) and 10 μm (right). (D) Alveolar epithelium stained for GSN (magenta) and HT2-280 (green). Scale bar: 100 μm. GSN staining intensity was quantified in HT2-280–negative, podoplanin-positive alveolar tissue. Mean grey intensities (per measured area) were normalized to the healthy control. Group data were split into high (Hi) and low (Lo) groups by using the mean value for the healthy control group.

Figure 7. GSN is released by airway epithelial cells due to CS exposure and levels are increased in the plasma of patients with COPD.

Primary human ALI cultures were exposed to gaseous CS by Vitrocell with n = 2 donors per group (biological replicates) and n = 2–9 inserts per donor (technical replicates). GSN quantity was determined by dot blot from (A) ALI cell lysates, (B) basal chamber media, and (C) apical chamber media. Data represent mean ± SEM. Statistically significant P values are noted. Statistics by 1-way ANOVA with Kruskal-Wallis post test. (D) GSN concentrations in human plasma were determined by ELISA for nonsmokers (n = 10), smokers without COPD (n = 38) and participants with COPD (n = 154). Data represent mean ± SEM. Statistically significant P values are noted. Statistics by 1-way ANOVA with Kruskal-Wallis post test. (E) Plasma GSN concentration compared with fibrinogen concentration in the same participant. Statistically significant P values are noted. Statistics by 1-way ANOVA with Kruskal-Wallis post test. Linear regression with P = 0.040. Grey, controls; Blue, participants with COPD.

Figure 8. Endogenous overexpression of GSN and exogenous GSN alter cellular migration, proliferation, and KRT8 dynamics.

Beas-2b cells were analyzed for wound healing, cellular proliferation, and cytokeratin expression in the context of GSN overexpression (GSN OE) or exogenous recombinant human GSN (rhGSN). (A) Wound healing assay with GSN OE. Wounds were measured at time 0 and 12 hours. Percent wound healing was calculated. Representative images are shown. Data represent n = 3–6 per group technical replicates and 2 separate experimental replicates. Statistical comparison by 2-tailed Student’s t test with Mann-Whitney test with P values noted. Images were acquired at 4X magnification. (B) Western blot for KRT8 and KRT18 was performed on Beas-2b cells with GSN OE compared with control, with or without 10% CSE exposure (24 hours). Data are representative of 2 separate experiments (n = 3 per group per experiment). (C) Quantification of Western blot band intensity for KRT8 and KRT18. Normalized data are representative of 2 separate experiments (n = 6 per group total). Statistics by 1-way ANOVA with Kruskal-Wallis post test with P values noted. (D) Cell proliferation was determined using the Cyquant assay in cells exposed to rhGSN at 10 or 30 μg/mL compared with media or BSA (30 μg/mL). Cells were assessed at 48, 72 and 96 hours. Simple linear regression analysis was performed. P values represent (*P = 0.0004, **P = 0.0145). (E) Wound healing assay was performed on Beas-2b cells treated with rhGSN at 10 or 30 μg/mL compared with cells treated with media or BSA (30 μg/ml). Wounds were measured at 0 and 12 hours. Percent wound healing was calculated. Representative images are shown. Data represent n = 6 wells per group. Statistical comparison by 1-way ANOVA with Kruskal-Wallis post test. (F) Beas-2b cells were exposed to rhGSN at 30 μg/mL for 24 hours with subsequent RT-PCR for KRT8, KRT18, and ACTA2. n = 3 biological replicates per group. Statistical comparison by a parametric 2-tailed Student’s t test.