Authors: Zhiwei Lin, Yueting Jiang, Huifang Liu, Juhua Yang, Bin Yang, Ke Zhang, Peiren Tang, Bo Xiang, Baoqing Sun
Categories: Research, Pulmonary disease, Chronic obstructive, Dysbiosis, Oxidative stress, Metagenomics, Transcriptome, Immune system, Precision medicine
Source: Journal of Translational Medicine
Authors: Zhiwei Lin, Yueting Jiang, Huifang Liu, Juhua Yang, Bin Yang, Ke Zhang, Peiren Tang, Bo Xiang, Baoqing Sun
Chronic Obstructive Pulmonary Disease (COPD) is characterized by progressive airflow limitation and chronic inflammation. Although airway microbes and host immunity are known contributors, the molecular mechanisms underlying disease severity remain unclear. This study explores microbial dysbiosis and host immune responses across varying COPD severities.
We conducted integrated metagenomic and transcriptomic analyses on bronchoalveolar lavage fluid from two a discovery cohort and a validation cohort. We investigated microbial diversity, pathogenic bacterial enrichment, and host gene expression patterns. Functional metagenomics was used to assess antibiotic resistance genes. Host-microbe network analyses explored correlations between pathogens and immune-metabolic pathways. Diagnostic models utilizing microbial-immune biomarkers were developed, trained on a subset of the discovery cohort, tested on remaining discovery samples, and validated by quantitative polymerase chain reaction (qPCR) in the validation cohort to distinguish COPD from controls and stratify disease severity.
Severe COPD exhibited reduced microbial diversity and an increased presence of pathogenic bacteria, including Moraxella osloensis and Streptococcus species. These pathogens were associated with dysregulated inflammatory signaling, and significant neutrophil activity, evidenced by the formation of Neutrophil Extracellular Traps (NETs), and oxidative stress, which correlated with airway remodeling and a decline in lung function. Functional metagenomics showed a significant increase in antibiotic resistance genes in severe cases, linked to chronic treatment pressures. Host-microbe network analyses revealed strong correlations between these pathogens and disrupted immune-metabolic pathways, such as altered energy metabolism and inflammatory cascades, consistent across both cohorts. Diagnostic models based on microbial-immune biomarkers demonstrated high accuracy in differentiating COPD patients from controls and in stratifying disease severity.
This study identifies microbial and immune signatures associated with COPD severity, providing mechanistic insights into its pathophysiology. The findings may inform precision medicine strategies by targeting airway dysbiosis and immune dysregulation. While causal relationships could not be established in this cross-sectional study, the findings provide a foundation for future mechanistic investigations using advanced in vitro and in vivo models.
The online version contains supplementary material available at 10.1186/s12967-025-06986-2.
Chronic Obstructive Pulmonary Disease (COPD) is a chronic and heterogeneous pulmonary disorder characterized by persistent airflow limitation, chronic inflammation, and progressive lung tissue destruction. The disease is commonly driven by exposure to noxious particles or gases, particularly cigarette smoke [1, 2]. Despite advancements in diagnostic tools, such as spirometry and high-resolution computed tomography (HRCT), and growing insights into airway pathophysiology, effective clinical management of COPD remains a major challenge [3]. The substantial heterogeneity among patients, including differences in disease severity, etiological factors, and immune responses, complicates accurate diagnosis and the development of individualized therapeutic strategies [4].
Historically, COPD was thought to result primarily from chronic exposure to irritants that induce airway inflammation and remodeling. Bacterial infections were considered a secondary factor, typically involving pathogens such as Pseudomonas aeruginosa and Haemophilus influenzae [5]. However, accumulating evidence now suggests that COPD is associated with a broader disruption of the airway microbiome, underscoring the role of microbial communities in influencing disease progression and clinical outcomes.
Advances in high-throughput sequencing technologies have enabled comprehensive profiling of the respiratory microbiome, allowing unbiased identification of bacterial, viral, and fungal communities within bronchoalveolar lavage fluid (BALF) samples [6]. These studies have revealed a marked reduction in microbial diversity and the expansion of pathogenic taxa in patients with COPD [7, 8]. Although airway dysbiosis is increasingly recognized as a key factor driving disease progression, the mechanisms by which dominant pathogens interact with the broader microbiota and the host immune system remain poorly defined [9]. A deeper understanding of these host–microbe interactions is essential for identifying critical mediators of COPD progression and for developing more targeted therapeutic approaches.
Microbial dysbiosis, oxidative stress, and immune activation represent core elements of COPD pathogenesis, together fueling a self-reinforcing cycle of chronic airway inflammation and structural damage [10]. Dysbiosis often presents as reduced microbial diversity and an increased dominance of pathogenic bacteria such as Pseudomonas aeruginosa and Haemophilus influenzae, which promote recurrent infections and contribute to lung injury [9]. These organisms release virulence factors, such as lipopolysaccharides, that stimulate host immune responses, leading to excessive production of reactive oxygen species (ROS) by both neutrophils and bacteria. In our previous research, we demonstrated that ROS induce lipid peroxidation, protein oxidation, and DNA damage. These processes compromise epithelial function, enhance mucus secretion, and impair mucociliary clearance [11]. Concurrently, immune activation driven by dysregulated signaling promotes elevated production of proinflammatory cytokines and the formation of neutrophil extracellular traps (NETs) [11, 12]. While NETs serve as an antimicrobial defense, their overproduction contributes to epithelial damage and supports bacterial biofilm formation, thereby exacerbating infection [13]. When integrated with metagenomic data, transcriptomic analyses can provide a more holistic view of how specific microbial signatures influence host immune responses. This multi-omics approach is increasingly recognized as critical for unraveling the complex interplay between airway pathogens and the host in COPD.
Despite recent advances, substantial knowledge gaps remain in elucidating the mechanisms underlying COPD pathogenesis, especially regarding the interactions among airway dysbiosis, immune dysregulation, and pulmonary function decline [2]. Many previous microbiome and transcriptomic studies have been limited by small sample sizes, inconsistent sampling protocols, and variability in DNA/RNA extraction and analytical methods, all of which restrict reproducibility and generalizability [14]. Moreover, the integration of molecular data with clinical characteristics and pulmonary function indices is still underdeveloped [15]. Although phenotypic clustering has been used to classify COPD into subtypes, linking these phenotypes to detailed microbial and immune signatures remains an underutilized avenue for precision medicine. Defining COPD endotypes, distinct subgroups characterized by specific microbial, molecular, and immunological features related to lung function impairment, holds great potential for therapeutic stratification. Realizing this potential, however, requires robust methodologies capable of integrating complex multi-omics datasets with comprehensive clinical and pulmonary function profiling.
Our study aims to address these challenges by integrating metagenomic and transcriptomic analyses of BALF samples with detailed clinical and pulmonary function data. Unlike prior studies that have focused on single-omics approaches or a narrow range of pathogens, our multi-omics design enables the identification of novel microbial and immune signatures associated with lung function decline [16]. Furthermore, we developed diagnostic models for COPD detection and for stratifying disease severity, thereby laying the foundation for precision medicine approaches tailored to endotypes defined by airway microbial–immune interactions and objective lung function measurements.
In this study, we systematically investigated the airway microenvironment in COPD patients compared to controls using a multi-omics framework. Two independent cohorts were analyzed to validate the relationship between lung function impairment and airway microbiota–host interaction networks. The discovery cohort comprised 88 COPD patients (53 with severe to very severe impairment and 35 with mild to moderate impairment) and 20 controls. The validation cohort included 43 COPD patients (22 with severe to very severe impairment and 21 with mild to moderate impairment) and 18 controls. BALF samples from both cohorts underwent metagenomic sequencing to comprehensively profile the airway microbiota. In parallel, transcriptomic analyses of host cells isolated from BALF were conducted to investigate gene expression profiles related to inflammation, immune regulation, and tissue remodeling. By integrating microbial and transcriptomic data and correlating these profiles with lung function parameters, this study sought to uncover key host–microbe interaction networks driving disease progression, identify biomarkers for clinical stratification, and elucidate the molecular pathways underlying COPD pathogenesis. Our findings are intended to enhance the mechanistic understanding of COPD and support the development of precision medicine strategies for its management.
This single-center, observational study was designed to investigate microbial and host transcriptomic alterations in COPD and their association with lung function severity. The study included both a discovery and a validation cohort. The discovery cohort comprised 88 patients with COPD and 20 control subjects, all recruited from The First Affiliated Hospital of Guangzhou Medical University between September 2023 and January 2025. The validation cohort included 43 COPD patients and 18 Controls, also recruited from the same institution. Ethical approval was obtained from the Ethics Committee of The First Affiliated Hospital of Guangzhou Medical University (approval 2022 No.121 for the discovery cohort and 2024 No. G-007 for the validation cohort). Written informed consent was obtained from all participants prior to enrollment. The study was conducted in accordance with the principles of the Declaration of Helsinki.
COPD was diagnosed in accordance with the Global Initiative for Chronic Obstructive Lung Disease (GOLD) 2023 guidelines. Eligible patients were required to demonstrate a post-bronchodilator forced expiratory volume in 1 s to forced vital capacity ratio (FEV1/FVC) of less than 0.70, confirming persistent airflow limitation, along with chronic respiratory symptoms such as dyspnea, chronic cough, or sputum production. Exclusion criteria (1) known respiratory disorders other than COPD; (2) history of lung surgery or tuberculosis; (3) cancer diagnosis; (4) recent blood transfusion (within four weeks of enrollment); (5) autoimmune diseases; (6) participation in a blinded drug trial; (7) antibiotic use within the past eight weeks; (8) endocrine disorders such as diabetes or connective tissue diseases; (9) other systemic inflammatory conditions; and (10) severe cardiovascular, cerebrovascular, hepatic, renal, or psychiatric disorders, or a history of pulmonary surgery. Only patients who were clinically stable and eligible for bronchoscopy were included.
In the discovery cohort, COPD patients were stratified according to lung function severity using GOLD classification based on post-bronchodilator FEV1% predicted. Of the 88 patients, 53 (60.2%) were classified as having severe to very severe lung function impairment (SLF group), corresponding to GOLD Stage III (FEV1 30–49% predicted) or Stage IV (FEV1 < 30% predicted). The remaining 35 patients (39.8%) were categorized as having mild to moderate impairment (NLF group), corresponding to GOLD Stage I (FEV1 ≥ 80% predicted) or Stage II (FEV1 50–79% predicted). In the validation cohort, 22 of 43 patients (51.2%) were classified into the SLF group, and 21 (48.8%) into the NLF group using the same criteria. All participants underwent comprehensive clinical phenotyping, including pulmonary function testing and BALF collection for subsequent metagenomic and transcriptomic analyses. Control subjects were age- and sex-matched individuals without a history of chronic respiratory disease, with normal spirometry values (FEV1/FVC > 0.70 and FEV1 ≥ 80% predicted), and no significant exposure to COPD-related risk factors.
BALF samples were collected from all participants following a standardized bronchoscopy protocol previously established by our team [11]. Under local anesthesia, a flexible bronchoscope was introduced into the airway, and bronchoalveolar lavage was performed in lung segments identified as representative based on HRCT imaging. To minimize contamination from the upper airways, the first 10 mL of lavage fluid was discarded. The remaining BALF was immediately placed on ice and transported to the laboratory within two hours.
In the laboratory, BALF samples were centrifuged at 300 × g for 10 min at 4 °C to separate the cellular fraction from the supernatant. The supernatant was then aliquoted and stored at –80 °C for nucleic acid extraction and downstream analysis. Samples exhibiting insufficient volume, visible contamination, or hemolysis were excluded. Rigorous quality control procedures were implemented throughout the sample collection and processing phases to ensure integrity and reproducibility.
Pulmonary function tests (PFTs) were conducted in accordance with the American Thoracic Society (ATS) guidelines to evaluate respiratory performance. Following bronchodilator administration, FEV1, forced vital capacity (FVC), and the FEV1/FVC ratio were measured. The percentage of predicted FEV1 (FEV1%) and FEV1/FVC values were calculated using standardized reference equations adjusted for age, sex, and height (Tables 1 and 2). These indices were chosen based on their established utility in assessing airflow obstruction in COPD.Table 1Basic information of discovery cohort (COPD and Controls)CharacteristicsCOPDControlsPNo8820Age, year66.83 ± 8.1765.50 ± 8.250.514BMI, kg/m^2^21.04 ± 3.5622.197 ± 2.090.063Current smoker38 (43.18%)8 (40.00%)0.795Lung FunctionFEV1 predicted2.55 ± 0.552.78 ± 0.590.101FEV1, L1.15 ± 0.553.00 ± 0.59 < 0.001FEV1 (A/Pd),%47.15 ± 21.82108.47 ± 8.08 < 0.001FVC predicted3.25 ± 0.593.33 ± 0.740.619FVC, L2.35 ± 0.743.683 ± 0.76 < 0.001FEV1/FVC0.48 ± 0.150.81 ± 0.06 < 0.001GOLD-stageI8 (9.10%)//II27 (30.68%)//III31 (35.22%)//IV22 (25.00%)//Inflammatory markersWBC,10^9^/L7.71 ± 2.786.32 ± 1.940.037NEU,10^9^/L5.41 ± 2.543.82 ± 1.600.009LYM,10^9^/L1.52 ± 0.591.69 ± 0.630.271MONO,10^9^/L0.58 ± 0.270.36 ± 0.10 < 0.001EOS,10^9^/L0.17 ± 0.190.40 ± 1.080.349FEV1 Forced expiratory volume in one second; FVC Forced vital capacity; WBC White blood cell; NEU Neutrophil; LYM Lymphocyte; MONO Monocyte; EOS EosinophilTable 2Basic information of discovery cohort (SLF and NLF)CharacteristicsSLFNLFPNo5335Age, year67.09 ± 7.9066.42 ± 8.670.711BMI, kg/m^2^20.65 ± 3.2321.65 ± 3.970.199Current smoker22 (41.50%)16 (45.71%)0.697Lung functionFEV1 predicted2.64 ± 0.482.41 ± 0.610.060FEV1, L0.83 ± 0.261.64 ± 0.51 < 0.001FEV1 (A/Pd),%32.15 ± 9.8269.88 ± 13.71 < 0.001FVC predicted3.32 ± 0.583.16 ± 0.600.216FVC, L2.11 ± 0.612.71 ± 0.78 < 0.001FEV1/FVC0.40 ± 0.120.61 ± 0.08 < 0.001GOLD-stageI31 (58.49%)//II22 (41.51%)//III/8 (22.85%)/IV/27 (77.15%)/Inflammatory markersWBC,10^9^/L8.20 ± 2.726.96 ± 2.760.040NEU,10^9^/L5.89 ± 2.444.67 ± 2.560.027LYM,10^9^/L1.52 ± 0.671.52 ± 0.450.967MONO,10^9^/L0.60 ± 0.270.56 ± 0.260.512EOS,10^9^/L0.18 ± 0.190.15 ± 0.190.586FEV1 Forced expiratory volume in one second; FVC Forced vital capacity; WBC White blood cell; NEU Neutrophil; LYM Lymphocyte; MONO Monocyte; EOS Eosinophil
Comparative analyses between cohorts included assessments of demographic variables and inflammatory markers such as white blood cells (WBC), neutrophils (NEU), lymphocytes (LYM), monocytes (MONO), and eosinophils (EOS), to further characterize the systemic immune profile across different disease stages.
Microbial DNA and RNA were extracted independently from the same samples to ensure optimal yield and integrity for downstream applications. Prior to nucleic acid extraction, host DNA was selectively depleted by treatment with 1 U of nuclease (Thermo Fisher Scientific, USA) and 0.5% Tween 20 (Sigma, USA), following established protocols [17]. DNA extraction was performed using the QIAamp® UCP Pathogen DNA Kit (Qiagen, Valencia, CA, USA) according to the manufacturer’s instructions. RNA extraction was carried out separately using the QIAamp® Viral RNA Kit (Qiagen). To enrich for microbial transcripts, human ribosomal RNA was removed from total RNA using the Ribo-Zero rRNA Removal Kit (Illumina, San Diego, CA, USA). Complementary DNA (cDNA) was synthesised from purified RNA using reverse transcriptase and deoxynucleotide triphosphates (dNTPs). All extracted nucleic acids were quantified with a Qubit Fluorometer (Thermo Fisher Scientific) and stored at − 80 °C until library preparation.
DNA and cDNA libraries were constructed using the Nextera XT DNA Library Prep Kit (Illumina, FC-131–1096). For metagenomic libraries, DNA was tagmented to fragment and attach adapters, using a minimum input of 1 ng. For transcriptomic libraries, total RNA was converted to cDNA using the SuperScript IV First-Strand Synthesis System (Thermo Fisher Scientific, 18,091,050) with random hexamer primers for first-strand synthesis, followed by second-strand synthesis with the NEBNext Ultra II Non-Directional RNA Second Strand Synthesis Module (New England Biolabs, E6111). Tagmented DNA and cDNA underwent PCR amplification (12 cycles) with Nextera XT Index Kit primers (Illumina, FC-131–2001) to add unique dual-index barcodes for multiplexing. Libraries were purified using AMPure XP beads (Beckman Coulter, A63881) at a 1.8:1 bead-to-sample ratio to remove primer dimers and fragments < 150 bp. Library quality was assessed using the Agilent 2100 Bioanalyzer with the High Sensitivity DNA Kit (Agilent, 5067–4626) to confirm fragment sizes of 250–300 bp and absence of adapter dimers. DNA concentration was measured with the Qubit dsDNA HS Assay Kit (Thermo Fisher Scientific, Q32854) on a Qubit 4 Fluorometer, targeting ≥ 2 ng/µL. Effective concentrations were quantified by quantitative PCR (qPCR) using the KAPA Library Quantification Kit (Roche, 07960140001) on a LightCycler 480 II (Roche), employing a tenfold serial dilution standard curve to ensure equimolar pooling at 4 nM per library.
Libraries were sequenced on the Illumina NextSeq 550 platform using the NextSeq 500/550 High Output Kit v2.5 (Illumina, 20,024,907) for 150-bp paired-end sequencing, targeting an average of 20 million reads per sample. This yielded approximately 6 Gb of raw data for metagenomic libraries (equivalent to 20 M paired-end reads at 150 bp) and 10 Gb for transcriptomic libraries to accommodate the higher complexity of host gene expression [18]. A PhiX Control v3 (Illumina, FC-110–3001) was spiked at 1% to enhance base calling accuracy. The run was configured for 400 million total reads, with samples randomized across flow cell lanes to minimize lane-specific biases. Sequencing was performed in a single run using the NextSeq System Suite (v1.4) to eliminate batch effects. Real-time quality metrics were monitored via Illumina BaseSpace Sequence Hub, ensuring cluster density of 180–220 K/mm^2^ and Q30 scores > 80%.
Raw FASTQ files were processed with fastp (v0.20.1) to remove low-quality reads, adapter sequences, and reads with > 10% ambiguous bases (N) or base quality scores ≤ 20 for over 50% of the read length. Reads were trimmed if the average quality score in a 4-bp sliding window fell below 15, and paired-end reads were discarded if either read failed quality thresholds. Adapters were automatically trimmed using fastp, and reads < 50 bp post-trimming were excluded. Quality metrics, including GC content, read duplication rates, adapter contamination, and per-base quality scores, were assessed with FastQC (v0.11.9), confirming mean Q30 scores > 85% and duplication rates < 10% after filtering. For metagenomic data, human-derived reads were removed by mapping to the human genome (GRCh38) using Bowtie2 (v2.4.2), retaining unmapped reads. For transcriptomic data, reads were aligned to GRCh38 using STAR (v2.7.9a) with parameters, yielding 80% mapped reads and an average of ~ 16 M high-quality reads per sample after filtering. Alignment statistics, including mapping rates and read distribution, were generated with SAMtools (v1.12).
All clean data were mapped to the human genome hg38 using HISAT2 with default parameters, to remove host data [19]. Microbial taxonomy was profiled using Metaphlan and Kraken 2 (https://ccb.jhu.edu/software/kraken2/), leveraging their respective reference databases to generate genus- and species-level abundance profiles [20]. Within-sample diversity was assessed using alpha-diversity indices (Shannon, Richness, Chao1, and ACE), while between-sample diversity was evaluated using beta-diversity metrics (Bray–Curtis dissimilarity). Beta-diversity results were visualized with principal coordinate analysis (PCoA), non-metric multidimensional scaling (NMDS), and Principal component analysis (PCA). Statistical significance was assessed using PERMANOVA (p < 0.05).
Differentially abundant taxa were identified using Wilcoxon rank-sum tests and LDA Effect Size (LEfSe), highlighting key microbial taxa associated with COPD and lung function clusters [21]. Functional pathway analysis was conducted using HUMAnN2 [22, 23], which mapped microbial reads to the Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG), and MetaCyc databases [24]. Enriched pathways were visualized to explore microbial functional adaptations. Antibiotic resistance genes (ARGs) were identified using DeepARG [25]. Comparisons of ARG diversity and abundance were performed across COPD and control groups, as well as between SLF and NLF clusters.
All clean data were mapped to the human genome hg38 using HISAT2 with default parameters, and then FeatureCounts. Normalized read counts were generated with DESeq2, and differentially expressed genes (DEGs) were identified with an adjusted p-value < 0.01 and absolute log2 fold-change ≥ 2. With the application of clusterProfiler package, GO, KEGG, and Reactome pathway enrichment analysis of DEGs, as well as Gene Set Enrichment Analysis (GSEA) were performed, with an adjusted P value < 0.05 showed significant enrichment [26]. Weighted Gene Co-expression Network Analysis (WGCNA) identified gene modules associated with clinical traits such as lung function indices. Enriched modules were analyzed for functional significance [27]. Overlapping hub genes from WGCNA and DEGs were identified and further analyzed for biological relevance. The immune-related genes list was downloaded from the ImmPortDB database (https://immport.org/shared/home) [28]. The CIBERSORT algorithm was applied to infer immune cell proportions using the LM22 gene signature matrix [29].
FastSpar was used to infer correlations between microbial taxa and host gene expression [30]. Significant correlations were visualized in heatmaps, while Sankey diagrams and network analyses were constructed to link microbial species, host genes, and functional pathways. Networks highlighted interactions relevant to immune dysregulation, oxidative stress, and COPD progression.
To evaluate the diagnostic potential of host transcriptomic features, two predictive models were one for discriminating COPD patients from healthy controls (HC), and another for stratifying COPD patients into subgroups with severe to very severe lung function impairment (SLF) versus mild to moderate impairment (NLF). A robust feature selection and modeling pipeline was implemented, combining least absolute shrinkage and selection operator (LASSO) regression and random forest (RF) algorithms for feature selection, followed by logistic regression as the primary classification algorithm. Both model development and independent validation were conducted using well-defined discovery and validation cohorts to ensure generalizability and reproducibility.
For the COPD versus HC classification task, the discovery cohort included 88 COPD patients (53 SLF and 35 NLF) and 20 Controls. This dataset was randomly partitioned into a training set (62 COPD 35 SLF and 27 NLF; 13 Controls) and an internal test set (26 COPD 18 SLF and 8 NLF; 7 Controls). The independent validation cohort comprised 43 COPD patients (22 SLF and 21 NLF) and 18 healthy controls. Gene expression markers identified in the discovery cohort were validated using quantitative polymerase chain reaction (qPCR) in the validation cohort.
For the SLF versus NLF classification model, the discovery cohort was similarly divided into a training set (40 SLF and 28 NLF patients) and an internal test set (13 SLF and 7 NLF patients). The same validation cohort (22 SLF and 21 NLF patients) was used for independent validation with qPCR verification of selected markers.
Feature selection was performed using a two-step strategy. First, LASSO regression with tenfold cross-validation was applied to minimize binomial deviance, yielding a sparse set of predictive features. Second, an RF model was trained using 500 decision trees, with the number of variables tried at each split set to the square root of the total number of input features. Features were ranked by mean decrease in Gini impurity, and those above the median were retained. The intersection of features selected by both LASSO and RF was then overlapped with differentially expressed genes to ensure biological relevance and reduce overfitting. The resulting refined feature set was used to construct logistic regression models for both classification tasks.
Model hyperparameters were tuned using fivefold cross-validation. Final models were trained on the discovery training set, evaluated on the internal test set, and validated on the independent cohort. Diagnostic performance was assessed using receiver operating characteristic (ROC) curves, with area under the curve (AUC) values calculated for both discovery and validation datasets. Additional metrics including sensitivity, specificity, positive predictive value (PPV), negative predictive value (NPV), and overall accuracy were also computed. All statistical analyses were performed using R software (version 4.3.2), as detailed in the Data and Software Availability section.
qPCR was performed to validate biomarker expression in BALF samples from the validation cohort. Total RNA was extracted from BALF cells using the RNeasy Mini Kit (Qiagen), followed by reverse transcription to cDNA with the SuperScript IV First-Strand Synthesis System (Thermo Fisher Scientific). qPCR was conducted using SYBR Green Master Mix (Applied Biosystems) on a QuantStudio 5 Real-Time PCR System, with primers specific to the selected biomarkers. Relative expression was quantified using the 2^−ΔΔCt^ method, normalized to GAPDH as the reference gene.
GAPDH (Glyceraldehyde-3-phosphate dehydrogenase).
Forward: 5'-GAAGGTGAAGGTCGGAGTCA-3'.
Reverse: 5'-GAAGATGGTGATGGGATTTC-3'.
GABARAP (GABA Type A Receptor-Associated Protein).
Forward: 5'-AAGAAGAGCATCCGTTCGAGA-3'.
Reverse: 5'-TTCTACTATCACCGGCACCC-3'.
IGSF6 (Immunoglobulin Supersfamily Member 6).
Forward: 5'-TCTATGTGACCGGTGTGTGC-3'.
Reverse: 5'-GCCGAGCACTCTTCTTCTTT-3'.
MARCKS (Myristoylated Alanine-Rich C-Kinase Substrate).
Forward: 5'-TTGAAGAAGCCAGCATGGGTG-3'.
Reverse: 5'-CTTCACGTGGCCATTCTCCTGT-3'.
STPG4 (Sperm-Tail PG-Rich Repeat Containing 4).
Forward: 5'-AGCTTTCTCCGGGGCAATAC-3'.
Reverse: 5'-CAGGGCCTTCATGGGGAATA-3'.
SPINDOC (SPIN1 Docking Protein).
Forward: 5'- GGACTCACCCAAAGACAGGG-3'.
Reverse: 5'- GGAGGGGGAGCTGGAGAG-3'.
Descriptive statistics were used to summarize clinical, microbiome, and transcriptomic data. Continuous variables were expressed as mean ± standard deviation (SD) for normally distributed data or as median with interquartile range (IQR) for non-normally distributed data. Normality was assessed using the Shapiro–Wilk test. Between-group comparisons were conducted using independent t-tests for normally distributed variables and Wilcoxon rank-sum tests for non-normally distributed variables. Categorical variables were analyzed using the chi-square test or Fisher’s exact test, as appropriate. All statistical analyses were performed using R software (version 4.3.2), and the complete dataset and corresponding analysis scripts, including all package dependencies and pipelines, are publicly available via the Zenodo repository (10.5281/zenodo.15632167), as detailed in the Data and Software Availability section.
The discovery cohort consisted of 88 patients diagnosed with COPD and 20 age-matched Controls. Clinical characteristics and pulmonary function were assessed, with results summarized in Table 1. No significant differences were observed in mean age or BMI between the COPD and control groups. Smoking status was comparable, with 43.18% of COPD patients and 40.00% of controls identified as current smokers. Predicted FEV1 (p = 0.101) and predicted FVC (p = 0.619) values did not differ significantly between COPD patients and controls, indicating comparable baseline demographic characteristics between the two groups. In contrast, measured FEV1 was significantly reduced in COPD patients compared to controls (p < 0.001), reflecting pronounced airflow limitation. Similarly, measured FVC was markedly lower in the COPD group (p < 0.001). The FEV1/FVC ratio, a key diagnostic indicator of airflow obstruction, was also significantly decreased in COPD patients (p < 0.001), consistent with the characteristic pathophysiology of the disease. Systemic inflammation was evident in COPD patients, characterized by significantly elevated WBC counts, neutrophil counts, and monocyte counts compared to controls (p < 0.001 for all). In contrast, lymphocyte and eosinophil counts did not differ significantly between groups. These findings underscore the pronounced airflow limitation and systemic inflammatory response in COPD, driven primarily by neutrophilic and monocytic activity, aligning with the established inflammatory and obstructive mechanisms of the disease.
PCoA based on Bray–Curtis distances demonstrated clear clustering patterns separating COPD samples from controls (Fig. 1A). Statistical testing via PERMANOVA confirmed significant differences in community structure between the two groups (R^2^ = 0.0417, p < 0.01). Alpha diversity metrics were consistently reduced in COPD patients. The median Shannon index was significantly lower in COPD samples (2.58, IQR: 2.09–3.01) compared to controls (3.21, IQR: 2.77–3.59, p < 0.05), with similar reductions observed in other indices, including the Simpson index (0.72 vs. 0.84, p < 0.05), Chao1, and Pielou’s Evenness (Fig. 1B).Fig. 1Distinct Microbiome Composition and Reduced Diversity in COPD. A PCoA based on Bray–Curtis distances demonstrate distinct clustering between COPD (n = 88) and Controls (n = 20). B Alpha diversity metrics, including Shannon, Simpson, Chao1, and Pielou’s Evenness indices, show significantly reduced microbial diversity in COPD compared to controls. , P < 0.05. , P < 0.01., P < 0.001.** C** Relative abundances of microbial genera reveal dominance of Streptococcus, Haemophilus, and Pseudomonas in COPD, whereas controls exhibit a more even distribution of taxa. D Species-level analysis identifies enrichment of Haemophilus influenzae, Pseudomonas aeruginosa, and Streptococcus pneumoniae in COPD. (E) LEfSe was used to identify taxa with significantly different relative abundances between groups. Species enriched in the COPD group included Rothia mucilaginosa, Porphyromonas gingivalis, Gemella sanguinis, Abiotrophia defectiva, Treponema denticola, and Tannerella forsythia. In contrast, Corynebacterium segmentosum, Streptococcus pneumoniae, Prevotella melaninogenica, and Haemophilus parainfluenzae were enriched in the Control group
Genus-level analysis revealed that COPD samples showed enrichment of pathogenic taxa, such as Streptococcus, Haemophilus, and Pseudomonas, whereas controls exhibited a more balanced distribution of genera (Fig. 1C). Notably, the top five genera in COPD accounted for 68.3% of total reads, compared to only 47.1% in controls, indicating reduced microbial evenness and dominance by a few pathogenic genera. At the species level, Haemophilus influenzae, Pseudomonas aeruginosa, and Streptococcus pneumoniae emerged as the dominant pathogens in COPD, whereas Controls harbored a wider range of low-abundance species (Fig. 1D).
LEfSe analysis of species-level differences highlighted the enrichment of key pathogenic species, Rothia mucilaginosa, Porphyromonas gingivalis, Gemella sanguinis, Abiotrophia defectiva, Treponema denticola, and Tannerella forsythia, in COPD patients compared to Controls (Fig. 1E). These taxa were significantly overrepresented in the COPD cohort (p < 0.05). Conversely, species enriched in Controls included Corynebacterium segmentosum, Streptococcus pneumoniae, Prevotella melaninogenica, and Haemophilus parainfluenzae, all of which demonstrated significantly higher relative abundances in the control group (p < 0.05).
Metagenomic analyses revealed significant alterations in functional pathways and ARGs profiles between COPD patients and healthy controls. In COPD samples, LEfSe-based Gene Ontology (LEfSeGO) analysis identified a marked enrichment of pathways associated with oxidative stress responses and adaptive microbial mechanisms. These included activities such as oxidoreductase activity and carbon–sulfur lyase activity. In contrast, pathways involved in maintaining cellular homeostasis, such as cytochrome-c oxidase activity and respirasome function, were significantly more abundant in control samples (Fig. 2A, Supplementary Data 1).Fig. 2Functional Pathway Shifts and Antibiotic Resistance in the COPD Microbiome. A LEfSeGO pathways show oxidative stress enrichment in COPD and homeostasis in controls. B LEfSeKO highlights microbial survival pathways in COPD and metabolic stability in controls. C LEfSeKEGG pathway analysis shows reduced pathways in COPD, except lysosome. D LEfSeMetaCyc pathways are consistently lower in COPD than controls. E LEfSeARGs are enriched in COPD, including resistance to tetracyclines and beta-lactams
A similar trend was observed in LEfSe-based KEGG Orthology (LEfSeKO) analysis, which showed that the COPD microbiota was enriched for pathways involved in energy-coupling factor (ECF) transport systems, iron acquisition mechanisms, and oxidative phosphorylation. Meanwhile, control samples exhibited greater abundance of pathways related to large subunit ribosomal proteins, indicative of preserved protein synthesis and translational machinery (Fig. 2B, Supplementary Data 2). At the broader pathway level, LEfSe-based KEGG Pathway analysis demonstrated that most functional pathways were significantly depleted in COPD samples compared to controls. Notably, lysosome-related pathways were among the few that were significantly enriched in COPD microbiota. In contrast, key pathways such as glycosaminoglycan degradation, bacterial secretion systems, and aminoglycoside biosynthesis were markedly reduced in COPD (Fig. 2C, Supplementary Data 3).
Further, LEfSeMetaCyc pathway analysis revealed a pronounced reduction in overall functional diversity within COPD-associated microbial communities. All identified pathways, including those involved in UDP-N-acetylglucosamine biosynthesis, peptidoglycan biosynthesis, and other central carbohydrate metabolism pathways, were significantly diminished in COPD compared to control samples (Fig. 2D; Supplementary Data 4). Finally, LEfSeARG analysis indicated an increased abundance of ARGs in COPD microbiota. Enriched resistance genes included those targeting beta-lactam antibiotics (e.g., penicillin-binding proteins PBP-1B and PBP-2X), macrolide–lincosamide-streptogramin (MLS) antibiotics (e.g., 23S rRNA methyltransferase RLMA), and fluoroquinolones (e.g., aminotransferase PATB) (Fig. 2E, Supplementary Data 5). These findings suggest that the airway microbiota in COPD is not only functionally impaired but also harbors a heightened reservoir of clinically relevant resistance determinants.
Transcriptomic profiling identified 1442 COPD-related differentially expressed genes (CDEGs), including 823 upregulated and 619 downregulated genes, defined by an p-value < 0.01 and |log2 fold change|≥ 2 (Fig. 3A, Supplementary Data 6). Functional enrichment of these CDEGs revealed significant biological processes, including “leukocyte chemotaxis,” “response to lipopolysaccharide,” and “regulation of inflammatory response” (Fig. 3B, Supplementary Data 7). KEGG pathway analysis further highlighted central pathways such as “NF-κB signaling,” “cytokine–cytokine receptor interaction,” and “complement activation,” demonstrating the overactivation of inflammatory and immune pathways in COPD (Fig. 3C, Supplementary Data 8). To explore gene co-expression patterns, WGCNA was conducted, identifying three key co-expression modules—brown, blue, and yellow—that were significantly correlated with COPD (Fig. 3D, E). The brown module was associated with inflammatory and metabolic processes, including “peptidyl-serine phosphorylation” and “regulation of mRNA metabolic process” (Fig. 3F, Supplementary Data 9), while the blue module was linked to innate immune activation, highlighting processes such as “leukocyte chemotaxis,” “response to lipopolysaccharide,” and “neutrophil extracellular trap formation” (Fig. 3G, Supplementary Data 10). The yellow module was primarily enriched for ribosomal and translational pathways, including “ribosome biogenesis” and “cytoplasmic translation,” suggesting alterations in protein synthesis in COPD (Fig. 3H, Supplementary Data 11).Fig. 3Immune and Inflammatory Gene Expression Alterations in COPD. A Volcano plot showing CDEGs in COPD. B, C GO and KEGG pathway enrichment analyses highlight pathways such as leukocyte chemotaxis, NF-κB signaling, and blood coagulation.** D** WGCNA cluster dendrogram identifies gene modules linked to COPD. E Module–trait heatmap shows strong correlations between gene modules and disease status. F–H Enrichment analysis of key modules reveals pathways WGCNA.** I** Venn diagram shows overlap between CDEGs and WGCNA hub genes. J, K Overlapping genes are enriched in processes. L, M GSEA of COPD
A subset of 159 COPD-related hub genes (CHGs) was identified as the intersection of CDEGs and WGCNA hub genes, representing a core gene set implicated in COPD pathophysiology (Fig. 3I, Supplementary Data 12). Enrichment analysis of these intersecting genes revealed significant GO terms such as “response to molecule of bacterial origin” and “leukocyte migration,” alongside KEGG pathways like “NF-κB signaling” and “TNF signaling,” reinforcing their roles in inflammation and immune regulation (Fig. 3J, K, Supplementary Data 13). GSEA further elucidated pathway dysregulation, identifying enrichment of GO terms such as “cytosolic ribosome” and “azurophil granule membrane” (Fig. 3L, Supplementary Data 14), while KEGG analysis highlighted pathways including “ferroptosis” and “antigen processing and presentation” (Fig. 3M, Supplementary Data 15).
Immune infiltration analysis demonstrated notable alterations in COPD patients compared to Controls, characterized by significantly elevated neutrophils and eosinophils, indicative of heightened innate inflammatory responses. In contrast, key adaptive immune cells, including regulatory T cells, CD4 + T cells, B cells, and plasma cells, were markedly reduced, suggesting suppressed adaptive immune regulation (Fig. 4A, B).Fig. 4Immune Cell Dysregulation and Microbial Correlations in COPD. A Immune infiltration analysis showing elevated neutrophils and eosinophils and reduced adaptive immune cells in COPD patients. B Boxplots of significantly altered immune cell subsets between COPD and controls. C Venn diagram identifying 40 COPD immune-related genes. D, E GO and KEGG enrichment analyses of CIGs, highlighting inflammatory responses, leukocyte migration, and NET formation. F Heatmap showing correlations between BIGs and microbial taxa (Prevotella melaninogenica, Streptococcus sanguinis). G Sankey diagram linking BIGs to microbial species and KEGG pathways. H Network diagram illustrating host-microbe interactions and associated pathways
To further investigate the immune dysregulation in COPD, CHGs from transcriptomic analysis were intersected with immune-related gene lists, yielding 40 COPD immune-related genes (CIGs) (Fig. 4C, Supplementary Data 16). Functional enrichment analysis of these CIGs revealed biological processes such as "regulation of inflammatory response," "response to molecules of bacterial origin," and "leukocyte migration," underscoring their involvement in innate immune activation. KEGG pathway analysis identified key pathways including "NF-kappa B signaling," "chemokine signaling," and "neutrophil extracellular trap (NET) formation" (Fig. 4D, E, Supplementary Data 17).
Correlation analysis linked specific microbial taxa to COPD immune-related genes (BIGs). Prevotella melaninogenica and Streptococcus sanguinis showed strong associations with genes involved in oxidative stress (SOD2), inflammation (SERPINE1), and neutrophil extracellular trap (NET) formation (AQP9) (Fig. 5F). Sankey diagrams mapped differentially expressed genes (e.g., AQP9, LDHA, and SOD2) to their respective microbial species and KEGG pathways, illustrating their functional connections (Fig. 4G). An integrative network diagram further visualized the intricate relationships between microbial species, host genes, and enriched KEGG pathways, highlighting robust correlations between microbial taxa, metabolic pathways, and host transcripts involved in inflammation, immune regulation, and tissue remodeling (Fig. 4H).Fig. 5Taxonomic Profiles in Severe and Mild Lung Function Clusters. A PCA of lung function indices demonstrates partial separation of SLF (n = 54) and NLF (n = 34) groups. B PCoA based on Bray–Curtis dissimilarity reveals no significant differences in overall microbiome composition between SLF and NLF. C Alpha-diversity indices, including Richness, ACE, Chao1, Shannon, Simpson, and Evenness, demonstrates no significant differences between the two groups, indicating similar taxonomic diversity. D, E Relative abundances of microbial genera and species show consistent dominance of Pseudomonas, Haemophilus, and Streptococcus in both SLF and NLF clusters.** F** Differential abundance analysis identifies taxa
Within the COPD cohort, patients were stratified into two subgroups based on lung function severe to very severe lung function impairment (SLF, n = 53) and mild to moderate lung function impairment (NLF, n = 35), as detailed in Table 2. No significant differences were observed in age, BMI, or smoking status between the SLF and NLF groups, ensuring comparability of baseline characteristics. Pulmonary function tests demonstrated significant differences between the SLF and NLF subgroups. FEV1 was significantly lower in the SLF group compared to the NLF group (p < 0.001), indicating greater airflow limitation. Similarly, FVC was significantly reduced in SLF patients (p < 0.001), and the FEV1/FVC ratio was markedly lower in the SLF group compared to the NLF group (p < 0.001), reflecting more severe obstructive impairment. Systemic inflammatory markers were also elevated in the SLF group. White blood cell counts and neutrophil counts were significantly higher in SLF patients compared to NLF patients (p < 0.001 for both), indicating a heightened inflammatory burden. No significant differences were observed in lymphocyte, monocyte, or eosinophil counts between the subgroups. These results suggest that patients with severe to very severe COPD exhibit greater airflow limitation and a more pronounced neutrophilic inflammatory profile compared to those with milder disease, highlighting a strong association between lung function impairment and systemic inflammation in COPD progression.
PCA revealed partial separation between SLF and NLF groups, reflecting the differences in spirometric measurements (Fig. 5A). Despite significant differences in lung function, a PCoA based on Bray–Curtis dissimilarity indicated no significant differences in overall microbiome composition between SLF and NLF groups (PERMANOVA: R^2^ = 0.0117, p = 0.39) (Fig. 5B). Alpha-diversity metrics, including Richness, ACE, Chao1, Shannon, Simpson, and Evenness, also showed no significant differences between the groups (Fig. 5C). Genus- and species-level analyses further demonstrated broadly similar microbial communities in both groups, with Pseudomonas, Haemophilus, and Streptococcus remaining the dominant genera across SLF and NLF patients (Fig. 5D, E). However, LEfSe analysis identified significant microbial differences between the SLF and NLF groups in COPD patients. Moraxella osloensis, Streptococcus sp. NPS 308, and Propionibacterium sp. oral taxon 193 were markedly enriched in the SLF group, indicating their potential role in exacerbating lung function decline (Fig. 5F).
Comparative metagenomic analysis between COPD patients with SLF and those with NLF revealed significant alterations in microbial functional pathways and resistome composition. LEfSeGO analysis indicated a higher enrichment of pathways associated with ribosomal structure and protein synthesis in NLF samples, including structural constituent of ribosome and magnesium ion binding. These findings suggest a relatively stable microbial metabolic profile in the NLF group. In contrast, SLF samples exhibited reduced enrichment of these biosynthetic pathways, accompanied by a mild elevation in stress response pathways such as iron binding, indicative of a more perturbed microbial environment under disease progression (Fig. 6A, Supplementary Data 18).Fig. 6Metabolic Adaptations and Antibiotic Resistance in Lung Function Clusters. A LEfSeGO pathways show higher ribosomal function in NLF, while SLF exhibits increased stress-related pathways like iron binding. B LEfSeKO pathways show ribosomal proteins enriched in NLF, while SLF has higher hydroxyethylthiazole kinase and NADPH dehydrogenase (quinone). C LEfSeKEGG pathways, including thiamine metabolism and glycolysis/gluconeogenesis, are reduced in SLF. D LEfSeMetaCyc pathways, such as pyruvate fermentation and adenosine salvage, are lower in SLF. E Resistome analysis reveals higher antibiotic resistance diversity in SLF, including resistance to macrolides, aminoglycosides, and β-lactams
LEfSeKO analysis further supported this distinction. Pathways involving large and small subunit ribosomal proteins were more prominent in NLF samples, while SLF samples showed higher representation of enzymes such as hydroxyethylthiazole kinase and NADPH dehydrogenase, which may reflect microbial adaptations to oxidative and nutritional stress (Fig. 6B, Supplementary Data 19). At the metabolic pathway level, LEfSe-based KEGG pathway analysis revealed that several key metabolic processes, thiamine metabolism, glycolysis/gluconeogenesis, and arginine and serine metabolism, were significantly diminished in SLF samples relative to NLF. Ribosome-related pathways were also markedly reduced in SLF (Fig. 6C, Supplementary Data 20).
Consistent with these trends, LEfSeMetaCyc analysis showed that microbial pathways such as pyruvate fermentation to isobutanol, adenosine salvage, and inositol phosphate metabolism were substantially lower in SLF samples compared to NLF, further indicating metabolic attenuation in advanced disease stages (Fig. 6D, Supplementary Data 21). Finally, resistome profiling (Fig. 6E, Supplementary Data 22) revealed a broader diversity and increased abundance of ARGs in the SLF group. Hierarchical clustering demonstrated significant enrichment of genes conferring resistance to MLS, aminoglycosides, and β-lactams in SLF samples, likely reflecting cumulative selective pressure from chronic antibiotic exposure. Although ARGs were also detectable in NLF samples, their abundance and diversity were substantially lower.
Transcriptomic profiling identified 32 significantly upregulated severe lung function-related differentially expressed genes (SLFGs), defined by an adjusted p-value < 0.01 and |log2 fold change|≥ 2 (Fig. 7A, Supplementary Data 23). Functional enrichment analysis of these SLFGs revealed pronounced activation of biological processes such as reactive oxygen species metabolic process, cellular oxidant detoxification, and blood microparticle formation (Fig. 7B, Supplementary Data 24). KEGG pathway analysis further highlighted key signaling cascades, including Spliceosome, Relaxin signaling pathway, and Carbohydrate digestion and absorption (Fig. 7C, Supplementary Data 25). These findings were supported by GSEA, which demonstrated significant upregulation of hallmark pathways such as cytosolic small ribosomal subunit and granulocyte activation in SLF samples (Fig. 7D), underscoring heightened immune and metabolic stress responses in advanced COPD.Fig. 7Transcriptomic Changes and Microbial Interactions in Lung Function Clusters. A Volcano plot of 32 upregulated SLFGs. B-C GO and KEGG enrichment highlight pathways related to oxidative stress and metabolic regulation. D GSEA shows significant enrichment of immune and metabolic pathways in SLF. E Immune cell proportions reveal no significant differences between SLF and NLF. F–H Correlation, Sankey, and network analyses link SLFGs with microbial taxa and KEGG pathways, emphasizing host–microbiota interactions in SLF
Immune cell deconvolution analysis indicated differences in immune cell composition between COPD patients and healthy controls; however, no statistically significant variation was observed in immune cell proportions between the SLF and NLF subgroups (Fig. 7E), suggesting that transcriptional reprogramming in severe COPD may be driven more by functional gene expression changes than by compositional shifts in immune infiltrates.
Microbiome–host correlation analysis uncovered distinct and biologically meaningful interactions between microbial species and SLFGs, offering insights into the host–microbiome interplay in severe disease (Fig. 7F). Notably, HBA2 and HBB exhibited strong positive correlations with Moraxella osloensis. A Sankey diagram further visualized these host–microbe interactions, linking SLFGs and enriched KEGG pathways, such as Glycolysis/Gluconeogenesis and Pyrimidine metabolism, to dominant microbial taxa including Haemophilus influenzae and Streptococcus gordonii (Fig. 7G). This integrative view supports the hypothesis that microbial metabolic activity contributes to epithelial inflammation, oxidative stress, and immune dysregulation in patients with severe airflow limitation. Network analysis identified Haemophilus influenzae as a central microbial node with strong associations to SLFGs involved in oxidative stress response, immune activation, and metabolic remodeling (Fig. 7H). These results suggest that specific microbial taxa may exert a direct influence on host transcriptional programs, thereby exacerbating lung function decline in advanced COPD.
To address the diagnostic and stratification challenges posed by the clinical heterogeneity of COPD, we developed two diagnostic one to distinguish COPD patients from Controls and another to stratify COPD patients into severe lung function and non-severe lung function groups.
In the validation cohort, COPD patients exhibited significant lung function impairments, including reduced FEV1, FVC, and FEV1/FVC ratio compared to controls (all p < 0.001), consistent with typical COPD pathophysiology (Table 3). Notably, there were no statistically significant differences in predicted FEV1 or predicted FVC values between the two groups (all p > 0.05), indicating comparable baseline physiological characteristics. Markers of systemic inflammation, including white blood cells, neutrophils, and monocytes, were significantly elevated in COPD patients (all p < 0.05), whereas lymphocytes and eosinophils showed no significant differences (all p > 0.05). Baseline demographic factors such as age, BMI, and smoking status were also comparable between groups (all p > 0.05), supporting the internal validity of the comparison. As summarized in Supplementary Tables 1 and 2, the clinical and demographic features of the validation cohort closely aligned with those of the discovery cohort, enhancing the reliability and generalizability of the multi-omics analyses and the robustness of the diagnostic models used to distinguish COPD from controls and stratify disease severity.Table 3Basic information of validation cohort (COPD and Controls)CharacteristicsCOPDControlsPNo4318Age, year66.186 ± 8.627965.611 ± 7.16310.804BMI, kg/m^2^22.008 ± 2.85321.845 ± 2.7510.837Current smoker18 (41.86%)8 (44.4%)0.848Lung FunctionFEV1 predicted2.5753 ± 0.535622.7739 ± 0.660580.223FEV1, L1.32 ± 0.657693.0028 ± 0.74631 < 0.001FEV1 (A/Pd),%52.167 ± 23.97108.41 ± 8.8527 < 0.001FVC predicted3.2237 ± 0.675553.3011 ± 0.818990.703FVC, L2.4628 ± 0.862583.6561 ± 0.90244 < 0.001FEV1/FVC0.52631 ± 0.139060.82154 ± 0.066609 < 0.001GOLD-stageI6 (13.95%)//II15 (34.89%)//III10 (23.25%)//IV12 (27.91%)//Inflammatory markersWBC,10^9^/L7.6347 ± 2.06085.8267 ± 1.96890.002NEU,10^9^/L5.2547 ± 1.38143.5783 ± 1.6303 < 0.001LYM,10^9^/L1.4672 ± 0.797431.69 ± 0.38480.147MONO,10^9^/L0.50209 ± 0.24640.37556 ± 0.142920.015EOS,10^9^/L0.17372 ± 0.187220.135 ± 0.106510.415FEV1 Forced expiratory volume in one second; FVC Forced vital capacity; WBC White blood cell; NEU Neutrophil; LYM Lymphocyte; MONO Monocyte; EOS Eosinophil
The first model aimed to differentiate COPD patients from Controls, and three biomarkers were selected by using LASSO and Random Forest, including GABARAP, IGSF6, MARCKS. All identified biomarkers overlapped with CHGs (Fig. 8A, Supplementary Data 26). In the training cohort (62 COPD vs. 13 controls), transcriptomic analysis revealed significantly elevated expression of GABARAP, IGSF6, and MARCKS in COPD patients compared to controls (p < 0.05). A logistic regression model based on these three genes achieved an AUC of 0.833 (95% CI: 0.739–0.926) (Fig. 8B). In the test cohort (26 COPD vs. 7 controls), transcriptomic data confirmed significantly higher expression of GABARAP, IGSF6, and MARCKS in COPD patients (p < 0.05), with the logistic regression model yielding an improved AUC of 0.885 (95% CI: 0.711–1.000) (Fig. 8C). To further validate the model’s robustness, we analyzed BALF samples from an independent validation cohort (43 COPD vs. 18 controls) using qPCR. Relative gene expression was normalized to GAPDH and calculated using the 2^−ΔΔCt^ method, with control samples set as the reference (1.0). Consistent with prior findings, GABARAP, IGSF6, and MARCKS showed significantly elevated expression in COPD patients (p < 0.05). The logistic regression model based on qPCR data achieved an AUC of 0.800 (95% CI: 0.673–0.926) (Fig. 8D). These results demonstrate that the three-gene model, based on GABARAP, IGSF6, and MARCKS, exhibits strong discriminatory power for distinguishing COPD patients from healthy individuals across both transcriptomic and qPCR datasets, confirming its diagnostic potential.Fig. 8Development and Validation of Diagnostic Models for COPD. A Venn diagram showing three biomarkers selected by LASSO and RF, all overlapping with CHGs. B–D ROC curves for the logistic regression model using three biomarkers, with AUC values of 0.833 (95% CI: 0.739–0.926) in the training cohort (B), 0.885 (95% CI: 0.711–1.000) in the internal test cohort (C), and 0.800 (95% CI: 0.673–0.926) in the validation cohort (D).** E** Venn diagram showing two biomarkers selected by LASSO and RF, all overlapping with SLFGs. F–H ROC curves for the logistic regression model using two biomarkers, with AUC values of 0.859 (95% CI: 0.771–0.947) in the training cohort (F), 0.912 (95% CI: 0.774–1.000) in the internal test cohort (G), and 0.937 (95% CI: 0.865–1.000) in the validation cohort (H)
In the validation cohort, COPD patients were stratified into SLF and NLF subgroups. Significant differences in pulmonary function were observed between FEV1, FVC, and FEV1/FVC ratio were markedly lower in the SLF subgroup compared to the NLF subgroup (all p < 0.001), indicating more severe airflow limitation in SLF patients (Table 4). In terms of systemic inflammation, SLF patients showed significantly elevated white blood cell and neutrophil counts (all p < 0.05), while lymphocytes, monocytes, and eosinophils were not significantly different between groups (all p > 0.05). Baseline clinical characteristics, including age, BMI, and smoking status, were statistically comparable (all p > 0.05), supporting the robustness of the subgroup comparisons. As shown in Supplementary Tables 3 and 4, the SLF and NLF subgroups in the validation cohort exhibited similar clinical and demographic profiles to those in the discovery cohort, reinforcing the consistency and generalizability of the stratification strategy and the diagnostic modeling.Table 4Basic information of validation cohort (SLF and NLF)CharacteristicsSLFNLFPNo2221Age, year65.864 ± 9.453166.524 ± 7.89060.805BMI, kg/m^2^21.974 ± 2.749122.043 ± 3.02570.938Current smoker9 (40.90%)9 (42.85%)0.897Lung FunctionFEV1 predicted2.6409 ± 0.543192.5067 ± 0.531950.418FEV1, L0.81864 ± 0.227581.8452 ± 0.53648 < 0.001FEV1 (A/Pd),%31.764 ± 9.532173.543 ± 13.172 < 0.001FVC predicted3.2705 ± 0.694513.1748 ± 0.668560.648FVC, L1.9414 ± 0.555833.009 ± 0.79101 < 0.001FEV1/FVC0.43931 ± 0.121430.61746 ± 0.089875 < 0.001GOLD-stageI10 (45.45%)//II12 (54.55%)//III/6 (28.57%)/IV/15 (71.43%)/Inflammatory markersWBC,10^9^/L8.36 ± 2.42056.8748 ± 1.25770.016NEU,10^9^/L5.8341 ± 1.27054.6476 ± 1.24750.004LYM,10^9^/L1.3977 ± 0.756461.54 ± 0.850680.565MONO,10^9^/L0.53591 ± 0.268460.46667 ± 0.221910.363EOS,10^9^/L0.20273 ± 0.228620.14333 ± 0.129740.304FEV1 Forced expiratory volume in one second; FVC Forced vital capacity; WBC White blood cell; NEU Neutrophil; LYM Lymphocyte; MONO Monocyte; EOS Eosinophil
Similarly, LASSO and RF regression were applied to identify biomarkers differentiating SLF from NLF clusters. Two biomarkers were selected, all of which overlapped with SLFGs, including STPG4, SPINDOC (Fig. 8E, Supplementary Data 27). In the training cohort (40 SLF vs 20 NLF), both STPG4 and SPINDOC were significantly upregulated in the SLF group compared to the NLF group (p < 0.05). A logistic regression model constructed with these two biomarkers demonstrated strong performance in distinguishing SLF from NLF patients, achieving an AUC of 0.859 (95% CI: 0.771–0.947) in the training cohort (Fig. 8F). The findings were further validated in the test cohort (13 SLF vs 7 NLF), where the same pattern of biomarker elevation was observed, and the model achieved an AUC of 0.912 (95% CI: 0.774–1.000) (Fig. 8G). To confirm the reliability and clinical utility of these biomarkers, we performed qPCR validation using BALF samples from an independent external validation cohort (22 SLF vs 21 NLF). Relative expression levels of STPG4 and SPINDOC were quantified using GAPDH as an internal reference, with expression values normalized to the NLF group. qPCR results were consistent with transcriptomic trends, with significantly higher expression of both genes in the SLF group (p < 0.05). The corresponding logistic regression model based on qPCR data yielded an AUC of 0.937 (95% CI: 0.865–1.000) (Fig. 8H), demonstrating excellent discriminatory power. These results collectively support the robustness of STPG4 and SPINDOC as predictive biomarkers for lung function stratification in COPD.
This integrative study combines metagenomic and transcriptomic analyses of BALF to elucidate the intricate host-microbe interactions driving COPD across varying severities. Our findings reveal diminished airway microbial diversity, pronounced enrichment of Moraxella osloensis and Streptococcus species, elevated antibiotic resistance, and dysregulated immune-metabolic pathways in severe COPD, complemented by high-accuracy diagnostic models for disease detection and severity stratification. These insights challenge existing paradigms and lay a robust foundation for precision medicine in COPD management.
The reduced airway microbial diversity in COPD, particularly in patients with severe lung function impairment, aligns with studies linking dysbiosis to chronic respiratory diseases [31–33]. This microbial imbalance, driven by chronic inflammation and stressors like cigarette smoke, fosters a niche for pathogenic bacteria [34]. Our identification of Moraxella osloensis as a dominant driver in severe COPD redefines the traditional focus on Haemophilus influenzae and Pseudomonas aeruginosa [14, 35, 36]. Unlike Pseudomonas, often tied to acute exacerbations, Moraxella osloensis likely sustains chronic inflammation via lipopolysaccharide-mediated activation of Toll-like receptor 4 (TLR4), triggering downstream MAPK and IRF3 pathways, which amplify neutrophil extracellular trap formation, a process where neutrophils release DNA to trap pathogens, and NF-κB signaling [31, 37–40]. In vitro studies confirm Moraxella species induce robust IL-8 production in airway epithelial cells via TLR4, potentially synergizing with Streptococcus to exacerbate tissue damage [37, 38, 41]. The lack of broad microbial community differences between SLF and NLF subgroups suggests specific pathogens drive disease progression, challenging assumptions of uniform dysbiosis. This finding, contrasting with studies emphasizing Pseudomonas dominance, may reflect our cohort’s stable COPD profile, warranting further exploration of patient-specific factors like smoking history or geographic variation [36].
Functionally, the SLF microbiota exhibits enriched oxidative stress and antibiotic resistance pathways. ARGs, targeting beta-lactams, macrolides, and fluoroquinolones, reflect chronic antimicrobial exposure, a hallmark of advanced COPD [42, 43]. For instance, PBP-1B, linked to persistent Streptococcus infections, underscores challenges in antibiotic efficacy [44–47]. Depleted metabolic pathways, such as glycolysis and thiamine metabolism, indicate a microbial community adapted to a nutrient-scarce, inflammatory environment, potentially impairing epithelial homeostasis and exacerbating airway remodeling, as seen in cystic fibrosis [48, 49]. In contrast, milder COPD cases show preserved ribosomal pathways, suggesting metabolic stability that may slow progression [50]. These findings challenge reports of consistent microbial functionality across COPD severities and highlight microbiome-targeted therapies, such as narrow-spectrum antimicrobials or probiotics, to restore balance and curb resistance [51–53].
Transcriptomic profiling reveals a suite of COPD-related differentially expressed genes, with a subset specific to severe cases, enriched in NF-κB signaling, cytokine interactions, and NET formation, underscoring neutrophilic inflammation and oxidative stress [54–56]. Strong correlations between Moraxella osloensis and Streptococcus species with genes like HBB, FCN3, and AQP9 indicate direct microbial modulation. HBB, associated with Streptococcus, supports redox homeostasis under oxidative stress, a mechanism implicated in COPD airway damage [57]. FCN3, linked to Streptococcus gordonii, encodes a ficolin that enhances complement activation, amplifying inflammation [58]. AQP9, tied to NET formation, aligns with reports of excessive NETosis driving epithelial damage and biofilm formation [59, 60]. Network analyses position Haemophilus influenzae and Streptococcus as central nodes in immune and metabolic dysregulation, consistent with chronic lung diseases [45, 57, 61]. These findings extend prior omics study, but contrast with Pseudomonas-centric reports, possibly due to differences in exacerbation frequency or sampling methods, necessitating further validation [11, 62]. Our diagnostic models, leveraging GABARAP, IGSF6, and MARCKS for COPD detection, and STPG4 and SPINDOC for severity stratification, achieve superior accuracy. GABARAP, a mediator of autophagosome formation, reflects cellular stress responses [63, 64]. MARCKS, involved in mucus secretion and epithelial repair, supports airway hyperresponsiveness [65]. STPG4 and SPINDOC, enriched in glycolysis and pyrimidine metabolism, indicate metabolic reprogramming in severe COPD [66, 67]. Validated via qPCR, these biomarkers can integrate with FEV1/FVC ratios in emergency settings for rapid stratification, offering cost-effective precision. Implementation via portable qPCR devices could enhance accessibility, though cost and scalability require evaluation.
Therapeutically, based on the comparison between SLF and NLF groups, the marked enrichment of Moraxella osloensis underscores its potential role in COPD lung function decline, supporting the development of targeted antimicrobials or vaccines. Elevated ARGs necessitate alternatives like NF-κB inhibitors or DNase I to disrupt NETs, effective in COPD models [54, 56]. Depleted microbial metabolism highlights prebiotic interventions to restore homeostasis, supported by gut-lung axis studies [68, 69]. The gut-lung axis, linking airway and intestinal microbiota, may amplify systemic inflammation, with gut dysbiosis noted in COPD and cardiovascular diseases [52]. Integrated therapies combining nebulized probiotics with anti-inflammatory agents could address both dysbiosis and immune dysregulation, though delivery challenges, like mucus barriers, require optimization [70].
While our study provides valuable insights into COPD pathophysiology, several limitations should be acknowledged. Establishing causality is limited by the airway microbiome’s complexity and the lack of robust experimental models for COPD. Unlike cystic fibrosis models [71], current animal models fail to fully replicate the chronic, heterogeneous pathology of COPD, including recurrent infections and airway remodeling. Cell models, such as air–liquid interface cultures, cannot capture complex microbial-immune interactions. Future advanced in vitro models, like airway organoids co-cultured with patient-derived microbiota [72], could validate causal mechanisms and test therapeutic targets, such as NET inhibitors or microbiome-modulating therapies. Moreover, the functional roles of specific taxa enriched in SLF, such as Moraxella osloensis, require further investigation. Experimental studies, including co-culture systems and in vivo models, are needed to determine their contributions to inflammation, immune modulation, and lung function decline. While our sample size of 131 COPD patients and 38 Controls is substantial for BALF-based multi-omics studies, the heterogeneity of COPD may warrant larger cohorts to further validate our findings. Our study focused on disease-associated microbial signatures; however, deeper sequencing could improve detection of low-abundance taxa and strain-level resolution in future work.
To address this critical limitation, longitudinal cohort studies should be conducted to track key microbial taxa, such as Moraxella osloensis and Streptococcus sp., alongside immune markers, and lung function metrics across multiple time points. This ongoing effort aims to correlate these changes with exacerbations and disease progression, providing a clearer picture of temporal relationships. Looking ahead, future research should focus on large-scale longitudinal studies to map the dynamic shifts in microbial-immune profiles and their impact on lung function decline, ensuring robust integration of multi-omics data with clinical outcomes to validate our diagnostic biomarkers. While qPCR validation aligns with our current translational objectives, complementary protein-level analyses such as Western blot will be valuable in future studies to further substantiate our transcriptomic findings. Future studies could also incorporate flow cytometry or single-cell RNA-seq to achieve higher-resolution immune profiling and complement deconvolution-based analyses. Additionally, establishing standardized protocols for serial BALF sampling will be essential to maintain data consistency across time points. This study will serve as a vital multi-omics framework, guiding hypothesis-driven longitudinal research to uncover the temporal mechanisms driving COPD progression and informing the development of targeted therapies, such as microbiome-modulating interventions or anti-inflammatory treatments.
This study provides a comprehensive multi-omics framework that elucidates the intricate interplay between airway microbial dysbiosis, host immune responses, and lung function decline in COPD. By integrating metagenomic and transcriptomic analyses of bronchoalveolar lavage fluid across well-characterized discovery and validation cohorts, we identified key microbial taxa, such as Moraxella osloensis and Streptococcus sp., enriched in severe lung function impairment patients. These taxa correlate with upregulated inflammatory pathways, including NF-κB signaling and neutrophil extracellular trap formation, which drive oxidative stress and airway damage. Notably, the enrichment of antibiotic resistance genes in severe lung function patients underscores the challenge of antimicrobial resistance in COPD management, highlighting the need for novel therapeutic strategies, such as microbiome-targeted interventions or NET inhibitors. Our diagnostic models, leveraging immune biomarkers, achieved robust performance, offering a pathway for precise patient stratification and personalized treatment. Despite limitations, including the challenge of establishing causality, our findings lay a foundation for future research to validate these insights and develop targeted therapies for COPD.
The Major clinical research project of Guangzhou Medical University (GMUCR2024-02009), Guangdong Provincial Clinical Research Center for Laboratory Medicine (2023B110008).
Additional file 1. Additional file 2.