The identification of DEGs in COPD

Through analysis with the GSE76925, 248 DEGs, including 204 down-regulated and 44 up-regulated genes, were found between the COPD sample and health control (adjusted P value ≤ 0.05 & |Log2FC|≥ 1, supplementary Table 1). The distribution and expression of TOP 10 up-regulated genes (CA3, FGG, HTR2B, CD19, DPYS, TNFRSF17, TM4SF19, COL10A1, COMP, CHIT1) and TOP 10 down-regulated genes (KLF5, SUMO1, ADI1, DNTTIP2, B3GNT5, ID4, PPIB, CD164, HMCN1, MED21) according to log2 (FC) were exhibited in Fig. 1A-B. Then, GO enrichment analysis of DEGs was performed to understand their function in COPD. As a result, 449 enriched GO terms were identified: 56 in MF, 25 in CC, and 368 in BP. The most significant terms (with the minimum P value) included leukocyte migration, protein processing, peptide hormone processing, phosphatidylinositol-5-phosphate binding, regulation of microtubule nucleation, signaling receptor ligand precursor processing, de novo protein folding, heat shock protein binding, collagen-containing extracellular matrix, and protein stabilization. Genes such as PPIB, LEP, BCHE, PLA2G7, SH3PXD2B, DNAJB4, DNAJC2, and SRPX2 are involved in these pathways. Moreover, 8 KEGG pathways were enriched, including Primary immunodeficiency, Basal transcription factors, B cell receptor signaling pathway, Cytokine-cytokine receptor interaction, Arachidonic acid metabolism, inflammatory mediator regulation of TRP channels, Ether lipid metabolism, and, alpha-Linolenic acid metabolism. Genes such as PLA2G4A, PLA2G2D, TNFRSF13B, CD79A, CD19, IL1RAP are involved in these pathways (Fig. 1C-D, Table 2, supplementary Table 2 and 3).

Fig. 1

Differential expression gene screening. A The volcano plot revealed DEGs between COPD and controls from the GSE76925 dataset (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE76925). Red plots represented upregulated DEGs, blue plots represented downregulated DEGs, and gray spots represented undifferentiated genes. B The expression of the TOP 20 differential expression genes. Red represented COPD samples, blue represented control samples. C annotation of differential expressed genes. Column colors represent different GO types, and lengths represent the number of enriched genes in term. D Chord diagrams of the KEGG pathway. The color of the left gene ribbon represents the logFC of the gene, and different ribbons on the right represent different pathways

Table 2 Top 10 GO enrichment and Top 5 KEGG enrichment of 248 DEGs in GSE76925The identification of DE-ORGs

Then, the intersection of DEGs and ORGs was performed, and a VENN plot was constructed to obtain DE-ORGs. 8 DE-ORGs including CA3, PPP1R15B, MAPT, MMP9, ECT2, PDE8A, HSPA1A, HIF1A, was obtained (Fig. 2A, supplementary Table 4). PPP1R15B was found on Chromosome 1, ECT2 on Chromosome 3, HSPA1A and CA3 on Chromosome 6, HIF1A on Chromosome 14, PDE8A on Chromosome 15, MAPT on Chromosome 17, and MMP9 on Chromosome 20, and no gene was distributed in XY Sex chromosomes (Fig. 2B). Subsequently, GeneMANIA was utilized to predict the hub-gene-associated genes and their functions. There were 20 genes with a total of 157 predicted links (Co-expression of 12.58%, Shared protein domains accounted for 10.16%, Physical Interactions accounted for 74.96%, and Genetic Interactions of 2.31%). These related genes (especially HSPA1A, MMP9, MAPT) were enriched in the regulation of pathways such as mitochondrion organization, apoptotic signaling pathway, mitochondrial outer membrane permeabilization, and apoptotic mitochondrial changes (Fig. 2C).

Fig. 2

Identification of differential expressed oxidative stress-related genes. A The intersection of DEGs and ORGs was performed, and a VENN plot was constructed to obtain differential expressed oxidative stress-related genes (DE-ORGs). B The chromosome localization of DE-ORG. C GeneMANIA Network. The middle is the key gene, the outer circle is the related genes with similar functions to the key gene, different colors connecting lines indicate different networks, and different colors in the pie chart indicate different functions of the genes

Screening of hub genes using machine learning algorithms

