Identification and Validation of Autophagy-Related Genes in Chronic Obstructive Pulmonary Disease

Introduction

Chronic obstructive pulmonary disease (COPD) is a chronic respiratory disease characterized by incompletely reversible airway obstruction, with a high mortality and disability rate.¹ Compared with healthy individuals, patients with COPD have an increased risk of other diseases, such as lung cancer.² Previous studies reported that the risk factors of COPD include genetic factors, smoking and airway inflammation.^3–5 Accumulating evidence has shown that several biological functions are involved in the pathogenesis of COPD, including cell proliferation, apoptosis and autophagy.^6–8 Among these biological functions, autophagy plays essential roles in the development of COPD.

Autophagy is the conserved mechanism that delivers endogenous or exogenous cytoplasmic materials to the lysosomes for degradation.⁹ Autophagy is related to various diseases including respiratory diseases. For instance, miR-93 regulates tumorigenicity and therapy response of glioblastoma by targeting autophagy.¹⁰ In addition, IL-4 induces autophagy in B cells leading to exacerbated asthma.¹¹ Some signaling pathways have been reported to affect the biological function of COPD through autophagy. PI3K/AKT/mTOR pathway regulates autophagy to induce apoptosis of alveolar epithelial cells in COPD.¹² However, autophagy-related genes of COPD remain largely unknown and need to be further explored. Exploring and revealing the potential autophagy-related genes of COPD will provide us potential biomarkers to treat COPD.

Ezzie et al completed one COPD-related dataset GSE38974, which analyzed the differentially expressed genes between COPD patients and healthy individuals.¹³ Their results showed that 70 miRNAs and 2667 mRNAs were differentially expressed in the two groups. In this study, the dataset was analyzed again from other perspectives. We explored the differentially expressed autophagy-related genes of COPD by analyzing the dataset GSE38974 from GEO database. Then, protein–protein interactions (PPI), correlation analysis, gene-ontology (GO) enrichment analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis were applied for the differentially expressed autophagy-related genes. Finally, the expression levels of key differentially expressed genes were further verified in COPD patients and healthy individuals.

Materials and Methods

Autophagy-Related Genes Datasets and Microarray Data

