Identification and Validation of CDKN1A and HDAC1 as Senescence-Related Hub Genes in Chronic Obstructive Pulmonary Disease.

Purpose: Cellular senescence participates in the occurrence and development of chronic obstructive pulmonary disease (COPD). This study aimed to identify senescence-related hub genes and explore effective diagnostic markers and therapeutic targets for COPD.
Methods: The microarray data from the GSE38974 dataset was downloaded from the Gene Expression Omnibus (GEO) database. The overlapping genes between genes from the GSE38974 dataset and CellAge database were considered differentially expressed senescence-related genes (DESRGs). Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed using R software. Protein-protein interaction (PPI), miRNA-mRNA network, and competitive endogenous RNA (ceRNA) network were constructed and visualized by Cytoscape software. GSE100281 and GSE103174 datasets were employed to validate the expression and diagnostic value of hub genes. Reverse transcription-quantitative polymerase chain reaction (RT-qPCR) was used to measure the mRNA levels of hub genes in peripheral blood mononuclear cells (PBMCs) from COPD and control samples.
Results: A total of 23 DESRGs were identified between COPD samples and healthy controls. Enrichment analysis revealed that DESRGs were mainly related to apoptosis and senescence. Moreover, four hub genes and two key clusters were acquired by Cytohubba and MCODE plugin, respectively. CDKN1A and HDAC1 were verified as final hub genes based on GSE100281 and GSE103174 datasets validation. The mRNA expression level of CDKN1A was negatively related to forced expiratory volume in 1 second/forced vital capacity (FEV1/FVC), and HDAC1 expression had the opposite correlation. Finally, an HDAC1-based ceRNA network, including 6 miRNAs and 11 lncRNAs, was constructed.
Conclusion: We identified two senescence-related hub genes, CDKN1A and HDAC1, which may be effective biomarkers for COPD diagnosis and treatment. An HDAC1-related ceRNA network was constructed to clarify the role of senescence in COPD pathogenesis.

Introduction

Chronic obstructive pulmonary disease (COPD), ranking the third leading cause of death all over the world, is the fifth disease burden around the world.1 A recent epidemiological survey has shown that 13.6% of Chinese adults suffer from COPD.2 Early diagnosis and effective treatment are crucial to reducing COPD mortality. Currently, the confirmed diagnostic method is the lung function test. The major pharmacological treatment is bronchodilators and inhaled corticosteroids (ICSs). However, the early diagnosis rate and effect of treatment are unsatisfactory. Thus, more effective diagnostic markers and treatment targets need to be further developed.

Cell senescence characterized by permanent cell cycle arrest can be induced by several signals, including oxidative stress, mitochondrial dysfunction, and telomere dysfunction.3,4 As an aging-related disease, COPD is closely related to cellular senescence. Recent studies have revealed that p16, a marker of cellular senescence, was up-regulated in cigarette smoke extract (CSE) treated endothelial progenitor cells (EPCs) and COPD patients.5,6 Senescence-associated secretory phenotype (SASP), characterized by secretion of proinflammatory cytokines, chemokines, and matrix metalloproteinases (MMPs), is increased in COPD.7 SIRT6, a member of the sirtuin family, class III histone deacetylases (HDACs), is an anti-aging factor.8 A study has indicated that the expression level of SIRT6 was decreased in CSE-stimulated human bronchial epithelial cells (HBECs), and overexpression of SIRT6 could reverse cellular senescence induced by CSE.9 Therefore, we predicted that senescence-related genes (SRGs) would be effective diagnostic markers or therapeutic targets in COPD.

Non-coding RNAs, including miRNAs and lncRNAs have been proved to be associated with COPD pathogenesis.10,11 LncRNAs act as miRNAs sponges that can regulate gene expression by competitively binding miRNAs, and this interaction among mRNAs, lncRNAs and miRNAs is also called competing endogenous RNA (ceRNA).12 Previous studies have revealed that ceRNA plays a pivotal role in COPD. The expression of lncRNA OIP5-AS1 is elevated in CSE-stimulated 16HBE cells and can enhance the inflammation and cell apoptosis by modulating the miR-410-3p/IL-13 signaling.13 LncRNA GAS5 functions as a ceRNA and regulates the miR-223-3p/NLRP3 axis to promote cell pyroptosis in COPD.14 However, the role of SRGs-based ceRNA in COPD has not been investigated.

