Analysis of Key Genes and miRNA-mRNA Networks Associated with Glucocorticoids Treatment in Chronic Obstructive Pulmonary Disease

Introduction

Chronic obstructive pulmonary disease (COPD) is a chronic inflammatory disease of the airway that is characterized by chronic inflammation of the airway, lung parenchyma, and pulmonary vessels. Alveolar macrophages, T lymphocytes (especially CD8+), and neutrophils are increased in different parts of the lung, and eosinophils are increased in some patients.¹

Among the drugs for the treatment of chronic airway inflammatory diseases, glucocorticoids (GCs) are recognized as the most rapid and potent anti-inflammatory drugs. They have achieved good clinical efficacy in the treatment of bronchial asthma^2,3 and COPD.^4–11 In patients with moderate to severe acute exacerbation of COPD (AECOPD), systemic use of GCs improves forced expiratory volume in one second (FEV₁) and oxygenation status and shortens rehabilitation and hospitalization time.¹² Compared with systemic GCs, inhaled corticosteroids (ICSs) have fewer adverse reactions and can replace or partially replace systemic GCs.^13–15

However, some studies have shown that patients with COPD do not benefit from inhaled corticosteroid (ICS) therapy alone. Specifically, ICS treatment cannot effectively delay the long-term decline in FEV₁ or reduce the mortality of patients with COPD.^16,17 GCs resistance caused by HDAC2 activity decreases, dysfunction of corticosteroid glucocorticoid receptor-α and GR-β increases in epithelial cells, etc., may be a reason why COPD patients cannot benefit from GC treatment.^18–20 Recent studies have shown that microRNAs are associated with GCs resistance in COPD.²¹

MicroRNAs (miRNAs) are a kind of noncoding single-stranded RNA that are conserved and widely distributed in eukaryotes. The size is approximately 19~25 nucleotides, which mainly plays a posttranscriptional regulatory role by inhibiting the expression or translation of target genes.^22,23 Studies have shown that multiple miRNAs are involved in the occurrence and development of COPD.^24–30 In patients with moderate to severe COPD, short-term and long-term treatment with ICS can affect the expression of miRNAs.³¹ However, there are few reports on the molecular mechanism of miRNA-mediated COPD, especially in the treatment of COPD by GCs. The regulatory mechanism and molecular network of GC-mediated miRNAs regulating mRNA expression and inhibiting the COPD inflammatory response are not clear. This study intended to mine key genes and construct a miRNA-mRNA regulatory network by analysing miRNA and mRNA expression microarray datasets related to GC treatment of COPD in the GEO database. To explore its mechanism of action and provide an important theoretical reference and scientific basis for GC treatment of COPD.

Materials and Methods

miRNA and mRNA Microarray

The microarray data of the miRNA and mRNA expression profiles related to GC treatment of COPD were retrieved from the GEO database of NCBI. The screening criteria were as follows: ① patients with COPD; ② the experimental group was treated with GCs (ICS or oral GCs), and the control group was treated with placebo. Finally, the miRNA microarray dataset GSE76774 and mRNA microarray dataset GSE36221 that met the requirements were downloaded. Both studies were derived from the GLUCOLD study in the Netherlands. Based on the GPL8786 platform (Affymetrix multispecies miRNA-1 array), GSE76774 included 69 patients with moderate and severe COPD (smokers or ex-smokers). Bronchopulmonary biopsy was performed before treatment (22 cases in the treatment group and 23 cases in the control group), 6 months (17 cases in the treatment group and 24 cases in the control group), and 30 months (17 cases in the treatment group and 20 cases in the control group) after treatment. Based on the GPL6244 platform (Affymetrix Human Gene 1.0 ST Array), GSE36221 included 90 patients with moderate and severe COPD (smokers or ex-smokers). Bronchopulmonary biopsy was performed before treatment (19 cases in the treatment group and 21 cases in the control group), 6 months (23 cases in the treatment group and 17 cases in the control group), and 30 months (19 cases in the treatment group and 17 cases in the control group) after treatment. The database or software involved in the research is shown in Table 1.

