From November 1, 2018, to May 31, 2023, we conducted a prospective screening of all patients who visited the outpatient clinic of Gangnam Severance Hospital. COPD was defined according to the recommendations of the Global Initiative for Chronic Obstructive Lung Disease (GOLD) 2018 as patients presenting symptoms of COPD or those with risk factors and a post-bronchodilator FEV₁/FVC ratio of less than 0.7. Patients were included if they met the diagnostic criteria for COPD. The exclusion criteria were as follows: patients with asthma, those unable to provide sputum samples, and those who had already received antibiotic treatment within the previous month. Stable Disease (SD) was defined as the absence of ECOPD episodes in the past year. ECOPD was defined as an acute worsening of respiratory symptoms requiring additional therapy. Sputum samples were collected from patients with SD and ECOPD on their enrollment dates (SD and AE-1, respectively). Additionally, patients with ECOPD provided a second sputum sample (AE-2) 2 weeks after initiating treatment with antibiotics, corticosteroids, or both.

Qualified spontaneous sputum samples, based on epithelial cell counts, were collected within 24 h of enrollment before any systemic therapy (SD and AE-1). We considered samples with epithelial cells <10/LPF in sputum Gram stain grading as uncontaminated by oropharyngeal secretions. AE-2 sputum samples were collected 2 weeks after medical treatment with antibiotics, corticosteroids, or both, either at outpatient clinics or upon hospital admission for ECOPD.

DNA was extracted using the DNeasy PowerSoil Pro Kit (Qiagen, Hilden, Germany) following the manufacturer’s protocol. Extracted DNA was quantified using Quant-IT PicoGreen (Invitrogen).

Sequencing libraries were prepared following the PacBio Amplicon Template Preparation and Sequencing protocols, targeting the 27F and 1492R regions. Input gDNA (2 ng) was polymerase chain reaction (PCR)-amplified with 10× LA PCR Buffer II (Mg²⁺-free), 2.5 mM dNTP mix, 2.5 mM MgCl₂, 500 nM each of forward and reverse PCR primers, and 5 U of TaKaRa LA Taq (Takara, Kusatsu, Japan). The PCR conditions were initial heat activation at 94°C for 5 min, followed by 25 cycles of 30 s at 94°C, 30 s at 53°C, and 90 s at 72°C, with a final extension at 72°C for 5 min. The primer pairs with asymmetric barcoded adapters used for the amplifications were: 27F-F: 5′- AGRGTTYGATYMTGGCTCAG −3′ and 1,492-R: 5′- RGYTACCTTGTTACGACTT -3′. PCR products were purified using SMRTbell cleanup beads, quantified using Quant-IT PicoGreen (Invitrogen), and qualified using the TapeStation D5000 Screen Tape (Agilent Technologies, Waldbronn, Germany). For PacBio Sequel IIe sequencing, 500 ng of pooled amplicon DNA was used for library preparation. A total of 10 μL of the library was prepared using the PacBio SMRTbell prep kit 3.0. SMRTbell templates were annealed using the Sequel II Bind Kit 3.1 and Int Ctrl 3.1. Sequencing was performed using the Sequel II Sequencing Kit 2.0 and SMRT cells 8 M Tray. Data acquisition involved 10-h movies captured for each SMRT cell using the PacBio Sequel IIe (Pacific Biosciences) sequencing platform, conducted by Macrogen (Seoul, Korea). Subsequent steps followed the PacBio Sample Net-Shared Protocol, which is available on the PacBio website.¹

We conducted a power analysis to evaluate the required sample size. With an effect size of 0.5, an alpha error of 0.05, and a power of 0.80, we set a target of including at least 35 participants.

Categorical variables were presented as frequencies (percentages). Continuous variables were presented as mean ± standard deviation (SD) for normally distributed data and as median (interquartile range [IQR]) for non-normally distributed data. Normality assumptions for continuous variables were confirmed using the Shapiro–Wilk test. Baseline characteristics among three or more distinct patient groups were compared using one-way analysis of variance (ANOVA) or Kruskal–Wallis tests. Tukey’s Honest Significant Difference test served was used post-hoc for ANOVA, while Dunn’s test followed the Kruskal–Wallis test. Fisher’s exact test was used to analyze categorical data.

