Bacterial c

Nature Microbiology (2023)Cite this article

34 Altmetric

Metrics details

Most microbes evolve faster than their hosts and should therefore drive evolution of host–microbe interactions. However, relatively little is known about the characteristics that define the adaptive path of microbes to host association. Here we identified microbial traits that mediate adaptation to hosts by experimentally evolving the free-living bacterium Pseudomonas lurida with the nematode Caenorhabditis elegans as its host. After ten passages, we repeatedly observed the evolution of beneficial host-specialist bacteria, with improved persistence in the nematode being associated with increased biofilm formation. Whole-genome sequencing revealed mutations that uniformly upregulate the bacterial second messenger, cyclic diguanylate (c-di-GMP). We subsequently generated mutants with upregulated c-di-GMP in different Pseudomonas strains and species, which consistently increased host association. Comparison of pseudomonad genomes from various environments revealed that c-di-GMP underlies adaptation to a variety of hosts, from plants to humans. This study indicates that c-di-GMP is fundamental for establishing host association.

Host-associated microorganisms have important effects on the physiological functioning and fitness of their plant and animal hosts1,2,3. These host–microbiota interactions are often studied using a host-centric view, with a focus on microbiota-mediated host functions. This view neglects the important fact that most microbes evolve faster than their hosts due to their shorter generation times and higher mutation rates, and thus that fitness improvements for the microbes may disproportionately drive the associations4. An important step in the evolution of a host–microbe association is the emergence of a more specialized interaction that allows free-living bacteria to reliably enter the host, persist and finally be released into the environment to colonize new hosts (Fig. 1a)4. Thus far, little is known about the traits and molecular processes that determine how bacteria adapt to such an association with the host.

a, Host-associated microbes transition from a free-living phase to host association, the latter comprising host entry, persistence and release. Six P. lurida populations were passaged ten times across these stages with the host C. elegans (EVOhost) or without the host as a control (EVOctrl). b, Host-adapted bacterial populations significantly increased fitness (given as c.f.u.s per worm) relative to the ancestor (two-sided t-tests and FDR-corrected Tukey post hoc comparisons, 6 replicates per treatment). c, Evolved bacteria remain beneficial to the host, determined by nematode population growth (two-sided t-test, 6 replicates per treatment). d, A wrinkly colony morphotype only emerged during host adaptation and dominates within worms (comparison of morphotype abundance within treatments: generalized linear models and FDR-corrected Tukey post hoc test, 6 replicates per treatment). e, Evolved host-specialist bacteria (tagged with red fluorescent dTomato) colonize the worm gut (autofluorescent vesicles in nematode intestinal cells in cyan). Scale bar, 10 µm. f, PCA of key traits of the host-associated life cycle (in a) reveals a distinct profile for evolved wrinkly host specialists compared with ancestral bacteria (see Supplementary Table 5 for individual traits measurements). b–d, Boxplots show median (centre line), upper and lower quartiles (box limits) and the interquartile range (whiskers).

We studied the evolutionary transition from free-living to host association through controlled experimental evolution, using the bacterium Pseudomonas lurida and the nematode host Caenorhabditis elegans as a model. This bacterium is occasionally found in the natural microbiota of C. elegans5,6. Under laboratory conditions, the presence of P. lurida is associated with increased population growth rates of C. elegans and can provide protection against pathogens, yet both host and bacterium can proliferate without each other and thus do not depend upon one another5,7,8. To select host-adapted bacteria, we serially passaged 6 P. lurida populations either with or without a host-associated phase (Fig. 1a; EVOhost or EVOctrl, respectively). All populations were inoculated from the same clonal ancestor. After 10 passages through hosts, the bacteria reached on average 5–10 times higher bacterial load in the host than their ancestor, a significant change not observed for the control that evolved without exposure to hosts but otherwise had identical conditions (Fig. 1b and Extended Data Fig. 1). The increased bacterial fitness did not come at a cost to the host, as nematode population growth (used as a proxy for nematode fitness5) did not change significantly, but rather increased in the presence of the adapted bacteria (Fig. 1c).

As a result of passaging, bacterial populations diversified in colony morphology. At the end of our experiment a ‘wrinkly’ morphotype was dominant in all host-associated experimental replicate populations and absent in the controls, whereas ‘fuzzy’ and ‘smooth’ (ancestral) morphotypes were present across treatments (Fig. 1d and Supplementary Table 1). Despite their significant advantage in hosts, the wrinkly morphotypes declined during growth on agar, while smooth and fuzzy types increased in abundance (Extended Data Fig. 1 and Supplementary Table 2). As the wrinkly types were unique to and reached very high abundance in worm-adapted bacteria, we considered them host specialists. These specialists can be found in clusters within the intestinal tract of the nematode, especially in the anterior and posterior parts (Fig. 1e and Extended Data Fig. 2). Notably, the evolved wrinkly morphotype is similar to wrinkly P. fluorescens that emerge at the air–liquid interface in static microcosms9 and to rugose variants of various pathogenic bacteria10,11,12. Our experiments suggest that this morphological change also occurs in beneficial bacteria adapting to host association. For a further characterization of these adaptations, we focused on 47 clones of the distinct and genetically stable morphotypes (Supplementary Table 3) isolated from the final populations of our evolution experiment.

An analysis of trait changes across the distinct stages of host association revealed specific adaptations of wrinkly morphotypes to the interaction with C. elegans. In detail, we characterized two traits of importance for the free-living stage and four traits for host association (as listed in Fig. 1a). We found that the wrinkly isolate profiles were significantly distinct from the ancestral trait profile (Fig. 1f and Supplementary Table 4). This was mainly due to significant increases in short-term persistence, release from the host and in vitro biofilm formation (Fig. 1f, Extended Data Fig. 3 and Supplementary Tables 4–6)—all traits that define late-phase interactions with the host. The overall pattern of improved host association was also recovered by analysing the genetically diverse populations from the end of the evolution experiment, where the host-associated populations similarly increased in persistence and release (Extended Data Fig. 4 and Supplementary Tables 7 and 8). In detail, biofilm formation can enable persistent contact with the host and increase stress tolerance13,14, as exemplified by many pathogens15, thereby improving survival in the nematode’s digestive tract. As a consequence of increased biofilm formation, aggregated cells may be expelled more easily16, thereby explaining the observed increase in release. Such shedding also enhances the chance for transmission to other hosts4, which restarts the cycle of host association. Notably, wrinkly isolates did not differ from ancestors in early colonization, yet showed a significant decrease in colony expansion and swarming on plates (Fig. 1f, Extended Data Fig. 3 and Supplementary Tables 4 and 5). The latter result is consistent with a decrease in motility described for E. coli that evolved to become a mutualist in stinkbugs17, but contrasts with findings that sufficient swarming is required for colonization initiation of zebrafish and bobtail squid18,19. These contrasts are probably due to differences in symbiont recruitment between the host systems, defined by either aquatic environments for zebrafish and squid, or terrestrial environments for C. elegans and stinkbug. Moreover, our observations of increased biofilm formation and reduced motility may indicate an evolved life-history trade-off between the traits defining host association and the free-living stage. We conclude that experimental evolution in the presence of the nematode host leads to the emergence and spread of a host-specialist type. We next asked whether the improved host association has a common genetic basis.

Whole-genome sequencing of the isolated morphotypes and the ancestor revealed several independent mutations in wrinkly host specialists that affect the bacterial second messenger cyclic diguanylate (c-di-GMP). In particular, a comparison of non-silent genomic variation identified variant genes specific to wrinkly host specialists (Fig. 2a and Supplementary Table 9). Two of the genes, wspE and wspF, code for a hybrid sensor histidine kinase and a methylesterase in the wrinkly spreader (wsp) operon, respectively20. These genes are part of a two-component system that regulates c-di-GMP levels (Fig. 2g) and wrinkly formation in beta- and gamma-proteobacteria, including pseudomonads20,21,22,23. We found additional mutations unique to the host specialists in the gene rph, encoding RNase PH that has not been linked to c-di-GMP signalling previously. Using both a fluorescence-based c-di-GMP sensor and liquid chromatography–mass spectrometry (LC–MS), we found a roughly twofold c-di-GMP increase in three wrinkly isolates, each with a single mutation in either wspE, wspF or rph, when compared with the ancestor (Fig. 2b, Extended Data Fig. 5 and Supplementary Table 10). This points to a loss-of-function mutation in wspF (which downregulates c-di-GMP) and alterations in active sites of WspE and Rph that all converge at upregulating c-di-GMP. We aligned evolved and ancestral amino acid sequences (Extended Data Fig. 6) and confirmed a disruption in WspF functional domains, as well as a disrupted receiver domain in WspE that probably prevents its de-autophosphorylation and thus constantly activates downstream WspR24. Amino acid substitutions in the exoribonuclease domain of Rph further link its ribonuclease activity to c-di-GMP metabolism. As we observed similar increases in c-di-GMP levels in other wrinkly, but not in smooth or fuzzy mutants (Extended Data Fig. 5), we subsequently asked whether the wrinkly-specific mutations indeed cause improved host association.

