Exploring the Intervention Mechanism of Effective-Component Combination of Bufei Yishen Formula III for Airway Epithelial Barrier Injury: A miRNA-mRNA Regulatory Network Analysis

Introduction

Chronic obstructive pulmonary disease (COPD) is a common heterogeneous pulmonary condition characterized by persistent (often progressive) airflow limitation and clinically manifests by chronic respiratory symptoms such as dyspnea, cough, and chest tightness.¹ The prevalence, morbidity, and mortality of COPD are substantial, making it a significant public health concern.^2,3 According to the Global Initiative for Chronic Obstructive Lung Disease (GOLD) definition, the global prevalence of COPD among individuals aged 30–79 years reached 10.3% in 2019.⁴ Bronchodilators, including long-acting β2 agonists, long-acting anticholinergics, and inhaled corticosteroids, alone or in combination, are common treatments for COPD and have been included in the GOLD guidelines.⁵ These drugs can effectively dilate the bronchial tubes and reduce respiratory symptoms, whereas their side effects remain noteworthy. For example, prolonged treatment with inhaled corticosteroids increases the risk of diabetes, pneumonia, and bone fractures in patients with COPD.^6,7 Traditional Chinese medicine (TCM) treatment is gaining attention as a complementary or alternative therapy.

Airway epithelial cells are integral to the innate immune system and constitute the primary defense against harmful substances and pathogens such as cigarette smoke, particulate matter, and microorganisms.⁸ Chronic exposure to cigarette smoke disrupts the barrier function of the airway epithelium, leading to a series of pathological changes, including airway inflammation, mucus hypersecretion, and small airway obstruction.^9,10 Airway epithelial barrier injury is one of the major causes of an excessive inflammatory response and airway remodeling in COPD. Clinical and experimental studies indicate that cigarette smoke decreases the expression of apical junctional proteins—zonula occludens-1 (ZO-1), E-cadherin (E-cad), and occludin (OCC)—in the lung tissues and airway epithelium of patients with COPD, leading to impaired barrier function.^11,12 Therefore, protecting the airway epithelial barrier and elucidating the underlying pathogenic mechanisms may be beneficial in the treatment of COPD.

TCM has certain advantages in the treatment of COPD, with clear clinical efficacy, fewer side effects, and lower costs.¹³ The Bufei Yishen formula (BYF; ZL.201110117578.1) is a Chinese herbal formula comprising 12 herbs indicated for the treatment of COPD. Clinical trials report that BYF can effectively reduce the frequency and duration of acute exacerbations, slow the rate of decline in lung function, and improve patients’ clinical symptoms, exercise tolerance, and quality of life.^14,15 However, due to its complex composition, elucidating its therapeutic mechanism remains challenging. Through a series of animal and in vitro investigations, five major active ingredients of BYF have been identified and validated: astragaloside IV, 20-S-ginsenoside Rh1, icariin, nobiletin, and paeonol, which constitute ECC-BYF III (patent application number, 201811115372.3).^16,17 Previous studies have shown that ECC-BYF III could improve COPD by regulating mitochondrial dysfunction,¹⁸ attenuating inflammation,¹⁹ reducing oxidative stress²⁰ and improving airway remodeling.²¹ Furthermore, it has also been demonstrated that ECC-BYF III can preserve airway epithelial barrier integrity by inhibiting activation of AHR/EGFR, which provides new perspectives for COPD therapy.²²

Although prior studies have highlighted the therapeutic potential of ECC-BYF III in COPD, its precise molecular mechanism—particularly its effect on the miRNA-mRNA regulatory network in epithelial barrier injury—remains not fully understood. miRNAs as a class of non-coding RNAs can negatively regulate gene expression by interacting with the 3’ untranslated regions (UTRs) of target mRNAs, thereby participating in diverse physiological and pathological processes.²³ Numerous miRNAs display aberrant expression patterns in COPD airways. For example, miR-29b modulates inflammatory cytokine expression in bronchial epithelial cells by targeting BRD4,²⁴ miR-34a regulates expression of anti-aging proteins SIRT1 and SIRT6 in bronchial epithelial cells,²⁵ and miR-146a-5p influences mucus hypersecretion in COPD rat models via inhibiting the EGFR/MEK/ERK pathway.²⁶ Both clinical and experimental evidence have demonstrated that dysregulated miRNAs play an important role in COPD, including cigarette smoke-associated immune and inflammatory responses.^27,28 Beyond analysis of individual miRNAs, several studies have used the negative regulatory relationship between miRNAs and mRNAs to construct COPD-related miRNA-mRNA interaction networks.^29,30 However, few studies have investigated the effects of ECC-BYF III on COPD-associated epithelial barrier dysfunction through modulation of the miRNA-mRNA regulatory network.

