Consistent processing of published and new microbiome datasets included in this study

To identify reproducible associations between airway microbiota and COPD, we performed 16S rRNA gene sequencing of sputum samples of 187 patients with COPD and HC recruited from the EndAECOPD cohort (see Supplementary Table 4 for subject characteristics), and analyzed them in the context of 1,555 additional sputum samples from four publicly available and geographically diverse datasets. In total, two case-control studies of stable COPD versus HC, two studies of COPD stability versus exacerbation, and one study with exacerbation-recovery time series data were included (Table 1 [6, 20,21,22]). All of these datasets had considerable sample size. Except for SRP066375 [20], other datasets including COPDMAP [6], BCCS & BCES [21], and BEAT-COPD [22] were published with detailed study design and subject demographics. Although unpublished, SRP066375 showed comparable quality control statistics with published datasets (Supplementary Table 5), indicative of its reasonable data quality for downstream analysis.

Table 1 Characteristics of the large-scale COPD microbiome datasets included in this study

All datasets were sequenced at high depth (mainly on Miseq platform) with similar primers targeting V3-V5 hypervariable region of the 16S rRNA gene (Table 1) and processed separately using DADA2 algorithm [14] for data cleaning and taxonomic profiling. Despite this, hypervariable region and sequencing platform potentially lead to ASV-level heterogeneity [23]. Downstream microbiome analysis was therefore performed at genus level according to the MiBioGen consortium initiative [24] and previous experiences [23, 25], which sacrificed taxa classification resolution in exchange for less data heterogeneity.

Microbial co-occurrence network in COPD stability reveals reduced negative interactions compared to HC

We first performed health-COPD comparative analysis using data of the two Asian cohorts: EndAECOPD and SRP066375 (Table 1). Beta diversity analysis demonstrated a clear separation of samples between COPD and HC, and alpha diversity measured by Shannon index was reduced in COPD samples in both cohorts (Fig. 1A, B). We next applied weighted co-occurrence analysis with an ensemble of similarity measures [16, 17] to characterize potential interactions between different taxa (Fig. 1C, D). Network parameters reflecting node (microbe) importance including degree, betweennes centrality, and stress centrality showed substantial inconsistency across cohorts, which were significantly different between COPD and HC in EndAECOPD but not in SRP066375 (Fig. 1E). Yet when focusing on the sign of edge (microbial interaction), we found that patients with COPD in both cohorts had significantly reduced percent of negative interactions, i.e. total number of negative interactions as a proportion of total interactions, in their microbial network compared with HC (Fig. 1E-G). This trend remained significant when microbial networks were reconstructed separately by using either GBLM, Spearman, or Pearson correlation that determines the sign of edge weights in the ensemble network inference method in both cohorts (see Methods and Supplementary Table 6). To validate geographical reproducibility of this finding, we analyzed data of Haldar et al [26] that compared COPD (N = 218) versus HC (N = 124) sputum microbiomes separately retrieved from two UK cohorts. Again we observed significantly lower percent of negative interactions in microbial network between COPD versus HC (47.0%[616/1,311] versus 34.4%[188/547], P < 0.0001, χ2 test) even though COPD samples had significantly higher alpha diversity indices (Supplementary Fig. 2). Further appraisal of the effects of edge weights (interaction strength) on the observed change of negative interactions showed that the absolute weights of negative interactions were also decreased in COPD versus HC in all above data albeit not reaching statistical significance in SRP066375 (Supplementary Fig. 3). Collectively, we demonstrate that bacterial network disturbance characterized by reduced antagonistic interactions, rather than diversity loss, is a reproducible feature of COPD dysbiosis.

Fig. 1