a, Overview of genes with non-silent changes in evolved bacterial isolates. Data points represent mutant isolates with one or multiple mutations in a given gene (the total number of isolates with a given morphology is specified in brackets). A cross indicates genes with variants lost in the evolved isolates as compared with the ancestor. b, Fluorescence sensor and LC–MS detected higher intracellular c-di-GMP concentrations in evolved wrinkly mutants compared with the ancestor (Welch’s ANOVA and Games–Howell post hoc comparisons, 5 replicates per treatment). Scale bars, 10 µm. c, Competitive fitness (c.f.u.s per worm relative to ancestor, dashed line) of evolved wspF, wspE and rph mutants (left), rescued mutants (middle) and reconstructed mutants in ancestral background (right) during persistence in C. elegans MY316 (3 < n < 5). d, Competitive fitness of evolved wspF, wspE and rph mutants during persistence in the non-native C. elegans strain N2 (8 < n < 10). e, Competitive fitness of wspF, wspE and rph mutants (5 < n < 6) with additional ΔwspR mutation compared to ancestral Pl_MYb11 (dashed line). f, Competitive fitness of wspF, wspE and rph expressing heterologous phosphodiesterase (PDE) PA2133 from plasmid (pJN2133) and ancestral MYb11 expressing constitutively active diguanylate cyclase (DGC) GCN4-WspR from plasmid (pJStrep-GCN4-WspR), each compared to their respective empty vector control (n = 4, dashed line). c–f, Persistence competition experiments were performed with ≥3 replicates per treatment and analysed with ANOVA (d) or LMM and FDR-corrected Dunnett post hoc tests; *P < 0.05, **P < 0.01, ***P < 0.001. Boxplots show median (centre line), upper and lower quartiles (box limits) and the interquartile range (whiskers). g, Graphical hypothesis of adaptive c-di-GMP manipulation via Rph and the Wsp system. Solid lines indicate previously established regulatory interactions, dashed lines emerging hypotheses. Red indicates the inferred consequences of the studied mutated gene from experimental evolution.

A functional genetic analysis of wspE, wspF and rph demonstrated their direct involvement in host adaptation. For this analysis, we assessed the competitive fitness of mutants relative to the ancestor during host colonization. First, we re-assessed the three selected wrinkly mutants and found them to be significantly more competitive than the ancestor (Fig. 2c, left panel, and Supplementary Table 11), alongside increased biofilm formation and decreased swarming in vitro (Extended Data Fig. 7 and Supplementary Table 12). Thereafter, we rescued these mutants with the corresponding ancestral alleles, which indeed abolished the mutants’ fitness increase (Fig. 2c, middle panel, and Supplementary Table 11). Thirdly, an experimental introduction of each mutation into the ancestral background resulted in a significantly higher competitiveness, at least for the wspF and rph mutations (Fig. 2c, right panel, and Supplementary Table 11). A similar fitness advantage was observed for the wspE and wspF mutants when either was subjected to quartet competition with the ancestor and the two other morphotypes (Extended Data Fig. 8 and Supplementary Table 13). Notably, fitness advantages of evolved mutants were consistently observed in a non-native host strain (the C. elegans laboratory strain N2) (Fig. 2d and Supplementary Data Table 14). While two of these genes are components of the Wsp system, which regulates c-di-GMP during surface sensing in other pseudomonads25,26, Pl_MYb11 in theory possesses a variety of c-di-GMP modifying enzymes. This includes 34 genes coding for GGDEF and 22 coding for EAL domains with putative diguanylate cyclase (DGC) and c-di-GMP-specific phosphodiesterase (PDE) functions, respectively. We validated the role of the Wsp system’s cognate DGC in host adaptation using wspR knockouts in our evolved host-specialist mutants. This change abolished the mutants’ competitive advantage in the host (Fig. 2e) and caused a change from wrinkly to smooth colony morphology (Extended Data Fig. 7 and Supplementary Table 15), thus linking the DGC wspR to wspE and wspF (as expected) and rph (previously unknown). In addition, we directly manipulated c-di-GMP levels by heterologous expression of a PDE and a DGC from P. aeruginosa23,27, which respectively resulted in either decreased or improved persistence in C. elegans, as expected (Fig. 2f, Extended Data Fig. 7 and Supplementary Table 16). We thus conclude that changes in wspE, wspF and rph that converge on increasing c-di-GMP levels via the Wsp system enhance bacterial fitness in the host (Fig. 2g). As upregulation of this second messenger mediates a fundamental life-history switch13, we next investigated whether it more generally mediates host association across pseudomonads.

Genetic manipulation of wspF and a bioinformatic analysis of Pseudomonas genomes revealed a general involvement of wsp genes in host association. For the former, we generated wspF deletion mutants for P. lurida strain MYb193 and the distantly related P. alkylphenolia MYb187 (both naturally associated with C. elegans), and further obtained mutant and wildtype P. fluorescens strain SBW25, a model for wrinkly formation21. We found that the mutants had significantly higher competitive fitness in the C. elegans host than their respective wildtypes (Fig. 3a and Supplementary Table 17). Furthermore, we correlated the presence of wsp and rph genes in 1,359 whole Pseudomonas genomes from NCBI with the bacterial isolation source, a proxy for lifestyle (Extended Data Fig. 9 and Supplementary Table 18). Pseudomonas isolates containing any of the wsp genes or the complete, highly syntenic (Supplementary Table 19)20 wsp operon were significantly more often isolated from a host than isolates lacking these genes (Fig. 3b and Supplementary Table 19). These findings may seem surprising for genes with opposite regulatory effects (for example, wspE versus wspF), yet are probably explained by the syntenic inheritance of the entire operon with its set of interacting genes (Supplementary Table 19; see also ref. 20). Further, rph was more prevalent in isolates from healthy/undiagnosed hosts than from diseased hosts. Across lifestyles, we additionally detected signatures of negative selection for wspE, wspF and rph, which additionally suggest that they are functionally stabilized by selection when present (Supplementary Table 20). We propose that the presence of these genes allows the finetuned regulation of c-di-GMP and thereby, adjustment to a host-associated lifestyle.

a, WspF deletion increases intra-host competitive fitness of different Pseudomonas species relative to their wild type (c.f.u.s per worm, wildtype indicated by dashed line; LMM and FDR-corrected Dunnett post hoc tests, 3 < n < 5,). Boxplots show median (centre line), upper and lower quartiles (box limits) and the interquartile range (whiskers). b, Presence/absence of wsp genes and rph, as well as completeness of the wsp operon (that is, presence of wspA-F and wspR) co-vary with isolation source of sequenced pseudomonads (Χ2 goodness-of-fit tests).

Together, our study demonstrates that bacteria can improve their association with a host by shifting their life history from a motile to a sessile, persisting lifestyle. This lifestyle shift results from correlated changes in a suite of life-history traits (Fig. 1f), which together represent a transition in life-history strategy. One way to interpret this transition is as a shift along the r–K life-history continuum, from an r-like strategy characterized by high reproductive rates to a K-like strategy characterized by persistence under high density conditions28,29. To demonstrate whether such a transition would generally lead to increased host association, we used an extension of a previously published mathematical model of microbial evolution towards host association30. Exploration of a broad parameter space with this model confirmed that increased within-host persistence is often the optimal strategy for microbial adaptation to hosts (Extended Data Fig. 10 and Supplementary Discussion), suggesting that the results from our study may be generally applicable.

In our experiments, the lifestyle shift from primarily free-living to host-associated is mediated by the Wsp system and subsequently, activity of the bacterial second messenger c-di-GMP. C-di-GMP is well known to regulate key physiological functions in bacteria, including the regulation of virulence in bacterial pathogens22,31. Our work demonstrates that this regulatory system promotes the adaptation of pseudomonads to diverse host systems, from plants to humans, not only in pathogens but extending to beneficial host–bacterial relationships. Given the importance of beneficial microorganisms in the functioning of their hosts, understanding the mechanisms that mediate non-pathogenic associations is crucial. Our study suggests that c-di-GMP plays an essential role in many such associations.

