Autophagy Induced by BCL2-Related ceRNA Network Participates in the Occurrence of COPD

Introduction

Chronic obstructive pulmonary disease (COPD) is major cause of morbidity and mortality in chronic diseases worldwide, and is one of the most significant causes of death in all disease at most countries.¹ The increasing morbidity of COPD imposes a huge economic and social burden.^2–5 Lung function is one of the diagnostic criteria of COPD that can also reflect the severity of COPD.¹

According to published studies, numerous biological functions take part in the occurrence of COPD, for instance, apoptosis, cell proliferation, and autophagy.^6–8 However, the mechanism of autophagy is not clear in COPD. Autophagy depends on a lysosomal degradation pathway, and under certain circumstances, this pathway is considered a non-apoptotic cell death pathway and is necessary for cell differentiation, homeostasis, survive, development, which also takes part in the inflammation, aging, cancer and so on.^9,10 It is related to lung,¹¹ liver,¹² and heart disease.¹³ In addition, autophagy plays a vital role in the process of infection via participating in host-pathogen interactions.^14,15 Studies have also demonstrated that autophagy participates in the occurrence and development of idiopathic pulmonary fibrosis (IPF),¹⁶ pulmonary arterial hypertension (PAH),¹⁷ and acute lung injury (ALI).¹⁸ Autophagy plays a significant role in cellular behavior in COPD.^19–21 However, there has been little research on autophagy participation in the occurrence of COPD.²² In this study, we identified the role of autophagy in COPD via bioinformatics analysis and further research.

Studies indicate that autophagy-related genes have number of biological functions²³ that take part in phagocytosis by LC3-associated phagocytosis (LAP),^24,25 secretion, and exocytosis.²⁶ BCL2 is one of the autophagy-related genes. Qin et al found that nur77 bound to BCL2 and decreased the affinity of Beclin-1 for BCL2 to promote autophagy in COPD.²⁷ Recent evidence suggests that the upregulation of BCL2 is related to the “abnormal” inflammation in COPD, and takes part in auto-maintenance.²⁸ BCL2 plays an important role in maintaining the internal environment of airway epithelial cells. Cigarette smoke (CS) blocks epithelial cells in chronic airway inflammation by inhibiting Bik expression.²² Research has found that BCL2 could sustain hyperplastic mucus cells, and CS could enhance the expression of BCL2 and help promote the survival of mucus cells.¹⁰ Noncoding RNAs, which cannot encode proteins, are classified into three groups according their nucleotide lengths: less than 50 nt, 50–500 nt, and more than 500 nt and are significant in life processes.²¹ MiRNA and lncRNA are noncoding RNAs, and miRNAs participate in the occurrence and development of most cancers, for instance malignant lung tumors,²⁹ liver carcinoma,³⁰ breast cancer,³¹ and gastric carcinoma,³² among others. In addition, studies have indicated that miRNAs have an important function in COPD³³ along with lncRNAs.³⁴ lncRNAs interact with miRNAs and regulate the expression of genes. Consequently, lncRNAs can act as competing endogenous RNAs (ceRNA) to compete with their target miRNAs and regulate gene expression. CeRNA also takes part in the development of COPD. Zhao reported that lncRNA RP11-86H7.1 promoted inflammatory reactions, thus promoting the development of COPD.³⁵

The reason there are still few effective treatments is that COPD pathogenesis remains unclear. In order to further understand this biological process, we constructed a BCL2-related ceRNA network to clarify the processes that induce autophagy in COPD. Meanwhile, it is vital as well to find potential biomarkers that can be used to guide treatment for COPD.

Methods

Data Sources

