Identifying miRNA Modules and Related Pathways of Chronic Obstructive Pulmonary Disease Associated Emphysema by Weighted Gene Co-Expression Network Analysis

Introduction

Chronic obstructive pulmonary disease (COPD) is a major health burdens worldwide causing high morbidity and mortality.¹ The disease is predicted to become one of the top four causes of death globally by 2030.² COPD is a chronic inflammatory disease characterized by progressive airflow limitation involving peripheral loss of lung parenchymal tissue (emphysema) and mucus hypersecretion with chronic bronchitis.¹ These pathophysiological changes result in a progressive decline of forced expiratory volume in the first second, after full inspiration (FEV1), inadequate lung emptying on expiration, and subsequent static and dynamic hyperinflation.³ COPD is a heterogeneous and complex disease caused by diverse cellular and pathophysiological changes under distinct genetic backgrounds.⁴ Studies have associated single genes, such as the gene encoding matrix metalloproteinase 12 (MMP-12) and glutathione S-transferase, to a decline in lung function⁵ and an increased COPD risk.⁶ Genome-wide association studies have linked genetic loci (including FAM13A, markers near the alpha-nicotinic acetylcholine receptor, hedgehog interacting protein (HHIP), and others) with COPD^7–11 These gene expression and polymorphism studies have shown that COPD is a polygenic disease involving multiple pathogenic signaling pathways.

MicroRNAs (miRNAs) are a class of small, single‑stranded noncoding RNAs. By binding to the 3′terminal noncoding region of target mRNAs, miRNAs regulate gene expression at the transcriptional level in the form of complete complementation or incomplete complementation.^12,13 miRNAs affect biological processes, including differentiation, proliferation, and apoptosis, through their regulation of gene expression under normal and pathological conditions.^14,15 Studies have suggested pathophysiological roles for miRNAs in COPD: For example, Cao et al detected increased expression levels of miR-183, miR-200b, and miR-200c in the blood of patients with COPD¹⁶ and, Donaldson et al found that plasma levels of muscle-specific miR-499, miR-133, and miR-206 were elevated in patients with COPD.¹⁷ By using differential gene expression analysis, these studies focused on one or several miRNAs and only revealed simple up- and down-regulation, which cannot build the relationships between miRNAs and disease.

The weighted gene co-expression network analysis (WGCNA) is a systems biology method for describing the correlation patterns among genes across microarray samples.¹⁸ The method is based on the theory that genes with closely functional linkages or involved in similar pathways may have similar expression profiles.^19,20 By using information from thousands of genes with greatest changes, WGCNA identifies interest gene clusters (which we call gene modules). Next, WGCNA identifies associations between modules and external clinical traits, providing an effective way to explore the mechanisms behind certain traits.¹⁸ WGCNA has been widely used in cancer research.^21–23 Qin et al used WGCNA to investigate and contrast the molecular processes differing between the bronchiolitis and emphysema phenotypes of COPD.²⁴

The aim of this study was to use the WGCNA method to construct co-expression modules for the miRNA expression data from human biological specimens and to identify correlations between different modules and clinical COPD traits. This methodology provides an effective way to explore the role of miRNAs from different perspectives and identifies promising diagnostic biomarkers and molecular COPD pathophysiological pathways.

Methods

Microarray Data Used

We downloaded the miRNA expression profile in the GSE33336 dataset from the Gene Expression Omnibus (GEO) public functional genomics data repository of the National Center of Biotechnology Information (www.ncbi.nlm.nih.gov/geo/), platform GPL6955. To obtain this dataset, researchers set up Agilent human miRNA oligo arrays (8x15K array, Part No. G4470A, Agilent Technologies) based on miRBase V9.1 on lung tissues from 29 emphysematous lungs obtained from patients undergoing curative resection for lung cancer.²⁵ Briefly, smokers with at least a 20-pack-year history who quit smoking at least 10 months prior to surgery were considered. All patients had FEV1/ (forced vital capacity) FVC ratio <0.70, indicating the presence of COPD. Thus, 29 former smokers with COPD were selected and classified according to gas transfer measurement results (single breath carbon monoxide diffusion coefficient, KCO) as having “mild” emphysema (>75% predicted KCO) or “moderate” emphysema (40–75% predicted KCO). Specimens from patients on inhaled or oral steroids, with lung pathologies that might confound spirometry measurements, or with alpha 1 antitrypsin (α1AT) enzyme deficiency were excluded. Detailed information can be found in other publications.^25,26