In this study, for the first time, we elucidated the relationship between SRGs and COPD by the bioinformatics method. We screened differentially expressed senescence-related genes (DESRGs) based on the Gene Expression Omnibus (GEO) and CellAge database. Then, we performed enrichment analysis, protein-protein interaction (PPI) network analysis, and correlation analysis. Furthermore, hub genes were acquired, and another two separate datasets were used for validation. The expression levels of hub genes in peripheral blood mononuclear cells (PBMCs) of COPD and control samples were further measured by reverse transcription-quantitative polymerase chain reaction (RT-qPCR). ceRNA among hub genes, miRNAs, and lncRNAs was also constructed. The flow chart of this study is shown in Figure 1. We believe that these results will provide a novel platform for exploring COPD pathogenesis and improving the prognosis of COPD patients.

Figure 1

The flowchart of the whole study.

Materials and Methods

Data Acquisition

In this study, three mRNA expression data of lung tissues were downloaded from the GEO datasets (). The keywords for searching were “chronic obstructive pulmonary disease”, “expression profiling by array”, “Homo sapiens”, and “tissue”. A GPL4133-platform dataset, GSE38974, which contained 9 healthy samples and 23 COPD samples, was identified as the test dataset.15 GSE103174 and GSE100281 were selected as validation datasets. GSE103174, which contained 11 healthy samples and 23 COPD samples, was acquired from the GPL13667 platform, and GSE100281, which included 16 healthy samples and 79 COPD samples, was acquired from the GPL11532 platform.16 A total of 279 senescence-related genes (SRGs) were obtained from the CellAge database ().17

Differentially Expressed SRGs in COPD and Correlation Analysis

GEO2R was employed to analyze the gene expression between the control and COPD groups in the GSE38974 dataset. SRGs between the control and COPD groups were acquired by intersecting genes from the GSE38974 dataset and CellAge database. The cut-off criteria of DESRGs was an adjusted p-value of <0.050 and log fold change (FC) of >1. Spearman correlation analysis was used to conduct the correlation analysis of DESRGs.

Functional Enrichment Analysis of DESRGs in COPD

We used the “cluster Profiler R” package () for Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis and the “GOplot R” package () for visualization. An adjusted p-value of < 0.050 indicated statistical significance.

PPI Network Analysis

The PPI network was predicted by the STRING database () by using an interaction score of > 0.4.18 Then, the interaction data were downloaded and visualized using Cytoscape software (version 3.9.1). The key clusters were obtained by Molecular Complex Detection (MCODE).19 Cytohubba was employed to screen important genes by several topological algorithms.20 In this study, the top six hub genes were obtained by five algorithms, including Degree, Maximal Clique Centrality (MCC), Maximum Neighborhood Component (MNC), Stress and EcCentricity. Then, the genes obtained from the above five algorithms were intersected to gain the final hub genes.

mRNA-miRNA Co-Expressed Network Construction

The Encyclopedia of RNA Interactomes (ENCORI, ) and miRWalk database () were used to predict the target miRNAs of 4 hub genes.21,22 A co-expressed network of mRNAs and miRNAs was constructed by Cytoscape software.

Clinical Samples Collecting and Preparation

Briefly, 5-mL venous blood samples were collected from 20 healthy controls and 24 COPD samples at Qilu Hospital of Shandong University. Thereafter, PBMCs were separated using isolation fluids (TBD sciences, Tianjin, China). The samples meeting the following criteria were included in the COPD group: the value of forced expiratory volume in 1 second/forced vital capacity (FEV1/ FVC) of < 0.7 after bronchodilator treatment and the age of 40 to 80 years. Healthy controls matched COPD samples in age and gender. The Medical Ethics Committee of Qilu Hospital of Shandong University approved the study procedure, which complied with the Declaration of Helsinki (as revised in 2013). All patients signed informed consent forms before the experiment.