We performed evolution experiments with P. lurida strain MYb11 (Pl_MYb11) and its natural host C. elegans strain MY316 (Ce_MY316) (ref. 5). In preparation for all experiments, we thawed frozen worm stocks (−80 °C) and raised worms on nematode growth medium agar (NGM32) seeded with E. coli OP50. In additional persistence colonization experiments, we used the standard laboratory strain C. elegans N2 as a non-native host for the evolved bacteria. A standard bleaching protocol was used to collect sterile and synchronized L1 larvae, which were then raised to L4 stage on E. coli OP50 (20 °C), unless stated otherwise.

P. lurida strains MYb11 and MYb193, and P. alkylphenolia MYb187 were isolated from Ce_MY316 (ref. 5), and P. fluorescens SBW25 from sugar beet leaves9. Bacteria were cultured on tryptic soy agar (20 °C, 48 h) and tryptic soy broth (28 °C, 150 r.p.m., overnight) unless stated otherwise.

Bacterial populations originating from a clone of Pl_MYb11 were serially passaged on NGM in the presence of Ce_MY316 (host treatment, 6 replicates) or without worms (negative control, 6 replicates). For each replicate, a lawn of Pl_MYb11 was seeded onto NGM and cultured for 3.5 d. For each cycle of the host treatment, 10 C. elegans L4 larvae were added per plate and incubated until the worms reached the F1 generation (3.5 d). In the negative controls, bacteria were maintained on NGM without worms. At the end of every cycle, bacteria were collected from either worms or plates in the host-associated and control treatments, respectively, 10% of the population (bottleneck) was transferred to the next cycle and a sample frozen (−80 °C). A similar number of colony-forming units (c.f.u.) was used to bottleneck the negative control. A total of 10 cycles were performed.

Frozen bacteria from cycle 10 were recovered and before further experiments were conducted, these were subjected to one more cycle of the evolution experiment to minimize any potential selective effects of freezing/thawing. To focus on evolved differences between populations of the host treatment and the negative control, rather than physiological responses to recent host exposure, bacteria were grown on NGM for 2 d as a common garden treatment and then used in subsequent assays.

Bacterial fitness during host association was quantified as c.f.u.s per worm. In preparation, bacterial lawns (125 µl, optical density (OD)600 = 2) were seeded on NGM and 5 synchronized L4 Ce_MY316 added. After 3.5 d at 20 °C, worms were collected with M9 buffer containing 0.025% Triton-100 and 25 mM of the paralyzing antihelminthic tetramisole. The worms were washed in buffer using a custom-made filter tip washing system33 and collected in M9 with Triton-100. Worm-free supernatant was collected as a background sample. Following homogenization by bead beating, serial dilution and plating were used to quantify c.f.u.s. C.f.u.s per worm was calculated as the difference in c.f.u. between worm and supernatant samples, divided by the number of worms per population. For diversified populations, colony morphologies were scored as smooth, fuzzy or wrinkly.

Worm population growth resulting from 5 L4 larvae over 3.5 d was quantified as a proxy for host fitness. Bacteria and worms were prepared as for colonization assays and washed worms frozen in 48-well plates. Photographs of worms were automatically scored in ImageJ2 (ref. 34): worms were detected as particles, approximated by ellipses, and those fitting C. elegans-like dimensions (major axis 0.18–1.3 mm, minor axis ≤0.1 mm (ref. 35)) were counted. Detection quality was validated by correlating automatic worm counts with counts of two independent experimenters (r(58) = 0.736, P = 2.106 × 10−11).

To quantify early colonization, persistence and release from L4 stage worms, bacterial lawns were prepared from ancestral Pl_MYb11 and evolved populations (post common garden) or clonal morphotypes (overnight cultures). In early colonization assays, we quantified bacteria that entered L4 Ce_MY316 that were previously raised on non-colonizing E. coli O50. Colonization levels were then assayed as above resulting in c.f.u.s per worm as a measure of early colonization.

For persistence and release assays, worms were raised on the respective assay bacteria (from L1 until L4 stage), mimicking the development of worms in the F1 generation of the evolution experiment. Worms were then collected, washed using the filter tip washing system and samples divided into supernatant (supernatant 1) and worm sample (100 µl each). Worms were then suspended in 200 µl M9 and incubated for 1 h, after which 100 µl supernatant containing released bacteria (supernatant 2) was collected. The c.f.u.s released per worm were determined by the difference in c.f.u.s between supernatant 2 and supernatant 1. Along with this, we quantified c.f.u.s maintained in worms of this sample as a measure of persistence.

To measure bacterial growth, bacterial populations (common garden treatment or overnight cultures) were adjusted to OD600 = 0.1 and 50 µl spotted on NGM. After incubation (24 h or 3 d at 20 °C), lawns were scraped off, homogenized and c.f.u.s determined by serial dilution.

Colony expansion and swarming were assayed on NGM containing 0.5% or 3.4% agar, respectively. In either case, 0.5 µl of cell suspension (OD600 = 1) was spotted on surface-dried agar plates. Colony diameter was measured after 24 h, 3 and 7 d.

In vitro biofilm formation was assayed in microtitre plates as described previously36. Notably, assays were performed in a randomized layout in Nunclon Delta surface-treated plates. Staining was performed after 48 h of incubation (20 °C, orbital shaking at 180 r.p.m.). Absorption of dyed biofilm solutions was measured at 550 nm using Gen5 microplate reader and Imager software (Biotek, v.3.08.01). To illustrate biofilm formation in liquid, glass test tubes were filled with 2 ml tryptic soy broth, inoculated with single colonies of ancestral Pl_MYb11 or evolved host-specialist mutants (wspE, wspF, rph) and incubated at 20 °C for 48 h until photographing.

Representative colonies with visually distinct morphologies were isolated from evolved cycle 10 populations. The evolved populations were thawed, serially diluted and plated (48 h, 20 °C). Unique morphotypes from all evolved populations were re-streaked and archived as frozen stocks (Supplementary Table 3). All morphotypes were thawed and re-streaked once, and showed stable colony morphology during 2 d of incubation.

Macrocolonies of Pl_MYb11 morphotypes and mutants were prepared as described previously37. Briefly, 5 µl of overnight culture were spotted on tryptic soy agar plates supplemented with 40 μg ml−1 Congo Red and incubated at 20 °C. After 24 h or 48 h, photographs were taken using a Leica fluorescence dissecting scope (LEICA M205 FA).

The wrinkly morphotype MT12 was labelled with red fluorescent dTomato (dT) using Tn7 transposon-based chromosomal insertion as previously described38,39. Insertion of the label did not affect the wrinkly morphology of the colonies.

Fluorescently labelled MT12 was used to localize colonization in Ce_MY316 using confocal laser scanning microscopy (ZEISS LSM 880). For this, synchronized L1 stage larvae were exposed to labelled bacteria for 72 h (20 °C), then collected using gravity washing and mounted for microscopy as previously described39. Overviews of complete worms were created using a ×25 LD LCI Plan-Apochromat multi-immersion objective (numerical aperture (NA) = 0.8) and details imaged using a ×40 C-Apochromat water immersion objective (NA = 1.2), in both cases using Immersol W (2010) with a refractive index of 1.334. Bacterial fluorescence and worm autofluorescence were sequentially excited (561 nm and 488 nm) and detected with an Airyscan detector (R-S sensitivity mode; longpass filter ≥570 nm; bandpass filter 495–550 nm). Data were processed with the automatic Airyscan processing function of ZEISS Efficient Navigation 2. For a list of the genetically modified bacteria used in this study, see Supplementary Table 21. After looking at the colonization of >10 worms in at least 3 biological replicate populations of consecutive weeks of experiments, a representative worm was imaged for Fig. 1e and Extended Data Fig. 2.

Total DNA was isolated using a cetyl-trimethylammonium-bromid-based protocol40. For Illumina MiSeq (paired-end, 300 bp) sequencing, libraries were prepared using the Nextera DNA Flex kit. Read quality was inspected using FastQC (v.0.11.8) (ref. 41) and reads trimmed using Trimmomatic (v.0.3.9) (ref. 42). Paired reads were aligned to the Pl_MYb11 reference genome (RefSeq: GCF_002966835.1; Bowtie2 v.2.3.3 (ref. 43)) and duplicate regions removed using Picardtools (v.2.22.2) (ref. 44). Variants were called using BCFtools (v.1.10.2) (ref. 45) and VarScan (v.2.3.9) (ref. 46), and then annotated (snpEff47,48). We filtered for non-synonymous variants not present in the ancestral control in R49,50. Gene ontology was inferred using Pseudomonas.com51. To infer genes coding for enzymes with putative DGC or PDE activity, we searched for proteins with GGDEF and EAL domains using the InterProScan of the conserved domains database (CDD) via Pseudomonas.com51.

