Construction of Potential miRNA-mRNA Regulatory Network in COPD Plasma by Bioinformatics Analysis.

Background: Chronic obstructive pulmonary disease (COPD) has become a major cause of morbidity and mortality worldwide. Increasing evidence indicates that aberrantly expressed microRNAs (miRNAs) are involved in the pathogenesis of COPD. However, an integrative exploration of miRNA-mRNA regulatory network in COPD plasma remains lacking.
Methods: The microarray datasets GSE24709, GSE61741, and GSE31568 were downloaded from the GEO database and analyzed using GEO2R tool to identify differentially expressed miRNAs (DEMs) between COPD and normal plasma. The consistently changing miRNAs in the three datasets were screened out as candidate DEMs. Potential upstream transcription factors and downstream target genes of candidate DEMs were predicted by FunRich and miRNet, respectively. Next, GO annotation and KEGG pathway enrichment analysis for target genes were performed using DAVID. Then, PPI and DEM-hub gene network were constructed using the STRING database and Cytoscape software. Finally, GSE56768 was used to evaluate the hub gene expressions.
Results: A total of nine (six upregulated and three downregulated) DEMs were screened out in the above three datasets. SP1 was predicted to potentially regulate most of the downregulated DEMs, while YY1 and E2F1 could regulate both upregulated and downregulated DEMs. 1139 target genes were then predicted, including 596 upregulated DEM target genes and 543 downregulated DEM target genes. Target genes of DEMs were mainly enriched in PI3K/Akt signaling pathway, mTOR signaling pathway, and autophagy. Through the DEM-hub gene network construction, most of the hub genes were found to be potentially modulated by miR-497-5p, miR-130b-5p, and miR-126-5p. Among the top 12 hub genes, MYC and FOXO1 expressions were consistent with that in the GSE56768 dataset.
Conclusion: In the study, potential miRNA-mRNA regulatory network was firstly constructed in COPD plasma, which may provide a new insight into the pathogenesis and treatment of COPD.

Introduction

Chronic obstructive pulmonary disease (COPD), characterized by persistent respiratory symptoms, progressive airflow limitation, and lack of effective treatment, is one of the leading cause of morbidity and mortality all over the world.1,2 In China, two national surveys of COPD showed that the overall prevalence of COPD in people aged 40 years or older has increased from 8.2% during 2002–20043 to 13.7% during 2014–2015,4 which has almost doubled, putting increasing economic pressure on the patients and society. As we all know, exposure to cigarette smoke (CS) is the largest risk factor for COPD. However, most studies have shown that no more than 80% of patients with COPD are caused by CS,5,6 that is to say, there are many other factors that contribute to COPD susceptibility.

MicroRNAs (miRNAs) are a class of endogenous small noncoding RNA molecules, consisting of 18–25 nucleotides. They can regulate gene expression by directly degrading target messenger RNAs (mRNAs) or blocking protein translation,7,8 and thus participate in many biological processes such as cell proliferation, differentiation, apoptosis, disease occurrence and progression.9,10 Studies have shown that dysregulated miRNAs are related to the onset and development of multiple diseases, including COPD.11–13 Since miRNAs can also be stably expressed in body fluids, miRNAs (especially circulating miRNAs) may be an integrated part of the COPD screening programs during future preventive treatment. Hu et al14 identified that plasma miR-125b could distinguish acute exacerbation of COPD from stable COPD and healthy controls, and correlated with inflammatory cytokine expressions, including tumor necrosis factor (TNF)-α, interleukin (IL)-6, IL-8, and leukotriene B4 (LTB-4). Tang et al15 showed that compared with healthy non-smokers, the plasma levels of miR-29b was reduced in smokers and COPD patients, with a positive correlation with decreased pulmonary function. Moreover, a recent study performed by Sundar et al16 for the first time reported that plasma-derived exosome miRNAs functioned as novel COPD biomarkers, which represents a greater translational potential for miRNAs in the diagnosis and treatment of COPD.

Despite many studies on the expression and function of miRNAs in COPD have been carried out, however, little study on the role of miRNA–mRNA regulatory network in COPD plasma has been conducted to date. In this study, we for the first time screened DEMs in COPD plasma compared with normal plasma by analyzing three datasets (GSE24709, GSE61741, and GSE31568) from the GEO database. Transcription factor-DEM, DEM-target gene, functional analysis of target genes, and DEM-hub gene network were then performed. In addition, the plasma levels of hub genes were further verified in the GSE56768 dataset. Herein, we aimed to construct potential miRNA–mRNA regulatory network in COPD plasma, and to provide promising targets for COPD diagnosis and treatment.