Feature genes were further selected using Boruta, LASSO regression, and SVM-RFE based on DE-ORGs in the training set (GSE76925), considering COPD as the outcome variable. First, the LASSO regression was applied. Eight DE-ORGs, including CA3, PPP1R15B, MAPT, MMP9, ECT2, PDE8A, HSPA1A, and HIF1A, were chosen because their regression coefficients remained above zero with a minimal penalty term (λ) of 0.000317374 (Fig. 3A-B, supplementary Table 5). SVM-RFE analysis was then used. The model prediction accuracy was highest for the first time when the number of characterized genes was 5. Additionally, the comparable genes were MMP9, ECT2, PPP1R15B, CA3, and MAPT (Fig. 3C, supplementary Table 6). Boruta analysis was performed on 8 DE-ORGs. The importance of seven genes (CA3, PPP1R15B, MAPT, MMP9, ECT2, HSPA1A, and HIF1A) was found to be Confirmed. The importance of PDE8A was defined as Tentative and therefore excluded (Fig. 3D, supplementary Table 7). Additionally, 8 LASSO-characterized genes, 5 SVM-RFE-characterized genes, and 7 Boruta-characterized genes were intersected to produce the hub genes, which were then further characterized. Thus, 5 hub genes (CA3, PPP1R15B, MAPT, MMP9, and ECT2) were identified (Fig. 3E, supplementary Table 8).

Fig. 3

Machine learning algorithms screen for hub genes. A Determination of the number of factors by the LASSO regression analysis. The X-axis was the logarithm of lambdas, and Y-axis was the coefficient of the variable. As lambdas increased, the coefficients of the variables converged to 0. When the optimal lambda was reached, variables with coefficients equal to 0 were excluded. B Ten-fold cross-validation of tuning parameters in LASSO analysis. The X-axis was logarithm of lambdas, The Y-axis was the partial likelihood Deviance. C SVM-RFE analysis. The horizontal coordinate indicated the number of characterized genes, and the vertical coordinate indicated the model prediction accuracy. D Boruta analysis was performed. The horizontal coordinate indicated the number of characterized genes, and the vertical coordinate indicated the importance of the characterized genes. E Venn diagram of characterized genes

Construction and validation of the diagnostic model

Next, a nomogram of five hub genes was created to forecast the likelihood of COPD development. Each gene received an individual score, with total points computed by aggregating all scores. The risk of developing COPD was predicted using total points; the greater the total points, the higher the probability of developing COPD and the higher the likelihood of developing COPD (Fig. 4A). The nomogram's prognostic value was assessed by the construction of the calibration curve. The outcome demonstrated that there was no discernible difference between the reference line and the nomogram's estimated probability (P = 0.0763, Fig. 4B). Additionally, the area under the curve was 0.8459459, indicating that the created nomogram had a high degree of prediction accuracy (Fig. 4C). DCA demonstrated a notably improved benefit rate of the constructed nomogram over the diagonal line (All) and the horizontal line (None), highlighting its additional advantages (Fig. 4D).

Fig. 4

Construction and validation of diagnostic model. A Construction of a nomogram of COPD risk based on hub gene expression in the training Set. B Calibration curves for nomogram. The horizontal coordinate was the predicted event rate, and the vertical coordinate was the actual observed event rate, both ranging from 0 to 1. The solid blue line represented the entire cohort, and the solid black line, bias-corrected by bootstrapping (1000 repetitions), represented the observed nomogram performance. C ROC curve. The horizontal coordinate was the rate of false positives (1-specificity), and the vertical coordinate was the rate of true positives (sensitivity or sensitivities). The AUC value was the area under the ROC curve. D DCA Curve. The diagonal line (All) represented all samples with all interventions; the horizontal line (None) represented all samples with no intervention

Relationship and regulation of hub genes