To prepare amino acid sequence alignments of ancestral and mutated WspE, WspF and RPH, nucleotide sequences were translated using EMBOSS Transeq52 (frame 1; bacterial codon table; forward for wspE and wspF, reverse for rph) and resulting amino acid sequences aligned using Clustal Omega52 (v.1.2.4; ClustalW with character counts and standard settings). For annotation and visualization of protein domains, domain predictions of the respective sequences were collected from Pfam/InterPro (sourced from Pseudomonas.com51) and visually highlighted in protein visualizations prepared with DOG (v.2.0)53.

To quantify intracellular concentrations of c-di-GMP in ancestral Pl_MYb11 and evolved wrinkly isolates (MT12: wspFEVO, MT14: wspEEVO and MT22: rphEVO), we used an established plasmid-based biosensor54. Bacterial strains carrying the plasmid were grown on gentamicin-selective plates (70 h, 20 °C). For microscopy, single colonies were resuspended in 1X PBS, spotted on 2% agarose patches on microscopy slides and sealed.

Bacterial fluorescence was visualized using confocal laser scanning microscopy (ZEISS LSM 700 with ×40 Plan-Apochromat oil immersion objective (NA = 1.4) and Immersol 518F with a refractive index of 1.518). Fluorescence of the sensor and normalizer were sequentially excited (555 nm and 488 nm) and detected with a photomultiplier tube detector and a variable secondary dichroic transmitting light with wavelengths ≤630 nm and ≤550 nm, respectively. The excitation and detection settings were kept identical across all measurements.

Fluorescence intensity per cell was measured in Image J34: all cells and five background areas were identified as regions of interest, and area, integrated density and mean grey values were measured. Data from the untransformed images were used to calculate the corrected total cell fluorescence55.

In addition to single-cell measurements, we quantified c-di-GMP at the population level. For this, colonies of evolved wrinkly (MT12, MT21, MT25, MT26), smooth (MT13, MT33) and fuzzy (MT11) isolates were grown as described above, resuspended in 1X PBS and adjusted to OD600 = 0.1. Cell suspensions (200 µl) were then transferred to black, flat-bottomed 96-well plates with transparent bottoms (Greiner Bio-One CELLSTAR 96-well, cell culture-treated) in triplicate. After shaking (10 s, orbital shaking at 1 mm amplitude), fluorescence was sequentially excited (454 nm and 460 nm, bandwidth 9 nm, 10 flashes) and emission detected (585 nm and 510 nm, bandwidth 20 nm; optimal gain, 20 µs integration) in a plate reader (Tecan, Infinite M200Pro), with 1X PBS serving as the background control.

To infer c-di-GMP concentration, we calculated the relative fluorescence intensity, or the ratio between TurboRFP and AmCyan fluorescence intensities, as previously described54, and compared average relative fluorescence intensities between ancestral Pl_MYb11 and evolved wrinkly, smooth and fuzzy morphotypes. For the images used in Fig. 2, linear LUT was used at full range. Brightness and contrast were applied equally to all images.

To quantify intracellular c-di-GMP using LC–MS in parallel reaction monitoring mode, ancestral and evolved Pl_MYb11 (MT12, MT14 and MT22) were grown in LB medium to an OD600 of 1.8 and pelleted by centrifugation. After washing with salt-free LB medium, pelleted cells were snap frozen and stored (−80 °C). Cells were mixed with 10 pmol of internal standard (cyclic-di-GMP-13C20,15N10, Toronto Research Chemicals) in 60 μl of water. Extraction of c-di-GMP was performed as previously described56 with the following modifications: extraction solution (240 μl of 1:1 acetonitrile (ACN)/methanol (MeOH)) was added and samples were vigorously vortexed. Following incubation on ice (15 min) and centrifugation (20,800 × g, 4 °C, 2 min), extract supernatant was collected and solvent extraction repeated twice (200 μl of 2:2:1 ACN/MeOH/water). Pooled extracts were dried, resuspended in 50 μl of water and centrifuged to remove insoluble compounds. Concentrations of solubilized protein precipitates were determined using the Pierce BCA protein assay kit (Thermo Fisher). For LC–MS/MS, 1 μl extract was injected onto an EASY-nLC 1000 UHPLC (Thermo Fisher) and separated on a 15-cm ReproSil-Pur C18-AQ nano LC column (0.1 mm i.d., 1.9 μm, 120 Å, Altmann Analytik) at 400 nl min−1. Eluent A was 10 mM NH4OAc with 0.1% HAc, eluent B was 100% MeOH. Chromatographic conditions were 5% eluent B (5 min), followed by a linear gradient from 5% to 20% B (15 min) and an increase to 70% B (1 min), followed by 70% B (5 min) and 5% B (5 min); higher-energy collisional dissociation of the m/z 691.1021 and m/z 721.0714 precursors was performed on a Q Exactive HF Orbitrap MS (Thermo Fisher). Peak areas for the qualifying57 product ions m/z 248.0778 (light) and m/z 263.0965 (heavy) determined in Skyline (v.21.1.0.146.3, MacCoss Lab software)58 were used to calculate total c-di-GMP amounts, which were normalized to total protein amount as obtained by the BCA assay.

A two-step allelic replacement method based on previously described protocols21,59 was used to introduce the evolved mutant alleles into an ancestral background and also to revert mutations by introducing ancestral alleles in the mutant background. We applied the following modifications: ~700 bp long PCR amplicons surrounding each mutation were cloned into pUISacB allowing for sucrose selection. The constructs were transformed into competent E. coli cells and transferred to Pseudomonas isolates via conjugative mating with an E. coli helper strain containing pRK2013 (ref. 60). Primers (see Supplementary Table 22) were designed using NCBI’s BLAST tool61 and NCBI Primer-BLAST62, NEBuilder v.2.3.0 (New England Biolabs) and Oligo Analyse Tool (Eurofins Genomics). BLASTn and alignments with Clustal Omega63 were performed using default settings.

To manipulate intracellular c-di-GMP levels and study the consequences for host association and colony morphology, we expressed a heterologous PDE and a heterologous DGC in our evolved host-specialist mutants (wspE, wspF and rph). The PDE PA2133 from P. aeruginosa was expressed from plasmid pJN2133 (ref. 23). A constitutively active GCN4-WspR fusion construct27 was synthesized (Eurofins) and then cloned into pJStrep to generate a C-terminal StrepII-tagged GCN4-WspR construct. Empty pJStrep (a modified pJN105 (ref. 64) vector containing the StrepII tag coding sequence) and empty pJN105 plasmids were used as controls. All plasmids were introduced into Pl_MYb11 and evolved mutants using a previously described electroporation protocol65.

Competition experiments were performed as described for the short-term persistence assays. Co-inoculated bacteria were OD-adjusted and mixed in equal volumes before seeding as lawns on NGM agar. A Pl_MYb11 labelled with dTomato39 was used, which is equivalent to the ancestral Pl_MYb11, as no differences were observed in short-term persistence (analysis of variance (ANOVA), F value = 0.99, d.f. = 1, P = 0.35). C.f.u.s per worm were determined by subtracting c.f.u.s in supernatants from those in worm samples. A competitive index was calculated as the ratio of c.f.u.s per worm of evolved or constructed mutants to c.f.u.s per worm of the ancestor.

Whole-genome sequences from NCBI were mined for c-di-GMP modulating genes (focus: wsp operon, rph) with bacterial lifestyle in members of the genus Pseudomonas. First, candidate genomes were obtained (NCBI Nucleotide’s command line search tool; size: 5–8 million bp). This retrieved 2,279 sequences, for which sample information from NCBI’s Biosample database was collected. When available, host, host disease status, isolation source and sample type were used to manually classify genomes as originating from free-living or host-associated isolates with or without/unknown disease (Supplementary Table 18). Next, we downloaded all available Pseudomonas reference sequences for rph and wsp genes from pseudomonas.com51. These were used to identify candidate sequences of rph, wspA, wspB, wspC, wspD, wspE, wspF and wspR. These target gene candidates were found in the selected genomes using BLAST (R package ‘rBLAST’) and filtered on the basis of sequence lengths and percent identities of the BLAST hits (Extended Data Fig. 9). Percent identity and sequence length were selected to maximize the chance that genes were correctly identified (red rectangles in Extended Data Fig. 9). If at least one candidate gene was identified during BLAST searches with the reference genes as query, this gene was considered present in the respective genome. We then used χ2 goodness-of-fit tests to infer whether isolates with and without the target genes differed in the relative proportions of host-associated lifestyles (Supplementary Table 19).