RNA Isolation and Reverse Transcription-Quantitative Polymerase Chain Reaction (RT-qPCR)

The total RNA of PBMCs was extracted using RNAfast200 Reagent (Fastegen, Shanghai, China). Then, cDNA was synthesized using the PrimeScript RT reagent kit (Takara, Japan). RT-qPCR was performed using the TB Green Premix Ex Taq II (TaKaRa). ACTB was considered an endogenous control. The gene-specific primer sequences are listed in . The 2−ΔCT method was employed to quantify the mRNA levels of the genes.

ceRNA Construction

An HDCA1-related ceRNA network was constructed. ENCORI was employed to predict the target lncRNA of six miRNAs. The predicted results were intersected using the “UpSetR” package.23 Then, the ceRNA network among mRNA, miRNA, and lncRNA was visualized by Cytoscape software.

Statistics

Statistical analysis was performed using the R language (version 3.6.3), GraphPad prism 7.0 (GraphPad Software, CA, USA) and IBM SPSS Statistics 28 (USA). Student’s t-test and Mann–Whitney U-test were used to compare gene expression levels between the control and COPD groups. Fisher-Freeman-Halton exact probability test was used to analyze the correlation between the Global Initiative for Chronic Obstructive Lung Disease (GOLD) stage and other clinical characteristics. A p value of < 0.050 was considered statistically significant.

Results

Gain of DESRGs and Correlation Analysis

In this study, GSE38974 was used to analyze gene expression levels between the control and COPD groups. The principal component analysis (PCA) revealed that data of GSE38974 can be well repeated (Figure 2A). After GEO2R analysis, 19,287 genes were acquired. Figure 2B shows that 272 genes were obtained after the interaction with SRGs from the CellAge database. With the thresholds of the adjusted p-value of < 0.050 and log fold change (FC) > 1, 23 DESRGs, including 17 upregulated genes and 6 downregulated genes, were obtained (). Volcano plots and heatmaps were created to visualize DESRGs (Figure 2C and D). The top 5 significantly different genes were SNAI1, MMP9, ATM, YAP1, and NOTCH3. Meanwhile, expression levels of DESRGs between the control and COPD groups were explored and visualized in the form of a boxplot. Figure 3A shows that all the DESRGs were significantly different between the two groups in the GSE38974 dataset. Furthermore, Spearman correlation analysis was used to explore the correlation among DESRGs. The specific results are shown on the heatmap in Figure 3B. The red color represents positive correlation, and the blue color represents negative correlation. The value of p of < 0.050 was considered statistically significant.

Figure 2

Identification of DESRGs between the COPD and control samples in the GSE38974 dataset. (A) PCA plot of the GSE38974 dataset. (B) Venn diagram of overlapping genes between the GSE38974 dataset and CellAge database. (C and D) Volcano blot and heatmap of DESRGs between the COPD and healthy control samples.

Figure 3

Expression and correlation analysis of 23 DESRGs. (A) Expression levels of DESRGs in the COPD and healthy control samples visualized by box plot. ** p<0.01, ***p<0.001. (B) Spearman correlation among 23 DESRGs.

GO and KEGG Enrichment Analyses

GO and KEGG analyses were performed to explore the biological functions of DESRGs. The top 5 GO biological processes of DESRGs were intrinsic apoptotic signaling pathway (p = 1.83e-06), intrinsic apoptotic signaling pathway in response to DNA damage (p = 1.15e-04), regulation of apoptotic signaling pathway (p = 1.79e-04), signal transduction in response to DNA damage (p = 1.79e-04) and response to gamma radiation (p = 1.79e-04) (Figure 4A and B; Table 1). According to the results of the KEGG analysis, DESRGs were mainly involved in transcriptional misregulation in cancer, microRNA in cancer, endocrine resistance, the hypoxia-inducible factor (HIF-1) signaling pathway and cell cycle (Figure 4C and D). The results indicated that DESRGs were closely related to senescence and apoptosis.

Figure 4