Identification of Differentially Expressed miRNAs (DEMs)

For the original array data, the context correction and normalization were performed using a robust multiarray average process. The data were subsequently converted into expressive measures using the R affy package. DEMs were identified using the Limma package. The cutoff values for DEM screening using the Benjamini & Hochberg procedure were a p ≤0.05 and an absolute Log₂Ratio ≥1.

WGCNA

The WGCNA package is a comprehensive collection of R functions for weighted correlation network analysis.^18,27 Datasets were constructed using a previously described method.¹⁸ First, we constructed a signed weighted correlation network by creating a matrix of pairwise correlations between all pairs of miRNAs chosen according to their variance. The resulting Pearson correlation matrix was transformed into a matrix of connection strengths (eg, an adjacency matrix). Next, we applied the R pickSoftThreshold function that performs an analysis of network topology and chose a proper soft-thresholding power of 4.²⁸ The topological overlap was calculated to measure the network interconnectedness. Average linkage hierarchical clustering was used to group genes based on the topological overlap dissimilarity measure (1-topological overlap) of their network connection strengths. Using a dynamic tree-cutting algorithm and merging the threshold function at 0.25, we identified 25 modules in the datasets. We generated heatmap plots to visualize the co-expression module structure of topological overlap in the gene network. The associations among modules were summarized by a hierarchical clustering dendrogram of the eigengenes and by a heatmap plot of the corresponding eigengene network.

Association of Modules to External Traits

To identify modules that were significantly associated with the traits of gender, age, COPD status (mild or moderate), predicted FEV1, predicted KCO, and smoking pack years we correlated the module eigengenes (MEs; ie, the first principal component of a module)²⁹ with external traits using Pearson’s correlation coefficients, modules with p-values <0.05 were identified as trait-related modules.

Identification of Hub miRNAs and Construction of Protein–Protein Interaction