Co-occurrence network analysis of sputum microbiome between COPD and HC. A, B. Comparison of beta diversity (weighted UniFrac distance) and alpha diversity (Shannon index, Mann–Whitney U-test) between COPD and HC in the EndAECOPD cohort (n = 187) (A) and the SRP066375 cohort (n = 190) (B). C, D. Co-occurrence network maps of COPD and HC illustrating identified microbial interactions in EndAECOPD (C) and SRP066375 (D). Each node represents a microbe and microbial interactions are illustrated by connecting edges. The number of interactions for a microbe is calculated as degree. Node size is proportional to the number of interactions for each microbe. Selected taxa of clinical relevance are noted with white color. E. Comparison of key network parameters illustrating node and edge characteristics respectively between COPD and HC in both cohorts. Degree, betweenness centrality, and stress centrality were compared between COPD and HC by Mann–Whitney U-test. Percent of negative interactions, i.e. total number of negative interactions as a proportion of total interactions, were compared by χ2 test. F, G. Visualization of the opposing network by displaying positive and negative interactions between the most abundant taxa in COPD and HC in EndAECOPD (F) and SRP066375 cohort (G). Edge transparency is proportional to interaction strength between microbes in terms of edge weights. Negative interactions are classified if the sign of the edge weights is negative, and vice versa. # denotes paraphyletic group; • denotes unclassified genus; •• denotes classified but unnamed genus

Interactome analysis identifies a set of taxa whose antagonistic networks rather than abundances confer consistent changes in COPD

In general, differential abundance analysis is used to identify disease-associated microbial biomarkers, yet is sensitive to both technical artifacts and individual variabilities in environmental exposure and lifestyle [27]. This has led to the identification of varying biomarker species in different discovery cohorts, as suggested by EndAECOPD and SRP066375 datasets that have only a few overlapping taxonomic biomarkers (Streptococcus, Catonella, Dialister, Pavimonas, and Fusobacterium) (Fig. 2A, B). In contrast, through network analysis, we identified 20 taxa that showed consistently altered degree of negative interactions (percent change of > 5%) between COPD versus HC in the two cohorts (Fig. 2C), 10 of which (Haemophilus, Campylobacter, Gemellaceae, Fusobacterium, Leptotrichia, Atopobium, TM7.3, Oribacterium, Actinomyces, and Rothia) were also recovered in the aforementioned UK study by Haldar et al [26] (Supplementary Fig. 2). Most of these taxa (18/20) showed fewer negative interactions with other taxa in COPD, corresponding to the reduced antagonistic interactions of the entire network.

Fig. 2

Network analysis of biomarker taxa in patients with COPD. A, B Enriched and depleted taxa in patients with COPD identified by Linear discriminant analysis effect size (LEfSe) analysis in EndAECOPD (A) and SRP066375 cohort (B). Genera significantly enriched or depleted in COPD in both cohorts were noted with red or green, respectively. (C) Venn diagram and barplot showing overlapped genera with altered degree of negative interactions (percent change of > 5%) between COPD versus HC in EndAECOPD and SRP066375 cohort. Genera that were further recovered in data of Haldar et al. were noted with red. (D-G). Visualization of the Campylobacter and Haemophilus interactomes between COPD and HC in EndAECOPD (D, F) and SRP066375 cohort (E, G). # denotes paraphyletic group; • denotes unclassified genus; •• denotes classified but unnamed genus

Campylobacter, recently found to be a strongest microbial biomarker of eosinophilic COPD, was the top contributor to the decreased network antagonism in COPD, while Haemophilus, an established neutrophilia-associated pathogen, acted in the opposite way [3] (Fig. 2C). Changes in relative abundance of these two genera in COPD versus HC were neither observed in most of previous microbiome studies (Supplementary Table 1) nor in EndAECOPD and SRP066375 cohorts (Fig. 2A, B; Supplementary Fig. 4). In comparison, our network analysis identified distinct Campylobacter and Haemophilus interactomes regarding proportion of negative interactions between COPD and HC in both cohorts (Fig. 2D-G), suggesting the superiority of interactomics over differential abundance analysis in the identification of disease-associated biomarkers.

Airway microbial network displays reduced negative interactions during COPD exacerbation, which recovers following treatment