GO and KEGG enrichment analyses of 23 DESRGs. (A) Bubble plot of enriched GO terms. (B) The circle plot showed the top 15 GO terms. The inner red circle represents the z-score values, and the outer circle represents the number of genes in the GO terms. Red indicates upregulated SRGs, and green indicates downregulated SRGs. (C) Chord plot of enriched KEGG pathways. (D) The top 10 enriched KEGG terms showed by the bubble plot. The size of the bubble represents the count of genes, and the color of the bubble represents different significance levels of KEGG terms.

Table 1

The Top 15 GO Enrichment Analysis of 23 DESRGs

Ontology	ID	Description	Gene Ratio	p.adjust
BP	GO:0097193	Intrinsic apoptotic signaling pathway	8/23	1.83e-06
BP	GO:0008630	Intrinsic apoptotic signaling pathway in response to DNA damage	5/23	1.15e-04
BP	GO:2001233	Regulation of apoptotic signaling pathway	7/23	1.79e-04
BP	GO:0042770	Signal transduction in response to DNA damage	5/23	1.79e-04
BP	GO:0010332	Response to gamma radiation	4/23	1.79e-04
CC	GO:0000792	Heterochromatin	3/22	0.005
CC	GO:0005811	Lipid droplet	3/22	0.005
CC	GO:0070603	SWI/SNF superfamily-type complex	2/22	0.081
CC	GO:0071013	Catalytic step 2 spliceosome	2/22	0.081
CC	GO:1904949	ATPase complex	2/22	0.081
MF	GO:0003714	Transcription corepressor activity	4/23	0.025
MF	GO:0047485	Protein N-terminus binding	3/23	0.025
MF	GO:0030332	Cyclin binding	2/23	0.031
MF	GO:0001046	Core promoter sequence-specific DNA binding	2/23	0.039
MF	GO:0070888	E-box binding	2/23	0.039

Abbreviations: BP, biological process; CC, cellular component; DESRGs, differentially expressed senescence-related genes; GO, gene ontology; MF, molecular function.

Identification of PPI Network and Hub Genes

The PPI network was constructed using the STRING database to identify the interactive relationship between DESRGs, and the results, including 19 nodes and 47 edges, were further visualized by Cytoscape software. Figure 5A shows the PPI network of DESRGs. The interaction number of each DESRG is presented in Figure 5B. GAPDH, SNAI1, CDKN1A, HDAC1, and ATM were the top 5 genes with the maximum interaction numbers. Next, two significant clusters were obtained by the MCODE app. Cluster 1, which had the highest score, contained 5 nodes and 8 edges, including ATM, CDKN1A, GAPDH, NOTCH3, and SNAI1 (Figure 5C). Cluster 2 included 3 nodes and 3 edges (Figure 5D). Then, Cytohubba was employed to explore hub genes. Four hub genes, including ATM, CDKN1A, GAPDH, and HDAC1, were gained by intersecting the top 6 genes obtained by Degree, MCC, MNC, Stress and EcCentricity algorithms (Table 2). Furthermore, the diagnostic value of four hub genes in COPD was examined by receiver operating characteristic (ROC) curves. The best diagnostic genes were HDAC1 (area under the curve (AUC): 1.000) and ATM (AUC: 1.000), followed by GAPDH (AUC: 0.971) and CDKN1A (AUC: 0.874) (Figure 6). These results revealed that ATM, CDKN1A, GAPDH, and HDAC1 may be hub genes of COPD, and the above-mentioned genes had good diagnostic value.

Figure 5

PPI network and key cluster of 23 DESRGs. (A) PPI network visualized by Cytoscape software. Red round rectangles represent upregulated DESRGs, and blue ellipses represent downregulated DESRGs. (B) The interaction number of DESRGs. (C) Cluster 1 (score: 4) includes 5 nodes and 8 edges. (D) Cluster 2 (score: 3) includes 3 nodes and 3 edges.

Table 2

Four Hub Genes Identified by Five Algorithms of CytoHubba

Gene Symbol	Description	LogFC	p.adjust	Changes
GAPDH	Glyceraldehyde-3-phosphate dehydrogenase	−1.062	2.79E-04	Down
HDAC1	Histone deacetylase 1	−1.036	1.08E-04	Down
CDKN1A	Cyclin-dependent kinase inhibitor 1A	1.472	4.10E-04	Up
ATM	ATM serine/threonine kinase	1.724	2.99E-09	Up