Module membership (MM) represents the correlation between the gene expression profile and the module eigengene. We considered the top 40 miRNAs with the highest module membership values in the modules as hub miRNAs and visualized them with Cytoscape. We obtained the target genes from the database “microWalk” (http://mirwalk.umm.uni-heidelberg.de). These target genes have been validated by existing experimental procedures. The target genes of hub microRNAs were used to construct a protein–protein interaction (PPI) network in the STRING database (https://string-db.org/) to identify associated functional proteins and important microRNAs.

Functional Annotation and Enrichment Analysis

We performed functional annotations with the ClusterProfiler based on the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases to assign biological functions to the miRNAs in these modules. Target genes for the miRNAs in each module were subjected to this analysis. We visualized the top 20 most enriched function terms (ordered by their p-values).

Results

Identification of DEMs in COPD

The miRNA-sequencing data from the GEO dataset were subjected to differential expression analysis based on a threshold of p ≤0.05 and an absolute Log₂Ratio ≥1. We found 422 DEMs that were significantly downregulated and 82 DEMs that were significantly upregulated in the samples from patients with moderate COPD samples when compared with the samples from patients with mild COPD (Table S2). The Additional file 1 (Figure S1) displays volcano plots of the miRNAs.

Identification of miRNA Coexpression Networks and Modules

We selected 9800 probes for network construction. miRNAs exhibiting similar patterns of expression were grouped into modules via hierarchical average linkage clustering. The sample dendrogram and trait heatmap grouped the selected samples into different clusters and provided the clinical trait data distribution map (Figure 1). Network topology was analyzed using various soft threshold powers at 4 (Additional file 2: Figure S2), according to the method by Zhang and Horvath.³⁰ After using a dynamic tree cutting algorithm (Figure 1), we obtained 25 distinct co-expression modules; 43 uncorrelated miRNAs were assigned into a grey module ignored in the following study. We used the Pearson correlation coefficient to analyze the interactions of these co-expression modules (Figure 1). The darker background indicates a higher module correlation.

Figure 1 Sample clustering and module detection. Notes: (A) Sample cluster tree and trait indicator. The leaves of the tree correspond to the sample. The red color represents female or male, and COPD status. The color intensity was proportional to older age, higher FEV1, higher KCO and higher smoke-pack-year. (B) Clustering dendrogram of miRNAs, with dissimilarity based on topological overlap, together with assigned module colors. (C) Visualizing the miRNA network using a heatmap plot. The heatmap depicts the Topological Overlap Matrix (TOM) among all miRNAs in the analysis, light color represents low overlap and progressively darker red color represents higher overlap. Blocks of darker colors along the diagonal are the modules. The miRNA dendrogram and module assignment are also shown along the left side and the top.

Association of miRNA Co-Expression Modules with Clinical Traits

We correlated the 25 modules with traits of interest and searched for the most significant associations to assess the physiological significance of the modules. In the heatmap of module-trait correlation (Figure 2), rows indicate modules and columns indicate external traits. miRNAs clustered in the royal blue module had the strongest negative correlation with the predicted FEV1 (correlation, −0.41; p-value, 0.02), miRNAs clustered in the light yellow module had the strongest positive correlation with the moderate COPD status (correlation, 0.46; p-value, 0.008), miRNAs clustered in the grey60 module had the strongest negative correlation with the smoking pack years (correlation, −0.61; p-value, 3e-04). Thus, we considered mainly modules with strong correlation to interesting traits to ensure biological significance. The royal blue, light yellow, and grey60 modules were analyzed further: The additional file 6 (Table S1) shows the miRNAs in these three modules associated to clinical traits of interest.

Figure 2 Module-trait associations. Each row corresponds to a module eigengene, and each column to a trait. Each cell contains the corresponding correlation and p-value. The table is color-coded by correlation according to the color bar.

Functional Annotation and Enrichment Analysis

We performed functional annotation and enrichment analyses to explore the biological functions associated with the genes clustered in each module of interest (royal blue, light yellow, and grey60). The GO database analysis resulted in three main annotated categories (biological processes, cellular components, and molecular functions). For each module and each function database, we generated bar plots to visualize the top 20 most significantly enriched terms. The additional files 3–5 (Figures S3–S5) show the summarized results of the GO enrichment analysis with important findings (p-value < 0.05).

For the royal blue module, the top five significantly enriched biological processes pathways included intracellular signal transduction, intracellular protein transport, neutrophil degranulation, positive regulation of GTPase activity, and endocytosis. The top five significantly enriched cellular component pathways were glutamatergic synapse, neuron projection, membrane raft, neuronal cell body, and axon. Moreover, the top five significantly enriched molecular function component pathways were DNA-binding transcription activator activity (RNA polymerase II-specific), protein kinase binding, protein serine/threonine kinase activity, transcription regulatory region DNA binding, and protein domain specific binding. As shown in Figure 3, the KEGG annotation results showed axon guidance, proteoglycans in cancer, MAPK signaling pathway, hepatocellular carcinoma, and sphingolipid signaling pathway as the top five significantly enriched pathways.

Figure 3 Kyoto Encyclopedia of Genes and Genomes pathways enrichment. (A) royal blue module (B) light yellow module (C) grey60 module.

For the light yellow module, the top five significantly enriched biological processes pathways were intracellular signal transduction, intracellular protein transport, positive regulation of GTPase activity, brain development, and the Wnt signaling pathway. The top five significantly enriched cellular component pathways were glutamatergic synapse, neuronal cell body, neuron projection, membrane raft, and dendrites. Further, the top five significantly enriched molecular function component pathways were DNA-binding transcription activator activity (RNA polymerase II-specific), GTPase activator activity, Rab GTPase binding, protein serine/threonine kinase activity, and protein kinase binding. As shown in Figure 3, the KEGG annotation results showed MAPK signaling pathway, proteoglycans in cancer, axon guidance, hippo signaling pathway, and inositol phosphate metabolism as the top five significantly enriched pathways.

For the grey60 module, the top five significantly enriched biological processes pathways were intracellular signal transduction, positive regulation of GTPase activity, intracellular protein transport, regulation of small GTPase mediated signal transduction, and the Wnt signaling pathway. The top five significantly enriched cellular component pathways were glutamatergic synapse, neuron projection, neuronal cell body, early endosome membrane, and axon. In addition, the top five significantly enriched molecular function component pathways were GTPase activator activity, DNA-binding transcription activator activity, RNA polymerase II-specific, protein serine/threonine kinase activity, DNA-binding transcription repressor activity, RNA polymerase II-specific, and transcription regulatory region DNA binding. As shown in Figure 3, the KEGG annotation results showed proteoglycans in cancer, axon guidance, endocytosis, signaling pathways regulating pluripotency of stem cells, and MAPK signaling pathway as the top five significantly enriched pathways.

Identification of Hub miRNAs and PPI

Module memberships (MMs) represent each microRNA’s Pearson correlation coefficient with its corresponding eigengene module. We considered the top 40 miRNAs with the highest module membership values as hub miRNAs and used Cytoscape to visualize them (Figure 4 and Table S3). The hub microRNAs miR-452, miR-149, miR-133a, miR-181a, and miR-421 may play a role in emphysema-associated COPD. We used the target genes of these hub microRNAs for protein–protein interaction analysis to get a systematic perspective of the molecular mechanisms of this disease. In PPI network, the degree means the connection between proteins. The higher the degree, the stronger the connection is. Protein with highest degree is the top protein. As shown in Figure 5, the top three proteins in the royal blue module were KRT18, insulin-like growth factors-1 (IGF1), and KRT84 (their degrees were 26, 21, and 19, respectively). The top three proteins in the light yellow module were KRT71, KRT84, and KRT2 (their degrees were 11, 10, and 10, respectively). The top three proteins in the grey60 module were KRT18, KRT2, and KRT76 (their degrees were 21, 17, and 17, respectively). miR-520b targets the KRT18 protein in both the royal blue and the grey60 modules.

Figure 4 Visualization of the network connections among the most connected miRNAs in the royal blue module (A) light yellow module (B) and grey60 module (C) generated by the Cytoscape software. Edge weight indicates the TOM similarity between two nodes.

Figure 5 The protein–protein interaction (PPI) analysis of target gene of hub microRNAs in royal blue module (A) light yellow module (B) and grey60 module (C). Degrees describe the importance of the protein nodes.

Discussion

For this study, we applied a systems biology approach (WGCNA) to identify miRNA modules and associated pathways of emphysema-associated COPD for the first time. We identified 25 modules by WGCNA analysis and found that the royal blue, light yellow, and grey60 modules were significantly correlated to our traits of interest (predicted FEV1, moderate COPD status, and smoking pack years). These three modules were analyzed further.

Axon guidance, proteoglycans in cancer, and the MAPK signaling pathway were common processes in the top 20 KEGG pathways of the three modules. Axon guidance plays an important role during the development of neural circuits. Several receptors and pathways have also been shown to play roles in axon guidance.³¹ Slit guidance ligand (SLIT) and roundabout (ROBO) were shown to be involved in neural development by inhibiting the migration of axons.³² Subsequently, SLIT was shown to be a potential anti-inflammatory molecule.³³ Lin et al used an integrated bioinformatics analysis and found that the expressions of SLIT2 and ROBO2 genes are downregulated in patients with COPD and that the expressions are significantly negatively correlated with the COPD stages.³⁴ The Wnt pathway is an evolutionarily conserved family that plays a crucial role in axon and dendrite development.³⁵ Wnt signaling modulates several cellular processes involved in COPD, including oxidative stress, apoptosis/proliferation, inflammation, mucus hypersecretion, protease/antiprotease imbalance, autophagy, senescence, metabolism reprogramming, mitochondrial dysfunction and stem-cell/progenitor-cell renewal.³⁶

Proteoglycans are a family of charged molecules that contain a core protein and at least one glycosaminoglycan. They are primarily found on the cell surface and in the extracellular matrix, where they serve as structural components. They are thought to influence a myriad of cellular processes, including those linked with cancer development.³⁷ Proteoglycans are found in the extracellular matrix, in plasma membrane of cells, and as intracellular structures in the lungs.³⁸ Changes in the proteoglycan composition in the lungs have been reported in human lung diseases, such as idiopathic pulmonary fibrosis, acute respiratory distress syndrome (ARDS), COPD, lymphangioleiomyomatosis, and asthma.³⁹ Studies have shown that proteoglycans and their associated glycosaminoglycans play important roles modulating pulmonary inflammation.^40–42 Takahashi et al found that proteoglycans maintain lung stability in an elastase-treated mouse model of emphysema, suggesting that the loss of proteoglycans in human emphysema contributes to disease progression.⁴³

Mitogen-activated protein kinases (MAPKs) are serine-threonine protein kinases including c-Jun NH2-terminal kinase (JNK), p38 MAPK, and extracellular signal- regulated kinase (ERK). MAPKs regulate several cellular activities such as proliferation, differentiation, apoptosis, survival, inflammation, and innate immunity.⁴⁴ MAPK signaling pathways have also been found to play key roles in the pathogenesis of many diseases.⁴⁵ The JNK and p38 MAPK signaling pathways are activated by various types of cellular stress such as oxidative stress, bacterial lipopolysaccharide, and proinflammatory cytokines (tumor necrosis factor-α and interleukin 1β);⁴⁴ these factors also play a role in the pathogenesis of COPD.⁴ Patients with COPD have increased p38 MAPK activation in macrophages and in other cells within the alveolar wall. Studies have shown that p38 inhibitors exhibit anti- inflammation effects both in vitro and in vivo.⁴⁶ However, many p38 inhibitors have failed in clinical trials due to their unacceptable safety profiles.⁴⁷

In addition to KEGG enrichment analysis of the key modules in our study, we used the STRING database to generate a protein–protein interaction (PPI) network of some hub genes in key modules. The PPI network analysis results showed that KRT18 had the highest degree value in all the module networks. Cytokeratins are major structural proteins of epithelia that are essential for tissue integrity.⁴⁸ Cytokeratin 18 (CK18, also known as KRT18) is an acidic cytokeratin member of the protein family of cytoskeletal intermediate filaments.⁴⁹ It is widely expressed in normal epithelia, including in pulmonary bronchial and alveolar epithelial cells.^50,51 CK18 is involved in important signaling pathways including apoptosis, cell cycle, and cancer progression.⁴⁹ During apoptosis, specific caspase-cleaved CK18 (ccCK18) fragments are generated. Hacker and Lichtenauer et al demonstrated that elevated levels of ccCK-18 in the serum of COPD patients correlate with disease severity, and that hypoxic conditions may trigger the release of the ccCK-18 fragments.^52,53

Studies have shown that IGF-1 signaling plays an important role during lung development and diseases (such as congenital disorders, cancers, inflammation, and fibrosis).⁵⁴ IGF-1 signaling is involved in the pathogenesis of COPD by promoting inflammation, muscle dysfunction metabolism, and the aging process.⁵⁴

In this study, we identified several hub miRNAs that may play roles in the development of emphysema-associated COPD. miR-452, miR-149, miR-133a, miR-181a, and miR-421 had been found to be abnormally expressed in patients with emphysema,⁵⁵ pulmonary fibrosis,⁵⁶ pulmonary arterial hypertension,⁵⁷ and lung cancer.⁵⁸ The expression of miR-452 in alveolar macrophages is reduced in smokers; miR-452 targets matrix metalloproteinase-12 (MMP-12), which is associated with the development of emphysema.⁵⁵ Shen et al reported that treatment with cigarette smoke extract (CSE) decreases the level of miR-149-3p in a murine monocytic cell line (THP-1): Transfection of miR-149-3p mimic into CSE-induced THP-1 cells decreased the levels of the TLR-4 and NF-κB p65 proteins, which play an important role in COPD.⁵⁹ Transforming growth factor (TGF)-beta1-induced miR-133a could inhibit myofibroblast differentiation and pulmonary fibrosis.⁵⁶ Moreover, myofibroblast differentiation has been observed in patients with COPD.⁶⁰ Liu et al found that the expression of miR-133 was lower in patients with acute COPD exacerbations than in stable patients and healthy controls.⁶¹ Hypercapnia, a common occurrence in patients with COPD, was found to increase the airway smooth muscle contractility via caspase-7-mediated miR-133a-RhoA signals.⁶² miR-181a was found to target endocan and thereby ameliorate the inflammatory response in monocrotaline-induced pulmonary arterial hypertension.⁶³ Studies have found that endocan is associated with disease severity, exacerbation, and lung function decline in patients with COPD.^63–65 Yuan et al demonstrated that miR-421 induces apoptosis in the lung and increases the inflammatory response in a bronchopulmonary dysplasia mouse model.⁶⁶ Knock-down of miR-421 causes levels of KEAP to rise, inhibiting Nrf2-dependent antioxidant expression, and ultimately increasing the intracellular ROS levels in non-small cell lung cancer.⁵⁸ Nrf2 and oxidative stress play important roles in COPD pathogenesis.⁶⁷ The role of miR-421 in COPD needs further exploration; but we hypothesize that the profiles of miR-452, miR-149, miR-133a, miR-181a, and miR-421 may be associated with the inflammation, myofibroblast differentiation, oxidative stress, and apoptosis evident in COPD. This hypothesis needs be confirmed by others.

We identified a few miRNAs related to in interesting clinical traits. Jing et al found that miR-558 could reduce the damage of HBE cells exposed to cigarette smoke extract (CSE) by targeting TNF receptor superfamily member 1A (TNFRSF1A) and inactivating TAK1/MAPK/NF-κB pathway.⁶⁸ miR-558 is demonstrated to participate in the lung cancer cell growth, tumorigenesis, and invasion.^69,70 Li et al demonstrated that miR −587 acts as an oncogene in non-small-cell lung carcinoma.⁷¹ Several studies showed that miR-621, miR −639, miR −566, miR −769, miR −597 play a role in proliferation, metastasis, migration, and invasion in a variety of cancer cells.^72–76

We are aware of the limitations in this study. First, normal lung tissues were not included in the study and this omission may have caused biases. Second, we analyzed a single platform dataset; therefore, the relatively small sample size may not fully represent all patients with COPD. Third, all the COPD lung samples were from patients who also had lung cancer, and the presence of the cancer may have affected our results. Finally, we did not provide biological experimental verification of the hub miRNAs and their associations with genes. The results need to be verified by functional experiments as a next step of the research.

Conclusions

In this study, we identified three microRNA modules related to emphysema-associated COPD through WGCNA. Axon guidance, proteoglycans in cancer, and the MAPK signaling pathway may play a role in the pathological process of emphysema-associated COPD. KRT18, miR-133a, miR-181a, and miR-421 may be associated with emphysema-associated COPD according to our combined KEGG enrichment and PPI network analyses. These findings should be verified and extended in future experimental work.