Analyses were performed using R (R Foundation for Statistical Computing, Vienna, Austria)². Phyloseq, microbiomeMarker, and NetCoMi packages were used for ASV sequencing, linear discriminant analysis (LDA) effect size (LEfSe), and network analyses, respectively. The MaAsLin2 package was used to apply a mixed-effects linear model for paired samples to conduct a within-patient analysis. A Spearman’s correlation plot was constructed using the Corrplot package. Permutational multivariate analysis of variance (PERMANOVA) was performed using the vegan package.

Thirty-eight patients were enrolled in the study, with 17 assigned to the SD group and 21 to the AE group (Table 1). Seventeen sputum samples were obtained once from the SD group and 38 samples were obtained twice from the AE group, taken at two different time points. Four AE-2 samples were excluded from longitudinal analysis in AE due to follow-up loss.

Although no significant differences were observed in age, sex, height, history of smoking, sputum collection season, underlying diseases, and initial chest radiographic findings between groups, weight and body mass index (BMI) were significantly lower in the AE than in SD (59.7 ± 9.6 kg vs. 67.3 ± 11.2 kg, p = 0.030; and 21.9 kg/m² [IQR 20.6–23.0] vs. 24.0 kg/m2 [22.1–25.3], p = 0.025, respectively, Table 1). The AE group showed a higher use of antibiotics (66.7%) compared to the SD group (17.6%), with a p-value of 0.007. Similarly, steroid use was more frequent in the AE group (95.2%) than in the SD group (11.8%), with a p-value of <0.001.

AE patients had more severe COPD grades, with a higher proportion in grade 3 (52.4%) and grade 4 (23.8%), while SD patients were mostly in grade 2 (76.5%), showing a significant difference in COPD grade distribution (p = 0.003).

Additionally, in the analysis for longitudinal assessment within the AE group, there were no significant differences in baseline characteristics between the included and excluded patients (Supplementary Table 1).

Pulmonary function tests revealed that forced expiratory volume in the first second (FEV₁) was significantly lower in AE (1.1 liters [0.8–1.5] vs. 1.7 liters [1.5–2.1], p = 0.001) as was z-score of FEV₁ by Global Lung Function Initiative-2012 (GLI-2012) reference (−3.3 [−3.8–-2.7] vs. −2.3 [−2.6–−2.0], p < 0.001). Forced expiratory flow (FEF, FEF_25–75%, FEF₂₅, FEF₅₀, and FEF₇₅) and peak expiratory flow were lower in the AE group than in the SD group (Table 2).

In 55 sputum samples, 3,470 amplicon sequence variants were obtained, with a median feature sequence frequency of 36,284.

In alpha diversity, the observed feature count was significantly higher in SD (133.0 [98.0–145.0]) compared to AE-1 (88.0 [72.0–125.0], p = 0.036) and AE-2 (92.0 [38.0–110.0], p = 0.002; Figure 1A). Shannon’s index was also significantly higher in SD (3.8 [3.3; 4.0]) than in AE-2 (3.1 [2.2; 3.9], p = 0.034; Figure 1B). Similarly, Fisher’s alpha index was higher in SD (17.1 [12.7; 19.6]) compared to AE-1 (10.9 [8.5; 16.1], p = 0.031) and AE-2 (11.2 [4.2; 14.0], p = 0.002; Figure 1C). The under-73 group also showed significantly higher diversity metrics compared to the 73-and-older group, with an observed feature count of 137.3 ± 32.8 vs. 82.9 ± 32.7, Shannon’s index of 3.9 [3.5–4.2] vs. 3.0 [2.8–3.6], and Fisher’s alpha index of 18.2 ± 4.9 vs. 10.4 ± 4.6 (all p < 0.001; Supplementary Figure 1A). After adjusting for age through stratified analysis, the SD group showed higher alpha diversity indices compared to the AE-1 group (Supplementary Figure 2A). In contrast, the alpha diversity indices did not show any significant differences based on smoking status, COPD severity, or the season in which the sample was collected (Supplementary Figures 1C–E). In stratified analysis conducted based on smoking status and sampling season, the SD group showed significantly higher alpha diversity compared to the AE group (Supplementary Figures 2C,D).