Spearman correlation analysis examined the relationships among the hub genes. ECT2 and PPP1R15B displayed a strong positive correlation (cor = 0.76, P < 0.05), whereas CA3 showed a notable negative correlation with ECT2 (cor = −0.41, P < 0.05). Additionally, MMP9 did not correlate lower with other genes (Fig. 5A). The chord plot (Fig. 5B) was then selected based on the results of |cor|> 0.3 and P < 0.05. Functional similarity analysis based on GO-MF, GO-BP, and GO-CC dimensions was carried out to acquire similarity scores between genes to learn about the functional similarities of the hub gene. MAPT, ECT2, and CA3 were the top three genes with the greatest functional similarity (Fig. 5C). To investigate whether oxidative stress indicators and the hub genes (including 8 biomarkers and 5 antioxidant biomarkers) are related, relevance scores between the hub gene and oxidative stress biomarkers were obtained using the GeneCards database. MMP9 was highly correlated (Relevance scores > 15) with TOS, ROS, and TOAC, and MPAT was highly correlated with TOS and TOAC (Fig. 5D, supplementary Table 10). A correlation network was constructed between the genes and relevant oxidative stress biomarkers (Relevance scores > 7). 11 nodes (including 2 hub genes, 6 oxidative biomarkers, and 3 Antioxidant biomarkers) and 15 edges were obtained (Fig. 5E). MMP9 was related to 8 oxidative stress biomarkers, and MPAT was 7 oxidative stress biomarkers.

Fig. 5

Relationship and regulation of hub genes. A Heat map of gene expression correlation. Red represented positive correlation; blue represented negative correlation; darker color meant higher correlation, and blank indicated no significance. B Chord diagram of gene expression correlation. Different colors represented different genes, the wider the line, the higher the correlation. C Hub gene functional similarity. Horizontal coordinates were hub gene functional similarity scores. D Hub gene correlation with oxidative stress biomarkers. E Correlation network between hub gene and oxidative stress biomarkers. Purple squares were hub genes. Yellow circles were oxidative biomarkers, and blue circles were antioxidant biomarkers; larger nodes indicated greater degrees; darker edges indicated higher Relevance scores

The hub gene-related biological pathways were examined using GSEA. CA3 was enriched in 403 GO_BP terms, including Cytoplasmic Translation, Microtubule Cytoskeleton Organization, Cytoplasmic Translation, Mitochondrial Gene Expression, and Mitochondrial Translation (Fig. 6A, supplementary Table 12), 77 GO_CC terms, including Collagen Containing Extracellular Matrix, Cytosolic Ribosome, External Encapsulating Structure, Ribosomal Subunit, Ribosome (Fig. 6B, supplementary Table 12), and 116 GO_MF terms, consisting of Extracellular Matrix Structural Constituent, Glycosaminoglycan Binding, Extracellular Matrix Structural Constituent Conferring Compression Resistance, Lyase Activity, and Structural Constituent of Ribosome (Fig. 6C, supplementary Table 12). As for KEGG pathways, CA3 was enriched in 29 pathways, including Valine Leucine, Type I Diabetes Mellitus, Drug Metabolism Cytochrome P450, Oxidative Phosphorylation, Ribosome, and Isoleucine Degradation (Fig. 6D, supplementary Table 12). PPP1R15B was enriched in 434 GO_BP terms including Mitochondrial Gene Expression, Mitochondrial Organization, Proton Transmembrane Transport, and Mitochondrial Translation (Fig. 6E, supplementary Table 13), 85 GO_CC terms included Mitochondrial Matrix, Inner Mitochondrial Membrane Protein Complex, Mitochondrial Protein Containing Complex, Organelle Inner Membrane, and Respirasome (Fig. 6F, supplementary Table 13), and 97 GO_MF terms included ATP Dependent Activity Acting on RNA and Proton Transmembrane Transporter Activity (Fig. 6G, supplementary Table 13). PPP1R15B was enriched in 33 KEGG pathways, such as Pyrimidine Metabolism, Parkinson's' Disease, Oxidative Phosphorylation, Lysosome, and Spliceosome (Fig. 6H, supplementary Table 13). MAPT was enriched in 391 GO_BP terms including Protein_N Linked Glycosylation via Asparagine and Positive Regulation of Protein Polyubiquitination (Fig. 6I, supplementary Table 14), 41 GO_CC terms including Axonemal Microtubule, Cornified Envelope, Manchette, Outer Membrane, Outer Membrane (Fig. 6J, supplementary Table 14), 79 GO_MF terms including N Methyltransferase Activity, Wnt Receptor Activity, Lysine N Methyltransferase Activity, Histone H3 Methyltransferase Activity, and Histone Lysine N Methyltransferase Activity (Fig. 6K, supplementary Table 14), and 19 KEGG pathways including Valine Leucine And_Isoleucine Degradation, Basal Cell Carcinoma, Pathogenic Escherichia Coli Infection, Peroxisome, Wnt Signaling Pathway (Fig. 6L, supplementary Table 14). MMP9 was enriched in 1304 GO_BP terms, including Antigen Receptor-Mediated Signaling Pathway, Adaptive Immune Response, Antigen Processing and Presentation, Antigen Processing and Presentation of Exogenous Antigen, and Positive Regulation of T Cell Proliferation (Fig. 6M, supplementary Table 15), 149 GO_CC terms including Vacuolar Membrane, Tertiary Granule, Ficolin 1 Rich Granule, External Side of Plasma Membrane, and Side of Membrane (Fig. 6N, supplementary Table 15), 161 GO_MF terms including Antigen Binding, Cytokine Receptor Activity, and, Cytokine Receptor Binding (Fig. 6O, supplementary Table 15), and 57 KEGG pathways containing Toll-Like Receptor Signaling Pathway, Leishmania Infection, Hematopoietic Cell Lineage, and Antigen Processing and Presentation, (Fig. 6P, supplementary Table 15). ECT2 was enriched in 363 GO_BP terms, including Response to Transforming Growth Factor Beta, Mitochondrion Organization, Nucleoside Triphosphate Biosynthetic Process, and Mitochondrial Respiratory Chain Complex Assembly (Fig. 6Q, supplementary Table 16), 74 GO_CC terms including Vesicle Tethering Complex, Respirasome, Organelle Inner Membrane, Mitochondrial Matrix, and Mitochondrial Protein Containing Complex (Fig. 6R, supplementary Table 16), 66 GO_MF terms including O Methyltransferase Activity, Catalytic Activity Acting on a Nucleic Acid, DNA Apurinic or Apyrimidinic Site Endonuclease Activity, Oxidoreduction Driven Active Transmembrane Transporter Activity, and Catalytic Activity Acting on RNA (Fig. 6S, supplementary Table 16), and 25 KEGG pathways comprising Oxidative Phosphorylation, Parkinson's Disease, Pyrimidine Metabolism, Huntington's Disease, and Spliceosome (Fig. 6T, supplementary Table 16).