To clarify how bacterial interactomes would change during the course of COPD exacerbation, we performed stability-exacerbation comparative analysis using data of COPDMAP and BCCS & BCES cohorts (Table 1). Beta diversity analysis clearly distinguished samples between stability and exacerbation in the two cohorts, however, the change in Shannon index was incongruent (Fig. 3A, B). In contrast, we observed significantly lower degree of antagonistic interactions in microbial network during exacerbations of COPD in both COPDMAP (P = 0.012) and BCCS & BCES (P = 0.011) (Fig. 3C, D). We further evaluated the dynamic change of COPD interactomes at four clinical timepoints: stability, exacerbation, post-therapy (2 weeks after treatment), and recovery (6 weeks post-exacerbation) in the BEAT-COPD cohort. Longitudinal analysis revealed the deviation of exacerbation microbiomes in beta-diversity to other timepoints (Fig. 3E, F), demonstrating the presence of dysbiosis at exacerbation. However, such dysbiosis was not captured by alpha diversity analysis (Fig. 3G). By contrast, co-occurrence network analysis highlighted major changes in interactomes, with a decrease in percent of negative microbial interactions during exacerbations compared with baseline (P = 0.001), post-therapy (P = 0.002), or at recovery (P = 0.021) (Fig. 3H). Further breakdown of interactomes identified eight taxa (Veillonellaceae, Streptococcus, Actinobacillus, Porphyromonas, Prevotella[paraphyletic group], Lactobacillus, Atopobium, and Leptotrichia) with decreased negative interactions during exacerbation across COPDMAP, BCCS & BCES, and BEAT-COPD (Fig. 3I). These findings highlight microbial interactome disturbance rather than diversity loss as a reproducible feature of COPD exacerbation.

Fig. 3

Longitudinal analysis of the COPD interactome during exacerbations. (A, B). Comparison of beta diversity (weighted UniFrac distance) and alpha diversity (Shannon index, Mann–Whitney U-test) between stability and exacerbation samples in COPDMAP (n = 715) (A) and BCCS & BCES cohort (n = 174) (B). (C, D). Visualization of the interactomes in terms of positive and negative interactions between the most abundant taxa at stability and during exacerbation in COPDMAP (C) and BCCS & BCES (D). # denotes paraphyletic group; • denotes unclassified genus; •• denotes classified but unnamed genus. (E, F). Microbial clustering (E) and one-way ANOVA analysis of weighted UniFrac distance for beta diversity (F) between stability, exacerbation, post-therapy, and recovery in the BEAT-COPD cohort (n = 476). *** P < 0.0001 (Exacerbation-stability distance versus Stability-stability distance). (G). One-way ANOVA analysis of alpha diversity (Shannon index) between stability, exacerbation, post therapy, and recovery in BEAT-COPD cohort. (H). Visualization of the interactomes in terms of positive and negative interactions between the most abundant taxa across timepoints in BEAT-COPD cohort. (I). Venn diagram and barplot showing overlapped genera with altered degree of negative interactions (percent change of > 5%) between exacerbation versus stability across the three cohorts

Airway microbial interactome is associated with key clinical determinants of COPD

Through unsupervised spectral clustering on the baseline sputum microbial network of COPD patients in EndAECOPD cohort, we identified two network clusters (Fig. 4A), with a cluster robustness of 99%. The cluster 2 (n = 18) had significantly lower proportion of negative bacterial interactions compared with the cluster 1 (n = 125) (Fig. 4B). Differences in microbial compositions were also evident between the two clusters as shown in Fig. 4C-E. Demographics including age, sex, smoking history, and medication use were similar between the two clusters except that patients in cluster 2 had lower body mass index (P = 0.035) (Supplementary Table 7). In assessment of COPD-related parameters, patients in cluster 2 were found to have lower lung function (FEV1% predicted) (P = 0.040) and higher mMRC score (P = 0.007) (Fig. 4F, G). Further evaluation of airway inflammatory status by profiling sputum cells and a panel of 19 cytokines revealed that patients in cluster 2 had exaggerated neutrophilic inflammation in terms of sputum total cells, neutrophil counts, interleukin (IL)-1β and IL-8 levels (Fig. 4H-L). During the 12-month follow-up, 71%[12/17] of patients in cluster 2 experienced an exacerbation while only 38%[46/120] of patients in cluster 1 had an exacerbation (P = 0.012). Patients in cluster 2 experienced significantly more moderate-to-severe exacerbations (P = 0.004) (Fig. 4M) and had higher proportion of frequent exacerbators (35%[6/17] versus 11%[13/120], P = 0.015). A life-table analysis showed that the exacerbation risk was significantly higher among patients in cluster 2 (Log-rank test, P = 0.004; Wilcoxon signed rank test, P = 0.006) (Fig. 4N).

