Evolutionary highways to persistent bacterial infection

Abstract

Persistent infections require bacteria to evolve from their naïve colonization state by optimizing fitness in the host via simultaneous adaptation of multiple traits, which can obscure evolutionary trends and complicate infection management. Accordingly, here we screen 8 infection-relevant phenotypes of 443 longitudinal Pseudomonas aeruginosa isolates from 39 young cystic fibrosis patients over 10 years. Using statistical modeling, we map evolutionary trajectories and identify trait correlations accounting for patient-specific influences. By integrating previous genetic analyses of 474 isolates, we provide a window into early adaptation to the host, finding: (1) a 2-3 year timeline of rapid adaptation after colonization, (2) variant "naïve" and "adapted" states reflecting discordance between phenotypic and genetic adaptation, (3) adaptive trajectories leading to persistent infection via three distinct evolutionary modes, and (4) new associations between phenotypes and pathoadaptive mutations. Ultimately, we effectively deconvolute complex trait adaptation, offering a framework for evolutionary studies and precision medicine in clinical microbiology.

Fig. 1

Study design. Upper panel: Every month, CF patients are seen at the CF clinic at Rigshospitalet in Copenhagen, Denmark. Here they deliver a sputum or endolaryngeal suction sample where selective microbiological culturing is performed. The longitudinally collected isolates have been genome sequenced and analyzed previously. Middle panel: Longitudinally collected isolates have been subjected to different phenotypic analyses for this study and are here (lower panel) analyzed using two data modelling approaches: Archetype analysis (AA) and Generalized Additive Mixed Model (GAMM). By integrating these approaches, we map dominant evolutionary trajectories and analyze mechanistic links between phenotypic and genetic adaptation

Fig. 2

Phenotypic characterization. We present summary statistics of our phenotype screen including a mean and mean standard deviation for each phenotype over all isolates as well as the P. aeruginosa PAO1 value and antibiotic breakpoint we use for normalization, respectively (above boxplot). We also show boxplots of continuous normalized variables (including the median as the center line, first and third quartile box bounds and whiskers representing 1.5× inter-quartile range). We also show the overall count of isolates with presence/absence of mucoidity, protease and hypermutator phenotypes and a maximum likelihood phylogeny (1000 bootstraps) of the DK19 clone type; nodes marked with white triangles have bootstrap values >=50. Blue stars represent late (>3 years) phenotypically naive (7–8 naive phenotypes out of 8) isolates from patient P7204. Circles with different shades of grey represent isolates from patients marked on the outer edge of the circle. We then compare the b expected adaptation over time based on field consensus versus c the measured raw adaptation of our isolate collection over time. The X-axis represents the time since colonization of a specific lineage or “colonization time”. Colors are linked with the expected change of the specific phenotype (b), so that blue denotes a “naive” phenotype and red denotes an “evolved” phenotype. For growth rate (in artificial sputum medium (ASM)), adhesion and aggregation, naive and evolved phenotypes are determined by comparison with the reference isolate PAO1 phenotype. For aztreonam and ciprofloxacin MIC, naive and evolved phenotypes are based on sensitivity or resistance as indicated by the EUCAST breakpoint values as of March 2017

Fig. 3

AA and GAMM models. We present a summary of the models underpinning our study of pathogen adaptation. a Screeplot showing the average residual sum of squares (RSS) for 25 iterations of each fit of a given number of archetypes. The “elbow” of the plot indicates that six archetypes are sufficient to model our dataset. b Characteristic trait profiles describing the five distinct phenotype levels that each of our 6 archetypes represents. We use the following abbreviations to represent our normalized data: grASM for growth rate in ASM, agg for aggregation, adh for adhesion, azt for aztreonam susceptibility, and cip for ciprofloxacin susceptibility. c Simplex plot of the AA showing the six archetypes (A1–A6) sorted by their characteristic growth rate (A3 and A5 vs A2 and A6), decreased sensitivity towards ciprofloxacin (A1 and A6), and increased aggregation and adhesion (A2 and A4). All further simplex visualizations are also sorted accordingly and can be interpreted using this key, which is annotated with the extreme phenotype values for each archetype. The complete analysis can be found in Supplementary Note 1. d P-values for GAMM models with multiple explanatory variables (columns) for the six predictor variables (rows) after model reduction. P-values are only shown for explanatory variables that showed a significant (p-value<0.01, GAMM with Wald-type tests) impact on the predictor in question. The complete analysis can be found in Supplementary Note 2

