Identification of Small Airway Epithelium-Related Hub Genes in Chronic Obstructive Pulmonary Disease

Introduction

Chronic obstructive pulmonary disease (COPD) is a chronic inflammatory disease characterized by progressive airflow obstruction, mucus hypersecretion, and chronic airway inflammation,^1,2 accounting for the fourth leading cause of death (126,000 deaths per year). Long-term exposure to cigarette smoke may induce an airway epithelium inflammatory response and alveolar wall destruction in COPD.³

Pulmonary small airway epithelia are the primary site of cellular and histological alterations in COPD,⁴ which are mainly lined by pseudostratified epithelial cells including secretory club cells, ciliated cells and basal cells.^5,6 It has been reported that airway epithelial cells exposed to tobacco might induce inflammation, oxidative stress, and protease/antiprotease imbalance.⁷ Emerging evidences have indicated that inflammation and imbalance of antioxidant defense systems were the major predisposing factors in the pathogenesis of COPD and acted as excellent targets for COPD therapies.^8,9 Therefore, essential biomarkers of pulmonary epithelium are urgent to be identified to elucidate profound alterations in the progression of the disease.

Previous studies have reported a correlation between susceptibility to COPD and gene expression.¹⁰ High-throughput microarray and in silico analysis may provide gene-expression profiling to identify differentially expressed genes (DEGs) associated with COPD patients.¹¹ Lin et al have ever identified candidate genetic biomarkers in COPD process via in silico analysis.¹² Herein, we applied an integrated analysis of GEO dataset containing COPD and control subjects by multiple bioinformatics approaches to investigate the potential biomarkers and molecular mechanisms for COPD.

In the current study, DEGs were preliminarily identified via “limma” package and WGCNA. Moreover, protein-protein interaction (PPI) network and Cytoscape were constructed to identify cluster modules related to COPD. The expressions and distributions of hub genes associated with clinical phenotype of COPD were identified via single-cell analysis and further confirmed in vitro. Our study may further promote the understanding of the potential molecular mechanisms of pulmonary small airway epithelium in COPD.

Materials and Methods

Dataset Preparation

The raw gene expression profile datasets were screened from Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo). The datasets search was performed using the keywords “airway epithelium or airway epithelia”, “COPD”, “smoking”. The selection criteria were as follows: 1) Pulmonary airway epithelium samples from healthy non-smoker and COPD patients were included; 2) All datasets contained clinical information of healthy and COPD patients, such as age, gender, ethnic, smoking status; 3) Datasets contained at least 30 normal and COPD cases. Based on the above criteria, the GSE11906 and GSE128708 datasets (containing 84 healthy non-smokers and 124 COPD samples) were obtained. The characteristic of the datasets was presented in Table S1. GSE11906 contained small airway epithelium obtained by fiberoptic bronchoscopy of 34 healthy non-smokers and 33 COPD patients. GSE11906 was used to screen DEGs and for WGCNA. GSE128708 was used for external validation of hub genes. In addition, single-cell RNA sequencing data of three COPD samples (GSE167295) and three control donor lungs (GSE173896) were collected to verify the distribution, expression and potential function of hub genes in COPD.

DEGs Screening and Enrichment Analysis