In the longitudinal comparison of AE samples, observed feature counts, Shannon’s index, and Fisher’s alpha index showed a decrease after 2 weeks (AE-1 vs. AE-2; 104.0 [74.0–136.0] vs. 92.0 [38.0–110.0], p = 0.011; 3.5 [3.0–3.9] vs. 3.1 [2.2–3.9], p = 0.023; and 13.2 [9.2–18.1] vs. 11.2 [4.2–14.0], p = 0.006, respectively, all by paired Wilcoxon test; Figures 1A–C). Antibiotic use during ECOPD did not result in a significant decrease in diversity (Supplementary Figure 1B). Additionally, AE-1 showed higher alpha diversity compared to AE-2 in stratified analysis by antibiotic use (Supplementary Figure 2B).

In the beta diversity analysis, un-weighted UniFrac showed a significant difference between the SD and AE-1 groups (R² = 0.031, p = 0.041 by PERMANOVA; Figure 1D). However, no significant differences were observed in weighted UniFrac and Bray–Curtis dissimilarity between the SD and AE-1 groups (R² = 0.024, p = 0.724 and R² = 0.025, p = 0.841, respectively; Figures 1E,F). In the longitudinal comparison of AE samples, the centroid distance of Bray–Curtis dissimilarity was significantly higher in AE-2 compared to AE-1 (0.656 [0.644–0.674] vs. 0.669 [0.653–0.691], p = 0.006 by paired Wilcoxon test; Figures 1G,H). Bray–Curtis dissimilarity, unweighted UniFrac, and weighted UniFrac analyses for beta diversity also showed significant differences according to COPD stage (R² = 0.029, p = 0.001 by PERMANOVA; R² = 0.030, p = 0.001; R² = 0.033, p = 0.006; Figures 1I–K, respectively).

Supplementary Figures 3, 4 illustrate the overall distribution of the microbiome. The predominant phyla in AE-1 and SD were Bacillota (77.5% [57.1–82.7] and 59.2% [40.2–75.3], p = 0.107, respectively), Pseudomonadota (7.4% [4.7–27.4] and 26.8% [16.4–45.6], p = 0.101, respectively), and Actinomycetota (3.5% [2.2–7.3] and 6.4% [2.1–9.4], p = 0.322, respectively). Fusobacteriota was significantly reduced in the AE-1 compared to the SD (0.2% [0.1–1.0] vs. 1.2% [0.4–1.7], p = 0.016).

At the genus level, significant differences were observed for Neisseria (AE-1 vs. SD; 2.6% [0.1–7.4] vs. 17.1% [3.6–38.4], p = 0.022), Fusobacterium (0.2% [0.1–0.8] vs. 0.7% [0.2–1.1], p = 0.049), Leptotrichia (0.1% [0.0–0.2] vs. 0.3% [0.1–0.6], p = 0.038), Bergeyella (0.0% [0.0–0.0] vs. 0.1% [0.0–0.2], p < 0.001), Peptostreptococcus (0.0% [0.0–0.0] vs. 0.0% [0.0–0.2], p = 0.041), Oribacterium (0.0% [0.0–0.1] vs. 0.2% [0.0–0.5], p = 0.048). At the species level, significant differences were found in Neisseria perflava (0.6% [0.0–2.4] vs. 6.1% [1.6–37.9], p = 0.011) and Haemophilus parainfluenzae (0.8% [0.0–2.1] vs. 6.0% [4.4–7.9], p = 0.011).