To assess whether our focal host-specialist genes (wspE, wspF, rph) were experiencing positive or purifying selection in the genus Pseudomonas, we performed MUSCLE codon-based multiple sequence alignments of nucleotide sequences (see dataset described above) using MEGA11 (ref. 66; default settings). Subsequently, we performed codon-based z-tests (default settings) to test for significant deviations from neutral selection. In addition, we analysed signatures of selection in Blast hits for a set of three Pseudomonas core genes (gyrB: PA0004, rpoD: PA0576 and a 16S rRNA methyltransferase: PA0419; see also ref. 67) in the set of genomes studied for wsp and rph presence/absence, also using multiple sequence alignments and codon-based tests of neutrality.

Before data collection, no statistical methods were used to pre-determine samples sizes, but our sample sizes are similar to those reported in previous publications. In all experiments, treatments and samples were blinded and randomized. Before data analysis, assumptions of parametric models (normality, homogeneity of variances) were checked by visual inspection (box-/qqplots) and with Shapiro–Wilk and Levene tests. When these were not met, non-parametric tests were applied. Boxplots show median (centre line), upper/lower quartiles (box limits) and 1.5× interquartile ranges (whiskers).

To check whether evolved populations differed from the ancestor in c.f.u.s per worm, we compared the shift in the evolved phenotype (ratio of c.f.u.s per worm of evolved populations to those of ancestral Pl_MYb11) to the ancestral phenotype using one-sample t-tests (alpha = 0.05, mu = 1) with false discovery rate (FDR68) correction for multiple testing. We applied this approach to analyse: bacterial colonization of individual worms, worm population growth, early colonization, persistence and release, colony expansion and swarming. To infer overall phenotypic shifts according to evolutionary treatment, a principal component analysis (PCA) including the assayed phenotypes was performed. We performed permutational analysis of variances (PERMANOVA, 1,000 permutations) followed by pairwise comparisons of groups (FDR-corrected) to test for differences in phenotype sets of ancestral and evolved groups, and plotted confidence ellipses (one standard deviation). Packages used included ggbiplot69, missMDA70, vegan71 and pairwise.adonis72.

Differences in proportions of the different colony morphologies (wrinkly, smooth and fuzzy) within worms were identified using generalized linear models (GLM; quasinormal distribution) with Tukey post hoc tests (using lme4 (ref. 73), lmtest74 and multcomp75).

Changes in morphotype proportions over time were tested using beta-regressions (using gamlss76).

Differences between morphotype phenotypes were detected using ANOVA or GLMs, followed by Tukey or Dunnett post hoc tests. To infer functional specializations across phenotypes, we used PCA and PERMANOVA.

Differences in biofilm formation and motility between evolved wrinkly host specialists (wspE, wspF and rph mutants) and ancestral Pl_MYb11 were analysed using nested ANOVAs followed by Tukey post hoc tests. When no batch effects (evolved population of origin) were detected, mutants were compared across populations. In the case of swarming diameters, however, the rph mutant was compared to the co-analysed ancestral Pl_MYb11 using a t-test.

Differences in c-di-GMP concentrations between evolved isolates were inferred using Welch’s ANOVA, nested ANOVA or ANOVA with Games–Howell or Dunnett post hoc comparisons.

We tested for differences in c.f.u.s per worm between morphotypes or mutants using GLMs and linear mixed models (LMMs), and Dunnett or Tukey post hoc comparisons.

All analyses and plotting were performed in R49,50,77,78.

We built a model to assess the selection gradient experienced by bacteria during the evolution experiment (Extended Data Fig. 10). We focused on the phase when bacteria are in contact with worms and considered a homogeneous population. The dynamics of the number of bacteria living in (any) host association n(t) can be described by the equation

where W(t) denotes the biomass of worms on the plate at time t. We consider that growing to saturation, bacteria on the plate are always in excess so that only the number of worms and the rate f at which they feed on bacteria limit the immigration of free-living bacteria to the host. We assume logistic growth of the bacterial population within the worms, with maximal rate r and a carrying capacity proportional to the biomass of worms W(t) and the per unit of worm biomass carrying capacity K. Finally, a fraction of the host-associated bacterial population is removed from the host at a rate δ, which encompasses bacterial death and expulsion to the environment. As in the evolution experiment, we assume that only host-associated bacteria are selected and continue to the next cycle, ignoring on-plate dynamics. We assume linear growth for the worm biomass, W(t) = g t + W0, encompassing both reproduction and development. We neglect the potential evolution of beneficial effects on worm growth and fix the parameters W0 = 10 and g = 711 d−1 to experimentally observed values.

We studied how the final number of host-associated bacteria, nf, is affected by changes in the parameters that describe the bacterial life cycle (r, δ, f, K). We defined a range of biologically plausible values for each of these parameters (that is, the trait space) that are informed by experimental data where possible:

10−1 d−1 < r < 101.25 d−1, that is, between a small fraction and around twice the maximum on-plate growth rate (~7 d−1).

10−0.5 d−1 < δ < 104 d−1, as the typical time for a worm to lose 50% of its microbiome (in the absence of feeding and replication) should range between seconds and days.

104 < K < 106.25, given the orders of magnitude from the maximal number of bacteria per worm measured experimentally (~105).

103 d−1 < f < 107.5 d−1, as the typical time for an empty worm to be colonized at 10% of its carrying capacity (K = 105) should vary between seconds and days, neglecting bacterial release and within-host replication.

For each point of the trait space, we numerically solved equation (1) to compute the expected final number of bacteria at tf = 3.5 d, nf = n(tf). Finally, we assessed the elasticity of nf along each direction of the trait space, which measures the expected relative change in nf with respect to a small relative change in one of the traits. We interpreted the vector of the elasticities as the selection gradient on the phenotypic traits79 and used the dominant element of this vector to define an ‘optimal evolution strategy’30 for each point of the trait space.

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.

Raw sequencing data are available at NCBI under Bioproject PRJNA862108. All other data are accessible at https://github.com/nobeng/c-di-GMP_host-association.

Custom code and associated data are available at https://github.com/nobeng/c-di-GMP_host-association.

McFall-Ngai, M. et al. Animals in a bacterial world, a new imperative for the life sciences. Proc. Natl Acad. Sci. USA 110, 3229–3236 (2013).

Article CAS PubMed PubMed Central Google Scholar

Toft, C. & Andersson, S. G. E. Evolutionary microbial genomics: insights into bacterial host adaptation. Nat. Rev. Genet. 11, 465–475 (2010).

Article CAS PubMed Google Scholar

Douglas, A. E. Fundamentals of Microbiome Science: How Microbes Shape Animal Biology (Princeton Univ. Press, 2018).

Obeng, N., Bansept, F., Sieber, M., Traulsen, A. & Schulenburg, H. Evolution of microbiota–host associations: the microbe’s perspective. Trends Microbiol. 29, 779–787 (2021).

Article CAS PubMed Google Scholar

Dirksen, P. et al. The native microbiome of the nematode Caenorhabditis elegans: gateway to a new host-microbiome model. BMC Biol. 14, 38 (2016).

Article PubMed PubMed Central Google Scholar

Johnke, J., Dirksen, P. & Schulenburg, H. Community assembly of the native C. elegans microbiome is influenced by time, substrate and individual bacterial taxa. Environ. Microbiol. 22, 1265–1279 (2020).

Article CAS PubMed Google Scholar

Kissoyan, K. A. B. et al. Natural C. elegans microbiota protects against infection via production of a cyclic lipopeptide of the viscosin group. Curr. Biol. 29, 1030–1037.e5 (2019).

Article CAS PubMed Google Scholar

Dirksen, P. et al. CeMbio—the C. elegans microbiome resource. G3 10, 3025–3039 (2020).

Article CAS PubMed PubMed Central Google Scholar

Rainey, P. B. & Travisano, M. Adaptive radiation in a heterogeneous environment. Nature 394, 69–72 (1998).

Article CAS PubMed Google Scholar

Starkey, M. et al. Pseudomonas aeruginosa rugose small-colony variants have adaptations that likely promote persistence in the cystic fibrosis lung. J. Bacteriol. 191, 3492–3503 (2009).

Article CAS PubMed PubMed Central Google Scholar

Anriany, Y. A., Weiner, R. M., Johnson, J. A., Rezende, C. E. D. & Joseph, S. W. Salmonella enterica serovar Typhimurium DT104 displays a rugose phenotype. Appl. Environ. Microbiol. 67, 4048–4056 (2001).

Article CAS PubMed PubMed Central Google Scholar