Methods

miRNA Microarray

According to the Gene Expression Omnibus (GEO) database (), the following words, that “COPD” (study keywords), “noncoding RNA profiles by array” (study type), “Home sapiens” (organism), and “plasma” (attribute name) were used to screen the miRNA datasets related to COPD plasma. Finally, three miRNA datasets (GSE24709, GSE6174117 and GSE31568) were chosen for subsequent analysis, which were all based on the GPL9040 platform (febit Homo Sapiens miRBase 13.0). Detailed information of the three datasets is shown in Table 1, and a flow chart of the study design is shown in Figure 1.

Table 1

Details for GEO COPD Data

Accession	Platform	Sample	Normal	COPD	Gene/microRNA
GSE24709	GPL9040	Blood	19	24	microRNA
GSE61741	GPL9040	Blood	94	47	microRNA
GSE31568	GPL9040	Blood	70	24	microRNA
GSE56768	GPL570	Blood	29	49	gene

Abbreviations: COPD, chronic obstructive pulmonary disease; GEO, Gene Expression Omnibus.

Figure 1

Flow chart of constructing the miRNA–mRNA regulatory network in COPD plasma.

Screening of DEMs

GEO2R (), a R-based web application that helps to analyze GEO data, was applied to compare and screen DEMs between COPD and normal plasma. The adjusted P value using the Benjamini and Hochberg false discovery rate (FDR) method was applied to correct for the occurrence of false positive results. Here, an adjusted P value <0.05 and |Log fold change (FC)|>1.5 were set as the threshold for identifying DEMs. Venn diagram was utilized to analyze the overlapping of DEMs in the three datasets (GSE24709, GSE61741, and GSE31568), and these consistently changed miRNAs were selected as candidate DEMs. Heatmaps of DEMs and candidate DEMs were generated using the TBtools, which is a toolkit integrating various biological data handling tools.

Prediction of Upstream Transcription Factors of DEMs

FunRich (), a stand-alone software tool mainly for analyzing functional enrichment and interaction network analysis of genes and proteins, was used to predict the potential upstream transcription factors of candidate DEMs. The P value < 0.05 was considered statistically significant.

Prediction of Downstream Target Genes of DEMs

miRNet (), an integrated tool kit that contains information about miRNA-target interactions and displays the association in a visual network, was used to predict the potential downstream target genes of candidate DEMs.

GO and KEGG Analysis of Target Genes

To investigate the biological function of DEM target genes, the database for annotation, visualization, and integrated discovery (DAVID) online tool was utilized to perform Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis. The adjusted P value < 0.05 was statistically significant.

Construction of PPI Network and Screening of Hub Genes

The protein-protein interaction (PPI) network of DEM target genes was performed by the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) online database () and then visualized in Cytoscape software (version 3.6.1). CytoHubba, a plugin of Cytoscape, was used to rank nodes in a preloaded PPI network by its features through several topological algorithms. Here, the hub genes were selected out as the top 30 nodes of the PPI network using the Maximal Clique Centrality (MCC) method, which has a better performance on the precision of predicting essential proteins. The score for a node v by MCC is defined as MCC(v) = ∑C∈S(v)(|C| - 1)!, where S(v) is the collection of maximal cliques which contain v, and (|C|-1)! is the product of all positive integers less than |C|.

Expression Analysis of Hub Genes by GSE56768 Dataset

GSE56768, downloaded from the GEO database, was introduced to analyze the expressions of hub genes due to no other data on mRNA expressions in COPD plasma. The dataset was based on the platform of GPL570 (Affymetrix Human Genome U133 Plus 2.0 Array), in which 49 COPD and 29 normal plasma were selected. Student’s t-test was used to identify the differentially expressed genes (DEGs) between COPD and normal plasma. The significantly hub genes needed to meet two criteria as follows: one was that upregulated DEM target genes were downregulated or downregulated DEM target genes were upregulated; the second was the P value should be less than 0.05.

Results

Identification of DEMs in COPD Plasma

