Microbial community organization designates distinct pulmonary exacerbation types and predicts treatment outcome in cystic fibrosis

Polymicrobial infection of the airways is a hallmark of obstructive lung diseases such as cystic fibrosis (CF), non-CF bronchiectasis, and chronic obstructive pulmonary disease. Pulmonary exacerbations (PEx) in these conditions are associated with accelerated lung function decline and higher mortality rates. An understanding of the microbial underpinnings of PEx is challenged by high inter-patient variability in airway microbial community profiles. We analyzed bacterial communities in 880 CF sputum samples and developed microbiome descriptors to model community reorganization prior to and during 18 PEx. We identified two microbial dysbiosis regimes with opposing ecology and dynamics. Pathogen-governed PEx showed hierarchical community reorganization and reduced diversity, whereas anaerobic bloom PEx displayed stochasticity and increased diversity. A simulation of antimicrobial treatment predicted better efficacy for hierarchically organized communities. This link between PEx type, microbiome organization, and treatment success advances the development of personalized clinical management in CF and, potentially, other obstructive lung diseases.


Introduction
Obstructive lung diseases, such as cystic brosis (CF), non-CF bronchiectasis, and chronic obstructive pulmonary disease (COPD), are characterized by chronic polymicrobial bacterial infection of the airways.
Intermittent increases in signs and symptoms of respiratory dysfunction, so-called pulmonary exacerbations (PEx), are associated with lung disease progression and mortality in these conditions [1][2][3] .Despite their importance, the pathophysiologic events underlying PEx are unclear but generally believed to involve transient perturbation of host-microbial dynamics in the airways.Management of these events typically involves frequent, often aggressive, antibiotic treatment, which is intended to decrease bacterial burden and blunt host in ammatory response that contributes to lung pathology.This care carries considerable cost and treatment burden and is limited by drug toxicity and ever-increasing antimicrobial drug resistance 4 .In CF, therapies that modulate the activity of the dysfunctional cystic brosis transmembrane conductance regulator (CFTR), the primary cellular defect in CF, have reduced the frequency of PEx for many, but not all, people with CF 5,6 .Thus, a better understanding of PEx remains a high priority in efforts to improve care and enhance quality of life for persons with obstructive lung conditions 7 .
Dysbiosis de nes disease-associated alterations of the microbiome that affects the taxonomic composition as well as the functional activity of the microbial community 8 .This serves as an umbrella term for a variety of non-exclusive community characteristics including diversity loss, symbiont loss, or pathobiont blooms.As such, the label dysbiosis may be of limited applicability in describing microbial dynamics in chronic obstructive lung diseases insofar as the pulmonary microbiome in these conditions displays a markedly different ecology from that in healthy lungs, and can be considered, by de nition, to represent a dysbiotic state 9,10 .Nevertheless, given that a pathologic microbiome persists even during periods of relative clinical stability, a better characterization of its reorganization patterns to classify (relative) pulmonary dysbiosis into distinct types could provide opportunities for improved management of PEx in chronic pulmonary conditions 8,10,11 .
The search, in cross-sectional studies, for common motifs in microbial community processes that drive PEx, particularly in CF, has been hampered by subject-speci c microbiome con gurations.Employing longitudinal sampling strategies revealed highly individual taxonomic pro les with context-dependent metabolic activities and signaling have been observed in numerous studies [12][13][14][15] .Unpredictable and illde ned onset of PEx, as well as personalized antimicrobial treatment schemes to manage PEx, further complicate analyses 16 .A strategy that is capable of consolidating process communalities against the background of natural case variability is therefore required 17 .
Recent studies on microbial community networks have found that the precise con guration of dependencies among members de ne their community role, as well as the dynamical behavior of the microbiome [18][19][20] .Moreover, the formation of network clusters (i.e. the coexistence of microbial subcommunities) modulates robustness to external perturbation including antimicrobial therapy 21,22 .
Recently, time evolution of gut, vaginal and oral microbiomes were modeled using alternative community descriptors anchored in information theory 23 .Switching patterns of microbiota that emerged as tradeoff between perturbation, accessible niches and internal forces were identi ed.The impact of complex community organization on medically relevant microbial behaviors such as pathogen virulence or resilience to antimicrobial therapy is understudied and remains largely unexplored for clinical applications.
In this study, we developed non-standard descriptors that aggregate ecological and compositional properties of the CF lung microbiome and used these to identify PEx types with communal patterns.We then analyzed the organization of the CF microbiome in these backgrounds and revealed two, fundamental dysbiosis states: a hierarchical community reorganization controlled by the dominant pathogen and a stochastic reorganization with blooming anaerobic taxa and high taxonomic turnover.Of note, the behavior of a focal pathogen was markedly different with different community hierarchy.Lastly, we modeled targeted antimicrobial treatment on data-inferred co-occurrence networks and observed that distinct community organizations signi cantly determined treatment outcomes.

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 study 15,[24][25][26][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 speci cally, 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][29][30][31][32][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 speci city of the lung microbiome, which typically overshadows potential communalities.Accordingly, we identi ed 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 speci city 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 symptoms 34 ), and zygosity of the cftr F508del allele.As expected, we found that individual subject was a strong predictor for ASV covariance ( , ), followed, to lesser effects, by age ( , ) and clinical state ( , ).