Table 1 Database or Software Involved in the Research Design

Methods

Data Processing and Differential Expression Analysis

The microarray data of miRNA and mRNA expression profiles before treatment and 6 months and 30 months after treatment were analysed by the limma package of R software. The screening criteria before treatment were set as p > 0.05; the screening criteria after 6 months and 30 months of treatment were set as p < 0.05 and | log2-fold change (FC) | > 0.3. The differentially expressed data were visualized by a volcano map and cluster map.

Functional Enrichment Analysis of DEGs

Using the clusterProfiler package in R software, the DEGs analysed by GSE36221 were analysed by Gene Ontology (GO) and pathway enrichment analysis of Kyoto Encyclopedia of Genes and Genomes (KEGG). Among them, GO analysis included biological process (BP), cellular component (CC), and molecular function (MF), with P < 0.05.

Construction of the DEG Protein Interaction Network

To further identify the relationship between different genes, the STRING (search tool for the retrieval of interacting genes) database was used for protein–protein interactions (PPIs), and the confidence score was set to > 0.4. The PPI network was visualized by Cytoscape software. The key genes were identified through the “cytohubba” module. The data were imported from the network result file obtained from the string database, and the degree algorithm was selected to obtain the top 25 key genes.

Prediction of miRNA Target Genes and Construction of the miRNA-mRNA Regulatory Network

The miRWalk Database was used to predict the target genes of differentially expressed miRNAs (DEMs), and the intersection between the predicted target genes and DEGs analysed by GSE36221 was obtained by Venny 2.1. According to the regulatory relationship between miRNAs and mRNAs, the relationship between miRNAs and mRNAs was clarified. The miRNA-mRNA regulatory network was visualized by Cytoscape software.

Cell Experiment

To verify the experimental results, miR-320d-2, miR-181a-2-3p, TFCP2L1, and SLAMF7 were randomly selected for validation using a COPD cell model.

Materials

BEAS-2B human bronchial epithelial cells were purchased from Beina Biotechnology Co., Ltd. (BNCC359274), Daqianmen cigarettes were purchased from Shanghai Tobacco Co., Ltd., TFCP2L1 antibody, SLAMF7 antibody were purchased from Proteintech Co., Ltd.

Method

Cell culture After the cell concentration was adjusted to 5×10⁴/mL by complete medium, the cells were uniformly inoculated on 96-well plates at 100ul/ well, and cultured at 37°C and 5%CO2 for 24h.

Cigarette Smoke Extract (CSE) Preparation

One unfiltered burning cigarette (consisting of 12 mg of tar, 0.9 mg of nicotine, 14 mg of carbon monoxide) was connected to a 50mL syringe through a rubber hose. 300mL of smoke was collected through continuous negative pressure suction. The smoke was then injected into an airtight container containing 25mL DMEM of high-sugar medium and thoroughly mixed. Then 1mol/L of NaOH was used to adjust the PH of the mixture to 7.4, and finally 0.22um filter was used to remove bacteria. The mixture at this point was defined as 100% CSE solution. The optical density value of CSE solution was detected with A UV spectrophotometer at the wavelength of 320nm, so that the optical density value of different batches was about 0.2.

CSE Concentration Screening

BEAS-2B cells were incubated with 3%, 5%, 8%, and 10% CSE solutions for 12, 24, and 48 hours, respectively. 0% CSE was used as the control group. Each group has 6 replicates. CCK8 solution was added with a dose of 10 μL/well and incubated for 1–3 hours. Finally, the OD value was measured at 450nm on the enzyme marker. According to the results, 5% CSE was selected for the next experiment.

Dexamethasone Concentration Screening

BEAS-2B cells were incubated with 10⁻¹⁰mol/L, 10⁻⁹mol/L, 10⁻⁸mol/L, 10⁻⁷mol/L and 10⁻⁶mol/L dexamethasone solutions for 12, 24, and 48 hours, respectively. 0% dexamethasone was used as the control group. The remaining experimental methods were the same as before. According to the results, 10⁻⁸mol/L dexamethasone was selected for the next experiment.

Cell Grouping and Medication