After screening with the threshold of an adjusted Pvalue <0.05 and |logFC|>1.5, 148 DEMs (78 upregulated and 70 downregulated) were identified in the GSE24709 dataset, 26 DEMs (19 upregulated and seven downregulated) were identified in GSE61741, and 71 DEMs (18 upregulated and 53 downregulated) were identified in GSE31568. The heatmaps of DEMs in the above three datasets are shown in Figure 2A–C, and the detailed information about DEMs is listed in –. Among these DEMs, six upregulated (miR-1246, miR-1258, miR-556-3p, miR-1468-5p, miR-126-5p, and miR-130b-5p) and three downregulated (miR-182-3p, miR-497-5p, and miR-492) DEMs were consistently changed in the three datasets (Figure 2D and E).

Figure 2

Identification of the candidate DEMs.

In the study, FunRich software was used to predict the upstream transcription factors of candidate DEMs. The transcription factors for upregulated and downregulated DEMs are presented in Figure 3A and B, respectively. For upregulated DEMs, the transcription factors were YY1 and E2F1. For downregulated DEMs, the transcription factors were SP1, YY1, EGR1, E2F1, ZFP161, ELF1, SP4, and MYC.

Figure 3

Potential transcription factors of DEMs predicted by FunRich.

A total of 1139 genes targeted by candidate DEMs were predicted by the miRNet database, among which there were 596 target genes for upregulated DEMs and 543 target genes for downregulated DEMs. For better visualization, upregulated and downregulated DEMs with their target genes were displayed in a DEM-target gene network in Figure 4A and B, respectively. In addition, the number of target genes for each DEM is listed in Table 2, and all predicted target genes are listed in .

Figure 4

Potential target genes of DEMs predicted by miRNet.

Table 2

Potential Target Genes of the Significantly Upregulated and Downregulated DEMs

Upregulated DEMs	Number	Downregulated DEMs	Number
hsa-mir-126-5p	119	hsa-mir-182-3p	40
hsa-mir-130b-5p	344	hsa-mir-492	45
hsa-mir-556-3p	54	hsa-mir-497-5p	461
hsa-mir-1246	44	Total	543
hsa-mir-1258	36
hsa-mir-1468-5p	17
Total	596

Abbreviations: DEMs, differentially expressed microRNAs; miRNA, microRNA.

Then, we performed GO and KEGG analysis on the 1139 target genes of DEMs. As shown in Figure 5A and B, biological process (BP) analysis showed that DEM target genes were particularly enriched in muscle tissue development, striated muscle tissue development, and activation of protein kinase activity; cellular component (CC) analysis revealed that DEM target genes were mainly enriched in nuclear speck, transferase complex, transferring phosphorus-containing groups, and protein kinase complex; molecular function (MF) analysis demonstrated that DEM target genes were significantly enriched in DNA-binding transcription activator activity, RNA polymerase II specific, cell adhesion molecule binding, and protein serine/threonine kinase activity.

Figure 5

GO annotation analysis for the target genes of DEMs in the biological process, cellular component, and molecular function.

Furtherly, KEGG pathway analysis showed that DEM target genes were significantly enriched in PI3K-Akt signaling pathway, human papillomavirus infection, focal adhesion, breast cancer, human T-cell leukemia virus 1 infection and cellular senescence (Figure 6A and B).

Figure 6

KEGG pathway analysis for the target genes of DEMs.

Construction of PPI and DEM-Hub Gene Network

The PPI network of genes targeted by upregulated and downregulated DEMs was established by the STRING database, respectively. Then we uploaded the above PPI network to Cytoscape and used its cytohubba plugin to screen for hub genes. The top 30 hub genes for upregulated and downregulated DEMs are shown in Figure 7A and B, respectively. For upregulated DEMs, the top 10 hub genes were MYC, VEGFA, MAPK8, STAT5A, FOXO3, FOXO1, IGF1, NR3C1, HSPA8, and EP300; For downregulated DEMs, the top 10 hub genes were FBXW7, CUL3, BTRC, ZBTB16, ANAPC13, CDC34, FBXL18, FBXL20, SMURF1, and UBE2V1 (Table 3).

Figure 7

Identification of the hub genes for DEMs in the PPI network.

Table 3

Top 10 Hub Genes of the Significantly Upregulated and Downregulated DEMs in the PPI Network Ranked by MCC