Fig. 6

GSEA result for key genes. A CA3 gene enriched GO-BP terms. B CA3 gene enriched GO-CC terms. C CA3 gene enriched GO-MF terms. D CA3 gene enriched KEGG pathway. E PPP1R15B gene enriched GO-BP terms. F PPP1R15B gene enriched GO-CC terms. G PPP1R15B gene enriched GO-MF terms. H PPP1R15B gene enriched KEGG pathway. I MAPT gene enriched GO-BP terms. J MAPT gene enriched GO-CC terms. K MAPT gene enriched GO-MF terms. L MAPT gene enriched KEGG pathway. M MMP9 gene enriched GO-BP terms. N MMP9 gene enriched GO-CC terms. O MMP9 gene enriched GO-MF terms. P MMP9 gene enriched KEGG pathway. Q ECT2 gene enriched GO-BP terms. R ECT2 gene enriched GO-CC terms. S ECT2 gene enriched GO-MF terms. T ECT2 gene enriched KEGG pathway

To assess immune status variation in COPD patients, immune cell infiltration was performed on the GSE76925 dataset using the CIBERSORT algorithm. Figure 7A displays the stacked abundance of 22 immune cell types. The Wilcoxon test revealed notable variations (P < 0.05) among follicular helper T cells, CD8 T cells, gamma delta T cells, and M0 macrophages (Fig. 7B) (supplementary Table 17). The relationships among immune cell types were investigated using Spearman correlation and visualized through heatmaps. As presented in Fig. 7C and supplementary Table 17, M1 correlated negatively with Dendritic cells activated (cor = −0.72, p < 0.05), but positively with Dendritic cells resting (cor = 0.45, p < 0.05). Spearman correlation analysis revealed that there was a positive correlation between MMP9 and M1 macrophages (p < 0.05, cor = 0.38) but a negative correlation between MMP9 and activated dendritic cells (p < 0.05, cor = −0.41) (Fig. 7D, supplementary Table 17). And, in the LASSO regression analysis, the parameter setting: family = ‘binomial’ results showed that the ideal lambda value (lambda.min) was 0.008925923. As shown in Fig. 7E, a total of four cells with regression coefficients that were not penalized to 0 were finally obtained (Macrophages M0, T cells CD8, T cells follicular helper, T cells gamma delta). In addition, these four differential immune cells showed a strong association with core genes and were notably expressed in the disease group (Fig. 7F-G).