To fill the gaps in existing research, this study utilized BEAS-2B cells and cigarette smoke extract (CSE) to establish an epithelial barrier injury model, which is widely used in other studies.^31,32 BEAS-2B cells, as a representative airway epithelial cell line, can well mimic the physiological and pathological features of airway epithelial cells in vivo.³³ CSE can simulate the direct effects of cigarette smoke on airway epithelial cells, including inflammation, oxidative stress, and barrier dysfunction.^32,34 Given that COPD is a highly complex disease, the miRNA-mRNA interaction network plays a crucial role in its pathogenesis.³⁰ Specifically, aberrant regulation of the miRNA-mRNA network may lead to pathological alterations such as airway epithelial cell dysfunction, inflammatory responses, and tissue remodeling.^35,36 Comprehensive analysis of the miRNA-mRNA network may provide critical insights into gene regulatory mechanisms underlying COPD and reveal potential therapeutic targets. Therefore, based on the CSE-induced epithelial barrier injury model, this study intends to investigate the impact of ECC-BYF III on the miRNA-mRNA network and elucidate its novel mechanisms for the treatment of COPD. These findings may advance understanding of the pharmacological effects of ECC-BYF III and provide new molecular targets and theoretical basis for the treatment of COPD.

Materials and Methods

Chemicals and Reagents

Cigarettes were obtained from Henan Tobacco Industry Co., Ltd (Hongqi Canal^® filter-tipped cigarettes, containing tar 10 mg, nicotine 1.0 mg, carbon monoxide 11 mg). Klebsiella pneumoniae (strain 46117) was provided by the National Medical Strain Preservation Center (Beijing, China). ECC-BYF III, comprising 20-S-ginsenoside Rh1, astragaloside IV, icariin, nobiletin, and paeonol, was purchased from Chengdu Must Biotechnology Co., Ltd. The purity of all compounds was determined by high-performance liquid chromatography to be greater than 99%.

Cigarette Smoke Extract Preparation

CSE was prepared by burning a cigarette and channeling the smoke into 2 mL of DMEM medium without fetal bovine serum (FBS) at a rate of 5 min per cigarette. Vigorous shaking dissolved the cigarette smoke into the medium, producing a coffee-colored solution. The optical density (OD) at 320 nm was measured using a UV spectrophotometer (Thermo Fisher, USA). An OD value of 1.8–2.2 indicated a 100% CSE solution, and this value was consistently maintained for each preparation. Subsequently, the solution was filtered through a 0.22 μm filter and diluted with FBS-free medium to the desired concentration. Fresh CSE was used within 30 min of preparation.

Cell Culture and Treatment

Human bronchial airway epithelial cells (BEAS-2B) were obtained from Shanghai ZiShi Biological Co., Ltd (Shanghai, China), and cultured in DMEM (Cat. 11965, Solarbio, Beijing, China) supplemented with 10% FBS (Cat. S711-011S, Lonsera, Shanghai, China). BEAS-2B cells were then divided into the following groups: control, model, and ECC-BYF III (70, 35, and 17.5 μg/mL).^21,22,37 The cells were seeded in 6-well plates for 12 h, followed by cultivation in FBS-free medium for 3 h. Based on previous studies, ECC-BYF III groups were first treated with ECC-BYF III for 3 h. Both model and ECC-BYF III groups were then exposed to 10% CSE for 24 h to induce epithelial barrier injury.^11,22

RNA-Seq Data of BEAS-2B Cells

For the samples obtained from the control, model, and ECC-BYF III groups (n=3), total RNA was isolated from BEAS-2B cells using the Trizol reagent (Cat. 15596026CN, Invitrogen, Carlsbad, CA, USA). Qualified total RNA samples were selected for RNA sequencing by Jingzhou Gene Technology Co., Ltd. Sequencing data were obtained for both miRNA and mRNA.

Identification of Differentially Expressed mRNAs and miRNAs

High-throughput RNA sequencing data were analyzed using the “edgeR” package within the R software (v4.3.1) to identify differentially expressed mRNAs (DEmRNAs) and miRNAs (DEmiRNAs), applying a significance threshold of P<0.05. Volcano plots were generated using the “ggplot2” package, while clustering analysis and heatmaps were created using the “pheatmap” package.

DEmiRNAs Target Gene Prediction