Yildiz, F. H. & Schoolnik, G. K. Vibrio cholerae O1 El Tor: identification of a gene cluster required for the rugose colony type, exopolysaccharide production, chlorine resistance, and biofilm formation. Proc. Natl Acad. Sci. USA 96, 4028–4033 (1999).

Article CAS PubMed PubMed Central Google Scholar

Hengge, R. Linking bacterial growth, survival, and multicellularity—small signaling molecules as triggers and drivers. Curr. Opin. Microbiol. 55, 57–66 (2020).

Article CAS PubMed Google Scholar

Pankey, M. S. et al. Host-selected mutations converging on a global regulator drive an adaptive leap towards symbiosis in bacteria. eLife 6, e24414 (2017).

Article Google Scholar

Hall-Stoodley, L., Costerton, J. W. & Stoodley, P. Bacterial biofilms: from the natural environment to infectious diseases. Nat. Rev. Microbiol. 2, 95–108 (2004).

Article CAS PubMed Google Scholar

Schlomann, B. H., Wiles, T. J., Wall, E. S., Guillemin, K. & Parthasarathy, R. Sublethal antibiotics collapse gut bacterial populations by enhancing aggregation and expulsion. Proc. Natl Acad. Sci. USA 116, 21392–21400 (2019).

Article CAS PubMed PubMed Central Google Scholar

Koga, R. et al. Single mutation makes Escherichia coli an insect mutualist. Nat. Microbiol. https://doi.org/10.1038/s41564-022-01179-9 (2022).

Robinson, C. D. et al. Host-emitted amino acid cues regulate bacterial chemokinesis to enhance colonization. Cell Host Microbe 29, 1221–1234.e8 (2021).

Article CAS PubMed PubMed Central Google Scholar

Isenberg, R. Y., Christensen, D. G., Visick, K. L. & Mandel, M. J. High levels of cyclic diguanylate interfere with beneficial bacterial colonization. mBio 0, e01671-22 (2022).

Google Scholar

Kessler, C., Mhatre, E., Cooper, V. & Kim, W. Evolutionary divergence of the Wsp signal transduction systems in Beta- and Gammaproteobacteria. Appl. Environ. Microbiol. https://doi.org/10.1128/AEM.01306-21 (2021).

Article PubMed PubMed Central Google Scholar

Bantinaki, E. et al. Adaptive divergence in experimental populations of Pseudomonas fluorescens. III. Mutational origins of wrinkly spreader diversity. Genetics 176, 441–453 (2007).

Article CAS PubMed PubMed Central Google Scholar

Jenal, U., Reinders, A. & Lori, C. Cyclic di-GMP: second messenger extraordinaire. Nat. Rev. Microbiol. 15, 271–284 (2017).

Article CAS PubMed Google Scholar

Hickman, J. W., Tifrea, D. F. & Harwood, C. S. A chemosensory system that regulates biofilm formation through modulation of cyclic diguanylate levels. Proc. Natl Acad. Sci. USA 102, 14422–14427 (2005).

Article CAS PubMed PubMed Central Google Scholar

Bourret, R. B. Receiver domain structure and function in response regulator proteins. Curr. Opin. Microbiol. 13, 142–149 (2010).

Article CAS PubMed PubMed Central Google Scholar

Laventie, B.-J. & Jenal, U. Surface sensing and adaptation in bacteria. Annu. Rev. Microbiol. 74, 735–760 (2020).

Article CAS PubMed Google Scholar

O’Neal, L. et al. The Wsp system of Pseudomonas aeruginosa links surface sensing and cell envelope stress. Proc. Natl Acad. Sci. USA 119, e2117633119 (2022).

Article PubMed PubMed Central Google Scholar

De, N., Navarro, M. V. A. S., Raghavan, R. V. & Sondermann, H. Determinants for the activation and autoinhibition of the diguanylate cyclase response regulator WspR. J. Mol. Biol. 393, 619–633 (2009).

Article CAS PubMed PubMed Central Google Scholar

Pianka, E. R. On r- and K-selection. Am. Nat. 104, 592–597 (1970).

Article Google Scholar

Andrews, J. H. Comparative Ecology of Microorganisms and Macroorganisms (Springer, 2017).

Bansept, F., Obeng, N., Schulenburg, H. & Traulsen, A. Modeling host-associating microbes under selection. ISME J. 15, 3648–3656 (2021).

Article PubMed PubMed Central Google Scholar

Valentini, M. & Filloux, A. Multiple roles of c-di-GMP signaling in bacterial pathogenesis. Annu. Rev. Microbiol. 73, 387–406 (2019).

Article CAS PubMed Google Scholar

Stiernagle, T. Maintenance of C. elegans. WormBook https://doi.org/10.1895/wormbook.1.101.1 (2006).

Papkou, A. et al. The genomic basis of Red Queen dynamics during rapid reciprocal host–pathogen coevolution. Proc. Natl Acad. Sci. USA 116, 923–928 (2019).

Article CAS PubMed Google Scholar

Rueden, C. T. et al. ImageJ2: ImageJ for the next generation of scientific image data. BMC Bioinformatics 18, 529 (2017).

Article PubMed PubMed Central Google Scholar

Mörck, C. & Pilon, M. C. elegans feeding defective mutants have shorter body lengths and increased autophagy. BMC Dev. Biol. 6, 39 (2006).

Article PubMed PubMed Central Google Scholar

O’Toole, G. A. Microtiter dish biofilm formation assay. J. Vis. Exp. https://doi.org/10.3791/2437 (2011).

Article PubMed PubMed Central Google Scholar

Serra, D. O., Richter, A. M. & Hengge, R. Cellulose as an architectural element in spatially structured Escherichia coli biofilms. J. Bacteriol. 195, 5540–5554 (2013).

Article CAS PubMed PubMed Central Google Scholar

Wiles, T. J. et al. Modernized tools for streamlined genetic manipulation and comparative study of wild and diverse proteobacterial lineages. mBio 9, e01877-18 (2018).

Article PubMed PubMed Central Google Scholar

Kissoyan, K. A. B. et al. Exploring effects of C. elegans protective natural microbiota on host physiology. Front. Cell. Infect. Microbiol. 12, 775728 (2022).

Article CAS PubMed PubMed Central Google Scholar

Schulenburg, V. D. et al. Extreme length and length variation in the first ribosomal internal transcribed spacer of ladybird beetles (Coleoptera: Coccinellidae). Mol. Biol. Evol. 18, 648–660 (2001).

Article PubMed Google Scholar

Andrews, S. FastQC: A Quality Control Tool for High Throughput Sequence Data (Babraham Institute, 2010).

Bolger, A. M., Lohse, M. & Usadel, B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics 30, 2114–2120 (2014).

Article CAS PubMed PubMed Central Google Scholar

Langmead, B. & Salzberg, S. L. Fast gapped-read alignment with Bowtie 2. Nat. Methods 9, 357–359 (2012).

Article CAS PubMed PubMed Central Google Scholar

Picard Toolkit (Broad Institute, 2019).

Li, H. A statistical framework for SNP calling, mutation discovery, association mapping and population genetical parameter estimation from sequencing data. Bioinformatics 27, 2987–2993 (2011).

Article CAS PubMed PubMed Central Google Scholar

Koboldt, D. C. et al. VarScan 2: somatic mutation and copy number alteration discovery in cancer by exome sequencing. Genome Res. 22, 568–576 (2012).

Article CAS PubMed PubMed Central Google Scholar

Cingolani, P. et al. A program for annotating and predicting the effects of single nucleotide polymorphisms, SnpEff. Fly 6, 80–92 (2012).

Article CAS PubMed PubMed Central Google Scholar

Cingolani, P. et al. Using Drosophila melanogaster as a model for genotoxic chemical mutational studies with a new program, SnpSift. Front. Genet. 3, 35 (2012).

RStudio Team. RStudio: Integrated Development for R (RStudio Inc., 2015).

R Core Team. R: A Language and Environment for Statistical Computing (R Foundation for Statistical Computing, 2016).

Winsor, G. L. et al. Enhanced annotations and features for comparing thousands of Pseudomonas genomes in the Pseudomonas genome database. Nucleic Acids Res. 44, D646–D653 (2016).

Article CAS PubMed Google Scholar

Madeira, F. et al. Search and sequence analysis tools services from EMBL-EBI in 2022. Nucleic Acids Res. 50, W276–W279 (2022).

Article CAS PubMed PubMed Central Google Scholar

Ren, J. et al. DOG 1.0: illustrator of protein domain structures. Cell Res. 19, 271–273 (2009).

Article CAS PubMed Google Scholar

Zamorano-Sánchez, D. et al. Functional specialization in Vibrio cholerae diguanylate cyclases: distinct modes of motility suppression and c-di-GMP production. mBio 10, e00670-19 (2019).