Fig. 4

Rapid early adaptation. GAMMs illustrate the significant impact of the explanatory variable colonization time on a growth rate in ASM, b ciprofloxacin sensitivity in ASM, and e the accumulation of all mutations (orange), nonsynonymous SNPs (blue) and indels (insertions and deletions). We also illustrate the proposed initial adaptation period by dark grey reference lines in a and b. (a/b/e notation) GAMMs are illustrated by solid smoothed trendlines, dashed two standard error bounds, and gray points as residuals. Y-axes are labelled by the predictor variable on which the effect of colonization time of the clone type has been estimated as well as the estimated degrees of freedom (edf) (for the e upper panel the edf is ordered as all mutations/NS SNPs). Residuals have not been plotted in the upper panel of e for clarity reasons. X-axes are the colonization time in years and patients are included as random smooths together with time. A rug plot is also visible on the x-axis to indicate the density of observations over time. Simplex visualizations of AA show (c) naive trait alignment of the first isolate of the twenty patients where we have analyzed the first P. aeruginosa isolate ever cultured at the CF clinic (blue circles) in contrast to evolved isolates that have been cultured at year 2–3 of colonization (red squares, all patients of the dataset). We contrast this trait-based ordination with (d) genetic adaptation, shown by a color overlay of the number of nonsynonymous mutations present in each isolate. Isolate 95 (purple circle) of the DK12 clone type has a very high number of mutations (>100) because one isolate in that lineage (isolate 96) is very different from the remaining 11 isolates. For the GAMM analysis shown in Fig. 4e, we filtered out the mutations from the errant DK12 96 single isolate that affected the whole lineage. Hypermutators are marked by purple triangles

Fig. 5

Mechanistic links, gyrA/gyrB/nfxB and retS/gacAS/rsmA. We use AA to illustrate phenotypic separation by isolates affected by distinct mutations in ciprofloxacin resistance genes gyrA, gyrB, and nfxB and the retS/gacAS/rsmA regulatory system. (a, b, left panel) As visualized by AA simplex plots, the diversity of trait profiles associated with isolates with mutations in DNA gyrase (gyrA/B) is in stark contrast to the constrained band of nfxB-mutated isolates. Mutations in DNA gyrase and nfxB do not co-occur in the same isolate but co-occur in different isolates of two lineages (patient P8804, genotype DK08 and patient P8203, genotype DK32). The differences in time of appearance during the colonization period and persistence of gyrA/B mutant isolates versus nfxB mutant isolates is shown in the lineage timelines plotted in the right column for gyrA/B (a, right panel) versus nfxB (b, right panel). Furthermore, gyrB-mutated isolates cluster more closely with A2 and A4 than gyrA mutated isolates, indicating a potential association with adhesion; GAMM predicts that gyrB mutation has a significant impact on adhesion (p-value « 0.01, GAMM with Wald-type tests). (c, left panel) Mutations in the retS/gacAS/rsmA system show a clear phenotypic change when retS is mutated alone (blue circles) or in combination with gacA or gacS (red squares and circles). The associated lineage plot (c, right panel) shows the appearance of double mutations (retS + gacA/S) after a colonization period by retS mutated isolates in three patient lineages. (a/b/c – lineage plot notation) Multiple isolates may be collected at the same sampling date based on differences in colony morphology or collected from different sinuses at sinus surgery, which explains the vertical overlap of isolates for some lineages. Lineage length is based on the span of time for which we have collected isolates and is indicated by gray bracketed lines, with only isolates affected by a mutation of the gene of interest plotted using shape and color (see legend of simplex plots). If multiple unique mutations occur in a lineage, this is specified by differential shading. (a/b only) Symbol size indicates the level of resistance to ciprofloxacin

Fig. 6

Evolutionary trajectories guided by different adaptation objectives. We present four different trajectories showing modes of evolution found in multiple patients: a Convergent evolution driven primarily by changes of a single phenotypic trait (decreased ciprofloxacin sensitivity). b Directed diversity with early/naive isolates showing a population moving in a broad and diverse plane from naive archetypes towards evolved archetypes. c General diversity where the population has no clear evolutionary trajectory. d A special case of convergent evolution with one outlier isolate (isolate 96 of DK12) but an otherwise clear trajectory first towards ciprofloxacin resistance and afterwards a gain in adhesive capabilities