Cell was divided into control group (Control), model group (Model) and dexamethasone group (DXM). The control group was cultured in complete medium (90% DMEM+10% FBS+1% dual antibody) for 24 hours. The model group was cultured with 5%CSE for 24h. The dexamethasone group was cultured with 5%CSE for 12h, and then cultured with 10⁻⁸mol/L dexamethasone for 12h.

Western Blot

The protein was extracted by the Protein Extraction Kit, and then quantified by the BCA method. The expression levels of TFCP2L1 and SLAMF7 in cells were detected by the WB method. The internal reference was GAPDH. The experiment was repeated three times.

RT-PCR

Total RNA of each group was extracted by Trizol reagent, then reverse transcribed into cDNA and amplified by PCR. Amplification program: 95 °C for 5 minutes, followed by 95 °C for 10 seconds, continued at 60 °C for 30 seconds, and finally 72 °C for 30 seconds, a total of 35 cycles. The expression levels of miR-320d-2, miR-181a-2-3p, TFCP2L1, and SLAMF7 were calculated by 2 ^−ΔΔCT method. The primer sequences are shown below: miR320d-2 Forward: AAAAGCTGGGTTGAGAGGA; miR181a-2-3P Forward: ACCACCGACCGTTGACTGTACC; U6-F Forward: CTCGCTTCGGCAGCACA, U6-R Reverse: AACGCTTCACGAATTTGCGT; SLAMF7 Forward: CAAGGTGTTCAAGCCGAAGG, SLAMF7 Reverse: GAGCACTCTGTGAGGATGGT; TFCP2L1 Forward: CTTGTCACCATACAGCCAGAAG; TFCP2L1 Reverse: ACATGCAGCACGTACTCCT; GAPDH Forward: CCTCTGACTTCAACAGCGACAC; GAPDH Reverse: TGGTCCAGGGGTCTTACTCC.

Statistical Analysis

The data was statistically analyzed by SPSS 25.0 and visualized by GraphPad Prism 8.0. P<0.05 indicated a statistically significant difference.

Results

DEMs

Before treatment, there were 816 miRNAs with no difference between the two groups. After six months of treatment, compared with the control group, there were 31 DEMs in the treatment group (10 upregulated and 21 downregulated). See Figure 1A and B for a volcano diagram and cluster diagram. After 30 months of treatment, compared with the control group, there were 23 DEMs in the treatment group (13 upregulated and 8 downregulated). See Figure 1C and D for volcano diagram and cluster diagram. Venny 2.1 was used to obtain DEMs with sustained, co directional,and differential expression. Four DEMs were ultimately obtained, including one upregulated expression (miR-320d) and three downregulated expression (miR-650, miR-155-5p, miR-181a-3p), as shown in Table 2.

Table 2 Continuously Express miRNAs in the Same Direction

Figure 1 Volcano map and Heatmap of DEMs (A) Volcano map of DEMs after six months of treatment (B) Heatmap of DEMs after six months of treatment (C) Volcano map of DEMs after 30 months of treatment (D) Heatmap of DEMs after 30 months of treatment. Red: up-regulation; green: downregulation; DEMs, differentially expressed miRNAs.

DEGs

Before treatment, there were 22,666 mRNAs with no difference between the two groups. After six months of treatment, compared with the control group, there were 387 DEGs in the treatment group (132 upregulated and 255 downregulated). See Figure 2A and B for a volcano diagram and cluster diagram. After 30 months of treatment, compared with the control group, there were 279 DEGs in the treatment group (55 upregulated and 224 downregulated). See Figure 2C and D for volcano diagram and cluster diagram. Venny 2.1 was used to obtain DEGs with sustained, co directional, and differential expression. 83 DEGs (20 were upregulated and 63 were downregulated) were ultimately obtained. The main DEGs are shown in Table 3.

Table 3 Main DEGs

Figure 2 Volcano map and Heatmap of DEGs (A) Volcano map of DEGs after six months of treatment (B) Heatmap of DEGs after six months of treatment (C) Volcano map of DEGs after 30 months of treatment (D) Heatmap of DEGs after 30 months of treatment. Red: up-regulation; green: downregulation. Abbreviation: DEGs, differentially expressed genes.