Article PubMed PubMed Central Google Scholar

Measuring Cell Fluorescence Using ImageJ (The Open Lab Book, 2014).

Bähre, H. & Kaever, V. in c-di-GMP Signaling: Methods and Protocols (ed. Sauer, K.) 45–58 (Springer, 2017).

Gao, X. et al. Functional characterization of core components of the Bacillus subtilis cyclic-di-GMP signaling pathway. J. Bacteriol. 195, 4782–4792 (2013).

Article CAS PubMed PubMed Central Google Scholar

Adams, K. J. et al. Skyline for small molecules: a unifying software package for quantitative metabolomics. J. Proteome Res. 19, 1447–1458 (2020).

Article CAS PubMed PubMed Central Google Scholar

Hmelo, L. R. et al. Precision-engineering the Pseudomonas aeruginosa genome with two-step allelic exchange. Nat. Protoc. 10, 1820–1841 (2015).

Article CAS PubMed PubMed Central Google Scholar

Figurski, D. H. & Helinski, D. R. Replication of an origin-containing derivative of plasmid RK2 dependent on a plasmid function provided in trans. Proc. Natl Acad. Sci. USA 76, 1648–1652 (1979).

Article CAS PubMed PubMed Central Google Scholar

Altschul, S. F., Gish, W., Miller, W., Myers, E. W. & Lipman, D. J. Basic local alignment search tool. J. Mol. Biol. 215, 403–410 (1990).

Article CAS PubMed Google Scholar

Ye, J. et al. Primer-BLAST: a tool to design target-specific primers for polymerase chain reaction. BMC Bioinformatics 13, 134 (2012).

Article CAS PubMed PubMed Central Google Scholar

Sievers, F. et al. Fast, scalable generation of high-quality protein multiple sequence alignments using Clustal Omega. Mol. Syst. Biol. 7, 539 (2011).

Article PubMed PubMed Central Google Scholar

Newman, J. R. & Fuqua, C. Broad-host-range expression vectors that carry the l-arabinose-inducible Escherichia coli araBAD promoter and the araC regulator. Gene 227, 197–203 (1999).

Article CAS PubMed Google Scholar

Choi, K.-H., Kumar, A. & Schweizer, H. P. A 10-min method for preparation of highly electrocompetent Pseudomonas aeruginosa cells: application for DNA fragment transfer between chromosomes and plasmid transformation. J. Microbiol. Methods 64, 391–397 (2006).

Article CAS PubMed Google Scholar

Tamura, K., Stecher, G. & Kumar, S. MEGA11: Molecular Evolutionary Genetics Analysis version 11. Mol. Biol. Evol. 38, 3022–3027 (2021).

Article CAS PubMed PubMed Central Google Scholar

Hesse, C. et al. Genome-based evolutionary history of Pseudomonas spp. Environ. Microbiol. 20, 2142–2159 (2018).

Article CAS PubMed Google Scholar

Benjamini, Y. & Hochberg, Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J. R. Stat. Soc. B 57, 289–300 (1995).

Google Scholar

Vu, V. et al. vqv/ggbiplot: A biplot based on gplot2. Github www.github.com/vqv/ggbiplot (2021).

Josse, J., & Husson, F. missMDA: a package for handling missing values in multivariate data analysis. J. Stat. Softw. https://doi.org/10.18637/jss.v070.i01 (2016).

Oksanen, J. et al. vegan: Community Ecology Package. Github https://github.com/vegandevs/vegan (2022).

Martinez Arbizu, P. pairwiseAdonis: Pairwise multilevel comparison using adonis. R package version 0.4 Github https://github.com/pmartinezarbizu/pairwiseAdonis (2020).

Bates, D., Mächler, M., Bolker, B. & Walker, S. Fitting linear mixed-effects models using lme4. J. Stat. Softw. 67, 1–48 (2015).

Article Google Scholar

Zeileis, A. & Hothorn, T. Diagnostic checking in regression relationships. R News 2, 7–10 (2002).

Hothorn, T., Bretz, F. & Westfall, P. Simultaneous inference in general parametric models. Biom. J. 50, 346–363 (2008).

Article PubMed Google Scholar

Rigby, R. A. & Stasinopoulos, D. M. Generalized additive models for location, scale and shape. J. R. Stat. Soc. C 54, 507–554 (2005).

Article Google Scholar

Wikham, H. ggplot: Elegant Graphics for Data Analysis (Springer, 2016).

Kassambara, A. ggpubr: ‘ggplot2’ based publication ready plots. R package version 0.6.0. https://rpkgs.datanovia.com/ggpubr/ (2023).

Caswell, H. Matrix Population Models (Sinauer, 2001).

Download references

We thank E. Stukenbrock (University of Kiel, Germany) for access to the LSM 880; K. Guillemin (University of Oregon, Eugene, United States), P. Rainey (Max-Planck Institute for Evolutionary Biology, Ploen, Germany), H. Schweizer (Northern Arizona University, United States), F. Yildiz (University of California Santa Cruz, United States) and G. O’Toole (Dartmouth Medical School, United States) for providing bacterial strains or plasmids; D. Rogers, J. Summers (both Max-Planck Institute for Evolutionary Biology, Ploen, Germany) for guidance in allelic exchange; J. Zimmermann (Schulenburg group, University of Kiel, Germany) for bioinformatic support; B. Pees (Schulenburg group, University of Kiel, Germany) for illustration support; S. Joel, J. Hofmann, J. Löwenstrom, J. Lorenzen, H. Griem-Krey, L. Bluhm and L. Rheindorf (all Schulenburg group, University of Kiel, Germany) for lab support; the Kiel BiMo/LMB for access to their core facilities; the Schulenburg lab for project feedback; and B. Bohannan (University of Oregon, Eugene, United States), R. Knight (University of California San Diego, United States) and P. Engel (Université de Lausanne, Switzerland) for advice on the manuscript. Funding was provided by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation), Project-ID 261376515 – SFB 1182, Projects A4 and Z3 (N.O., A.C., F.B., J.L., A. Tholey, A. Traulsen, H. Schulenburg); the DFG Research Infrastructure NGS_CC project 407495230 (J.F.) as part of the Next Generation Sequencing Competence Network project 423957469; the International Max-Planck Research School for Evolutionary Biology (N.O., A.C.); the Max-Planck Society (Fellowship to H. Schulenburg); and NIH project R01AI168017 (M.J.G.G., H. Sondermann).

Open access funding provided by Christian-Albrechts-Universität zu Kiel.

Thekla Schultheiß

Present address: Institute of Toxicology and Pharmacology, University of Kiel, Kiel, Germany

Department of Evolutionary Ecology and Genetics, University of Kiel, Kiel, Germany

Nancy Obeng, Anna Czerwinski, Daniel Schütz, Jan Michels, Thekla Schultheiß, Melinda Kemlein & Hinrich Schulenburg

Department of Systematic Proteome Research and Bioanalytics, University of Kiel, Kiel, Germany

Jan Leipert & Andreas Tholey

Max Planck Institute for Evolutionary Biology, Plön, Germany

Florence Bansept, Arne Traulsen & Hinrich Schulenburg

CSSB Centre for Structural Systems Biology, Deutsches Elektronen-Synchrotron DESY, Hamburg, Germany

María J. García García & Holger Sondermann

Institute of Clinical Molecular Biology, University of Kiel, Kiel, Germany

Janina Fuß

Section of Biology, University of Kiel, Kiel, Germany

Holger Sondermann

You can also search for this author in PubMed Google Scholar

N.O., A.C., F.B., A. Traulsen and H. Schulenburg conceptualized the project. N.O., A.C., M.J.G.G. and H. Sondermann developed the methodology. N.O., A.C., J.M., J.L., T.S., M.K. and J.F. conducted investigations. N.O., A.C., D.S. and F.B. analysed data. N.O., A.C., D.S., J.M., J.L., F.B., M.J.G.G., T.S., M.K., J.F., A. Tholey, A. Traulsen, H. Sondermann and H. Schulenburg contributed to the writing of the manuscript. N.O., A. Tholey, A. Traulsen, H. Sondermann and H. Schulenburg supervised the project. F.B. performed modelling.

Correspondence to Hinrich Schulenburg.

The authors declare no competing interests.

Nature Microbiology thanks Hassan Salem and the other, anonymous, reviewer(s) for their contribution to the peer review of this work.

Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