Upregulated DEMs		Downregulated DEMs
Gene Symbol	Score	Gene Symbol	Score
MYC	22,916	FBXW7	87,178,818,926
VEGFA	22,125	CUL3	87,178,299,448
MAPK8	18,605	BTRC	87,178,293,375
STAT5A	15,337	ZBTB16	87,178,292,663
FOXO3	14,769	ANAPC13	87,178,292,097
FOXO1	14,200	CDC34	87,178,292,085
IGF1	13,828	FBXL18	87,178,291,926
NR3C1	12,845	FBXL20	87,178,291,922
HSPA8	11,147	SMURF1	87,178,291,408
EP300	8381	UBE2V1	87,178,291,229

Abbreviations: DEMs, differentially expressed microRNAs; MCC, maximal clique centrality; miRNA, microRNA; PPI, protein-protein interaction.

To better investigate the molecular mechanisms of these DEMs in COPD plasma, the DEM-hub gene network was constructed by Cytoscape software. The interactions between upregulated DEMs and hub genes were as follows: miR-130b-5p interacted with five hub genes, including EP300, HSPA8, STAT5A, IGF1, and FOXO1; miR-126-5p interacted with five hub genes, including MYC, MAPK8, VEGFA, NR3C1, and FOXO3; The interactions between downregulated DEMs and hub genes were as follows: miR-497-5p interacted with eight hub genes, including SMURF1, BTRC, ZBTB16, FBXL20, CUL3, UBE2V1, FBXL18 and ANAPC13; miR-182-3p interacted with FBXW7; and miR-492 interacted with CDC34 (Figure 8).

Figure 8

The miRNA-hub gene regulatory network.

Validation of Hub Genes Expression

According to the DEM-hub gene network, plasma levels of the top 12 hub genes (MYC, VEGFA, MAPK8, STAT5A, FOXO3, FOXO1, FBXW7, CUL3, BTRC, ZBTB16, ANAPC13, and CDC34), were then identified using the GSE56768 dataset. For upregulated DEMs, compared with normal plasma, only the expressions of MYC and FOXO1 were consistently decreased in COPD plasma, while VEGFA expression was increased on the contrary (Figure 9A–F). For downregulated DEMs, no miRNAs were significantly changed (Figure 9G–L). Therefore, miR-126-5p-MYC and miR-130b-5p-FOXO1 were identified as two potential regulatory pathways in COPD plasma.

Figure 9

The mRNA expression of the top 12 hub genes was determined from the GSE56768 dataset.

Discussion

During the past few years, increasing number of studies have been conducted in COPD diagnosis and treatment, however, the prognosis of these patients is still poor partly due to the limited understanding of COPD pathogenesis. Recently, with microarray technology, it is easier to reveal thousands of genetic alterations in the progression of various diseases. In the present study, data from three datasets, namely GSE24709, GSE61741, and GSE31568, were extracted to identify DEMs between COPD and normal plasma. A total of six upregulated DEMs (miR-1246, miR-1258, miR-556-3p, miR-1468-5p, miR-126-5p, and miR-130b-5p) and three downregulated DEMs (miR-182-3p, miR-497-5p, and miR-492), which were consistently changed in the three datasets, were identified as candidate DEMs to be further analyzed. Among upregulated DEMs, detection of plasma and plasma exosome miR-1246 and miR-130b-5p can be used as diagnosis biomarkers for various cancers. Despite that, anti-miR-1246 alleviated lipopolysaccharide (LPS)-induced inflammation and permeability in vivo and vitro.18 Anti-miR-130b-5p prevented hepatic lipid accumulation and insulin resistance in nonalcoholic fatty liver disease.19 For downregulated DEMs, overexpression of miR-497-5p inhibited IL-1β-induced cartilage matrix degradation in chondrocytes, providing a potential therapeutic target for osteoarthritis.20 Notably, many of these miRNAs were reported in cancer biology, and need further studies for their roles in COPD.