We used “COPD,” “noncoding RNA profiles by array,” and “Homo sapiens” as key words and collected three miRNA datasets of blood samples, GSE31568,³⁶ GSE24709, and GSE61741,³⁷ from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov) from the GPL9040 platform. The Human Autophagy Database (HADb) network was developed in the Laboratory of Experimental Cancer Research headed by Dr. Guy Berchem, and integrates human genes and proteins that directly or indirectly take part in autophagy, and was up-to-date as described in literature. Autophagy genes were obtained from the “clustering” module of HADb (http://www.autophagy.lu),^38–40 which contained 231 autophagy genes. Then a mRNA dataset that belongs to a lung tissue sample, GSE151052, was acquired from the GEO database GPL17556 platform via searching the key words “COPD,” “expression profiling by array,” “Homo sapiens,” and “tissue.” The lung function datasets GSE71220⁴¹ also came from the GEO database by searching “COPD,” “expression profiling by high throughput sequencing,” “Homo sapiens,” and “tissue,” which belong to the GPL11532 platform. The detail of the five datasets was shown in Table 1. The flow charts of this study were presented in Figures 1 and 2.

Table 1 The Detail of Datasets

Figure 1 The detailed process of identifying autophagy participated in the occurrence of COPD.

Figure 2 A flow diagram of BCL2-related ceRNA network.

Differential Expression Analysis

GEO2R is a GEO web analysis tool to identify differentially expressed genes via experimental conditions between various groups in a GEO series. This tool was applied to analyze the differentially expressed genes in the GSE31568, GSE24709, and GSE61741 datasets. The screening criteria were adj.p<0.05 and log (fold change [FC])>1. The NetworkAnalyst network (https://www.networkanalyst.ca)^42–46 was used to analyze 38 differentially expressed autophagy-related genes in GSE151052. NetworkAnalyst is a comprehensive gene expression profiling and visual analytics network. We used the “gene expression table” module, input the expression data, and chose the “limma” statistical method; adj.p<0.05 was the screening criterion.

Prediction of Transcription Factors and Target Genes for Differentially Expressed miRNAs (DEmiRNAs)

Transcription factors of 41 DEmiRNAs were obtained from FunRich (version 3.1.3),⁴⁷ which was used mainly for functional enrichment analysis and interaction network analysis of genes and proteins. This stand-alone software tool was designed to handle multiple gene and protein datasets irrespective of the organism. The function of “miRNA enrichment” in FunRich was used in our study to obtain the transcription factors of 41 DEmiRNAs. The top 10 transcription factors were identified in our study, and P < 0.05 was considered to be statistically significant. Then the target genes of the 41 DEmiRNAs were predicted via DIANA microT-CDS (http://diana.imis.athena-innovation.gr/DianaTools/index.php?r=microT_CDS/index),^48,49 miRDB (http://mirdb.org),^50,51 miRTarBase (http://mirtarbase.cuhk.edu.cn),⁵² and TargetScan (http://www.targetscan.org).⁵³ These tools use differential algorithms to predict target genes. The overlap of the four databases was used in further analysis.

Functional Enrichment Analysis and Protein–Protein Interaction (PPI) Network Analysis

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis can enrich many single genes and determine the common function and pathway for those genes. GO can classify the genes, proteins, ncRNAs and complexes, and then describe their biological roles.^54,55 KEGG can link the genetic information and acquire further information via computerization. Consequently, GO and KEGG methods were used to do a functional enrichment and pathway analysis. Moreover, we constructed a PPI network to help us identify the interactions within proteins, so as to furthermore research the molecular mechanism of diseases. As a consequence, a PPI network was constructed of autophagy-related differentially expressed genes in GSE151052. Metascape (https://metascape.org),⁵⁶ a gene annotation and analysis resource, which was involved in gene enrichment analysis and PPI analysis. When autophagy-related differentially expressed genes were input into the network, a document containing the result of enrichment and PPI can be exported. String (https://www.string-db.org)⁵⁷ predicted protein-protein interactions with a direct and indirect association via computational prediction and integration with other databases. Metascape and the String network were used for enrichment and protein-protein interaction (PPI) network analysis in our study.

Construction of the BCL2-Related ceRNA Network

In our study, the ceRNA network of BCL2 was constructed with miRNA and lncRNA to reveal the interaction between RNAs and to induce autophagy in COPD. MiRNAs were obtained by using DIANA microT-CDS, miRTarBase, starBase (http://starbase.sysu.edu.cn),⁵⁸ TargetScan and the intersection with DEmiRNAs. LncRNAs were found via starBase and DIANA Lncbase (http://carolina.imis.athena-innovation.gr/diana_tools/web/index.php?r=lncbasev2%2Findex),⁵⁹ respectively. Finally, miRNAs and lncRNAs were input using Cytoscape (version 3.8.2)⁶⁰ as a tool, and the BCL2-related ceRNA was structured via setting up nodes and edges in Cytoscape.

Statistical Analysis

The 41 DEmiRNAs were first collected in GEO databases and statistical analysis was performed using the GEO2R tool on the GEO database, and the statistically significant standards were adj.p<0.05 and logFC>1. Then acquired autophagy-related differentially expressed miRNAs in GEO151052 were analyzed by “limma” via the NetworkAnalyst network. The expression information of BCL2, BECN1, MAPK8, ITPR1 as well as the lung function data of COPD patients in GSE71220 were collected. Spearman analysis was applied to analyze the correlation between gene expression and lung function using SPSS 24.0.

Results

Differentially Expressed miRNAs in COPD

We first obtained GSE31568, GSE24709, and GSE61741 datasets from blood samples. Those datasets came from the GPL9040 platform. GSE31568 contained 70 control samples and 24 COPD samples, GSE24709 contained 19 control samples and 24 COPD samples, and GSE61741 contained 94 control samples and 47 COPD samples. More details are shown in Table 1. Then we obtained the differentially expressed miRNAs by using the GEO2R tool. There were 453 upregulated DEmiRNAs and 410 downregulated DEmiRNAs in GSE31568, 464 upregulated DEmiRNAs and 399 downregulated DEmiRNAs in GSE24709, 435 upregulated DEmiRNAs and 413 downregulated DEmiRNAs in GSE61741. Details are shown in Figure 3A–C); a total of 863, 863, and 848 DEmiRNAs were obtained. According to the screening criteria adj.p<0.05 and logFC>1, we obtained 113, 186, and 110 DEmiRNAs, respectively. Finally, 41 overlapping DEmiRNAs were obtained (Figure 3D), including 30 upregulated DEmiRNAs and 11 downregulated DEmiRNAs as shown in Supplementary Table 1. In addition, we also predicted the transcription factors of 41 DEmiRNAs via the FunRich tool and identified the top 10 transcription factors with the most significant differences, EGR1, SP1, POU2F1, SP4, NKX6-1, RORA, MEF2A, SOX1, ZFP161, and YY1. More details are presented in Table 2 and Figure 4.

Table 2 The Top10 Trancription Factor for DEMs

Figure 3 DEmiRNAs between COPD blood samples and control blood samples. (A) The volcano plot of DEmiRNAs in GSE31568 by GEO2R tool. (B) The volcano plot of DEmiRNAs in GSE24709 by GEO2R tool. (C) The volcano plot of DEmiRNAs in GSE61741 by GEO2R tool. (D) Venn diagram of DEniRNAs in the three datasets GSE31568, GSE24709 and GSE61741.

Figure 4 The column chart of transcription factor for DEmiRNAs, the top10 was showed.

Screening Autophagy-Related Genes According to the Predicted Target Genes of DEmiRNAs

To explore the relationship between COPD and autophagy and to clarify whether autophagy participates in the occurrence of COPD, the target genes of DEmiRNAs were predicted and autophagy genes were collected. A total of 11,447 target genes were predicted by the DIANA microT-CDS database, 7368 target genes were predicted by MiRDB database, 3255 target genes were predicted by the miRTarBase database, and 15,035 target genes were predicted via the TargetScan database, the overlap were 1596 target molecules, as shown in Figure 5A. In addition, 231 autophagy genes were downloaded from HADb as shown in Supplementary Table 2. The overlap between 231 autophagy genes and 1596 target genes was 38 autophagy-related genes. The results were shown in Figure 5B and the autophagy-related gene list was shown in Supplementary Table 3.

Figure 5 Target genes were predicted and autophagy genes were obtained. (A) The combination of target genes. 1596 target genes were obtained via the four kinds of databases. (B) The union of predicted target genes and autophagy from HADb (Human Autophagy Database).

Validation That Autophagy-Related Genes Were Differentially Expressed in Tissue Samples in COPD Patients

In order to acquire data on the differential expression of autophagy-related genes in COPD and reduce the sample differences, we used tissue samples for validation. GSE151052, a lung tissue sample dataset, was downloaded from GEO, and we acquired the expression data on the 38 autophagy-related genes in GSE151052. Then we identified 19 autophagy-related differentially expressed genes that were significantly different (P < 0.05) via the NetworkAnalyst network: EIF4EBP1, DNAJB9, RAB1A, CHMP2B, BECN1, FKBP1A, PTEN, EDEM1, ITGA6, SAR1A, WIPI2, MAPK8, ITPR1, UVRAG, CDKN1B, TP53INP2, SIRT1, VEGFA, and BCL2. The box plots were made via the NetworkAnalyst network and the heatmap was displayed via http://www.bioinformatics.com.cn, a free online platform for data analysis and visualization. More details were presented in Figures 6 and 7.

Figure 6 Identified the differentially expressed of 38 autophagy-related genes between COPD lung tissue samples and control lung tissue samples in GSE151052. (A–S) The box plot of differentially expressed genes between COPD and control, the 19 autophagy-related differentially expressed genes are remarkable difference (p<0.05).

Figure 7 The hot map of 19 autophagy-related remarkable difference genes (p<0.05).

Functional Enrichment and PPI Network Analysis of the 19 Differentially Expressed Autophagy-Related Genes

In this study, we used GO and KEGG enrichment analysis to clarify the potential biological function of these genes via the Metascape network. The results showed that the 19 genes were mainly involved in the autophagy process in COPD. The top five enrichment items were ko04140 (autophagy), GO:0009991 (response to extracellular stimulus), GO:0001666 (response to hypoxia), CORUM:6472 (Bcl-2-Beclin1-UVRAG-PI (3) KC3 complex), and GO:0010506 (regulation of autophagy). They related to autophagy almost. The first item obviously confirmed our theory that autophagy plays an important role in the occurrence of COPD. In addition, the article 5 explained genes took part in the regulation of autophagy in COPD. So, KEGG and GO enrichment analysis showed that autophagy-related differentially expressed genes play significant role in autophagy and regulation of autophagy in the occurrence of COPD. The enriched ontology cluster results were presented in Figure 8A, The enrichment factor dot plot by http://www.bioinformatics.com.cn is shown in Figure 8B The enrichment network by cluster ID and the enrichment network by P-value are displayed in Figure 8C and D.

Figure 8 Enrichment analysis of 19 autophagy-related differentially expressed genes by metascape. (A and B) The top 14 pathways of 19 autophagy-related differentially expressed genes in barplot and dotplot. (C) The enrichment network of colored by cluster ID, nodes of the same color belong to the same cluster. (D) The enrichment network of colored by P-value. The dark the color, the more statistically significant the node is.

PPI network analysis was used to determine the interaction between 19 autophagy-related differentially expressed genes by Metascape and the STRING network. We found that there were close interactions among 12 genes, for instance, BCL2, BECN1, MAPK8, ITPR1, and so on. In addition, we also found ITGA6, CHMP2B, and TP53INP2 were independent of each other at 19 autophagy-related differentially expressed genes. The network is shown in Figure 9B. In addition, we obtained the relationship of BCL2, BECN1, MAPK8 and ITPR1 genes from the PPI network analysis via Metascape. It was shown that BCL2 with BECN1, MAPK8, and ITPR1 were densely connected, as well as with BECN1, but MAPK8 and ITPR1 had no obvious connection. The network was shown in Figure 9A.

Figure 9 PPI network of 19 autophagy-related differentially expressed genes. (A) The PPI mcode components was acquired by metascape. MCODE algorithm was then applied to this network to identify neighborhoods, BCL2, BECN1, MAPK8, ITPR1, where proteins are densely connected. (B) The PPI network of 19 autophagy-related differentially expressed genes was obtained by STRING.

Associations of BCL2, BECN1, MAPK8, and ITPR1 Expression Levels with Lung Function

Lung function is not only used to diagnose COPD but also as a tool to evaluate its severity.¹ In order to explore the effect of autophagy in COPD, we analyzed the relationship between BCL2, BECN1, MAPK8, and ITPR1 gene expression and lung function in COPD. We obtained the GSE71220 dataset from the GEO database and acquired the expression data of BCL2, BECN1, MAPK8, ITPR1 genes and FEV1/FVC and FEV1 predictive data. Then by Spearman relationship analysis we found that BCL2 and FEV1/FVC, BCL2 and FEV1 predict, ITPR1 and FEV1 predict all had significant correlation in COPD via SPSS 24.0. Spearman analysis results suggested that the expression level of BCL2 was related to FEV1/FVC (P = 0.008) and FEV1 predict (P < 0.001). The expression level of ITPR1 was related to FEV1 predict (P = 0.001) in COPD, while the expression levels of BECN1 and MAPK8 were not related to FEV1/FVC (P = 0.31, P = 0.151) and FEV1 predict (P = 0.472, P = 0.331). The expression level of ITPR1 was not related to FEV1/FVC (P = 0.185). Detailed information was presented in Tables 3 and 4.

Table 3 The SPEARMAN Relationship Between FEV1/FVC and BCL2, BECN1, ITPR1, MAPK8 in COPD

Table 4 The SPEARMAN Relationship Between FEV1 Predict and BCL2, BECN1, ITPR1, MAPK8 in COPD

Construction of the BCL2-Related ceRNA Network

We found there was a strong association between expression of BCL2 and FEV1/FEV and FEV1 predict in COPD, while ITPR1 was associated only with FEV1 predict. In addition, BCL2 is one of the most important oncogenes in apoptosis research, and the relationship between BCL2 and COPD is an active research topic. In order to explain the interaction between RNAs when autophagy is induced in COPD, we chose the BCL2 gene and structured a BCL2-related ceRNA network. We obtained 284, 68, 319, and 1255 predicted miRNAs via the DIANA (microT-CDS), miRTarBase, starBase, and TargetScan databases. The overlap of the four database-predicted results were 18 miRNAs: hsa-miR-448, hsa-miR-449a, hsa-miR-96-5p, hsa-miR-20a-5p, hsa-miR-153-3p, hsa-miR-30b-5p, hsa-miR-15b-5p, hsa-miR-365a-3p, hsa-miR-125a-5p, hsa-miR-497-5p,hsa-miR-182-5p, hsa-miR-15a-5p,hsa-miR-503-5p, hsa-miR-125b-5p, hsa-miR-195-5p, hsa-miR-34a-5p, hsa-miR-21-5p, and hsa-miR-16-5p. Finally, we found hsa-miR-153-3p and hsa-miR-497-5p not only were included in those 18 miRNAs but also belong to the DEmiRNAs. The intersection graph is shown in Figure 10A. As a result, hsa-miR-153-3p and hsa-miR-497-5p were used for further prediction. A total of 341 lncRNAs were obtained via the starBase database, and 1958 lncRNAs were obtained as shown in Figure 10B and C). Eventually, the BCL2-related ceRNA network was constructed by Cytoscape and consisted of 65 lncRNAs and 2 miRNAs, as well as BCL2, which clarified the relationship between RNAs when autophagy is induced in COPD. Details are shown in Figure 11.

Figure 10 BCL2-related miRNAs and LncRNAs. (A) The venn of BCL2-related miRNAs. The 18 BCL2-related miRNAs were obtained by four databases, the intersection with DEmiRNAs is 2 miRNA, hsa-miR-153-3p, hsa-miR-497-5p. (B) The venn plot of hsa-miR-153-3p miRNA-related LncRNA by Lncbase and Starbase databases. (C) The venn plot of hsa-miR-497-5p miRNA-related LncRNA by Lncbase and Starbase databases.

Figure 11 The ceRNA network of BCL2-miRNAs-LncRNAs. (A–C) The BCL2-related ceRNA was constructed via cytoscape. The red represents the BCL2, the yellow is miRNAs and the blue is lncRNAs.

Discussion

COPD is the third leading cause of disease-related death worldwide.⁶¹ Researchers have found that COPD is an age-associated disease,⁶² and with the increase in aging populations, the morbidity of age-related diseases such as COPD is increasing. However, because there are no effective therapies, mortality and disability in COPD are increasing,^63–66 along with its financial and social burdens. As a consequence, more research on the pathogenesis of COPD to guide treatment is needed.

Increasing evidence shows that autophagy may take part in the occurrence of COPD,⁶⁷ but there are few studies of this phenomenon. An increasing number of studies have suggested that miRNAs and lncRNAs play a role and act as biomarkers to diagnose and evaluate the prognosis.^68,69 They also have found a relationship with autophagy in cancer. Fu found that miR-9, miR-31, miR-34c, miR-196, and miR-221/222 are involved in autophagy.⁷⁰ Autophagy and proliferation could be suppressed by overexpression of miR-502 in colon cancer cells,^71–74 but its value in COPD is not clear. In this study, we first studied the DEmiRNAs and determined that autophagy plays a role in in the occurrence of COPD; furthermore, we obtained the autophagy-related gene BCL2 and constructed a BCL2-related ceRNA to clarify the mechanism of autophagy in COPD.

We obtained 1596 predicted genes via four predictive databases, 231 autophagy genes via HADb. We found 38 autophagy-related genes in blood samples of COPD. For further confirmation of the role of autophagy in COPD, we further undertook a differential expression analysis for 38 autophagy-related genes in lung tissue samples and found a total of 19 autophagy-related genes that were differentially expressed. Blood samples and tissue displayed a consistently differential expression trend, indicating that COPD is a systemic disease with not only pulmonary features but also extrapulmonary features.^75,76 For instance, cigarette smoking is the predominant etiology of COPD.¹ Cigarette smoke is rich in many compounds, such as high concentrations of acrolein and other carcinogenic organic compounds. Studies have reported an inextricable link between acrolein and oxidant stress,⁷⁷ and multiple lines of evidence indicate plasma acrolein and lung tissue acrolein are increased synchronously, which further supports COPD as a systematic disease.^78–80 Yasuo found that the longer the period of smoking, the higher acrolein in the plasma, inducing oxidant stress and small airway obstruction and resulting in the progression of COPD.⁸⁰ Our study could identify COPD is a systemic disease from the gene perspective.

Then we found the major enrichment item was autophagy, showing that autophagy takes part in the occurrence of COPD. Several studies have suggested that autophagy in part of the pathogenesis in the occurrence of COPD. CS is the leading risk factor in COPD. Some studies have shown that autophagy and selective autophagy create a harmful reaction in epithelial cells in COPD,^81–85 promoting airway inflammation and mucus hyperproduction,^86,87 clarifying the role of autophagy in occurrence of COPD. We also constructed a PPI network in order to understand the mechanism of autophagy-related differentially expressed genes and found a strong interaction among BCL2, BECN1, MAPK8, and ITPR1. In this study, we researched the relationship between BCL2, BECN1, MAPK8, and ITPR1 gene expression and lung function in COPD to find a biomarker to help further research on the influence of autophagy in COPD. We found BCL2 and FEV1/FVC, BCL2 and FEV1 predict, ITPR1 and FEV1 predict in COPD and had a significant correlation via Spearman relationship analysis. BCL2 encodes an integral outer mitochondrial membrane protein that blocks the apoptotic death of some cells such as lymphocytes, and as a consequence, BCL2 is associated with autophagy.^88,89 Studies have shown that BCL2 is correlated with lung function and may take part in the occurrence and progression of COPD.⁹⁰ BCL2 is a mitochondrial outer membrane permeabilization protein that mediates cell apoptosis, which induces the occurrence and progression of COPD.^91,92 As a result, autophagy is associated with the occurrence of COPD through the differential expression of BCL2 in mediating cell apoptosis.

In recent years, the ceRNA network has been widely studied and it has been shown that ceRNA affects the development of numerous cancers, for example, lung cancer,^93,94 breast cancer,⁹⁵ and gastric cancer.³² It was also studied in COPD,³⁵ but it was deficient. Last, we choose the autophagy-related gene BCL2 for further research. A BCL2-related ceRNA network was constructed to clarify the mechanism of autophagy and find biomarkers to guide treatment in COPD. The ceRNA network may also clarify the mechanism of RNAs in inducing autophagy and participating in the occurrence and development of COPD.

There are several limitations in our research. First, we used bioinformatics analysis throughout the whole research process and did not verify experimental studies. Second, we did not demonstrate the process of RNAs interaction in ceRNA network, and it needs to be further researched and experimentally validated. Third, we acquired autophagy-related genes in blood samples and in order to decrease the sample differences, we choose tissue samples to verify differential gene expression, which could increase errors.

Conclusion

In conclusion, 19 autophagy-related differentially expressed genes were found and were identified as participating in the occurrence of COPD. There was a correlation between BCL2 and lung function, and a BCL2-related ceRNA network was constructed to explain the interaction between RNAs that induce autophagy in COPD.