Among these genera, the proportions of Fusobacterium (under 73 vs. 73 and older: 0.7% [0.2–1.4] vs. 0.2% [0.0–0.7], p = 0.031), Leptotrichia (0.2% [0.0–0.6] vs. 0.0% [0.0–0.2], p = 0.019), and Oribacterium (0.1% [0.0–0.5] vs. 0.0% [0.0–0.1], p = 0.012) were significantly higher in patients under 73 years compared to those aged 73 and older (Supplementary Figure 5A). Similarly, among COPD GOLD grades, the proportions of Fusobacterium (grades 1–2 vs. 3–4: 1.0% [0.2–1.7] vs. 0.2% [0.0–0.5], p = 0.001), Bergeyella (0.1% [0.0–0.1] vs. 0.0% [0.0–0.0], p < 0.001), Peptostreptococcus (0.0% [0.0–0.1] vs. 0.0% [0.0–0.0], p = 0.006), and Oribacterium (0.1% [0.0–0.4] vs. 0.0% [0.0–0.2], p = 0.046) were significantly higher in patients with GOLD grades 1–2 compared to those with grades 3–4 (Supplementary Figure 5B).

LEfSe analysis revealed that the SD group showed significantly higher abundances of Neisseria (LDA 4.996, adjusted P [adj.P] = 0.021), Fusobacterium (LDA 3.688, adj.p = 0.047), and Peptostreptococcus (LDA 3.379, adj.p = 0.039) at the genus level compared to the AE-1 group. At the species level, N. perflava (LDA 5.074, adj.p = 0.010) and H. parainfluenzae (LDA 4.467, adj.p = 0.011) were significantly more abundant in the SD group compared to the AE-1 group (Figure 2; Supplementary Figure 6). No genera or species in the AE-1 group showed significantly higher abundance compared to the SD group (Figure 2).

Network analysis revealed that the SD group had a more complex and stronger microbial network than the AE-1 group (0.155 vs. 0.108 for edge density and 0.034 vs. 0.023 for natural connectivity, respectively; Figure 3A). Hub genera identified based on empirical quantiles of centralities included Haemophilus, Granulicatella, Neisseria, Lactobacillus, and Butyrivibrio in the SD group and Streptococcus, Gemella, Actinomyces, Klebsiella, and Staphylococcus in the AE-1 group.

During the recovery process from ECOPD, the AE-2 group showed a significant increase in the abundance of the genera Filifactor (coefficient 1.042, p < 0.001) and Oribacterium (coefficient 0.811, p = 0.008) compared to the AE-1 group, as analyzed by a mixed-effects linear model for longitudinal analysis (Figure 4). In contrast, Prevotellamassilia (coefficient − 1.154, p = 0.003) and Peptostreptococcus (coefficient − 2.438, p = 0.005) showed a significant decrease in AE-2. At the species level, the AE-2 group showed a significant increase in the abundance of the species C. sputigena (coefficient 0.821, p < 0.001) and L. trevisanii (coefficient 0.063, p < 0.001) compared to the AE-1 group. In contrast, A. oris (coefficient − 0.087, p < 0.001) and L. massiliensis (coefficient − 0.422, p < 0.001) showed a significant decrease in AE-2.

Network analysis revealed that AE-2 had a more complex microbial network than AE-1 (natural connectivity: 0.031 vs. 0.022; Figure 3B). The hub genera identified in the AE-1 group included Abiotrophia, Gemella, Neisseria, Streptococcus, and Staphylococcus. In contrast, the AE-2 group featured hubs such as Oribacterium, Stomatobaculum, Fusobacterium, Prevotella, and Parvimonas.

Spearman correlation analysis revealed that alpha diversity was significantly negatively correlated with age (Shannon’s index: rho = −0.413, p = 0.002) and COPD stage (rho = −0.414, p = 0.002; Figure 5A; Supplementary Table 1). Conversely, alpha diversity positively correlated with serum eosinophil levels (rho = 0.413, p = 0.003) and the z-score of FEV₁ (rho = 0.416, p = 0.002). Fusobacterium, Peptostreptococcus, Oribacterium, and Bergeyella were positively correlated with FEV₁ and alpha diversity indices (Figure 5; Supplementary Table 2). Supplementary Figure 7 and Supplementary Table 3 describe additional correlations with other species.