Fig. 7

Immune infiltration analysis. A Proportion of immune cells in high and low risk groups. B Differences in immune cells between high and low risk groups. C Heatmap of immune cell correlation. Purple indicated a positive correlation, green indicated a negative correlation, darker color indicated a higher correlation; blank indicated no significance. D hub Heatmap of gene-immune cell correlation. Red represented positive correlation, blue represented negative correlation, darker color meant higher correlation; significance was indicated by an asterisk. E The results of the LASSO regression analysis. Above panel: the horizontal coordinate is log Lambda and the vertical coordinate represents the local likelihood deviation; below panel: the horizontal coordinate is log Lambda and the vertical coordinate is the coefficient of the feature. F Results of correlation analysis between differential immune cells and core genes. Red represents positive correlation and blue represents negative correlation. G Expression of differential immune cells in disease and control groups. Red represents the COPD group and blue represents the normal group. *, P < 0.05, **, P < 0.01, ***, P < 0.001, ****, P < 0.0001

Regulation network prediction

The direct relationship between transcription factors (TFs) or noncoding (nc) RNAs and hub genes was evaluated to further predict the mechanism by which the hub gene was regulated in COPD. 31 TFs were predicted, and a network included 36 nodes (5 hub genes and 31 TFs) and 46 edges. Of these, FOXC1 was predicted to be related to 4 hub genes, i.e., PPP1R15B, CA3, MMP9, and MAPT (Fig. 8A, supplementary Table 18). There were 229 miRNAs predicted for the 5 hub genes in the Tarbase database. 199 miRNAs were predicted for the 3 hub genes (ECT2, MMP9, PPP1R15B) in the mirTarbase database. miRNA-hub gene regulatory pairs were obtained by retaining the predictions consistent with those of both databases and consisted of the 3 hub genes (ECT2, MMP9, PPP1R15B) and 50 miRNAs. Next, a total of 228 lncRNAs related to 2 hub genes (ECT2, PPP1R15B) were predicted using the miRNet database. A network including 242 nodes (2 hub genes, 12 miRNAs, and 228 lncRNAs) was constructed (Fig. 8B, supplementary Table 19). Then, to understand the associations between the hub gene and respiratory diseases and chronic complications, an analysis was performed with the CTD database. 7 diseases (chronic obstructive pulmonary disease, respiratory insufficiency, non-small cell lung carcinoma, asthma, pulmonary fibrosis, pneumoconiosis, pulmonary emphysema) related to 3 hub genes (MAPT, MMP9, CA3) were predicted. A network including 10 nodes (7 diseases and 3 hub genes) and 7 edges was constructed (Fig. 8C, supplementary Table 20). Also, 3 hub genes (MAPT, MMP9, CA3) were predicted as related to 37 drugs, and a network including 47 nodes (7 diseases, 3 hub genes, and 37 drugs) was constructed (Fig. 8D, supplementary Table 21).

Fig. 8

hub gene regulatory network. A TF-hub gene regulatory network. purple squares were hub genes, yellow circles were TFs, larger nodes indicated greater degree. B ceRNA regulatory network. purple squares were hub genes, green diamonds were miRNAs, blue circles were lncRNAs, larger nodes indicated larger degree; due to the large number of lncRNA nodes, only lncRNAs with degree > 5 were labeled. C hub Gene-Disease Relationship Network. Purple squares were hub genes, green circles were disease names. D Drug-Gene-Disease Network. White squares were hub genes, blue octagons were diseases, yellow circles were drugs, and dark yellow indicated that the drug was approved and light yellow indicated that the drug was unapproved

Hub gene expression in COPD