Target genes of DEmiRNAs were predicted using four online databases: DIANA-microT-CDS (http://www.microrna.gr/microT-CDS), miRWalk (http://mirwalk.umm.uni-heidelberg.de), miRDB (http://www.mirdb.org), and TargetScan v7.2 (www.targetscan.org). Since each database employs distinct algorithms for target prediction and ranks interactions differently, only miRNA-mRNA interactions identified in at least three databases were considered. An mRNA was defined as the target gene when a miRNA was found to target it in at least three databases.^38,39

Enrichment Analysis for DEmRNAs and Target Genes

Gene Ontology (GO) (including biological process (BP), cellular component (CC), and molecular function (MF)) and KEGG pathway enrichment analyses were performed for both DEmRNAs and predicted target genes, using the “ClusterProfiler” package and a hypergeometric distribution model. Pathways with P<0.05 were considered significantly enriched. The “GOCircle” package was used for circle diagrams, and the “ggplot2” package was used for bar plots.

The Screening of miRNA-mRNA Regulatory Pairs

Predicted target genes were intersected with DEmRNAs to identify mRNAs that are both differentially expressed and regulated by miRNAs. Based on the negative regulatory relationship between miRNAs and mRNAs, upregulated (downregulated) miRNAs were paired with downregulated (upregulated) mRNAs, to screen miRNA-mRNA pairs. Subsequently, miRNA-mRNA pairs whose dysregulation directions reversed after ECC-BYF III intervention were identified based on differential expression analysis (ECC-BYF III group vs model group). Moreover, Pearson’s correlation analysis was performed to explore the correlations between miRNAs and mRNAs.

COPD Rat Model Construction

Eighteen Sprague-Dawley rats (200 ± 20 g) were obtained from Beijing Vitonglihua Laboratory Animal Co., Ltd (Beijing, China) and randomly divided into control, model, and ECC-BYF III groups (n=6). From weeks 1 to 8, the COPD model was established by cigarette smoke exposure (3000 ± 500 ppm for 40 min, twice daily) combined with repeated Klebsiella pneumoniae infection (6×10⁸ CFU/mL, 0.1 mL nasal instillation every 5 days).⁴⁰ From weeks 9 to 16, control and model group rats received 0.5% CMC-Na by gavage, whereas ECC-BYF III group rats received ECC-BYF III (5.5 mg/kg/d) administration. The dosage of ECC-BYF III was determined based on our previous studies.^19,21 Lung tissue samples were collected at week 16 and subsequently used for qRT-PCR experiments. The animal experiments were approved by the Ethics Committee of Laboratory Animal Welfare of Henan University of Chinese Medicine (DWLLGZR202303082).

Quantitative Real-Time PCR Assay

Total RNA was extracted from lung tissue and BEAS-2B cells using Trizol reagent and assessed for quality and purity. The RNA was reverse transcribed to cDNA using HiScript II Q RT SuperMix (Cat. R223, Vazyme, Nanjing, China), and gene expression was quantified by qRT-PCR. The primers for mRNA and miRNA were synthesized by Azenta Life Sciences and are listed in Tables 1–3. GAPDH and RNU6B (U6) served as internal reference genes, and the 2^−ΔΔCt method was used to evaluate gene expression.

Table 1 miRNA Primer Sequences for qRT-PCR Experiments in BEAS-2B Cells

Table 2 mRNA Primer Sequences for qRT-PCR Experiments in BEAS-2B Cells

Table 3 mRNA Primer Sequences for qRT-PCR Experiments in Rats

Independent Dataset Validation

The scRNA-seq data GSE173896 (containing 5 COPD and 5 control lung tissue samples) and the bulk RNA-seq data GSE11906 (containing airway epithelium obtained by fiberoptic bronchoscopy from 33 COPD patients and 35 healthy non-smokers) were downloaded from the GEO database. For the GSE173896 dataset, the “Seurat” package was utilized for quality control and normalization processes. The top 2000 highly variable genes were selected for principal component analysis (PCA) using the FindVariableFeatures function. Cellular subtypes were annotated using the CellMarker database and the “singleR” package. The FindMarkers function was employed to identify differentially expressed genes between COPD and control groups within the epithelial cell population. For the GSE11906 dataset, the “limma” package was used to verify the reliability of hub genes between COPD and control groups. The “ggplot2” package was utilized to visualize the expression of hub genes.

Single-Sample Gene Set Enrichment Analysis (ssGSEA) Analysis

The ssGSEA analysis was performed using the “GSVA” package to evaluate the infiltration levels of 28 immune cell types. Spearman correlation analysis was performed to explore the associations between hub genes and the abundances of infiltrating immune cells. Statistical significance was set at P<0.05.

Molecular Docking Analysis

The two-dimensional structures of the five components of ECC-BYF III were downloaded from the PubChem database (https://pubchem.ncbi.nlm.nih.gov). The energies of the compounds were minimized using Chem3D v15.1 software, and three-dimensional (3D) structures were saved in mol2 format. The 3D structures of the hub genes (protein receptors) were obtained from the PDB (http://www.rcsb.org) or Uniprot (http://www.uniprot.org) databases. The protein receptors were processed to remove water molecules and other small molecules utilizing PyMOL software, followed by hydrogenation and charge adjustment with AutoDockTools v1.5.6 software. Molecular docking was performed using AutoDock Vina, and the results were visualized using PyMOL. A larger absolute value of the docking score indicated a more stable interaction between the compound and the protein receptor.

Statistical Analysis

Statistical analyses were performed using R (v4.3.1) or SPSS v26.0 (IBM Corporation, Armonk, NY, USA). The qRT-PCR data are presented as mean ± standard deviation (SD). Statistical significance of the qRT-PCR results for the miRNAs and genes was evaluated using one-way ANOVA. P<0.05 was considered statistically significant. Differential expression analyses of mRNAs and miRNAs were conducted using the “edgeR” package. GO functional enrichment analyses for DEmRNAs and target genes were performed using the “ClusterProfiler” package.

Results

Identification of Differentially Expressed mRNAs

A total of 2997 DEmRNAs were identified between the model and control groups, with 1483 upregulated and 1514 downregulated (P<0.05), as shown in Figure 1A. Furthermore, a heatmap was constructed to display these DEmRNAs (Figure 1B).

Figure 1 Volcano plot and heatmap of DEmRNAs between the model and control groups. (A) Volcano plot of 2997 DEmRNAs. (B) Heatmap of DEmRNAs. Red indicates upregulated genes, and blue indicates downregulated genes.

Identification of DEmiRNAs and Prediction of Target Genes

Compared with the control group, four DEmiRNAs were identified in the model group, of which three were downregulated and one was upregulated (P<0.05), as shown in Figure 2 and Table 4. Moreover, based on four databases (DIANA-microT-CDS, miRWalk, miRDB, and TargetScan v7.2), a total of 2430 genes were predicted as targets of these DEmiRNAs in at least three databases (Table 5).

Table 4 Four DEmiRNAs Between the Model and Control Groups

Table 5 The Number of Target Genes of DEmiRNAs

Figure 2 Volcano plot of four DEmiRNAs between the model and control groups.

Pathway Enrichment Analysis

Enrichment analysis of the 2997 DEmRNAs identified 96 significantly enriched KEGG pathways, including NF-κB, AMPK, TNF, and IL-17 signaling pathways (Figure 3A). For the 2430 target genes of the DEmiRNAs, 112 pathways, such as MAPK, FoxO, p53, and Wnt signaling pathways, were enriched (Figure 3C). Specifically, 53 pathways showed significant overlap based on the hypergeometric distribution model (P<0.05). Among these overlapping pathways, MAPK, PI3K-Akt, AMPK, and TGF-beta signaling pathways have been reported to be associated with COPD.

Figure 3 GO and KEGG enrichment analyses of DEmRNAs and target genes of DEmiRNAs. (A) KEGG enrichment analysis of DEmRNAs. The more yellow the color, the higher the significance of enrichment. (B) GO enrichment analysis of DEmRNAs. The bar height in the inner circle indicates the importance of GO terms, with colors corresponding to the Z-scores. The outer circle scatter plot shows the expression of genes (logFC) in each term. (C) KEGG and (D) GO enrichment analyses of target genes of DEmiRNAs.

GO analysis identified 651 significantly enriched GO terms based on DEmRNAs: 577 for BP, 65 for CC, and 9 for MF. The top 5 enrichment terms in each category are shown in the circle diagram in Figure 3B. In BP, the genes were mainly enriched in epithelial cell proliferation, tissue remodeling, positive regulation of epithelial to mesenchymal transition, and response to oxidative stress. CC terms included cell-substrate junction and protein complex involved in cell adhesion. MF terms included extracellular matrix binding, growth factor receptor binding, and fibronectin binding. Similarly, for the DEmiRNA target genes, 504 GO terms were enriched (397 BP, 61 CC, and 46 MF). As shown in Figure 3D, BP terms were primarily enriched in Wnt signaling, epithelial cell development, and epithelial cell proliferation; CC terms included adherens junction and focal adhesion; and MF terms included Wnt receptor activity, cAMP response element binding, and MAP kinase activity.

Construction of the miRNA-mRNA Regulatory Network

Among the 2430 predicted target genes and 2997 DEmRNAs, 411 overlapped genes were obtained, of which 131 were upregulated and 280 were downregulated. According to the negative regulatory relationship between miRNAs and mRNAs, 149 miRNA-mRNA pairs, comprising 4 miRNAs and 139 mRNAs, were identified (Figure 4A and B). The miRNA-mRNA regulatory network was visualized with Cytoscape software (Figure 4C).

Figure 4 Construction of the miRNA-mRNA regulatory network. (A) Venn diagram of target genes of downregulated miRNAs overlapping with upregulated DEmRNAs. (B) Venn diagram of target genes of upregulated miRNAs overlapping with downregulated DEmRNAs. (C) Regulatory network of miRNA-mRNA interactions. Red indicates upregulation; blue indicates downregulation. Triangles represent miRNAs; the symbol “V” represents mRNAs.

Identification of miRNAs and mRNAs Regulated by ECC-BYF III

For the 149 identified miRNA-mRNA pairs (comprising 4 miRNAs and 139 mRNAs), differential expression analysis (ECC-BYF III group vs model group) revealed that 41 mRNAs were significantly reversed after ECC-BYF III intervention (P<0.05). Similarly, the expression patterns of 3 DEmiRNAs were reversed following treatment. This resulted in 42 miRNA-mRNA pairs, consisting of 3 miRNAs and 41 mRNAs. Notably, 12 of these mRNAs—RNF44, YBX1, PGAM5, RELT, RFX1, NDRG1, ESM1, CTTN, BCL2L1, MAP2K7, ADAM19, and SMYD5—have been reported to be closely associated with COPD pathological processes or inflammatory responses.^41–48 Based on these findings, these 12 mRNAs were designated as candidate target genes, and they, together with their corresponding 3 miRNAs, formed 13 miRNA-mRNA pairs for subsequent qRT-PCR validation (Table 6).

Table 6 The Target Genes of the miRNAs Regulated by ECC-BYF III

Pearson’s correlation analysis was conducted for these miRNA-mRNA pairs. As shown in Figure 5, the expression of all 13 pairs was negatively correlated, consistent with the expected regulatory relationship. A significant negative correlation was observed between hsa-miR-3658 and several mRNAs, namely RNF44, YBX1, PGAM5, RELT, and RFX1. Moreover, hsa-miR-3685 expression was significantly negatively associated with ESM1, while SMYD5 expression showed a significant inverse correlation with hsa-miR-3936.

Figure 5 Pearson’s correlation analysis of the identified miRNA-mRNA pairs.

Validation of the Expression of miRNA and mRNA in BEAS-2B Cells

qRT-PCR results of BEAS-2B cells confirmed that the expression trends of hsa-miR-3685 and hsa-miR-3936 were consistent with the data analysis results (Figure 6). The expression of ESM1, the target gene of hsa-miR-3685, was significantly upregulated in the model group, and markedly decreased after ECC-BYF III treatment. The expression levels of the target genes (RNF44, BCL2L1, MAP2K7, ADAM19, and SMYD5) of hsa-miR-3936 were also elevated compared with controls, and then significantly decreased after ECC-BYF III intervention. Ultimately, two miRNAs (hsa-miR-3685, hsa-miR-3936) and six target genes (ESM1, RNF44, BCL2L1, MAP2K7, ADAM19, and SMYD5) associated with airway epithelial barrier injury were identified. Subsequently, we converted the obtained human miRNAs and mRNAs into their homologous rat RNA molecules using the bioDBnet database. Six rat-derived gene homologs were obtained, whereas no rat-derived miRNA homologs were identified. Thus, we conducted qRT-PCR experiments to validate these six hub genes in the lung tissues of COPD rat models.

Figure 6 qRT-PCR validation of identified miRNAs and mRNAs in BEAS-2B cells. The miRNA expression levels of hsa-miR-3685 (A) and hsa-miR-3936 (B). The mRNA expression levels of ESM1 (C), RNF44 (D), BCL2L1 (E), MAP2K7 (F), ADAM19 (G), and SMYD5 (H). Comparison with the control group: *P<0.05, **P<0.01. Comparison with the model group: #P<0.05, ##P<0.01.

Validation of the Expression in Rat Lung Tissues

qRT-PCR results showed that the expression levels of five hub genes were upregulated in lung tissues from COPD rats compared with controls, with ESM1, BCL2L1, ADAM19, and SMYD5 reaching statistical significance (Figure 7). Following ECC-BYF III treatment, the expression levels of all five genes were significantly reduced, consistent with the data analysis results. We found that the dysregulation of ESM1, RNF44, BCL2L1, ADAM19, and SMYD5 was reversed after ECC-BYF III treatment in both rat lung tissues and BEAS-2B cells, which deserves further study.

Figure 7 qRT-PCR validation results in lung tissues from COPD rat models. The mRNA expression levels of ESM1 (A), RNF44 (B), BCL2L1 (C), ADAM19 (D), and SMYD5 (E). Comparison with the control group: *P<0.05, **P<0.01. Comparison with the model group: ##P<0.01.

Validation of Hub Genes Using Independent Datasets

After quality control and normalization, the GSE173896 dataset yielded 38,735 cells and 26,685 genes. Subsequently, 25 clusters were identified through dimensionality reduction and clustering analysis, which were visualized by UMAP (Figure 8A). Based on marker genes (Figure 8C), nine cell types were classified, including T cells, epithelial cells, endothelial cells, macrophages, B cells, NK cells, monocytes, smooth muscle cells, and neutrophils, as shown in Figure 8B. Among epithelial cells, four hub genes—RNF44, BCL2L1, ADAM19, and SMYD5—were detected. The dysregulation direction of BCL2L1, ADAM19, and SMYD5 matched expectations, with SMYD5 showing statistical significance (Figure 8D).

Figure 8 Validation of the five hub genes in dataset GSE173896. (A) UMAP visualization of 25 clusters. (B) Visualization of nine identified cell types. (C) Marker genes for nine cell types. (D) Violin plots showing expression of three hub genes in epithelial cells.

Analysis of the independent dataset GSE11906 showed upregulation of all five hub genes (ESM1, RNF44, BCL2L1, ADAM19, and SMYD5) in the COPD group compared with controls, consistent with the qRT-PCR results. Notably, RNF44 and ADAM19 exhibited statistical significance, further demonstrating the reliability of the identified hub genes (Figure 9).

Figure 9 Validation of the five hub genes in dataset GSE11906. The expression levels of ESM1 (A), RNF44 (B), BCL2L1 (C), ADAM19 (D), and SMYD5 (E) between the COPD group and controls.

Characterization of Immune Cell Infiltration

The infiltration of 28 immune cell types was evaluated using the ssGSEA algorithm. Ten immune cell types displayed significant differences between the model and control groups. As shown in Figure 10A, central memory CD8 T cells, T follicular helper cells, CD56dim natural killer (NK) cells, and eosinophils were significantly increased in the model group, whereas activated CD4 T cells, gamma delta T cells, type 2 T helper (Th2) cells, regulatory T cells, NK cells, and macrophages were significantly decreased. ECC-BYF III intervention reversed the infiltration patterns of eight immune cell types, and the changes of activated CD4 T cells and Th2 cells were statistically significant (Figure 10B).

Figure 10 Immune cell infiltration analysis based on BEAS-2B cells sequencing data. (A) Box plots showing the proportion of 28 immune cell types in the model and control groups. (B) Box plots showing the immune cell types infiltration differences between the ECC-BYF III and model groups. *P<0.05.

In addition, correlation analysis between immune cells and hub genes revealed that RNF44 and BCL2L1 were positively correlated with CD56dim NK cells and eosinophils, and negatively correlated with activated CD4 T cells and neutrophils. ADAM19, ESM1, and SMYD5 showed significant positive correlations with CD56dim NK cells, and negative correlations with regulatory T cells, NK cells, and gamma delta T cells (Figure 11A). Moreover, correlation analysis among immune cell types indicated interactions between certain cell populations (Figure 11B). For instance, effector memory CD8 T cells were positively associated with immature B cells and neutrophils, but inversely correlated with activated B cells and plasmacytoid dendritic cells. Activated CD4 T cells were negatively associated with immature dendritic cells, while positively associated with monocytes.

Figure 11 Results of immune cell correlation analysis in BEAS-2B cells sequencing data. (A) Correlation analysis between hub genes and immune infiltrating cells. (B) Correlation analysis among the 28 immune cell types. *P<0.05, **P<0.01, ***P<0.001.

Molecular Docking Analysis

To further validate the intervention mechanism of ECC-BYF III in airway epithelial barrier injury, molecular docking was conducted to explore potential interactions between the five identified hub genes and the components of ECC-BYF III. The results showed that nearly all genes exhibited strong binding affinities with astragaloside IV, 20-S-ginsenoside Rh1, icariin, nobiletin, and paeonol, as presented in Table 7 and Figure 12.

Table 7 Binding Energy Between the Five Hub Genes and the Five Components of ECC-BYF III

Figure 12 Molecular docking results of the five hub genes with the five components of ECC-BYF III.

Discussion

COPD is a complex disorder that poses a major threat to human health. Airway epithelial cells serve as the first line of defense against external stimuli and play a critical role in lung diseases including COPD. Various environmental factors, such as tobacco smoke, particulate matter, ozone, and vehicular emissions, can directly or indirectly impair epithelial barrier function. Dysfunction of epithelial barrier is closely associated with the development of COPD. For example, CSE downregulates E-cad protein expression, thereby damaging epithelial barrier integrity in patients with COPD.⁴⁹ Cigarette smoke also increases bronchial epithelial permeability and promotes the release of inflammatory factors, contributing to sustained airway inflammation in COPD.⁵⁰ Persistent epithelial barrier injury and aberrant repair drive chronic airway inflammation and airway remodeling in COPD. Therefore, preserving epithelial integrity and barrier function is crucial to prevent or delay disease progression.

TCM offers an alternative therapeutic strategy for COPD, with proven efficacy and few adverse effects. BYF has been shown to be effective in the treatment of COPD in both clinical and experimental studies. Multicenter, randomized, double-blind trials have demonstrated that BYF could improve lung function, reduce the frequency and severity of acute exacerbations, and alleviate symptoms in patients with COPD.¹⁵ However, due to the complex composition of BYF, elucidating its precise mechanism of action has been challenging. Based on the principle of active ingredient compatibility, ten active components of BYF (ECC-BYF I) were screened in various cell models, with their bioactivities further validated in COPD rat models.¹⁷ The ECC-BYF II formulation was obtained by comprehensively evaluating the effects of each component of ECC-BYF I on lung function, lung histopathology, Inflammatory responses, and oxidative stress based on R-value analysis.¹⁶ Further refinement yielded five active ingredients, namely 20-S-ginsenoside Rh1, astragaloside IV, icariin, nobiletin, and paeonol, which were combined in fixed proportions to prepare the ECC-BYF III formulation. Qin et al reported that ECC-BYF III could suppress inflammatory response in COPD by inhibiting NLRP3 inflammasome activation through modulation of the GSK3β pathway.¹⁹ Moreover, ECC-BYF III could also maintain airway epithelial barrier integrity by suppressing AHR/EGFR activation.²² However, its regulatory mechanism on the miRNA-mRNA network has not been fully clarified.

In this study, expression levels of two miRNAs were markedly decreased in the model group, accompanied by significant upregulation of five target genes (ESM1, RNF44, BCL2L1, ADAM19, and SMYD5). Notably, ECC-BYF III treatment reversed these aberrant expression patterns. Several of these regulatory genes were reported to be related to COPD pathogenesis. Recent study has highlighted the protective role of ESM1 in acute pulmonary inflammation, demonstrating its potential value as an inflammatory biomarker.⁵¹ Another study found that serum ESM1 levels correlate with pulmonary function and asthma severity.⁵² Consistent with our findings, previous studies have reported elevated levels of serum and plasma ESM1 in COPD patients.^48,53 RNF44, an E3 ubiquitin ligase family member, mediates ubiquitination of multiple protein substrates.⁵⁴ Given the close association between the ubiquitination process and COPD, investigating RNF44 may reveal novel therapeutic avenues for COPD.^55,56 Elevated expression of BCL2L1 in alveolar macrophages of smokers has been associated with decreased apoptosis,⁵⁷ and alterations in BCL2L1 expression levels are correlated with the delayed apoptosis of neutrophils in COPD patients.⁴⁵ Studies have demonstrated that ADAM19 is associated with lung function,⁵⁸ airflow obstruction,⁵⁹ and genetic susceptibility to COPD.^47,60 As a negative regulator of inflammatory gene expression in macrophages, SMYD5 may influence the COPD-related inflammatory process.⁴⁶ The specific functions of hsa-miR-3685 and hsa-miR-3936 in COPD remain largely unknown, necessitating further molecular investigations to clarify their underlying mechanisms. Our findings suggest that ECC-BYF III may alleviate epithelial barrier injury in COPD by regulating these network targets and thereby modulating inflammation, apoptosis, and other biological processes.

KEGG pathway enrichment analysis of DEmRNAs and target genes of DEmiRNAs revealed 53 significantly overlapping pathways (hypergeometric distribution model, P<0.05). This high degree of overlap suggests that miRNAs and mRNAs are involved in the same or interconnected biological processes. Several of these pathways are known to regulate inflammation and epithelial barrier function in COPD. For example, reactive oxygen species can induce hypermethylation of the Rab26 promoter through MAPK pathway activation, thereby exacerbating airway epithelial inflammation.⁶¹ The MAPK pathway also regulates the expression of tight junction proteins (ZO-1, Claudin-1, OCC) and alters epithelial permeability, compromising barrier integrity.^62,63 Specifically, the activation of the PI3K/Akt pathway promotes oxidative stress and the release of pro-inflammatory cytokines.⁶⁴ Song et al⁶⁵ reported that PI3K signaling downregulates epithelial barrier proteins under particulate matter-induced airway inflammation, further impairing epithelial integrity. A previous study found that elevated TGF-β levels were associated with decreased lung function and may modulate the inflammatory response by affecting immune cell recruitment.⁶⁶ Furthermore, the TGF-β pathway plays a key role in airway remodeling during COPD progression.^67,68

Given the crucial role of immune cell infiltration in COPD progression and recovery, we employed the ssGSEA algorithm to quantify immune cell infiltration levels. The analysis revealed significant differences in ten immune cell subsets between the model and control groups. Following ECC-BYF III treatment, infiltration patterns of eight immune cell subtypes were altered, with activated CD4 T cells and Th2 cells showing statistically significant changes. Most of these immune cell populations have been previously reported to be associated with COPD pathogenesis. For example, activated CD4 T cells, particularly the Th1 and Th17 subsets, are key contributors to the chronic inflammatory response in COPD.⁶⁹ These cells secrete pro-inflammatory cytokines that exacerbate the inflammatory cascade in the lungs.⁷⁰ In addition, activated CD4 T cells release cytokines that promote fibroblasts’ proliferation and activation, contributing to airway remodeling in COPD.⁷¹ Th2 cells are known to secrete cytokines including interleukin (IL)-4, IL-5, and IL-13, which drive type 2 inflammation.⁷² This inflammatory response constitutes a fundamental pathophysiological process in chronic airway diseases such as asthma, COPD, and allergic rhinitis.⁷³ Furthermore, Th2 cells can exacerbate airway inflammation and airway hyperresponsiveness in COPD by recruiting and activating eosinophils through cytokine secretion.⁷⁴ Our results also showed that BCL2L1 was positively correlated with CD56dim NK cells and negatively correlated with Th2 cells. Previous studies have demonstrated that BCL2L1 promotes the survival of CD56dim NK cells by exerting anti-apoptotic effects,⁷⁵ whereas spontaneous apoptosis of Th2 cells is associated with downregulation of BCL2L1.⁷⁶ ESM1 exhibited a positive correlation with T follicular helper cells and a negative correlation with activated CD4 T cells. Moreover, ESM1 has been reported to inhibit the cytotoxic function of tumor-infiltrating CD4+ T cells and mediate immunosuppression in T lymphocytes.⁷⁷ Collectively, these findings suggest that the identified hub genes may regulate immune responses in COPD, which could influence disease progression and therapeutic outcomes.

This study has several limitations. First, COPD is a multifactorial disease, and epithelial barrier injury represents only one aspect of its complex pathogenesis. Future studies should integrate multiple models, such as inflammatory or oxidative stress models, to achieve a more comprehensive understanding of disease mechanisms. Second, the sample size used in this study was relatively small. To improve data reliability, rigorous preprocessing steps were applied, including the exclusion of lowly expressed genes, to enhance the accuracy of the high-throughput analysis. Moreover, compared with tissue samples, cellular samples exhibit lower heterogeneity, and the use of three biological replicates per group is a widely accepted strategy for high-throughput analysis.^78,79 Importantly, our findings were well validated in independent datasets, BEAS-2B cells experiments, COPD rat models and molecular docking, which further prove the reliability of our results. Nonetheless, in terms of clinical applicability, results derived from small sample groups may not be applicable to larger populations of COPD patients with diverse demographic characteristics, disease severities, and comorbidities. Therefore, expanding the sample size and incorporating clinical samples will be prioritized in future studies to enhance the applicability and translational relevance of our findings. Third, a relatively lenient threshold (P<0.05) was applied for high-throughput analysis, potentially increasing the risk of false positives. However, the reliability of our findings was robustly validated across multiple models. Moreover, differential expression analysis was performed with a strict threshold (FDR<0.05), identifying 1010 DEmRNAs but no DEmiRNAs between the model and control groups, precluding subsequent miRNA-mRNA network analysis. For the 1010 DEmRNAs, following ECC-BYF III intervention and target gene prediction analysis, only 11 reversed target DEmRNAs were obtained. In contrast, under the P<0.05 threshold, 2997 DEmRNAs were identified and 41 reversed target DEmRNAs were obtained. Notably, all 11 reversed target DEmRNAs obtained under FDR<0.05 were included within the 41 reversed DEmRNAs identified with P<0.05, indicating that P<0.05 possesses greater statistical power. Similarly, KEGG pathway enrichment analysis identified 25 significantly disturbed pathways under FDR<0.05, based on the 1010 DEmRNAs. Among these 25 pathways, 23 were replicated at P<0.05 (96 pathways were identified), further demonstrating the statistical power and reliability under the threshold of P<0.05. Therefore, the threshold of P<0.05 for high-throughput analysis was considered appropriate in this study. Overall, the identified miRNA-mRNA regulatory network provides valuable insights into the molecular mechanisms underlying epithelial barrier injury in COPD, which might identify new therapeutic strategies. Targeting the miRNA-mRNA network in combination with traditional therapeutic approaches may facilitate a more comprehensive and effective management of COPD.

In summary, based on RNA-seq data of BEAS-2B cells, miRNA databases, and multiple bioinformatics algorithms, we discovered that ECC-BYF III might improve airway epithelial barrier injury by regulating the hsa-miR-3685-ESM1 and hsa-miR-3936-RNF44/BCL2L1/ADAM19/SMYD5 networks. And our findings were well validated in independent datasets (single-cell and bulk RNA-seq data), BEAS-2B cells experiments, COPD rat models, and molecular docking, confirming the reliability of the results. This study provides a foundation for further mechanistic investigation and research into the prevention and treatment of COPD using TCM.

Conclusions

Our findings suggest that ECC-BYF III might improve airway epithelial barrier injury in COPD by regulating the hsa-miR-3685-ESM1 and hsa-miR-3936-RNF44/BCL2L1/ADAM19/SMYD5 networks (via downregulating miRNAs and upregulating mRNAs). However, the precise mechanisms of these regulatory networks and the related pathways of ECC-BYF III require further validation through clinical studies. Moreover, this study provides new ideas for exploring the pharmacological mechanisms of other diseases treated with TCM.