Compositional characterization of PEx time series
We aggregated a collection of 880 sputum samples from 11 adults with CF, comprising 18 PEx time series. The characteristics of the study subjects and sputum samples are provided in Table 1; sample inclusion criteria are provided in Tables S1 and S2. Subjects and sputum samples were chosen from a larger dataset that had been generated during the course of a long-term observational study15,24–27. Subjects were selected from this larger dataset based on the availability of near-daily sputum samples (i.e., a sample available from at least 60% of days) that spanned periods of clinical stability culminating with a PEx that prompted antibiotic treatment by the subject’s care team. More specifically, samples in each PEx time series were collected from 60 days prior to one day prior to the initiation of treatment for PEx (Figure S1). The time frame of 60 days prior to the start of PEx antibiotic treatment was selected to accommodate potential changes in the lung microbiome preceding symptom onset together with changes occurring during the acute PEx phase. As we have done in previous studies 15,24,28–33, samples were further characterized based on the subject’s clinical state at the time of collection: baseline samples were those collected between 60 and 15 days prior to the start of antibiotic treatment; exacerbation samples were collected between 14 and one day prior to the start of antibiotic treatment. Neither samples obtained during acute antibiotic treatment for PEx (treatment) nor within three weeks after PEx treatment stopped (recovery) were included in this analysis. Chronic (maintenance) antibiotic therapies such as inhaled tobramycin and aztreonam, and oral azithromycin used on each day were recorded. A mean of 49 (SD, 7.3) sputum samples were analyzed per PEx time series.
Identifying common microbiologic features of PEx in CF is challenged by the pronounced subject specificity of the lung microbiome, which typically overshadows potential communalities. Accordingly, we identified 1,949 amplicon sequence variants (ASVs) among the 880 samples, with only eight ASVs present in every subject. For further analyses, the data set was denoised to 194 core ASVs by removing taxa present with an average relative abundance below 0.0075%. To quantify the degree of subject specificity in the data set, a PERMANOVA test was performed to calculate the effect sizes of clinical and demographic covariates on data variance. Covariates included subject, subject age, subject sex, clinical state (baseline health or exacerbation of symptoms34), and zygosity of the cftr F508del allele. As expected, we found that individual subject was a strong predictor for ASV covariance (\({\omega }^{2}=0.51\),\(p<0.001\)), followed, to lesser effects, by age (\({\omega }^{2}=0.023\),\(p<0.001\)) and clinical state (\({\omega }^{2}=0.003\),\(p<0.003\)).
Identifying distinct PEx types using non-standard descriptors
To reduce subject-specific microbiome bias, we abandoned ASV composition as the sole sample descriptor, assembling the 194 core ASVs into five higher-order groups. The first group comprised conventional CF pathogens (Pseudomonas, Staphylococcus, Burkholderia, Haemophilus, Achromobacter, and Stenotrophomonas) based on the prominent role these species are believed to play in CF lung disease35,36. Three groups were categorized reflecting species oxygen requirement for growth37, considering that the CF lung microbiome is strongly conditioned by local oxygen gradients: strictly aerobic, strictly anaerobic, and facultatively anaerobic. The fifth group comprised uncultivated taxa with unknown oxygen requirements.
Building on these five ASV categories, we assembled the following non-standard descriptors for every sputum sample: i) the ratio of CF pathogens to strict anaerobes, ii) the relative abundance of the most abundant CF pathogen, iii) the Shannon diversity index of the core ASVs, iv) the Chao1 richness of the core ASVs, and v) a community classification using Dirichlet multinomial mixtures (DMM). The DMM model was implemented using the five ASV groups as input and identified six community classes. Two DMM community classes were dominated by CF pathogens, three by anaerobes, and one by facultatively anaerobic organisms. Community classes, selection of Dirichlet components, class distribution over the cohort and class-wise sample compositions are presented in Figure S2.
Using this suite of descriptors, we implemented a second PERMANOVA model and found that the variance explained by this model was reduced compared to that based solely on ASVs (R2 = 0.62 and R2 = 0.78, respectively). Most importantly, the effect size of subject bias decreased by 51% (Fig. 1A). Data independent of clinical state, were ordinated using principal component analysis (Fig. 1B), and the first three principal components were used to group similar samples. K-mer clustering identified three distinct PEx clusters or types using \({\chi }^{2}\) statistics (Figs. 1C and 1D).
Having identified three robust clusters, hereafter referred to as PEx types, we next analyzed the distribution of DMM communities among these. PEx Type 1 (hereafter called PAT) comprised communities dominated by conventional CF pathogens, including Pseudomonas, Burkholderia, Achromobacter, Haemophilus, Staphylococcus, and Stenotrophomonas. PEx Type 2 (AN1) and Type 3 (AN2), on the other hand, were driven by three distinct anaerobic community configurations. These results suggested that species sorting occurred in subjects’ lungs according to oxygen requirements38. Importantly, PEx proceeded in both aerobic and anaerobic communities.
To assign subjects and their PEx time series to a single PEx type, we performed Spearman’s rank association (Figures S3A and S3B). Two time series were excluded from further analyses due to inconclusive association to a single type (time series 9, 12). The reduced number of 789 samples distributed as 286, 254, and 249; the number of subjects as 4, 3, and 4; and the number of PEx as 6, 5, and 5 to PEx types PAT, AN1, and AN2, respectively (Figures S3C and S3D). We found that sample association with PEx type PAT was remarkably stable both at the level of individual PEx time series (60 days), as well as with subjects over time. On the contrary, more transition events were observed between PEx types AN1 and AN2 (Figures S4A and S4B). Overall, subjects showed a tendency to persist either in PAT or in AN1 or AN2 despite recurrent antibiotic treatment between time series (treatment samples excluded, Figure S4C).
In summary, aggregated measures of sample diversity, ecology and function were used to reduce the organism-driven subject bias and group PEx trajectories with similar properties. We identified three communal PEx types among subjects, termed PAT, AN1 and AN2, that displayed distinguishable microbiomes.
Temporal behavior of microbiomes in distinct exacerbation regimes
We studied the configuration of the lung microbiota in and between the identified PEx types and modelled common re-organization patterns over time as the community proceeded towards the start of PEx treatment. To elucidate underlying ecological processes, we first asked whether PEx types could be simply explained by the DMM community classes, i.e. different community compositions28,39 and whether oxygen availability could motivate shifts in microbiome configurations40. We assessed the distribution and temporal change of DMM community classes previously modeled from coarse-grained ASV groups (Figs. 2A, S2). Unexpectedly, no significant temporal evolution of DMM communities was observed within PEx types, indicating that the overall proportions of pathogens, anaerobes, facultative anaerobes, and aerobes persisted over most of the PEx cycle with few exceptions. These sporadic shifts occurred only between comparable community classes, i.e., due to continuous transitions (increase or decrease) of taxonomic groups.
Several studies have investigated microbiome structure and rearrangement prior to PEx with inconsistent results29,30,41,42. Neither pathogen load nor other recurrent organisms were consistent predictors for imminent PEx across larger patient cohorts. Here, we stratified the microbiota by the identified PEx types and analyzed diversity and richness over time in trajectories with similar properties. Mixed effect models were implemented to test time dependencies of Shannon and Chao1 for the three PEx types and corrected for confounders (subject and PEx cycle) (Figs. 2B and 2C). All PEx types displayed significant diversity evolution across samples culminating in antibiotic treatment (pPAT = 0.00126, pAN1 = 0.04462, pAN2 = 0.03672). The analogue analysis for Chao1 identified PAT and AN1 to exhibit significant dependency with time (pPAT = 4.333e-05, pAN1 = 0.025).
Interestingly, richness and diversity decreased towards treatment for the pathogen-dominated communities (PAT) and increased for anaerobic PEx types (AN1 and AN2). Furthermore, it is important to note that the time dependency of richness and diversity were consistently small and ranged between \({\eta }^{2}=\{0.02, 0.06\}\) for all PEx types. In short, we revealed that diversity evolves opposingly, as the microbial communities approached PEx treatment. Of note, despite the modest effect size, these results have the potential to explain the inconclusive reports of previous studies that were conducted without consideration of PEx regimes7.
Species turnover displays antagonistic patterns in pathogen or anaerobe communities
Evidence suggests that changes in airway microbial community structures may precede the onset of clinical symptoms of PEx by days or even weeks15,41,43. To determine the most likely time interval for such changes, samples were systematically grouped by collection time (days before the initiation of antibiotic treatment for PEx) and tested for significant differences of Shannon diversity and Chao1 richness (Figure S5). A split into 1–23 and 24–60 days before PEx treatment showed statistically significant relative changes in all three PEx types, in accordance with the previous result indicating diversity and time dependency.
Community turnover \(T\) describes the rate of species compositional change over time as defined by Ontiveros and colleagues44. We employed Aitchison distance to quantify community dissimilarity over time and assessed turnover as the slope of a fitted, linear model. We analyzed turnover \(T\) during onset of (1–23 days prior to treatment) and prior to (24–60 days) PEx by evaluating Aitchison distance between any two samples collected in an interval of one to 20 days in a subject-wise manner (Fig. 2D). Overall, dissimilarity was smaller in the pathogen-dominated PEx type PAT (ANCOVA p < 0.001) and turnover \(T\) was reduced during the 23 days compared to 24–60 days prior to PEx treatment (\({T}_{<24}=0.17; {T}_{24-60}=0.26; p<0.001\) for both tests). Interestingly, the anaerobic PEx types AN1 and AN2 again exhibited antagonistic patterns, with increased species turnover shortly before PEx treatment (ANCOVA p < 0.001 for both tests). Together, the previous results suggested two PEx regimes (PAT vs AN) with antagonistic temporal behavior.
Characteristic community reorganizations stratify pulmonary dysbiosis types
The detailed organization of interactions and dependencies throughout an ecological community predefines its emergent, dynamical capabilities45. In particular, resilience to perturbations such as antimicrobial treatment, community robustness, and the stabilizing effect of keystone organisms were previously attributed to properties of dependency networks18,19,22. Therefore, it is not only important to identify the most relevant CF pathogen in the airway microbiome, but to understand how the background community organization impacts the focal driver organism, modulates its virulence, and contributes to stability.
To study community organization, we inferred co-occurrence networks from sample subsets of individual PEx time series. In detail, for every network, 20 consecutively collected samples were used for robust inference18 and a sliding window was employed to work across the individual PEx time series (with a step size of one sample). This approach yielded 589 co-occurrence networks, where topology changes between successive networks were caused by the substitution of a single sample. The resulting graphs were subsequently analyzed by PEx type (nPAT = 222, nAN1 = 192, nAN2 = 175; detailed description in Table S2).
We studied the topology of the largest network components, defined as the ensemble of nodes belonging to the biggest connected subgraph of the network and therefore expected to be the most impactful for microbiome dynamics46. For PEx type PAT, a reduced number of organisms and associations were observed in the largest component, as well as increased betweenness centrality (Wilcoxon p < 0.001 for each pairwise test; Figs. 3A, 3B and 3C) in contrast to PEx types AN1 and AN2. Graph betweenness centrality measures the extent of centralized organization reinforcing effective communication patterns47.
In analogy to interaction networks, we furthermore examined network hierarchy of the microbial co-occurrence networks across the PEx types. In the seminal work of Barabasi and Oltvai on biological interactions networks, a “quantifiable signature of network hierarchy” was defined as “the dependency of the clustering coefficient \(C\) on the degree \(k\) of a node, which follows \(C\left(k\right)\sim{k}^{-1}\)”48. As a result, highly connected hub nodes should display low clustering, if located on top of the network hierarchy. To test this in co-occurrence networks, degree distributions p(k) together with clustering distributions Ci(k) were inferred and compared by the slope of a power law fit across all co-occurrence networks within the same PEx type (Figs. 3D and 3E). While the fit to degree distributions showed similar slopes (\({\alpha }_{PAT}=-0.71\), \({\alpha }_{AN1}=-0.6\), \({\alpha }_{AN2}=0.-69\), \({R}_{PAT}=-0.64\), \({R}_{AN1}=-0.67\), \({R}_{AN2}=-0.73\), \(p<0.001\)), the slope of the node clustering distribution differed significantly between PEx types (\({\alpha }_{PAT}=-0.97\), \({\alpha }_{AN1}=-0.27\), \({\alpha }_{AN2}=-0.28\)). The pathogen-driven PEx type PAT displayed the strongest descent of clustering with degree k, was indeed approximating − 1, and hence indicated a clear microbial hierarchy. These compelling results supported the hypothesis of a pronounced, hierarchical dysbiosis type. Moreover, the anaerobic PEx types AN1 and AN2 also exhibited weaker correlation between fits and data (\({R}_{PAT}=-0.61\), \({R}_{AN1}=-0.3\), \({R}_{AN2}=-0.29\), \(p<0.001\)) suggesting flat organization and unpronounced community structure.
To confirm that these results were driven by PEx types rather than sample diversity, we implemented independent linear mixed effect models for every graph readout, corrected for subject and calculated effect sizes of PEx types and covariates Shannon diversity and Chao1 richness. We found that betweenness centrality, clustering, number of vertices, and number of edges significantly depended on PEx types with effect sizes being 2.9 times, 5.5 times, 8.9 times and 3.4 larger than the most effective diversity measures, respectively (\({p}_{between}<0.004, {p}_{cluster}<0.048, {p}_{\#vertex}<2.7e-7, {p}_{\#edge}<7.1e-11\)). Clustering and betweenness centrality were statistically independent of tested diversity measures, while both diversities influenced edge numbers and richness graph size to a minor extent (Figure S6).
Next, we investigated which organisms preferentially occupied the most hierarchical positions in the communities and therefore likely controlled the overall microbiome dynamics. For each network, we identified the most hierarchical nodes, which were defined as nodes with a degree k > 90% and a clustering coefficient Ci < 10% of all nodes in the graph. ASV frequencies were then assessed on these positions (Figs. 3F). In PAT communities, Pseudomonas, Staphylococcus and Streptococcus were not only masters in hierarchy, but also belonged to the most abundant ASV group in the samples. In AN2, a stronger variation of taxa was observed in the most hierarchical nodes. Of note, one third of the ranking was occupied by various ASVs belonging to the genus Prevotella in these communities. In Figure S7 we visualized three representative PEx type communities with exemplary hierarchies using the Sugiyama algorithm for hierarchical graphs49, which corroborated these findings. We concluded that in anaerobic communities, individual species were less relevant for overall microbiome dynamics, and microbiome organization was increasingly stochastic and less well picked up by co-occurrence analysis. To the contrary, in the pathogen-driven PEx type, the network hierarchies were conserved, occupied by a few key organisms and well supported by a simple, centralized community organization.
To contextualize these observations, we examined the identified configurations in contrast to microbiota in lung homeostasis. The definition of dysbiosis employed in this work defines a reorganization of the microbiota in the disease microenvironment. In healthy lungs, the pulmonary microbiome shows neutral community dynamics50–52. After Hubbell, diversity and abundance distributions of neutral communities can be explained by stochastic immigration and extinction events alone53. To the contrary, in chronic lung disease, microbial interactions, local replication, and environmental adaptation become key for diversity and community dynamics, and the impact of dispersal diminishes13,50. Consequently, if neutrality is a property of microbial eubiosis in the lung, then microbial interactions can be considered a hallmark of dysbiosis of the pulmonary microbiota. Together with metabolic adaptations, such interactions promote outgrowth of certain taxa to high relative abundances. Accordingly, we propose that both observed community states resemble different fundamental kinds of dysbiosis: the first, a structured, interacting community under the governance of an abundant, conventional CF pathogen, and the second, a globally successful functional guild that gains abundance by adapting to selective environmental pressures. Of note, similar community archetypes characterized by species-sorting or mass effects were described in metacommunity theory, a framework for ecological community assembly and dynamics54. The transition between the two metacommunity archetypes was explained by changes in dispersal due to altered spatial arrangements38. Here, we speculate that subject-specific mucus accumulation and decreasing oxygen availability in the lung microenvironment determine CF dysbiosis states in equivalent ways.
Importantly, both community states are robust maladaptations to the disease conditions of the lung, which raises the question whether negative loops exist in the system that enable their dynamical stability13,55. We hypothesize that in the first state, functional adaptations of the dominant pathogens together with antimicrobial defense against microbial competitors provide important negative feedback, whereas limitations of available niche space stabilize the second regime.
CF pathogens drive PEx dynamics in hierarchical, but not in flat community organization.
The virulence of pathogenic bacteria depends on microbial interactions and the biochemistry of the microenvironment among other factors. For example, Pseudomonas aeruginosa tightly regulates biofilm formation, as well as the production of siderophores and exotoxins based on iron availability and oxygen levels56,57. Moreover, the fermentation products 2,3-butanediol and lactic acid produced by anaerobic members of the CF microbiome were reported to trigger quorum-sensing and further virulence58,59. Conversely, synergistic interactions such as metabolic cross-feeding affects pathogen growth and lowers the tolerance of P. aeruginosa to antimicrobial treatment independent of intrinsic antibiotic resistance profiles60,61.
The insight that bacterial organization appeared markedly distinct in the identified PEx types raised the important question whether microbiome organization could modulate pathogen importance or interfere with treatment outcomes in a foreseeable manner. As a first step, page rank was used as a statistical descriptor for network importance to compare the importance of conventional pathogens, strictly anaerobic, facultatively anaerobic, and strictly aerobic taxa in the community. We found that CF pathogens were differentially important for the CF community, displaying significantly higher page rank in hierarchic than in flat community organization (Fig. 4A). Next, pathogen dynamics in different community organizations were assessed using time series information. Previously, we demonstrated that stochastic ecological processes can be distinguished from interaction-driven processes by Fourier spectra inferred from the abundance changes of microbiota62. Here, the spectrum of every ASV per PEx time series was determined and their noise color was inferred. We observed that pathogens exhibited interaction-driven dynamics (pink noise) in steep hierarchies only (Fig. 4B). In flat anaerobe-dominated organization, pathogens instead showed stochastic behavior (white noise, AN2) or anaerobic taxa dominated community dynamics (AN1). These results supported the hypothesis that CF pathogen activity depends on the community background and its positioning within (Figure S7). Indeed, in graph theory it was demonstrated repeatedly that the topology of interaction networks was intimately linked with its dynamics63–65. Moreover, the particular configuration of microbial interactions has also been identified as a key element for robust composition forecasting66. In complex networks, interaction topology not only determined the importance of hubs for systems dynamics20, but also controllability of dynamics supporting our previous conclusions67.
To clarify the relevance of community organization for PEx treatment strategies, we asked whether community (network) organization could influence pathogen importance and modulate treatment outcomes. Recently, it was shown that dynamical fluctuation in response to external damage depended on the local network architecture of single layered and multiplex networks63–65. In an analogue, simplified approach, we modeled the response of our empirical co-occurrence networks to focal depletion of the dominant pathogen by antibiotic treatment. We hypothesized that community organization affected the degree of network disruption and consequently community dynamics and likely treatment outcomes. In the model, the most abundant pathogen was removed from the major component of 313 networks (PAT = 152, AN1 = 53, AN2 = 108) and resulting network disruption was assessed by monitoring modularity change, breakup into subcomponents, and size reduction after single node removal (Figs. 4C, 4D and 4E). While total pathogen removal may not always be achieved in practice, microbial interactions were expected to abate together with strain abundance. Here we modeled the extreme case for the purpose of hypothesis testing. We found that pathogen elimination from steeper background hierarchies resulted in significantly stronger topology disruption supported by all three topological parameters (Kolmogorov-Smirnov statistic in Table S3). We observed stronger change in modularity, increased disruption into unconnected subcomponents, as well as more pronounced size loss of the biggest connected component. Indeed, the depletion of the same organism resulted in divergent effects for the overall community architectures. The markedly different outcomes might serve as indicators for the degree of niche rearrangement after antibiotic treatment. We hypothesized that maintained niche accessibility should benefit repopulation after depletion, while niche reorganization could instigate the establishment of different community configurations. Although these data-derived hypotheses call for rigorous experimental testing, they are in line with previous clinical and experimental observations reporting altered responses of focal organism to antibiotic treatment in different background communities68,69. In fact, response to antimicrobials may be recognized as an emergent property of the entire microbiome70.
We concluded that the relevance of CF pathogens for microbial community dynamics and, by extension, likely also clinical course, was crucially shaped by community organization of the CF microbiome. Moreover, targeted treatment of pathogens resulted in distinct responses as a function of microbiome hierarchy, i.e., steep or flat community background.