a, b, Bacterial fitness during the evolution experiment. a, Bacterial fitness in host across cycles of the evolution experiment measured as colony forming units (CFU) per worm population after 3.5 days of exposure to Ce_MY316. b, In the negative control, bacterial fitness was assessed on nematode growth agar in absence of the host. For each data point, bacteria were collected at the bottleneck time point of the noted cycle. Replicate populations (n = 6) are shown as separate thin lines, with the mean shown as a thick line. c, Mean CFU per individual host in a worm population for the evolved bacterial populations of cycle 10. Five L4 C. elegans larvae proliferated on evolved or ancestral bacterial lawns for 3.5 days (reaching F1 generation) and CFUs were extracted from the whole worm population. CFUs per population were divided by the number of worms in the population. Overall, results are shown as boxplots, with boxes indicating 25% above and below the median, which is given as the thick line within boxes; replicate populations (n=6) are indicated as individual data points. d, e, Dynamic changes in morphotype composition during the free-living phase of the host-associated life cycle for bacterial populations from the end of the evolution experiment. c, Results for the replicate populations from the host-associated evolution treatment. Box plots show median (center line), upper and lower quartiles (box limits) and the interquartile range (whiskers). d, Results for the replicate populations from the control treatment. Proportions of the different colony morphotypes (see graphical legend) is shown across time of the host-associated life cycle. Time point 0 is at the end of the host-associated phase, when bacteria are transferred to the free-living phase, which itself lasts 168 hours. Fdr-corrected beta-regressions were used to predict proportions and test for a change in proportions over time (see Supplementary Table 2).

Confocal laser scanning micrographs (upper panel: longitudinal optical section; lower panels: maximum intensity projections showing longitudinal optical sections) revealing intact bacterial cells (red) within the intestinal system of a young adult Ce_MY316 (cyan). The upper micrograph shows an overview of the complete worm, and the lower micrographs show detailed views of the worm sections indicated by the dashed frames above. These include the posterior pharynx with the worm grinder and the first intestinal ring (left), a central intestinal (middle) and the anal region (right). The bottom left panel is identical to the micrograph shown in the main text (Fig. 1e). Scale bars = 50 µm (overview) and 10 µm (detailed views).

Phenotypes of morphotype clones isolated from independent host-evolved populations (left) and control populations (right), including smooth, fuzzy and wrinkly morphotypes, are shown. Results for each morphology are summarized as boxplots (1 < n < 4). Box plots show median (center line), upper and lower quartiles (box limits) and the interquartile range (whiskers). Dashed lines and grey shaded areas indicate the mean and standard error of ancestral traits, respectively. Differences between evolved morphologies and the ancestor were assessed with generalized linear models and fdr-corrected Tukey post-hoc tests. Letters indicate statistical differences between morphologies, asterisks indicate deviation from the ancestor (see Supplementary Table 6).

a, Principal component analysis on characteristic stages of host-association for ancestral, host evolved and control evolved bacterial populations. Individual data points refer to replicate populations colored according to evolution treatment (Supplementary Table 5). b–f, Shifts in phenotypes from the bacterial ancestor in the evolved populations for (b) early colonization, determined by CFU extracted from L4 larvae exposed to bacteria for 1.5 hour; (c) persistence in L4 larvae kept in M9 buffer for 1h (raised on bacteria); (d) CFU of bacteria released from L4 larvae into buffer within 1h (previously raised on bacteria from L1 to L4), (e) swarming distance on 0.5% agar within 24; and (f) colony expansion on 3.4% agar within 72h. All panels show ratios of evolved over ancestral populations for five replicates shown as individual data points. The dashed line indicates the mean values obtained for the ancestral population. The difference between evolved and ancestral phenotypes were assessed using one-sided t-tests (fdr-corrected; Supplementary Table 8). In all box plots, median (center line), upper and lower quartiles (box limits) and the interquartile range (whiskers) are shown.

a, Amount of intracellular c-di-GMP measured with a fluorescence sensor. Raw fluorescence intensity (RFI) is the ratio of TurboRFP (c-di-GMP-dependent) and AmCyan (plasmid copy number-dependent) and, thus normalized for copy number of the sensor plasmid. b, Total c-di-GMP determined by isotope dilution PRM analysis. c, C-di-GMP determined by isotope dilution PRM analysis and normalized by total protein amount. In a, b and c, c-di-GMPs levels are studied for the ancestor and three wrinkly isolates from the end of the evolution experiment, each with a single mutation in either wspF, wspE, or rph. We compared c-di-GMP levels replicate cell populations (n = 5). d, Intracellular c-di-GMP of isolates from end the of the evolution experiment with wrinkly, smooth, and fuzzy colony morphology measured with a fluorescence sensor (nested ANOVA and fdr-corrected Dunnett post hoc test; n = 5). e, RFI of the different isolates normalized by ancestral RFI to correct for replication-dependent effects (ANOVA and fdr-corrected Dunnett post hoc test; n = 5; Supplementary Table 10. For respective mutations see Supplementary Table 9). Box plots show median (center line), upper and lower quartiles (box limits) and the interquartile range (whiskers).

On top, illustrations of protein domains show the organization of WspF, WspE and ribonuclease PH with sites affected by Pl_MYb11 mutation highlighted by red arrows. Below, protein sequence alignments of ancestral and evolved proteins are shown with domains highlighted in the same colors as above.

a, Biofilm formation of wspF, wspE, and rph mutants compared to Pl_MYb11 (ancestral median = dashed line) after two days shaking incubation microtiter plates. Illustration of biofilms with photographs of biofilms from test tubes for after 48h incubation. Scale bars = 4mm. b, Swarming motility and c, colony expansion of isolates wspF, wspE, and rph mutants were compared to Pl_MYb11 (ancestral median = dashed line) after 24 and 72h, respectively. b-c, Scale bars = 0.5 cm. a-c, Experiments were performed with min. 3 replicates/treatment and analyzed with ANOVA and fdr-corrected Dunnett post-hoc tests or a one-sided t-test. Box plots show median (center line), upper and lower quartiles (box limits) and the interquartile range (whiskers). d, Macrocolonies of ancestral Pl_MYb11 and three wrinkly isolates from the end of the evolution experiment (with single mutations in either wspF, wspE, or rph). e, Macrocolonies of ΔwspR mutants in Pl_MYb11 wspF, wspE, or rph background. f, Macrocolonies of wspF, wspE, or rph isolates expressing the phosphodiesterase PA2133 from plasmid, Pl_ MYb11 expressing the constitutively active diguanylate cyclase GCN4-WspR from plasmid, and empty vector controls. d-f, Macrocolonies were grown on tryptic soy agar supplemented with Congo Red for one (d,e) to three (f) days. Scale bars = 1mm.

Three co-evolved morphotypes isolated from host-evolved replicate population T3 were paired with the ancestor. In one quartet the wrinkly wspF mutant (MT12) was present and in the other the wrinkly wspE mutant (MT14). Data points represent independent replicates, min. 3 replicates per treatment. Differences between morphotypes and ancestor were assessed using a linear mixed model and subsequent fdr-corrected Dunnett post-hoc comparisons; see Supplementary Table 13. Box plots show median (center line), upper and lower quartiles (box limits) and the interquartile range (whiskers).

Distribution of sequence lengths and percent identities of the BLAST results for individual genes. The proportion of BLAST results belonging to particular sequence length and percent identity classes are shown as blue shades with varying intensity (cf. scale). Red rectangles show the areas for which the presence of the considered gene is assumed, and were set to include the largest BLAST hit values at both maximum sequence length and percent identities.

a, Definition of the rates for the model of a microbial lineage being taken up by, replicating within and being expulsed from worms on a plate. b, Distribution of the optimal strategies across the whole traits space. c, Projection of the trait space over the axes (r, δ). For each point of the map, the color represents the proportions of times each of the 4 possible strategies are optimal, integrating over the values of f and K. The color scheme uses the CMYK color code: a purely cyan (respectively magenta, yellow) pixel indicates that the only optimal strategy for the considered values (r, δ) is ↓ δ (respectively ↑ f, ↑K). A color of a darker shade indicates that ↑ r is also optimal in a small proportion at that point, as shown on the additional color scales for each edge.

Supplementary Tables 1–22.

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

Reprints and Permissions

Obeng, N., Czerwinski, A., Schütz, D. et al. Bacterial c-di-GMP has a key role in establishing host–microbe symbiosis. Nat Microbiol (2023). https://doi.org/10.1038/s41564-023-01468-x

Download citation

Received: 21 September 2022

Accepted: 10 August 2023

Published: 31 August 2023

DOI: https://doi.org/10.1038/s41564-023-01468-x

Anyone you share the following link with will be able to read this content:

Sorry, a shareable link is not currently available for this article.

Provided by the Springer Nature SharedIt content-sharing initiative

Blog