Figure 6

Diagnostic value of 4 hub genes. ROC curves of ATM (A), CDKN1A (B), GAPDH (C) and HDAC1 (D).

Prediction of Target miRNAs of Hub Genes and Construction of miRNA-mRNA Network

In this study, we used ENCORI and miRWalk databases to predict target miRNAs. After the intersection of the results predicted by the ENCORI and miRWalk databases, 31 miRNAs of ATM, 25 miRNAs of CDKN1A, 6 miRNAs of HDAC1, and 4 miRNAs of GAPDH were obtained. Then, a network between miRNAs and 4 hub genes, including 67 nodes and 66 edges, was visualized by Cytoscape software (Figure 7). The common target miRNA of ATM, CDKN1A, and GAPDH was has-miR-1224-5p.

Figure 7

Co-expressed network of hub genes and target miRNAs. Red ellipses represent hub genes, and blue diamonds represent target miRNAs.

Validation of Differentially Expressed Hub Genes in Two GSE Datasets

To further confirm differential expression levels of hub genes between the control and COPD groups, GSE103174, containing 11 healthy samples and 23 COPD samples, and GSE100281, including 16 healthy samples and 79 COPD samples, were employed. Figure 8A shows that the expression levels of CDKN1A and HDAC1 were significantly different between the two groups in the GSE100281 dataset. In the GSE103174 dataset, CDKN1A, GAPDH, and HDAC1 were differentially expressed between the control and COPD groups (Figure 8B). These results indicated that CDKN1A and HDAC1 might be the final hub genes of COPD.

Figure 8

Validation of hub genes in another two GSE datasets. (A) Expression levels of ATM, CDKN1A, GAPDH and HDAC1 in the GSE100281 dataset. (B) Differential expression of 4 hub genes between the COPD samples and healthy control samples in the GSE103174 dataset. *p<0.050, ** p<0.01, ***p<0.001, ns: no significant difference.

Diagnostic Value Analysis of CDKN1A and HDAC1 in COPD Patients

In this part, the ROC curve was used to verify the diagnostic value of hub genes. The AUC can reflect both specificity and sensitivity. The criteria of AUC were 0.5–0.7 meaning moderate indicator, 0.7–0.9 on behalf of better indicator, and > 0.9 representing high indicator. The AUC of CDKN1A and HDAC1 in the GSE100281 dataset was 0.888 and 0.672, respectively (Figure 9A). In the GSE103174 dataset, the AUC of CDKN1A and HDAC1 was 0.678 and 0.839, respectively (Figure 9B). These findings suggested that CDKN1A and HDAC1 may be diagnostic markers for COPD.

Figure 9

Verification of the diagnostic value of CDKN1A and HDAC1 in verified datasets. ROC curves of CDKN1A and HDAC1 in the GSE100281 dataset (A) and GSE103174 dataset (B).

The Expression of CDKN1A and HDAC1 in PBMCs of COPD Patients and Correlation Analysis with Lung Function

We measured the mRNA expression levels of CDKN1A and HDAC1 in PBMCs of COPD and healthy control samples. A total of 20 healthy controls and 24 COPD samples were screened. The clinical characteristics of cases are listed in Table 3. According to the correlation analysis of the GOLD stage and other clinical characteristics of the patients, we found that smoking exposure (p < 0.001) and body mass index (BMI) (p = 0.022) were related to the GOLD stage (). The mRNA expression level of CDKN1A in PBMCs of COPD samples was higher than that in the healthy controls (p < 0.050) (Figure 10A), which was consistent with the results of the GSE100281 and GSE103174 datasets. The CDKN1A level was negatively related to FEV1/FVC value (r = −0.310, p = 0.041) (Figure 10B), while it was not related to FEV1% pred (Figure 10C). Compared to the healthy controls, the expression level of HDAC1 was decreased in PBMCs of COPD samples (p < 0.001) (Figure 10D). Correlation analysis revealed that the expression of HDAC1 was positively associated with FEV1/FVC value (r = 0.520, p < 0.001) and FEV1% pred (r = 0.457, p = 0.002) (Figure 10E and F).