The R package “limma” was used to identify DEGs between healthy non-smoker samples and COPD airway epithelium samples. Genes with an adjusted P value <0.05 and log₂ fold-change >2 were considered as being differentially expressed. The volcano plot was generated to visualize the distribution of the identified DEGs. The enrichment analysis of DEGs was executed by the KEGG and on the basis of the Hallmark gene sets collection in MsigDB databases (http://software.broadinstitute.org/gsea/msigdb/).

WGCNA

To identify whether the hub genes were related to the clinical features of healthy subjects and COPD patients, WGCNA of GSE11906 dataset was performed by using R package.¹³ Firstly, Pearson’s correlation analysis of all DEGs was conducted to construct an adjacency matrix. Subsequently, the adjacency matrix was transformed into a topological overlap matrix (TOM) when the power of β = 7. Genes with similar gene expression profiles were placed in the same module. A cluster dendrogram was created based on the TOM with a minimum cluster size of 30.

PPI Network Construction and Hub Gene Confirmation

After screening for genes associated with the response to GOLD classification based on module connectivity, the genes were uploaded to the Search Tool for the Retrieval of Interacting Genes (STRING) database to study the interactions among DEGs. The interaction data for the above-mentioned genes were downloaded from STRING and analyzed using Cytoscape (version 3.5.1) to construct a PPI network. Furthermore, CytoHubba software was used to identify the hub genes in the PPIs and the analysis was performed using the default parameters.¹⁴

scRNA-Seq Data Processing and Analysis

Single-cell analysis was performed to validate the expression and distribution of hub genes between healthy and COPD patients. In GSE167295 and GSE173896 datasets, Seurat R package v4.0 was applied to perform quality control, normalization, integration of samples, batch correction and cell clustering. Combined with the CellMarker database, single R package and relevant literature, cell clusters were annotated.^15,16 Analysis of differentially expressed genes and Partition-based graph abstraction (PAGA) were performed using R package and OmniAnalyzer software.

Cigarette Smoke Extract Preparation

Cigarette smoke extract (CSE) was prepared as follows: Two non-filtered cigarettes (containing 0.8 mg of nicotine and 10 mg of tar per cigarette, Shishi, China) were burned with a vacuum pump. Smoke was slowly bubbled through the 20 mL serum-free medium to collect CSE. The media was regarded as the 100% CSE solution. Subsequently, the pH of CSE was adjusted to 7.4 and filtered through a 0.22-μm filter (MilliporeSigma, Burlington, USA). The CSE solutions were stored at -80 °C until use.

Cell Culture

BEAS-2B cell line was purchased from the American Type Culture Collection (ATCC, www.atcc.org). Cells were cultured in DMEM medium containing 10% fetal bovine serum (FBS, Invitrogen, USA) and 1% penicillin-streptomycin (Gibco, CA) at 37 °C in the atmosphere containing 5% CO₂. The ethical approval of cell line application is not needed because the cell lines are available commercially.

CCK-8 Assay

CCK-8 assays were applied to assess cell viability (Meilun Biotechnology, Dalian, China). BEAS-2B cells were seeded in 96-well plates at a density of 4.0×10³ cells/well overnight and treated with different concentrations of CSE (0, 2.5%, 5%, 10%, 20%, 40%) for 24h. After incubation, 100µL mixture of CCK-8 and serum-free medium at a volume ratio of 1:10 was added and incubated for additional 1h at 37°C. The optical density was measured at 450 nm using a microplate reader (Thermo Fisher Scientific). Cell inhibition rate (%) = [(Ac−As)/(Ac−Ab)] × 100; As=absorbance of the experimental well, Ab=blank well absorbance and Ac=control well absorbance.

Reverse Transcription‑quantitative PCR (RT‑qPCR)

RNA was extracted using TRIzol^® reagent (Invitrogen, Thermo Fisher Scientific, Inc.) according to the manufacturer’s protocol. After quantification and RNA purity were confirmed, the extracted RNA was reverse transcribed into cDNA using Reverse Transcription Kit (Takara). Then, RT-qPCR was conducted on ABI 7300 Real-Time PCR system (Applied Biosystems, Foster City, CA, USA) with a SYBR Green PCR Kit (Thermo Fisher Scientific, Inc.) according to the following primer sequences. The corresponding gene primer sequences were listed in Table 1. The relative gene expression level was measured using the 2^−ΔΔCt method.

Table 1 Primer Sequences Information

Statistical Analysis

The non-parametric Wilcoxon rank-sum test was used to examine the relationship of continuous variables between the two groups. Statistical difference between control and CSE groups was conducted with Student’s t-test. Pearson’s correlation analysis was performed to identify the correlation between the two variables. All analyses were conducted in R software package (version 4.1.1). Statistical significance was defined when the P value was less than 0.05.

Results

Identification of DEGs in GSE11906

The overall study design was presented in Figure 1. The microarray dataset was obtained from GEO database to analyze the DEGs of COPD epithelium samples with a threshold of P <0.05 and log2 fold change >2. As shown in the volcano plot, 1130 DEGs (1099 up-regulated and 31 down-regulated) were identified in the GSE11906 dataset (Figure 2A). Heatmap of the top 20 upregulated and downregulated DEGs was shown in Figure 2B. GO analysis indicated that the DEGs were mainly enriched in intrinsic/integral component of membrane, endomembrane system, quinone metabolic process, alcohol dehydrogenase (NADP+) activity, polyketide metabolic process (Figure 2C). KEGG analysis demonstrated that the DEGs were mainly enriched in metabolic pathways, cytokine-cytokine receptor interaction, arachidonic acid metabolism, metabolism of xenobiotics by cytochrome P450 (Figure 2D). Therefore, we hypothesize that the epithelial DEGs might participate in the pathogenesis of COPD by regulating the above biological functions and pathways.

Figure 1 Flow chart of the study design.

Figure 2 Identification of DEGs in GSE11906. (A and B) Volcano plot and heat map of DEGs. (C and D) GO functional annotation and KEGG enrichment analysis of DEGs.

WGCNA Construction and Gene Module Identification

WGCNA was applied to analyze gene co-expression modules based on the clinical characteristics and gene expression profiles. As shown in Figures 3A, soft threshold power of β = 7 was selected in our study and the corresponding R² was 0.86, which satisfied scale-free topology of network. After removing outliers, the sample dendrogram and trait heatmap of COPD and healthy non-smoker subjects were drawn (Figure 3B). A total of nine modules were identified via hierarchical clustering analysis of adjacency value difference, and each module was presented using a different color in Figure 3C. As shown in Figure 3D, the midnightblue module was found to have a strong connection with COPD GOLD classification. Figure 3E illustrated the correlation between gene significance of COPD, GOLD and smoking pack-year and module membership in midnightblue module, which can be seen as highly correlated.

Figure 3 WGCNA network and module detection. (A) Scale Independence and mean network connectivity of soft-threshold power. A soft threshold of 7 is the optimal value. (B) Hierarchical clustering tree based on WGCNA. (C) Cluster dendrogram of differential genes represented the co-expression modules identified by WGCNA. (D) Module-trait weighted correlations and corresponding P values for the identified gene module and clinical features. (E) Correlation between gene significance of COPD, GOLD and smoking pack-year and module membership.

Construction of PPI Network and Screening of Hub Genes

The 98 key genes in midnightblue module were uploaded to STRING to analyze their interactions and the PPI network was presented in Figure 4A. Subsequently, we downloaded the genes interaction data and uploaded it to Cytoscape tool - Cytohubba plugin for screening hub genes. As shown in Figure 4B, we eventually obtained the top eight hub genes (AKR1C3, ALDH3A1, AKR1C1, CYP1A1, GPX2, CBR3, AKR1B1, GSR). Furthermore, the Pearson correlation coefficient analysis was calculated across these hub genes. The results indicated that there was significant correlation among these key genes (Figure 4C). According to the above results and literature reported, all of these hub genes were involved in antioxidant related actions. Therefore, we investigated the relationship between these hub genes and oxidative stress-related genes. As shown in Figure 4D, we found that the eight hub genes were significantly correlated with NQO1, RPE, ANG, ATM, CYP2E1, G6PD. Consistently, the results of Hallmark gene sets indicated that the eight genes were significantly correlated with reactive oxygen species pathway, apoptosis and fatty acid metabolism (Figure 4E).

Figure 4 Identification of eight hub genes. (A) Protein-protein interaction network of 98 genes correlated with COPD clinical features were obtained from STRING database. (B) Top eight hub genes identified by Cytohubba. (C) Pearson correlation coefficient analysis of eight hub genes. (D) The relationship between hub gene and oxidative stress related genes. (E) Enrichment plot of Hallmark gene sets with eight hub genes.

Single-Cell Transcriptomic Analysis

To further verify the relationship between the hub genes and airway epithelium in COPD, we obtained single-cell RNA-seq data from the GSE173896 and GSE167295 dataset, which comprised the data from three normal and three COPD tissues. The clinical characteristics of the datasets were shown in Tables S1 and S2. The dataset contained 12,692 cells from normal and 12,502 cells from COPD lung tissues with 13 cell types identified based on representative marker genes (Figure 5A and B). As shown in Figure 5C, we found that the percentage of T lymphocytes subset in the COPD group (42.67%) was significantly higher than that in the control group (19.19%), while there was no significant difference of epithelial cells between normal (9.25%) and COPD groups (7.99%). We further annotated and identified the epithelial cells into six cell sub-clusters, including ciliated cells, club cells, basal cells, goblet cells, alveolar epithelial cell type 1 and alveolar epithelial cell type 2 (Figure 5D). We then used the violin plot to illustrate the expressions of AGER, LAMP3, KRT5, FOXJ1, SCGB1A2, MUC5B in each subset, which further validated the correctness of identification in these epithelial cell subsets (Figure 5E).

Figure 5 Profiling of cell types in COPD using scRNA seq data. (A) Uniform Manifold Approximation and Projection (UMAP) representation of 25,194 cells from 3 COPD and 3 control donor lungs. (B) Visualization of clustering by UMAP plot with 13 distinct cell types identified. (C) Relative proportion of each cell subset between normal and COPD tissues. (D) UMAP plot of epithelial cells colored according to identified clusters. (E) Violin plots of normalized expression values for canonical cell specific marker genes for alveolar epithelial cell type 1, alveolar epithelial cell type 2, basal cells, ciliated cells and goblet cells.

Subsequently, we explored the enrichment and expressions of the above eight hub genes in each subtype of epithelial cells. However, only six hub genes, including AKR1C3, ALDH3A1, GPX2, CBR3, AKR1B1 and GSR, were expressed in the single-cell dataset. In the Figure 6A, we found that the expressions of AKRAC3, ALDH3A1, GPX2, CBR3 and AKR1B1 were significantly increased in the COPD group when compared with the normal group, while GSR expression showed the inverse trend. Further analysis of epithelial cell subpopulations demonstrated the different expression levels of these hub genes in different epithelial sub-clusters. Interestingly, we found that only ALDH3A1 was higher expressed in all subtypes of epithelial cells (Figure 6B). Gene-wise correlation heatmap indicated that ALDH3A1 was closely related to GSR and AKR1C3 (Figure 6C). The above evidence indicated the hub genes were differentially expressed in epithelial cells of the COPD group. Hence, we further illustrated the distribution and expressions of these hub genes in epithelial sub-clusters of the COPD group. Notably, we found that the expression of ALDH3A1 was enriched in ciliated cells (Figure 6D), which indicated that ALDH3A1 may exert important function on ciliated cells in COPD patients. To further explore the topology of epithelial cells, we applied PAGA, an algorithm that maps discrete connected and continuous connected cell-to-cell variation.¹⁷ In the PAGA plots, we found that ciliated cells were most closely related to basal cells when ALDH3A1 was expressed (Figure 7A and B). Moreover, according to the expression of ALDH3A1, ciliated cells were further subdivided into two subtypes. We conducted differential gene expression analysis between ALDH3A1⁺ ciliated cells and ALDH3A1⁻ ciliated cells (Figure 7C). KEGG pathway enrichment analysis showed that these DEGs were significantly enriched in P53 signaling pathway (Figure 7D).

Figure 6 Hub genes expressions in COPD epithelial cells. (A) Hub genes expressions in epithelium were assessed by scRNA-seq from COPD and non-smokers. (B) Hub genes expressions in pulmonary epithelial cell subpopulations. (C) Gene-wise correlation heatmap of six hub genes. (D) Distribution and expression levels of hub genes in epithelial cells of COPD patients.

Figure 7 Single-cell RNA transcriptomic analysis of ALDH3A1 in epithelium. (A and B) Partition-based graph abstraction (PAGA) indicated ALDH3A1 was related to ciliated cells. (C) volcano plots of DEGs between ALDH3A1+ ciliated cells and ALDH3A1- ciliated cells. (D) KEGG pathway enrichment analysis of DEGs in ALDH3A1 ± ciliated cells.

Verification of Hub Genes

To assess the robustness of hub genes, we downloaded another dataset (GSE128708) for verification. The results demonstrated that these eight hub genes were markedly increased in epithelium of the COPD group when compared with the normal group (Figure 8A). Moreover, we detected their mRNA expression levels in human pulmonary epithelial cells lines (BEAS-2B cells) after CSE incubation for 24h. CSE IC50 was further calculated, and the result was presented in Figure 8B. The IC50 was approximately 3% CSE. Subsequently, we applied the concentration to BEAS-2B and detected the mRNA expression of the hub genes. The results indicated that the expressions of AKR1C3, ALDH3A1, AKR1C1, CYP1A1, GPX2, CBR3, GSR were significantly elevated when the cells were exposed to CSE, while the expression of AKR1B1 was downregulated (Figure 8C).

Figure 8 Verification of hub genes. (A) External validation of eight hub genes in GSE128708. (B) IC50 of cigarette smoke extract. (C) mRNA expression levels of eight hub genes were detected via RT-qPCR (*P<0.05, **P<0.01, ***P<0.001, ****P<0.0001).

Discussion

Vascular inflammation and oxidative stress are implicated as essential factors in the pathogenesis and progression of COPD.^18,19 Previous studies have reported that cigarette smoke exposure may induce redox imbalance by shifting the balance between oxidants and antioxidants.²⁰ The airway epithelium is the first protective barrier between the pulmonary tissues and air, making it extremely susceptible to the stimulation of harmful substances, such as PM2.5 and smoke.^21,22 Therefore, the present in silico analysis was conducted to identify hub genes in pulmonary airway epithelium exposed to cigarettes in COPD patients.

In the present study, we preliminarily identified DEGs of GSE11906 and eventually confirmed eight hub genes associated with airway epithelial cells via WGCNA and PPI network, including ALDH3A1, AKR1C3, AKR1C1, CYP1A1, GPX2, CBR3, AKR1B1, GSR. The above molecules were regarded as metabolic regulators that participate in catabolic redox reactions, especially ALDH3A1.^23–26 The KEGG enrichment result also supported the evidence that the eight hub genes were enriched in reactive oxygen species (ROS) pathway and fatty acid metabolism. Consistently, the correlation coefficient map demonstrated that the hub genes were correlated with essential molecules in the regulation of oxidative stress, including SOD1, NQO1, G6PD.^27–29

The underlying mechanism of the above candidate genes in COPD remains to be investigated. Previous studies have reported that AKRs, including AKR1C3, AKR1C1 and AKR1B1, might catalyze the reduction of lipid peroxidation products, such as 4-hydroxynonenal (4-HNE), and their glutathione (GSH) conjugates.^30–32 Moreover, ALDH metabolized other aldehydes generated during the process of oxidative stress, including 4-HNE.³³ Reactive 4-HNE was the toxic products of lipid peroxidation,³⁴ which acted as the oxidative degradation of polyunsaturated fatty acids and may further induce cellular apoptosis in COPD pathogenesis.³⁵ GPX2 is believed to regulate cancer progression by regulating the level of hydroperoxides inside cells.³⁶ In addition, its antioxidant effect may protect epithelial tissue from harmful substances, such as hydrogen peroxide, free radicals, tumor necrosis factor.³⁷ Vecchio et al reported that cigarette smoke-exposed mouse was more prone to show a higher baseline production of GPX2 in alveolar macrophages,³⁸ which was consistent with our results. Additionally, CYP1A1 induced ROS production and was also correlated with the generation of lipid peroxidation products.^39,40

The specific expression distribution and underlying mechanism of ALDH3A1 have not been fully elucidated. Herein, we found that ALDH3A1 served as an essential role in the oxidative stress of COPD. ALDH3A1 is a metabolic enzyme responsible for the oxidation of various aldehydes.⁴¹ Previous studies have indicated that ALDH3A1 was associated with external environmental stress, and the metabolic role of ALDH3A1 was related to resistance to oxidative stress. Additionally, the expression of ALDH3A1 robustly increased in the fluid lining the epithelial cells of smokers or primary human bronchial epithelial cells after exposing to cigarette smoke.^42,43 In addition, KEGG enrichment analysis demonstrated that ALDH3A1-associated DEGs in ciliated cells were significantly enriched in P53 pathway. Similarly, previous studies also found that ALDH3A1 was regulated by P53 protein in HBE cell lines.⁴⁴ Therefore, we speculated that ALDH3A1 might exert underlying function of metabolism and oxidative stress on ciliated cells through P53 pathway. However, whether ALDH3A1 played a potential role via P53 pathway needed further experimental verifications.

The present study revealed candidate genes on the metabolism and oxidative stress of COPD. While we have to acknowledge that our study has several limitations. Firstly, the sample size was relatively small and the clinical information was not sufficient, resulting in a less profound investigation of WGCNA. Secondly, although we performed RT-qPCR to validate the expressions of the eight hub genes in BEAS-2B cell stimulated by CSE, further biological function and potential mechanism have not been verified. Therefore, the specific molecular regulatory mechanism and how these genes interact to regulate corticosteroid resistance still need to be clarified in future experiments.

Conclusion

In the present study, we identified eight hub genes related to the GOLD classification by integrating the microarray data of COPD, including AKR1C3, ALDH3A1, AKR1C1, CYP1A1, GPX2, CBR3, AKR1B1, GSR. Moreover, functional enrichment analysis indicated that the expressions of the above genes were significantly correlated with metabolism and antioxidant stress in COPD process. Single-cell transcriptomic analysis indicated that five hub genes, including AKR1C3, ALDH3A1, GPX2, CBR3 and AKR1B1, were elevated in pulmonary airway epithelium. Moreover, ALDH3A1 was significantly elevated in all subtypes of epithelial cells, especially in ciliated cells.