DEGs Function Analysis Results

The GO and KEGG pathway enrichment analyses of 83 DEGs were carried out by the clusterProfiler package of R software. A total of 218 GO enrichment entries were obtained, including 66 BP entries, 5 CC entries, and 11 MF entries. BP was mainly enriched in the production of molecular mediator of immune response; positive regulation of cell activation; immunoglobulin production; positive regulation of leukocyte activation; B-cell activation; lymphocyte differentiation; immune response-regulating cell surface receptor signalling pathway; mononuclear cell differentiation, etc. CC was mainly enriched in immunoglobulin complex; external side of plasma membrane; immunoglobulin complex, circulating; blood microparticle, etc. MF was mainly enriched in antigen binding; immunoglobulin receptor binding; immunoglobulin binding; glycosaminoglycan binding; integrin binding, etc. See Figure 3 for the first five BP, CC, and MF entries. KEGG pathway enrichment analysis identified four signalling pathways: PI3K/Akt signalling pathway; ECM receptor interaction; primary immunodeficiency; and adhesive spots (Figure 4).

Figure 3 Bar chart of GO enrichment analysis of DEGs.

Figure 4 Bar chart of KEGG enrichment analysis of DEGs.

Construction of the DEG PPI Network and Screening of Key Genes

The string platform was used to construct PPI networks for 83 DEGs, as shown in Figure 5. The results show that 58 node proteins (nodes) play a direct or indirect role in the treatment of COPD with GC, involving a total of 59 interactions (edges). Cytoscape was used for network topology analysis. According to the screening conditions, key genes such as CD2, KIT, FCER1A, SLAMF1, SLAMF 7, GPX3, and PRELP were obtained, as shown in Figure 6.

Figure 5 DEGs PPI network.

Figure 6 Key genes.

MiRNA Target Gene Prediction and Regulatory Network Construction

The miRWalk database was used to predict the downstream targets of mRNA, and a total of 6378 genes were obtained. A total of 18 overlapping genes were obtained from the intersection with the above 83 DEGs, as shown in Table 4. According to the negative regulatory relationship between miRNA and mRNA, nine pairs of miRNA-mRNA pairs were finally selected. The miRNA-mRNA regulation network was constructed and visualized by Cytoscape software, as shown in Figure 7. The flow chart of the whole research method and results are shown in Figure 8.

Table 4 Intersection Genes of miRNA Target Genes and DEGs

Figure 7 miRNA-mRNA regulatory networks.

Figure 8 Flow chart.

Cell Experiment

We determined the expression of SLAMF7 and TFCP2L1 using WB, and the expression of SLAMF7, TFCP2L1, miR-320d-2 and miR-181a-2-3p using PCR The study showed that the expression of miR-320d-2 and TFCP2L1 were low in the COPD group, and the expression of miR-181a-2-3p and SLAMF7 were high in the COPD group. After dexamethasone treatment, the expression of miR-320d-2 and TFCP2L1 were upregulated, while the expression of miR-181a-2-3p and SLAMF7 were downregulated (Figures 9 and 10). The results of the cell experiment were consistent with previous studies.

Figure 9 WB of Partial miRNA-mRNA. **, ##P value<0.01, ***, ###P value<0.001.

Figure 10 PCR of Partial miRNA-mRNA (A) miR-181a-2-3p (B) miR-320d-2 (C) SLAMF7 (D) TFCP2L1. **, ##P value<0.01, ***, ###P value<0.001.

Discussion

In this study, 83 DEGs (such as CD2, KIT, FCER1A, SLAMF1, SLAMF7, and NMI), 4 DEMs (miR-320d, miR-155-5p, miR-650, and miR-181a-3p) and 9 miRNA-mRNA regulatory pairs (such as miR-320d-SLAMF1, miR-320d-SLAMF7, miR-181a-2-3p-TFCP2L1) related to the treatment of COPD by GCs were screened, and a related biological pathway, the PI3K/Akt signalling pathway, was obtained (Figure 11).