Recent studies have shown that the expression of miRNA can be modulated by transcription factors,21,22 and thus we predicted transcription factors that can potentially regulate these candidate DEMs. Specificity protein 1 (SP1), a C2H2-type zinc-finger transcription factor,23 was predicted to account for the highest percentage of downregulated DEMs, and has been widely documented in COPD. In the year of 2012, Di et al24 showed that CS induced MUC5AC gene transcription predominantly through increased SP1 nuclear protein levels and increased SP1 binding to its promoter region. Following that, more experiments suggested that ERK-SP1 may be a promising target for inhibiting COPD-related inflammation and mucus hypersecretion.25–27 Yin Yang 1 (YY1) and E2F1 were predicted to regulate the expressions of both upregulated and downregulated DEMs. In the study of YY1, a transcriptional repressor that inhibits muscle gene expression and myogenesis, Natanek et al28 found that YY1 expression was related to smaller quadriceps fiber cross-sectional area (CSA), and its nuclear localization was in muscle of COPD patients but not controls, suggesting an important role of YY1 in COPD-related muscle atrophy. E2F1, a key proliferative factor,29,30 was found to be targeted by miR-197 and participate in the vascular remodeling of COPD.12 In fact, despite of ZFP161 and SP4, the remaining transcription factors have been reported in COPD,31–33 which in turn supports the importance of these candidate DEMs in the pathogenesis of COPD.

KEGG analysis showed that the DEMs were mainly enriched in PI3K-Akt signaling pathway, mTOR signaling pathway, and autophagy. In the previous studies, the activation of PI3K-Akt signaling pathway has been amply documented to have an association with COPD-related inflammation. A study performed by Horiguchi et al34 showed that pulmonary administration of Wortmannin, an inhibitor of PI3K, was capable of inducing the repair of alveolar injury and a significant recovery of the respiratory function. Conversely, the activation of mTOR signaling pathway was reported to suppresses CS-induced inflammation and emphysema development, likely through modulation of autophagy, apoptosis, and necroptosis.35 Wu et al36 identified that increased levels of autophagy in peripheral blood mononuclear cells (PBMCs) of COPD patients were negatively correlated with FEV1% predicted, and positively correlated with circulating expressions of inflammatory cytokines. Therefore, inhibition of PI3K-Akt, activation of mTOR and inhibition of autophagy may represent novel therapeutic strategies for COPD. Although the numbers of genes targeted by either upregulated or downregulated DEMs were similar, however, many of the target genes of upregulated DEMs were not involved in GO and KEGG pathways (results were not shown), including the transcription factors described above. This may hint that downregulated DEMs may play a more critical role in the pathogenesis of COPD.

By constructing DEM-hub gene network, we found that most of hub genes could be potentially targeted by miR-497-5p, miR-130b-5p, and miR-126-5p. Among the top 12 hub genes, only two genes (MYC, FOXO1) expressions were in accordance with that in the GSE56768 dataset, which may result from different sample sources. MYC is well known to play a vital role in the pathogenesis of squamous cell lung cancer (SCLC). In COPD, highly expressed c-MYC in lung tissue was reported to induce the pulmonary vessels structural remodeling.37 However, in our results, the expression level of MYC in plasma was opposite to that in lung tissue. Volckaert et al38 showed that conditional deletion of c-Myc from parabronchial smooth muscle cells (PSMCs) inhibited airway epithelial repair, and thus plasma MYC downregulation may explain the dysregulation of lung repair in COPD. FOXO1, a member of the forkhead-box O family, can act as a homeostasis regulator and coordinate responses to various environmental changes, including oxidative stress.39,40 Liu et al41 identified that AMPK alleviated endoplasmic reticulum stress in COPD through activation of FOXO1. In agreement with the study of Hueso et al,42 FOXO1 expression of our analysis was decreased in COPD plasma compared with normal plasma, suggesting that miR-130b-5p-FOXO1 may participate in the pathogenesis of COPD through excessive oxidative stress.

Although we for the first time investigated the potential miRNA–mRNA regulatory network in COPD plasma by integrating multiple microarray datasets, however, there were several limitations existing in our present study. Firstly, we only focused on miRNAs and mRNAs between COPD and normal plasma, however, some of these may vary in different stages of COPD, which need to be further identified. Secondly, the sample size of each dataset used here was not enough, compared with that usually required for biomarker analysis. Thirdly, the miRNA–mRNA interactions were only based on predictions from public databases, and further studies with experimental validations in vivo and in vitro are still needed to validate our analysis.

Conclusion

In summary, based on the GEO database and bioinformatics analysis, we firstly constructed two potential miRNA–mRNA pathways (miR-126-5p-MYC, miR-130b-5p-FOXO1) in COPD plasma as potential biomarkers for the diagnosis and treatment of COPD patients. We hope that these findings can contribute to future in-depth studies and achieve the goal of improving the prognosis of patients with COPD.