Table 3

The Clinical Characteristics of Patients

Variables	Non-Smokers	COPD
Number	20	24
Gender (Male/Female)	20/0	24/0
Age (Years)	62.3 ± 2.22	63.92 ± 1.81
FEV1% pred	107.4 ± 4.19	82.92 ± 3.74***
FEV1/FVC (%)	79.25 ± 1.28	60.83 ±1.90***

Note: ***p< 0.001.

Abbreviations: COPD, chronic obstructive pulmonary disease; FEV1%pred, forced expiratory volume in the first second of expiration for predicted values; FEV1/FVC, forced expiratory volume in 1 second/forced vital capacity.

Figure 10

Verification of CDKN1A and HDAC1 in PBMCs and correlation analysis between hub genes and lung function. (A) The mRNA levels of CDKN1A in the PBMCs of COPD and healthy control samples. Correlation analysis between CDKN1A and FEV1/FVC (B) and FEV1% pred (C). (D) Differential expression of HDAC1 between the PBMCs of COPD and healthy control samples. Spearman analysis between HDAC1 and FEV1/FVC (E) and FEV1% pred (F). *p<0.050, ***p<0.001.

Construction of ceRNA Between HDAC1 and Its Target Non-Coding RNAs

The expression level of HDAC1 was decreased in COPD samples. Meanwhile, the expression of HDAC1 was related to FEV1/FVC and FEV1% pred value. Increasing HDAC1 expression levels has been proved to alleviate the release of inflammatory cytokines.24 We proposed that HDAC1 was a promising therapeutic target for COPD patients. Therefore, we constructed an HDAC1-related ceRNA network. The target miRNAs of HDAC1 were has-miR-1271-5p, has-miR-3614-5p, has-miR-28-3p, has-miR-7-5p, has-miR-214-3p, and has-miR-330-3p. We predicted 7 target lncRNAs of has-miR-1271-5p, 10 target lncRNAs of has-miR-3614-5p, 11 target lncRNAs of has-miR-28-3p, 13 target lncRNAs of has-miR-7-5p, 19 target lncRNAs of has-miR-214-3p, and 14 target lncRNAs of has-miR-330-3p using the ENCORI database. LncRNA, which was the predicted target of at least two miRNAs, was selected as the final target lncRNA. Finally, 11 lncRNAs were obtained based on the above criteria (Figure 11A). The HDAC1-related ceRNA network was constructed and visualized by Cytoscape software (Figure 11B).

Figure 11

The HDAC1-related ceRNA network. (A) Upset blot of interacted lncRNAs of six HDAC1-targeted miRNAs. (B) The ceRNA network of lncRNA-miRNA-HDAC1 constructed by Cytoscape software. Red ellipses represent hub genes, and blue diamonds represent target miRNAs, and green triangles represent target lncRNAs.

Discussion

COPD is one of the leading mortal diseases worldwide.25 In the past decades, various biological processes, including inflammation, autophagy, and pyroptosis, have been considered to be involved in COPD pathogenesis.26–28 Recently, numerous studies have indicated that cellular senescence was associated with COPD. However, the value of SRGs and their mechanism of action in COPD have remained unclear. In this paper, for the first time, we identified key SRGs and explored the function of SRGs in COPD pathogenesis by bioinformatics analysis.

We screened 23 DESRGs of COPD from the GEO dataset and CellAge database. GO analysis revealed that DESRGs were enriched in apoptosis and response to damage. KEGG pathway analysis showed that DESRGs were involved in cell rest and HIF-1 signaling pathway. Previous studies have shown that DNA damage induces cellular senescence via several signaling pathways.4 HIF-1α prevents cellular senescence and is involved in the pathogenesis of aging-related chronic diseases. HIF-1 signaling pathway can interact with multiple pathways involved in cell senescence, including the SIRT1 pathway.29 Hence, we predict that SRGs may participate COPD pathogenesis via the HIF-1 signaling pathway.