Identifying distinct PEx types using non-standard descriptors
To reduce subject-speci c microbiome bias, we abandoned ASV composition as the sole sample descriptor, assembling the 194 core ASVs into ve higher-order groups.The rst 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 disease 35,36 .Three groups were categorized re ecting species oxygen requirement for growth 37 , considering that the CF lung microbiome is strongly conditioned by local oxygen gradients: strictly aerobic, strictly anaerobic, and facultatively anaerobic.The fth group comprised uncultivated taxa with unknown oxygen requirements.
Building on these ve 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 classi cation using Dirichlet multinomial mixtures (DMM).The DMM 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 (R 2 = 0.62 and R 2 = 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 rst three principal components were used to group similar samples.K-mer clustering identi ed three distinct PEx clusters or types using statistics (Figs.1C and 1D).
Having identi ed 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 con gurations.These results suggested that species sorting occurred in subjects' lungs according to oxygen requirements 38 .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 identi ed 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 con guration of the lung microbiota in and between the identi ed 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 rst asked whether PEx types could be simply explained by the DMM community classes, i.e. different community compositions 28,39 and whether oxygen availability could motivate shifts in microbiome con gurations 40 .We assessed the distribution and temporal change of DMM community classes previously modeled from coarse-grained ASV groups (Figs.2A, S2).Unexpectedly, no signi cant 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 results 29,30,41,42 .Neither pathogen load nor other recurrent organisms were consistent predictors for imminent PEx across larger patient cohorts.Here, we strati ed the microbiota by the identi ed 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 signi cant diversity evolution across samples culminating in antibiotic treatment (p PAT = 0.00126, p AN1 = 0.04462, p AN2 = 0.03672).The analogue analysis for Chao1 identi ed PAT and AN1 to exhibit signi cant dependency with time (p PAT = 4.333e-05, p AN1 = 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 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 regimes 7 .

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 weeks 15,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 signi cant differences of Shannon diversity and Chao1 richness (Figure S5).A split into 1-23 and 24-60 days before PEx treatment showed statistically signi cant relative changes in all three PEx types, in accordance with the previous result indicating diversity and time dependency.
Community turnover describes the rate of species compositional change over time as de ned by Ontiveros and colleagues 44 .We employed Aitchison distance to quantify community dissimilarity over time and assessed turnover as the slope of a tted, linear model.We analyzed turnover 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 was reduced during the 23 days compared to 24-60 days prior to PEx treatment ( 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 prede nes its emergent, dynamical capabilities 45 .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 networks 18,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 inference 18 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 (n PAT = 222, n AN1 = 192, n AN2 = 175; detailed description in Table S2).
We studied the topology of the largest network components, de ned as the ensemble of nodes belonging to the biggest connected subgraph of the network and therefore expected to be the most impactful for microbiome dynamics 46 .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 patterns 47 .
In analogy to interaction networks, we furthermore examined network hierarchy of the microbial cooccurrence networks across the PEx types.In the seminal work of Barabasi and Oltvai on biological interactions networks, a "quanti able signature of network hierarchy" was de ned as "the dependency of the clustering coe cient on the degree of a node, which follows " 48 .As a result, highly T T <24 = 0.17; T 24−60 = 0.26; p < 0.001 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 C i (k) were inferred and compared by the slope of a power law t across all co-occurrence networks within the same PEx type (Figs.3D and 3E).While the t to degree distributions showed similar slopes ( , , , , , , ), the slope of the node clustering distribution differed signi cantly between PEx types ( , , ).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 ts and data ( , , , ) suggesting at organization and unpronounced community structure.
To con rm 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 signi cantly 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 ( ). Clustering and betweenness centrality were statistically independent of tested diversity measures, while both diversities in uenced 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 identi ed the most hierarchical nodes, which were de ned as nodes with a degree k > 90% and a clustering coe cient C i < 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 graphs 49 , which corroborated these ndings.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 identi ed con gurations in contrast to microbiota in lung homeostasis.The de nition of dysbiosis employed in this work de nes a reorganization of the microbiota in the disease microenvironment.In healthy lungs, the pulmonary microbiome shows neutral community dynamics [50][51][52] .After Hubbell, diversity and abundance distributions of neutral communities can be explained by stochastic immigration and extinction events alone 53 .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 diminishes 13,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 rst, 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 dynamics 54 .The transition between the two metacommunity archetypes was explained by changes in dispersal due to altered spatial arrangements 38 .Here, we speculate that subject-speci c 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 stability 13,55 .We hypothesize that in the rst 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 at 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 bio lm formation, as well as the production of siderophores and exotoxins based on iron availability and oxygen levels 56,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 virulence 58,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 pro les 60,61 .
The insight that bacterial organization appeared markedly distinct in the identi ed PEx types raised the important question whether microbiome organization could modulate pathogen importance or interfere with treatment outcomes in a foreseeable manner.As a rst 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 signi cantly higher page rank in hierarchic than in at 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 microbiota 62 .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 at 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 dynamics [63][64][65] .Moreover, the particular con guration of microbial interactions has also been identi ed as a key element for robust composition forecasting 66 .In complex networks, interaction topology not only determined the importance of hubs for systems dynamics 20 , but also controllability of dynamics supporting our previous conclusions 67 .
To clarify the relevance of community organization for PEx treatment strategies, we asked whether community (network) organization could in uence pathogen importance and modulate treatment outcomes.Recently, it was shown that dynamical uctuation in response to external damage depended on the local network architecture of single layered and multiplex networks [63][64][65] .In an analogue, simpli ed 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 signi cantly 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 bene t repopulation after depletion, while niche reorganization could instigate the establishment of different community con gurations.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 communities 68,69 .In fact, response to antimicrobials may be recognized as an emergent property of the entire microbiome 70 .
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 at community background.

Conclusions
Our study on the human lung microbiome showcases the importance of community organization for understanding microbiome dynamics in homeostasis and dysbiosis.Using CF as a model disease, we employed a functional/ecological coarse graining to analyze both temporal and organizational aspects of pulmonary infection and its relevance for therapy outcomes.
We identi ed two archetypes of dysbiosis in the CF lung (driven by pathogens or anaerobes), characterized their community structures, and discussed stabilizing factors in the context of current graph and ecological theory.It is important to realize that the identi ed ecological features can cancel out if analyzed cross-sectionally due to their antagonistic nature.This might explain the di culties with establishing robust predictors for PEx thus far 7 .
We modeled the focal depletion of the most abundant pathogen in empirical co-occurrence networks and recognized that distinct background communities shaped the outcome of this treatment simulation.We concluded that the relevance of pathogenic taxa for microbiome dynamics, disease progression, and treatment effect is systematically linked to the organization of the background microbiome.
Our insights are limited to the ecological dynamics in airway microbiota observed in 11 subjects.Despite the large dataset investigated (880 samples which comprise a tidy and comparable subset of a collection incorporating > 21 patient years of near-daily sampling), this study cannot possibly reveal the full complexity of airway microbial ecology in CF.Furthermore, none of the subjects in our study were receiving CFTR modulator therapy at the time of sample collection.How this therapy will impact CF airway microbiology is the subject of on-going studies.Nevertheless, we believe our model and the observations made in this study contribute to generating novel hypotheses regarding CF lung pathology, building theory for targeted dysbiosis management, enhancing antibiotic stewardship, and advancing personalized medicine.

Study cohort, sample collection, and sample inclusion criteria
Expectorated sputum was collected from a cohort of people with CF as part of a long-term prospective study of CF airway microbiota.This study was approved by the Institutional Review Board of the University of Michigan Medical School (HUM00037056), and informed written consent was obtained from all participants.Subjects collected daily sputum samples at home, which were stored at either 4 o C or -20 o C before shipment to the University of Michigan for immediate storage at -80 o C. Electronic medical records were reviewed for subject demographic and clinical data.Among the 880 sputum samples and DNA sequences included in this study, 283 were reported previously (NCBI BioProjects PRJNA520924 and PRJNA611611) 24,25 .Sample inclusion criteria applied for downstream analyses are detailed in Tables S1,  S2.

DNA extraction
Sputum samples were thawed on ice and homogenized with 10% Sputolysin (MilliporeSigma, Burlington, MA, USA).Samples were treated with bacterial lysis buffer (Roche Diagnostics Corp., Indianapolis, IN, USA), lysozyme (MilliporeSigma), and lysostaphin (MilliporeSigma) as previously described 31 , followed by mechanical disruption by glass bead beating and digestion in proteinase K (Qiagen Sciences, Germantown, MD, USA).DNA was extracted and puri ed using the MagNA Pure nucleic acid puri cation platform (Roche Diagnostics Corp., Indianapolis, IN, USA) according to the manufacturer's protocol.
Sequencing controls, protocol and taxonomic annotation DNA libraries were prepared by the University of Michigan Microbial Systems Molecular Biology Laboratory as described previously 71 .Human Microbiome Project (HMP) or Zymo (Zymo Research, Irvine, CA, USA) mock community standards were included on each sequencing plate.In brief, the V4 region of the bacterial 16S rRNA gene was ampli ed using touchdown PCR with barcoded dual-index primers.
Touchdown PCR was performed consisting of 2 min at 95°C, followed by 20 cycles of 95°C for 20 sec, 60°C (starting from 60°C, the annealing temperature decreased 0.3°C each cycle) for 15 sec, and 72°C for 5 min, followed by 20 cycles of 95°C for 20 sec, 55°C for 15 sec, and 72°C for 5 min and a nal 72°C for 10 min.The resulting amplicon libraries were normalized and sequenced on an Illumina sequencing platform using a MiSeq Reagent Kit V2 (Illumina Inc., San Diego, CA, USA).The nal load concentration was 4.0-5.5 pM with a 15% PhiX spike to add diversity.
Samples were processed separately by subject through sample inference to avoid batch effects, then merged for the remaining processing steps.The data set was denoised removing all ASVs with < 0.0075% average abundance across all samples as previously described 73 and subsequently rari ed using R package vegan 74 .Sequencing data, taxonomic information and clinical meta data were organized in phyloseq objects for further analysis 75 .Sequencing error rates were determined by comparing 43 mock community pro les to reference sequences in mothur (v1.48) 76 using the seq.error command, which measures error as the sum of mismatches to the reference divided by the sum of bases in the query.In 24 sequencing runs, the median mock community error rate was 0.037% (range 0.012% − 0.690%).

Community typing with Dirichlet multinomial mixture models
To stratify representative community classes across subjects, Dirchlet multinomial mixtures (DMM) were inferred employing ASV groups as the taxonomic level 39 .Counts from ASVs were summarized into ve groups: CF pathogens, strict anaerobes, facultative anaerobes, strict aerobes and unknown according to the oxygen requirements of the respective taxon 37 .Next, the total data set was subject-strati ed to avoid bias.Thirty-six random data subsets with 650 samples each were generated by sampling with replacement from the total data collection.Subsequently, models with 1-25 DMMs were inferred stepwise for each subset using the R package DirichletMultinomial 77 .Laplace approximation, BIC, and AIC were queried independently to identify the optimal number of DMMs.

Variance testing and ordination
To quantify the impact of covariates on the lung microbiome, PERMANOVA was performed on rari ed ASV data and the identical data were remodeled by non-standard sample descriptors 74 .Bray-Curtis distance was employed for ASVs and Euclidian distance for scaled and centered non-standard descriptors.The model was designed to test the marginal effects of the individual covariates (function adonis2, parameter setting by = margin).For comparison, effect sizes ( ) of the covariates were calculated using the adonis_omegaSq function from the MicEco package 78 .Principal component analysis was performed on scaled and centralized sample descriptors in R. Sample descriptors included Shannon diversity, Chao1 richness, relative abundance of the most abundant CF pathogen, the ratio of CF pathogen counts to counts from anaerobic taxa and the classi cation to a particular DMM community class.

PEx type clustering and sample classi cation
Hierarchical k-mer clustering of samples was conducted on the rst three principal components of the ordinated sample data using the R package pheatmap 79 .Pearson correlation was employed as a similarity measure.Next, an contingency test was used to identify the best k cluster number for the classi ed sputum samples.Subsequently, entire PEx time series and networks were assigned to a single k-mer cluster by majority vote of the included samples.Two PEx time series were excluded from further analysis, because frequent type transitions prevented a conclusive association to a single PEx type.
Detailed inclusion criteria are explained in Tables S1, S2.

Co-occurrence network inference and network statistics
Co-occurrence networks were inferred from 20 samples collected at consecutive days with SparCC in a sliding window along each PEx time series 85 .Missing samples were imputed using R package seqtime 62 and jump size for the sample window was set to 1.
For downstream analysis on the largest network components, only strong ( ) and statistically signi cant (p < 0.01 after FDR correction) co-occurrence edges were included, as well as networks with < 10 imputed samples.Topological properties were assessed using R package igraph 49 .
Node degree and node clustering distributions were calculated across all networks classi ed in the same PEx type.To identify power laws and their respective slopes, linear regressions were performed on logtransformed data using R package ggplot2 83 .
For comparing the effect size of PEx clusters with covariates Shannon and Chao1 on graph topology, we rst calculated clustering coe cient, betweenness centrality, node counts, and edge counts for each largest component and averaged Shannon and Chao1 of all samples included for inference of the respective network.Next, we implemented independent LMMs, corrected for subject and assessed the effects sizes (partial ) as described previously.
To quantify the presence of certain ASVs in top hierarchical network positions, we rst selected positions with the relative highest degree (> 90% of all degree values) and relative lowest clustering (< 10% of all clustering coe cients).Subsequently, the frequency of ASVs on these positions was counted, normalized, and ranked.All calculations were performed in R.

Noise analysis of ASV time behavior
For each PEx time series, noise colors of participating ASVs were inferred using seqtime as previously described 62 .In short, missing samples were interpolated, and rounded to counts, negative interpolation values were set to 0 counts.Subsequently, the wrapper function identifyNoisetypes() performed a spectral density estimate and calculated a linear t to the resulting periodogram (log frequency vs log spectral density) of the ASV time series.According to the slope of the t, ASV time series were classi ed into categorical noise color groups.Noise colors were plotted as ratios of pathogens vs anaerobic ASVs with similar color (relative abundances) in the same sample.

Pathogen removal
To identify the dominant pathogen by network, the subset of samples used to infer the individual cooccurrence network was queried and the pathogen with highest cumulative abundance was selected.
Only networks with the dominant pathogen locating to the major component were used for further analysis.We calculated modularity, the number of unconnected components and node size of the largest    See image above for gure legend.
using the ve ASV groups as input and identi ed 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.

Figures Figure 1
Figures

Figure 2 See
Figure 2

Table Table 1
component independently for each co-occurrence network before and after pathogen removal.All parameters were normalized to the corresponding value before node removal.Density distributions of the normalized parameters were scaled and plotted with ggplot2 function geom_density() 83 .To test for difference of the cumulative parameter distributions, two-sided Kolmogorov-Smirnov tests were performed (ks.test(), stats package)86. is available in the Supplementary Files section.