Figure 11 Key genes and miRNA-mRNA networks associated with glucocorticoids treatment in COPD.

This study found that miRNAs such as miR-320d, miR-155-5p, miR-650, and miR-181a-2-3p changed continuously during the treatment of COPD by GCs, in which the expression of miR-650, miR-155, and miR-181a-3p were downregulated and the expression of miR-320d was upregulated. MiR-155 and miR-320d are associated with the inflammatory response^24,25 and oxidative stress³² in COPD. Clinical studies have shown that compared with nonsmokers, the expression of miR-155 in lung tissue²⁵ and induced sputum²⁴ of patients with COPD. The expression of miR-155 increased in lung tissue,²⁵ alveolar macrophages,²⁵ and bronchoalveolar lavage fluid²⁴ of mice exposed to cigarette smoke (CS). However, the inflammatory cells, cytokines, and chemokines in the bronchoalveolar lavage (BAL) of miR-155 knockout mice exposed to CS decreased, and lung inflammation induced by CS decreased significantly. Moreover, this inflammation was alleviated by a nasal drip of a specific miR-155 inhibitor.²⁵ miR-155 is also associated with oxidative stress in COPD.³³ In model rats of COPD, the effective components of the Bufei Yishen recipe can inhibit oxidative stress induced by PM2.5 through the miR-155/FOXO3a pathway.³³ Clinical studies have also confirmed that miR-320d is related to the inflammatory response in COPD. Compared with nonsmokers, mir-320d is upregulated in plasma-derived exosomal RNA in patients with COPD.³² MiR-320d is associated with the inflammatory response of the airway epithelium in COPD.³¹ It has been demonstrated that GCs play an anti-inflammatory role in COPD³¹ and asthma³⁴ through miR-155 and miR-320d. MiR-320d is a novel ICS mediator that regulates the proinflammatory response of the COPD airway epithelium.³¹ MiR-155 is the target of GCs in the treatment of allergic asthma. The anti-inflammatory effect of GCs can be reversed by inhibiting the expression of miR-155.³⁴ However, whether miR-155 has the same effect in theanti-inflammatory treatment of COPD remains to be confirmed. GCs may affect the inflammatory response and oxidative stress of COPD by regulating expression of multiple miRNAs, such as miR-155 and miR-320d, but its specific intervention mechanism is unknown.

This study found that DIO2, TFCP2L1, DCLK1, GALNT15, PER1, SLAMF7, PRELP, SLAMF1, TSC22D3, NMI, METTL7A, MS4A2, FAM107A, IL33, CD96, GPX3, TIMP4, APOD, and other key genes were continuously changed during COPD treatment by GCs, in which the expression of DIO2, DCLK1, SLAMF7, SLAMF1, NMI, MS4A2, IL33, and CD96 were downregulated, and the expression of TFCP2L1, GALNT15, PER1, PRELP, TSC22D3, METTL7A, FAM107A, GPX3, TIMP4, and APOD were upregulated. Studies have shown that GPX3,^35,36 TSC22D3³⁷, DIO2,³⁸ Per1,^39,40 IL-33,⁴¹ and other genes are related to the inflammatory response and oxidative stress of COPD. The expression of GPX3 decreased in HBE cells of COPD,³⁵ lung tissue of patients with severe emphysema,³⁶ and lung tissue of patients with advanced COPD.³⁵ In the treatment of COPD, PPARγ agonists may play an antioxidant role by blocking the downregulation of GPX3. Glucocorticoid-induced leucine zipper (TSC22D3) is a regulator of GC anti-inflammatory effects.³⁷ Beclomethasone-17-monopropionate (17-BMP) can upregulate the expression of TSC22D3 in alveolar macrophages of COPD patients and smokers.⁴² TSC22D3 can interact with multiple pathways, including nuclear factor-κB (NF-κB), activator protein 1 (AP-1), Raf-1, and Ras.³⁷ GC-induced leucine zipper (TSC22D3) can inhibit the transcriptional response of NF-κB and AP-1,⁴³ and its mechanism is related to the increase in TSC22D3 expression induced by GCs through transactivation and inhibition of NF-κB signalling, resulting in proinflammatory cytokine release from airway epithelial cells.⁴⁴ Per1 is an important part of the biological molecular oscillation system. Compared with nonsmokers, the level of Per1 in lung tissue and PBMCs in patients with COPD increased. The oscillation and dysregulated expression of biological clock genes such as BMAL1, Per1, and Per2 are related to the aetiology, chronic inflammation, and imbalance of autophagy levels in COPD.^39,40 IL-33 may become a new target for COPD therapy.⁴¹ Elevated levels of IL-33 in lung biopsy samples, epithelial and endothelial cells, serum or plasma, and sputum in patients with COPD are associated with decreased lung function.⁴¹ Nonhaematopoietic lung cells are not only an important cell source ofIL-33 but also express ST2 and respond to IL-33 stimulation (regardless of source), thus participating in lung inflammation and tissue remodelling.⁴⁵ The risk of COPD associated with the loss of function of IL-33 (rs146597587) is reduced. In contrast, the increased function of IL-33 and IL1RL1 variants increases the risk of COPD.⁴⁶ Compared with placebo, itepekimab targeting IL-33 can provide COPD patients with previous smoking potential benefits in reducing the frequency of exacerbation and improving lung function.⁴¹