A total of 222 genes were obtained from The Human Autophagy Database (http://www.autophagy.lu/index.html). The mRNA expression profile dataset of GSE38974 was downloaded from GEO (http://www.ncbi.nlm.nih.gov/geo/). GSE38974 is in GPL4133 platform (Agilent-014850 Whole Human Genome Microarray 4x44K G4112F), which included 23 COPD and 9 normal lung tissue samples.

Differentially Expressed Analysis of Autophagy-Related Genes

The normalized expression matrix of microarray data was downloaded from the GSE38974 dataset. Then the probes was annotated with the annotation files from the dataset. The repeatability of data in GSE38974 was verified by principal component analysis (PCA). The “limma” package of R software was used to identify the differentially expressed autophagy-related genes. Genes with an adjusted P-value <0.05 and absolute fold-change value >1.5 were considered as differentially expressed genes. The heatmap, volcano plot and box plot were conducted using “heatmap” and “ggplot2” packages of R software.

PPI Analysis and Correlation Analysis of the Differentially Expressed Autophagy-Related Genes

PPI analysis of differentially expressed autophagy-related genes was analyzed using STRING database (https://string-db.org/) and Cytoscape software (version 3.8.1). The correlation analysis of the differentially expressed autophagy-related genes was identified using Spearman correlation in the “corrplot” package of R software.

GO and KEGG Pathway Enrichment Analysis of Autophagy-Related Genes

GO and KEGG pathway enrichment analysis were conducted in R software using the package “GO plot”. The GO analysis consisted of cellular component (CC), biological process (BP) and molecular function (MF).

COPD Patients and Healthy Individuals

A total of 20 COPD patients (cases) and 20 age-matched healthy individuals (controls) were obtained from the Tianjin Medical University General Hospital between July 2020 and October 2020. The diagnosis of COPD was followed according to the Global Initiative for Chronic Obstructive Pulmonary Disease (GOLD) criteria: patients with post-bronchodilator FEV₁/FVC <70%, and <80% predicted FEV_1. Patients were excluded if they had bronchial asthma, bronchiectasis, pulmonary fibrosis, lung tumor and tuberculosis. The 20 healthy individuals were recruited from the hospital’s health checkup center. This study was conducted in accordance with the Declaration of Helsinki and approved by the Medical Ethics Committee of the hospital. Written informed consent was obtained from all the participants. Venous blood was collected from all cases and controls who participated in the study.

RNA Extraction and Quantitative Real-Time Polymerase Chain Reaction (qRT-PCR)

The blood samples from each participant were processed to isolate peripheral blood mononuclear cells (PBMCs) using Ficoll solution (Solarbio Life Sciences, Beijing, China). Total RNA was extracted from PBMCs with RNA Extraction Kit (Omega, Guangzhou, China). Reverse transcription was conducted using PrimeScript RT Master Mix Kit (Takara, Dalian, China). The mRNA level was assessed using TB Green Premix Ex Taq Kit (Takara, Dalian, China) following the instructions. Primers are available in Supplementary Table S1. The relative expression of mRNA was calculated by 2^−ΔΔCt method with the normalization to ACTB.

Statistical Analysis

The statistical analyses were performed using R software (version 3.6.2). Gene expression levels of our clinical samples were compared using Student’s t-test. P < 0.05 was considered statistically significant.

Results

Differentially Expressed Autophagy-Related Genes in COPD-Retrospective Analysis of Autophagy-Related Genes

To evaluate the intra-group data repeatability, we performed principal component analysis (PCA), and the results showed that the repeatability of data in GSE38974 is fine (Figure 1A). We next analyzed the expression of 222 autophagy-related genes in 23 COPD patients and 9 healthy individuals, and 40 autophagy-related genes were identified using the criteria of adjusted P-value <0.05 and absolute fold-change value >1.5, including 14 up-regulated genes and 26 down-regulated genes (Table 1). Following the analysis of the GSE38974 dataset with R software, the 40 differentially expressed autophagy-related genes between COPD and normal groups were presented in heatmap and volcano plot (Figure 1B and C). Moreover, box plots showed the expression patterns of 40 differentially expressed autophagy-related genes between COPD and normal samples (Figure 2A and B). The top five up-regulated genes included NAF1, HIF1A, CDKN1A, BAG3 and CHMP4B, and the top five down-regulated genes included ERBB2, GAPDH, FAM48A, ATG16L1, and ITPR1. (Figure 2A and B; Table 1).

Table 1 The 40 Differentially Expressed Autophagy-Related Genes in COPD Samples Compared to Healthy Samples

Figure 1 Differentially expressed autophagy-related genes in COPD and healthy samples. (A) Principal component analysis for GSE38974. (B) Volcano plot of the 222 differentially expressed autophagy-related genes. The red dots represent the significantly up-regulated genes and the blue dots indicate the significantly down-regulated genes. (C) Heatmap of the 40 differentially expressed autophagy-related genes in COPD and healthy samples.

Figure 2 The boxplot of 40 differentially expressed autophagy-related genes in COPD and healthy samples. (A) The boxplot of top 20 differentially expressed autophagy-related genes in COPD and healthy samples. (B) The boxplot of last 20 differentially expressed autophagy-related genes in COPD and healthy samples.

PPI Network and Correlation Analysis of the Differentially Expressed Autophagy-Related Genes

To determine the interactions among differentially expressed autophagy-related genes, we performed PPI analysis. The results demonstrated that these autophagy-related genes interacted with each other (Figure 3A) and showed the interaction number of each gene (Figure 3B). To explore the expression correlation of these autophagy-related genes, correlation analysis was performed. The results showed the relationship of the 40 differentially expressed autophagy-related genes in GSE38974 dataset (Figure 4).

Figure 3 Protein–protein interactions (PPI) analysis the 40 differentially expressed autophagy-related genes. (A) The PPI among 40 differentially expressed autophagy-related genes. (B) The interaction number of each differentially expressed autophagy-related gene.

Figure 4 Spearman correlation analysis of the 40 differentially expressed autophagy-related genes.

GO and KEGG Enrichment Analysis of the Differentially Expressed Autophagy-Related Genes

To analyze the potential biological functions of these differentially expressed autophagy-related genes, we conducted GO and KEGG enrichment analysis by using R software. The results revealed that the most significant GO enriched terms involved in autophagy, process utilizing autophagic mechanism, macroautophagy (biological process); autophagosome, membrane raft, membrane microdomain (cellular component); ubiquitin-like protein ligase binding, ubiquitin protein ligase binding, protein serine/threonine kinase activity (molecular function) (Figure 5A and B; Supplementary Table S2). In KEGG enrichment analysis, the differentially expressed autophagy-related genes mainly involved in the process of autophagy and mitophagy (Figure 6; Supplementary Table S3).

Figure 5 Gene Ontology (GO) enrichment analysis of 40 differentially expressed autophagy-related genes. (A) and (B) Bubble plot of enriched GO terms. Abbreviations: BP, biological process; CC, cellular component; MF, molecular function.

Figure 6 Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis of 40 differentially expressed autophagy-related genes.

Validation the Differentially Expressed Autophagy-Related Genes in COPD Patients – Prospective Analysis in New Population

To validate the reliability of the GSE38974 dataset, the expression levels of top five differentially expressed autophagy-related genes were further identified by qRT-PCR in our clinical samples. The clinicopathological variables of cases and controls are summarized in Table 2. Similar to the results of mRNA microarray in lung tissue samples, the expression levels of HIF1A, CDKN1A and BAG3 were significantly higher in COPD blood samples than in normal blood samples (Figure 7A). In addition, the expression levels of ERBB2 and ATG16L1 were significantly decreased (Figure 7B). However, the expression levels of NAF1, CHMP4B, GAPDH, FAM48A and ITPR1 showed no significant difference between the two groups (Figure 7A and B).

Table 2 Clinicopathological Variables of Cases and Controls in This Study

Figure 7 RNA expression of top five autophagy-related genes were measured in COPD and healthy samples. (A) RNA expression of NAF1, HIF1A, CDKN1A, BAG3 and CHMP4B were measured in blood samples using qRT-PCR. (B) RNA expression of ERBB2, GAPDH, FAM48A, ATG16L1 and ITPR1 were measured in blood samples using qRT-PCR. P-values were calculated using a two-sided unpaired Student’s t-test. **P<0.01; ***P<0.001; ****P<0.0001. Abbreviation: ns, non significant.

Discussion

The development of COPD is influenced by smoking, inflammation, airway remodeling and other factors. Increasing evidences suggested that autophagy may be involved in the pathogenesis of COPD. Smoking can lead to oxidative stress of lung cells, which further leads to autophagy activation of lung epithelial cells. Autophagy causes emphysema by inducing autophagic death of lung cells.¹⁴ Another evidence suggested that autophagy can promote the release of inflammatory cytokines, which contributes to the inflammatory response of COPD.¹⁵ However, extensive validations are needed to improve the understanding of autophagy in pathogenesis of COPD.

A series of recent studies have explored the link between PBMCs and lung tissue RNAs and COPD. Dang et al performed RNA profiling of PBMCs from smokers and COPD patients by microarray. In their study, some differentially expressed miRNAs and mRNAs including miRNA-320b, miR-24-3p, CD177 and IL6 were identified.¹⁶ In addition, a recent study found that circRNA0001859 was down-regulated in lung tissue of mice and may be a potential biomarker for the treatment of COPD.¹⁷ Moreover, one previous study demonstrated that autophagy‑associated protein levels of PBMCs in COPD patients were increased and were correlated with FEV₁% predicted values and circulating levels of cytokines.¹⁵ However, there are still few related studies in this field, and further investigations will be required to better understand this field.

To the best of our knowledge, there are several published cancer-related articles exploring autophagy-related genes.^18–20 For instance, a recent study reported 30 autophagy-related genes in lung adenocarcinoma, some of which are associated with the prognosis of patients.¹⁹ However, bioinformatics analysis of autophagy-related genes has not been explored in COPD. In this study, we were the first identified 40 potential autophagy-related genes in COPD through bioinformatics analysis. Some of these autophagy-related genes of COPD has been previously studied. For example, Ding et al demonstrated that the serum VEGFA is one promising diagnostic biomarker of asthma-COPD overlap syndrome.²¹ In addition, evidence suggested that SIRT1 and FOXO1 mRNA expression levels in PBMCs correlate to physical activity in COPD patients.²² We intend to explore more potential autophagy-related genes of COPD in the future.

The potential biological functions of these differentially expressed autophagy-related genes were also conducted through GO and KEGG enrichment analysis. GO and KEGG enrichment analysis of differentially expressed autophagy-related genes indicated several enriched terms related to autophagy and mitophagy. Several published articles have confirmed that autophagy can affect the progress of COPD. One in vitro study revealed that miR-21 could increase autophagy and promote the apoptosis of 16HBE cells in COPD.²³ Another recent study showed that FUNDC1 silencing could suppress the progress of COPD by inhibiting mitochondrial autophagy through interaction with DRP1.²⁴ Future experiments will be necessary to explore the potential biological functions of these differentially expressed autophagy-related genes.

Based on bioinformatics analysis results, the expression levels of top five differentially expressed autophagy-related genes were further identified by qRT-PCR in our clinical samples. The results of qRT-PCR showed that the expression levels of HIF1A, CDKN1A, BAG3, ERBB2 and ATG16L1 were consistent with the bioinformatics analysis results from mRNA microarray. Some of these verified genes have been reported to be associated with respiratory diseases. Wang et al revealed that the HIF1A gene rs10873142 polymorphism increases the risk of COPD in a Chinese Han population.²⁵ Another study also found that the expression level of CDKN1A is significantly higher in COPD samples when compared with normal samples, in agreement with the results from our clinical samples.²⁶ BAG3 has not been reported in COPD, but BAG3 could promote resistance to apoptosis through Bcl-2 family members in non-small cell lung cancer (NSCLC).²⁷ One evidence demonstrated that the rs4719839 SNP in ATG16L1 can influence the expression of miR-148 and negatively regulate ATG16L1 expression to increase the risk of ventilator-associated pneumonia.²⁸ However, the explicit mechanism of these verified genes in COPD remains largely unclear and needs to be further explored.

Several studies have reported that autophagy is involved in different respiratory diseases. For instance, one research demonstrated that YBX1 (Y-box binding protein 1) could mediate autophagy by targeting p110β and decreasing the sensitivity to cisplatin in NSCLC.²⁹ Moreover, miR-30a targets ATG5 (autophagy-related 5) and attenuates airway fibrosis in asthma by suppressing autophagy.³⁰ Recent study reported a novel role of CX3CR1 (C-X3-C motif chemokine receptor 1) in regulation of macrophage autophagy and promotion of pulmonary fibrosis in hyperoxic lung injured mice.³¹ In summary, autophagy plays an important role in respiratory diseases.

There are still some limitations in our research. First, the bioinformatics results were obtained from the lung tissues of smokers and COPD patients. However, we performed experimental verification in blood samples of the included individuals. Second, the number of clinical samples included in our study is limited, and we need to confirm our conclusions in a larger COPD cohort. Third, we only verified the expression level of the differentially expressed autophagy-related genes in clinical samples, but did not explore the potential mechanism of these genes in COPD cells and mouse models. Therefore, further research needs to be explored in the future.

Conclusion

In conclusion, we identified 40 potential autophagy-related genes of COPD through bioinformatics analysis. Moreover, the key genes HIF1A, CDKN1A, BAG3, ERBB2 and ATG16L1 may affect the development of COPD by regulating autophagy. These results expanded our understanding of COPD and might be useful in treatment of COPD.