The hub gene expression was analyzed in GSE76925 and GSE38974. In the training set, GSE76925, CA3, MAPT, and MMP9 expression was significantly higher, while ECT2 and PPP1R15B were notably reduced in COPD patients (Fig. 9A). While in the validation set GSE38974, 3 hub genes (CA3, ECT2, MMP9) showed consistent expression trends between the validation and training sets, and 2 of them (ECT2, MMP9) showed significantly differential expression (Fig. 9C). All hub genes in GSE76925, as well as the ECT2 and MMP9 in the validation set GSE38974, had area under the curve values greater than 0.7 (Fig. 9B, D), indicating that ECT2 and MMP9 had high discriminatory ability in COPD samples. Then, the mRNA expression in peripheral blood was detected. Compared with healthy controls, COPD patients had significantly increased MMP9 but significantly decreased ECT2 (Fig. 10A-B). Since MMP9 was significantly expressed in GSE76925 and GSE38974, as well as in the in-house cohort. As shown in the above results, MMP9 was highly related to oxidative stress biomarkers and respiratory tract disease; MMP9 was selected for the subsequent experiments. The CS-induced COPD model was constructed, and HE was performed to evaluate lung tissue histopathology. As a result, the tissues from the control group showed intact alveolar lumens, no detachment of epithelial cells from the mucosa and necrosis, no hyperplasia of submucosal glands, and few inflammatory cells in the airway periphery were observed. There were no secretions in the tracheal lumen, basement membranes, or smooth muscle hyperplasia. While in the CS-COPD group, the alveolar septa were thickened, the basement membranes and smooth muscles of the bronchi were hyperplastic, the bronchial walls were thickened, there was a large number of inflammatory cells infiltrating, epithelial cell detachment or secretion was observed in the lumen of the trachea, and there were varying degrees of stenosis in the trachea (Fig. 11A). The lung function was detected. The COPD group had lower minute ventilation volume, peak expiratory flow, and peak inspiratory volume (Fig. 11B-D) as well as declined cytokines (Fig. 11E-G). As presented in Fig. 12H-K, the SOD, CAT, and GSH were reduced, while MDA was increased in the COPD group, suggesting the occurrence of OS injury in the CS-COPD model. Also, MMP9 expression was detected by immunohistochemistry and qPCR. MMP9 was located in the extracellular matrix and cytoplasm, and MMP9 was highly expressed in the COPD group (Fig. 11L-N). Since MMP9 was highly positively related to M1 macrophages, we detected the M1 and M2 cell surface markers. The M1 macrophage markers, including CD86, iNOS, and MCP1, were all increased. In contrast, the markers of M2 macrophages, such as IL-10, ARG1, CD206, and CD163 were all reduced in the COPD group (Fig. 11O-P), indicating the significant difference of M1 macrophage exerted in COPD.

Fig. 9

hub gene expression and ROC curve. A Hub gene expression in GSE76925 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc = GSE38974). B ROC curves of hub gene in GSE76925. C Hub gene expression in GSE38974. D ROC curves of hub gene in GSE38974. *, P < 0.05, **, P < 0.01, ***, P < 0.001

Fig. 10

Hub gene expression in peripheral blood. A MMP9 expression was detected in peripheral blood from control and COPD Patients. B ECT2 expression was detected in peripheral blood from control and COPD Patients. ***, P < 0.001, ****, P < 0.0001

Fig. 11

MMP9 was highly expressed in CS-induced COP. A HE staining was used to detect the pathology changes in alveoli and trachea in CS-COPD model. (scale bar = 50 μm). B-D Minute ventilation volume, peak expriratory flow and peak inspiratory volume was detected to reflect lung function in the model. E–G The content of IL-6, IL-8, and TNF-α was detected. H–K The content of SOD, MDA, CAT and GSH was detected. L MMP9 expression was detected by qPCR. M–N MMP9 expression was detected by IHC (scale bar = 50 μm). O-P The M1 macrophage surface markers, including CD86, iNOS, MCP1, and M2 macrophage surface markers, including ARG1, CD163, CD206, and IL-10 were detected.**, P < 0.01, ***, P < 0.001, ****, P < 0.0001

Fig. 12

Deprivation of MMP9 alleviated the CSC-induced injure in BEAS-2B. A The screening of CSC concentration in BEAS-2B. B qPCR was used to detect the transfection efficiency of lentivirus. C The mRNA expression of MMP9 was detected in BEAS-2B. D The protein expression of MMP9 was detected in BEAS-2B. E Cell viability was detected by CCK-8. E–H The concentration of IL-6, IL-8 and TNF-α was detected. H–L The content of MDA, SOD, CAT, and GSH-Px was detected. M–N The content of ROS was detectedO-P) O-P Cell apoptosis rate was detected by flow cytometer. *, P < 0.05, **, P < 0.01, ***, P < 0.001, ****, P < 0.0001