Fig. 4

Clinical determinants correlate with microbial network. (A) Unsupervised stratification of microbial networks with spectral clustering for the 143 COPD patients in EndAECOPD cohort into two clusters (cluster 1, n = 125; cluster 2, n = 18). (B). Visualization of positive and negative interactions between the most abundant taxa in each cluster. Cluster 2 had significantly lower proportion of negative bacterial interactions compared with cluster 1 (P < 0.0001, χ2 test). # denotes paraphyletic group; • denotes unclassified genus; •• denotes classified but unnamed genus. (C). Principal Coordinates Analysis (PCoA) illustrating beta diversity (Weighted UniFrac distance) between network clusters. (D). Comparison of Shannon index and bacterial load (16S rRNA gene copies/uL sputum plug) betweeen clusters (Mann–Whitney U-test). (E). LEfSe analysis of differentially abundant taxa between clusters. (F-L). Comparison of baseline lung function (FEV1% predicted) (F), breathlessness (mMRC) score (G), and sputum total cells (H), neutrophils (I), eosinophils (J), IL-1β (K), and IL-8 (L) between the two identified patient clusters. (M). Comparison of follow-up moderate-to-severe exacerbations between the two clusters. (N). Survival analysis of differences between the two clusters in the time to first exacerbation

Frequent exacerbators have reduced antagonistic interactions in microbiota while exhibiting subtle compositional changes

Finally, we focused specifically on the association of microbial network with the frequent exacerbator phenotype (≥ 2 moderate-to-severe exacerbations during 1-year follow-up). We classified COPD patients according to exacerbator phenotype (n = 137; 6 patients lost to follow-up) in EndAECOPD. Characteristics between frequent and infrequent exacerbators were shown in Supplementary Table 8. Frequent exacerbators (n = 19) had significantly greater breathlessness (mMRC) score and higher level of sputum IL-1β (P = 0.016) and S100A8 (P = 0.029) compared to infrequent exacerbators (n = 118) (Fig. 5A–E). Conventional microbiome analysis indicated a broad similarity of microbiome signatures, with no significant differences observed in microbial composition, beta diversity, alpha diversity or bacterial load (Fig. 5F–I). Yet interestingly, the frequent exacerbator phenotype presented the typical network disturbance, i.e., reduced antagonistic interactions (Fig. 5J). Besides, 7 of the 10 stability-related taxa and 5 of the 8 exacerbation-related taxa identified in above network analysis also showed reduced negative interactions in frequent exacerators (Fig. 5K). These findings suggest that network disturbance, instead of alterations in microbial abundance or load, may be a key driver of the frequent exacerbator phenotype in patients with COPD.

Fig. 5

Network analysis of the interactome in frequent exacerbators. Classification of the COPD patients in the EndAECOPD cohort according to exacerbator phenotype. (A-E). Comparison of breathlessness (mMRC) score (A), lung function (FEV1% predicted) (B), sputum total cells (C) and cytokines (D, E) between exacerbator phenotypes. Of the 19 detected cytokines, IL-1β (D) and S100A8 (E) were significantly elevated in frequent exacerbators. Box plots reflect median and IQRs, with whiskers (5–95 percentile) bounding non-outlier values. ns, not significant (Mann–Whitney U-test). IFE = infrequent exacerbator, FE = frequent exacerbator. F. Barplots of microbial abundance displaying similar taxonomic composition between exacerbator phenotypes which was confirmed by LEfSe analysis revealing no differentially abundant taxa. (G). PCoA illustrating similar beta diversity (Weighted UniFrac distance) between exacerbator phenotypes. (H, I). Comparison of Shannon index (H) and bacterial load (I) between exacerbator phenotypes with Mann–Whitney U-test. (J). Opposing network analysis of exacerbator phenotypes. The frequent exacerbators had significantly lower proportion of negative interactions in microbial network compared with infrequent exacerbators (P < 0.0001, χ2 test). # denotes paraphyletic group; • denotes unclassified genus; •• denotes classified but unnamed genus. K. Change of negative interactions (%) for the stability-related taxa (blue, up) and exacerbation-related taxa (purple, down) identified in above network analysis in frequent versus infrequent exacerbators. The dashed line represents ± 5%