After the construction of the PPI network, four hub genes (ATM, CDKN1A, GAPDH, and HDAC1) were acquired. CDKN1A and HDAC1 were identified as final hub genes based on validation by another two independent GEO datasets. CDKN1A, also known as P21 or CIP1/WAF1, a cyclin-dependent kinase (CDK) inhibitor, has been found to play a key role in controlling cell cycle progression.30 Several studies have indicated that the expression of CDKN1A is upregulated in the COPD group.31,32 In three GEO datasets and the PBMCs of our clinical cases, the expression level of CDKN1A in COPD samples was higher than that in control samples, which is consistent with the above-mentioned findings. ROC curve analysis implied that CDKN1A had a better prognostic value. HDAC1, a member of the class I HDAC family, plays an important role in controlling inflammation in COPD in vivo and in vitro.24,33 Furthermore, HADC1 is a key regulatory factor in cellular senescence.34,35 Previous studies had shown that the expression level of HADC1 in COPD is decreased compared to the control group, which is consistent with our findings.36 A study has observed that erythromycin reverses the declined expression of HDAC1 induced by cigarette smoke extract (CSE) to alleviate the release of inflammatory cytokines in human macrophages.24 These results suggested that CDKN1A may be an effective diagnostic marker, and an HADC1 inducer may be a promising therapeutic agent for COPD patients.

Protein-coding RNA and non-coding RNA, acting as endogenous miRNA sponges, have been found to interact with each other, and the crosstalk between RNAs is called ceRNA.12 Studies have demonstrated that ceRNA is an important regular in hypoxia-induced pulmonary hypertension (HPH),37 interstitial lung disease (ILD),38 and COPD.39,40 However, the role of the senescence-related ceRNA network in COPD pathogenesis has not been explored yet. In this paper, we confirmed the downregulated HDAC1 in COPD samples. Meanwhile, we found that the expression of HDAC1 was positively correlated with FEV1/FVC and FEV1% pred. Thus, HDAC1, a senescence-related gene, was selected for analysis. An HDAC1-related ceRNA network, including 11 lncRNAs (XIST, NEAT1, ZFAS1, OIP5-AS1, LINC02381, MIR17HG, FGD5-AS1, AL035425.3, AL078639.1, AC079781.5 and AC016876.2) and 6 miRNAs was constructed. Among the target lncRNAs of miRNAs, several lncRNAs have been reported. A previous study has implied that NEAT1 is upregulated in aged bone marrow mesenchymal stem cells (BMSCs) and regulates mitochondrial function by sponging miR-27b-3p.41 Additionally, NEAT1 is a marker of COPD susceptibility, and NEAT1/miR-193a affects the development of COPD.42 In CSE-exposed 16HBE cells, XIST/miR-200c-3p/EGR3 can intensify inflammation action and apoptosis.43 OIP5-AS1 is a diagnostic marker for COPD, and the OIP5-AS1/miR-410-3p/IL-13 axis plays a crucial role in regulating COPD.13 Hence, we predicted that the HDAC1-related ceRNA might contribute to clarifying the effect of HDAC1 on senescence and the role of senescence in the development of COPD. The specific role and mechanism of ceRNA in COPD need to be further verified.

Although important discoveries were revealed by this study, there are also limitations. First, due to the lack of demographic and clinical information about patients in the three datasets, the correlation analysis between this information and hub genes was not performed. Second, we only investigated the mRNA expression levels of hub genes in peripheral blood, and further verification, including protein levels of genes, and the mechanism of SRGs in COPD pathogenesis needs to be conducted both in vivo and in vitro. Additionally, we did not confirm the interplay between RNAs in the ceRNA network, and it needs to be further studied and experimentally validated.

Conclusions

In the present study, we identified 23 DESRGs between the control samples and COPD samples via bioinformatics analysis. Based on our analysis and validation, CDKN1A and HDAC1 were verified as senescence-related hub genes of COPD. CDKN1A may be a promising diagnostic marker for COPD. The HDAC1-related ceRNA may contribute to elucidating the role of cellular senescence in COPD pathogenesis. HDAC1 seemed to be a crucial therapeutic target for COPD patients.