Deprivation of MMP9 alleviated the CSC-induced injury and promoted the macrophage polarization to M1 macrophage

To evaluate the appropriate concentration of CSC to stimulate the BEAS-2B, the freshly prepared 100% CSC was diluted to concentrations at 1%, 2%, 4%, 8%, and 15% in serum-free medium, and BEAS-2B was stimulated by different concentrations of CSC for 24 h. CCK-8 assessed cell viability. When the concentration of CSC reached 4%, the percentage of BEAS-2B viability decreased to 50%, and IC50 was calculated as 4.079%. Therefore, we chose 4% CSC to stimulate the BEAS-2B in subsequent research (Fig. 12A). And the lentivirus was constructed and transfected to deprive the MMP9 expression in BEAS-2B. The transfection efficiency was detected with qPCR. As a result, sh-MMP9#3 showed better expression suppression and was selected for further detection (Fig. 12B). The MMP9 level was up-regulated in CSC-BEAS-2B than in BEAS-2B but down-regulated when sh-MMP9 was transfected (Fig. 12C-D). The cell viability was decreased in CSC-BEAS-2B compared to BEAS-2B, which was increased after the deprivation of MMP9 expression (Fig. 12E). The cytokine levels were elevated in CSC-BEAS-2B than in BEAS-2B but decreased after the deprivation of MMP9 expression (Fig. 12F-H). ROS content and oxidative stress marker levels were determined to assess the presence of oxidative stress injury induced by CSC. MDA and ROS levels were increased, whereas CAT, SOD, and GSH-Px levels exhibited a decrease following CSC stimulation. Upon depletion of MMP9 expression, a decrease was observed in MDA and ROS levels, whereas CAT, SOD, and GSH-Px levels showed an elevation (Fig. 12I-N). The apoptosis rate of BEAS-2B was increased in CSC-BEAS-2B compared to BEAS-2B, which was decreased when sh-MMP9 was transfected (Fig. 12O-P). This data suggests the key role of MMP9 in CSC-induced oxidative stress injury.

Then, to detect the role of MMP9 in macrophage polarization, the co-culture of PMA-induced THP-1 cells with CSC-BEAS-2B was performed. THP-1 cells were suspension cells, showing a small round shape with regular morphology, and PMA-induced THP-1 cells were transformed into primary macrophages (M0 macrophages), which were irregularly polygonal in shape and grew adherently to the wall (Fig. 13A). The positive rate of CD68 was increased after the induction of PMA (Fig. 13B-C), which indicated the successful induction of M0 macrophage. Also, the expression of MMP9 in co-cultured macrophages was detected. As a result, MMP9 mRNA expression was increased in M0 that co-cultured with CSC-BEAS-2B compared to M0 or M0 that co-cultured with BEAS-2B, and decreased when MMP9 lentivirus was transfected in CSC-BEAS-2B (Fig. 13D). The surface markers of M1 and M2 were all increased in M0 co-cultured with CSC-BEAS-2B, compared to M0 or M0 that co-cultured with BEAS-2B. In contrast, transfection of MMP9 lentivirus in CSC-BEAS-2B led to a decrease in M1 surface markers like CD86, iNOS, and MCP-1, while M2 markers such as IL-10, ARG1, CD206, and CD163 increased (Fig. 13E-K, Supplementary Fig. 1). These results showed that deprivation of MMP9 alleviated the CSC-induced injury and promoted the macrophage polarization to M2 macrophage (Fig. 14).

Fig. 13

Deprivation of MMP9 promoted the macrophage polarization to M1 macrophage. A The morphology ofTHP-1 and PMA induced THP-1 (scale bar = 20 μm). B-C The surface marker CD68 was detected. D The mRNA expression of MMP9 was detected in Macrophage. E-K The M1 macrophage surface markers, including CD86, iNOS, MCP1, and M2 macrophage surface markers, including ARG1, CD163, CD206, and IL-10 were detected. *, P < 0.05, **, P < 0.01, ***, P < 0.001, ****, P < 0.0001

Fig. 14