In this study, GO and KEGG analyses were carried out for the screened DEGs. The results showed that the effect of GCs on COPD was mainly related to the production of molecular mediators of the immune response, positive regulation of cell activation, the PI3K/Akt signalling pathway, and so on. The therapeutic effect of GCs on COPD may be related to the intervention of the PI3K/Akt pathway and regulation of its downstream target genes through miR-155 and miR-320d. Although this hypothesis has not been confirmed in COPD, miRNAs have been widely confirmed to target and regulate their downstream genes by intervening in the PI3K/Akt pathway. Multiple miRNAs, such as miR-155,³⁰ miR-320d,⁴⁷ miR-650,⁴⁸ and miR-181,^49,50 can be used as upstream factors to regulate the PI3K/Akt pathway. The PI3K/Akt pathway is related to SLAMF7, NMI, TIMP-4, and other genes.^48–50 MiR-155 plays a role through the PI3K/Akt pathway in an LPS intraperitoneal injection-induced sepsis mouse model.³⁰ Curcumin inhibits the inflammatory response of macrophages in an LPS-induced sepsis mouse model, and its mechanism is related to curcumin interfering with the PI3K/Akt pathway through miR-155.³⁰ The PI3K/Akt signalling pathway is related to SLAMF7, NMI, TIMP-4, and ETC.^48–50 Studies have shown that SLAMF7 is expressed in LPS-stimulated monocytes via PI3K and NF-κB-dependent pathways and inhibits the production of proinflammatory cytokines, TNF-α and IL-12p70.⁵¹ SLAMF7 regulates the activation of NK cells by recruiting EAT-2 and activating the PI3K and phospholipase Cgamma signalling pathways in human NK cells.^52,53

This study found that there were nine pairs of miRNA-mRNA related to the treatment of COPD by GCs: miR-320d-SLAMF1, miR-320d-SLAMF7, miR-320d-NMI, miR-181a-2-3p-TFCP2L1, miR-650-FAM107A, miR-650-GPX3, miR-650-PRELP, miR-650-TIMP4, and miR-650-APOD. The in vitro experiments also confirmed the above experimental results. These DEMs and DEGs are the key nodes of GCs in the treatment of COPD.

Although this study constructed a potential miRNA-mRNA regulatory network and key genes for GC treatment of COPD based on bioinformatics, there are certain limitations: ①The reports on miR-650 and miR-181a-2-3p mainly focus on tumours, but there are no reports in COPD, which needs further research. ②Patients with stable COPD mostly use ICS/LABA, and some patients can benefit from it. The effect of ICS/LABA on miRNA and mRNA needs to be further studied.

Conclusion

GCs may intervene the expression of PI3K/Akt signal pathway through mir-320d, miR-320d and mir-155-5p, and target its downstream signal factors, to play an anti-inflammatory role in the treatment of COPD.