Identification of Cyclic-di-GMP-Modulating Protein Residues by Bidirectionally Evolving a Social Behavior in Pseudomonas fluorescens

Microbial communities comprise densely packed cells where competition for space and resources is fierce. Aging colonies of Pseudomonas fluorescens are known to repeatedly produce mutants with two distinct phenotypes that physically work together to spread away from the overcrowded population. ABSTRACT Modulation of the intracellular cyclic di-GMP (c-di-GMP) pool is central to the formation of structured bacterial communities. Genome annotations predict the presence of dozens of conserved c-di-GMP catalytic enzymes in many bacterial species, but the functionality and regulatory control of the vast majority remain underexplored. Here, we begin to fill this gap by utilizing an experimental evolution system in Pseudomonas fluorescens Pf0-1, which repeatedly produces a unique social behavior through bidirectional transitions between two distinct phenotypes converging on c-di-GMP modulation. Parallel evolution of 33 lineages captured 147 unique mutations among 191 evolved isolates in genes that are empirically demonstrated, bioinformatically predicted, or previously unknown to impact the intracellular pool of c-di-GMP. Quantitative chemistry confirmed that each mutation causing the phenotypic shift either amplifies or reduces c-di-GMP production. We identify missense or in-frame deletion mutations in numerous diguanylate cyclase genes that largely fall outside the conserved catalytic domain. We also describe a novel relationship between a regulatory component of branched-chain amino acid biosynthesis and c-di-GMP production, and predict functions of several other unexpected proteins that clearly impact c-di-GMP production. Sequential mutations that continuously disrupt or recover c-di-GMP production across discrete functional elements suggest a complex and underappreciated interconnectivity within the c-di-GMP regulome of P. fluorescens. IMPORTANCE Microbial communities comprise densely packed cells where competition for space and resources is fierce. Aging colonies of Pseudomonas fluorescens are known to repeatedly produce mutants with two distinct phenotypes that physically work together to spread away from the overcrowded population. We demonstrate that the mutants with one phenotype produce high levels of cyclic di-GMP (c-di-GMP) and those with the second phenotype produce low levels. C-di-GMP is an intracellular signaling molecule which regulates many bacterial traits that cause tremendous clinical and environmental problems. Here, we analyze 147 experimentally selected mutations, which manifest either of the two phenotypes, to identify key residues in diverse proteins that force or shut down c-di-GMP production. Our data indicate that the intracellular pool of c-di-GMP is modulated through the catalytic activities of many independent c-di-GMP enzymes, which appear to be in tune with several proteins with no known links to c-di-GMP modulation.

is synthesized from two molecules of GTP by a diguanylate cyclase (DGC) and hydrolyzed by a phosphodiesterase (PDE) (1). DGCs and PDEs are identified by the GGDEF and EAL/HD-GYP domains, respectively, and there are also hybrid proteins with both DGC and PDE functions (7,8). However, relatively little is known about other functional domains, if any at all, in most DGCs/PDEs.
Pseudomonas spp. typically possess greater than 50 predicted DGCs or PDEs (9), but systematically knocking out independent genes appears to have relatively little impact on the intracellular c-di-GMP pool (9). Each enzyme could collectively influence the total intracellular pool (10) or independently carry out a subcellularly localized role (11). Specific signaling pathways or regulators have been characterized to independently regulate several DGCs and PDEs (12)(13)(14)(15)(16)(17)(18)(19), and some have been demonstrated to be under allosteric control, where a certain range of c-di-GMP levels either initiate or prevent enzymatic activity (20)(21)(22)(23)(24). However, the function of the vast majority of bioinformatically predicted DGCs and PDEs has not been empirically validated, and their regulatory mechanisms remain entirely unknown (25).
Experimental evolution studies in Pseudomonas and Burkholderia readily capture diverse mutations that impact biofilm formation, including those that alter c-di-GMP levels (26)(27)(28). We have previously described the bidirectional evolution of mucoid (M) and dry (D) colony phenotypes in P. fluorescens Pf0-1, which appears to occur exclusively through c-di-GMP oscillation (29). An aging colony of M cells repeatedly produces spreading fans at the perimeter, which comprises the original M cells along with new D mutant cells that produce a dry and wrinkly colony morphology on their own. Additionally, an aging colony of D cells also repeatedly generates spreading fans, which comprise the original D cells along with a new mutant that appears phenotypically identical to the ancestral M. This bidirectional selection of M and D phenotypes occurs continuously because M and D cells self-organize in space and physically work together to spread out from the nutrient-poor and crowded conditions of dense colony growth, and neither type is capable of spreading out on its own. Genetic analyses of five evolved isolates revealed characteristic mutations in wsp genes to suggest that c-di-GMP production is likely elevated in D and reduced in M (29).
Here, we carry out a large-scale parallel evolution experiment using the bidirectional M-D system to identify 147 unique mutations in genes that are empirically demonstrated, bioinformatically predicted, or previously unknown to impact c-di-GMP levels. We expand the current model of c-di-GMP production by the Wsp signal transduction system and identify unique missense mutations outside the conserved catalytic domain of numerous DGCs that stimulate or reduce enzymatic activity. We also describe proteins with no previous associations to c-di-GMP production that appear to influence the activity of several DGCs, including a regulatory component of the branched-chain amino acid (BCAA) biosynthesis pathway. Lastly, we model the interconnectivity of the c-di-GMP regulome by mapping sequential mutations that disrupt or recover c-di-GMP production across functionally discrete systems.

RESULTS AND DISCUSSION
Bidirectional evolution of mucoid and dry phenotypes occurs through intracellular c-di-GMP pool oscillation. The bidirectional selection of mucoid (M) and dry (D) phenotypes occurs continuously in our experimental system (29) because neither type is capable of spreading out on its own, but a colony of initially mixed M and D cells readily spreads out in a radial pattern without requiring additional mutations (Fig. 1A). With the expectation that a large-scale utilization of this experimental evolution system could identify diverse mechanisms of c-di-GMP modulation, we initially set up 25 parallelly evolving lineages of the common ancestral M isolate, which is also referred to as M0. We randomly picked evolved isolates of D from M0, and each D isolate (D0) represents the first branching point and the ancestor of its respective lineage (Fig. S1). We then randomly picked several evolved isolates of M from independent D0 isolates, and they are referred to as M1; D1 evolved from M1, M2 evolved from D1, D2 evolved from M2, and so on (see Materials and Methods for a detailed description of the experimental procedures). We also established additional lineages with genetically engineered ancestral strains, as described below, to represent 33 lineages in total. As anticipated, we observed bidirectional transitions between M and D phenotypes across all lineages. In sum, we evolved over 600 isolates, analyzed genome sequences of 191 isolates, and quantified c-di-GMP in a subset of 56 isolates. Only the sequenced isolates within in each lineage are depicted in Fig. S1 in the supplemental material and summarized in Table S1. We also employed a numerical nomenclature to distinguish isolates from independent lineages. The first number represents the lineage number, and the second number represents the chronology of sequential transition steps. For example, 5.1 represents the first D0 isolate that evolved from M0 in lineage 5, and 5.2. represents the M1 isolate that evolved from D0, and so on.
Remarkably, we primarily detected a single mutation associated with each phenotypic transition, with few additional silent mutations across all genome sequences. Genome sequencing of successive isolates within each lineage provided an internal control, which confirmed the exclusive presence of the respective preceding mutations. Among the 147 unique mutations detected, 81 mutations were in wsp genes and 66 mutations were in other genes including those that are hypothetical. Motility assays showed that all D isolates were nonmotile and all M isolates were motile as expected, and liquid chromatography tandem mass spectrometry (LC-MS/MS) quantification confirmed that c-di-GMP was elevated in all D isolates relative to their respective parent M isolate (Table S1). Notably, many mutations in non-wsp genes sequentially flanked specific mutations in the Wsp signal transduction system, the most extensively studied c-di-GMP regulatory system in Pseudomonas spp. (19). Several missense mutations observed in Wsp proteins identified in this study overlap comparable mutations from various Pseudomonas and Burkholderia studies that we had previously compiled from the literature (19) (Fig. 2). However, the vast majority of mutations identified in this study had never been observed, and many occur outside functionally annotated domains. The data presented in Fig. 2 represent the most comprehensive record of functionally important residues of the Wsp system to date and empower our predictions of mutation consequences. We first focus on sequential mutations in the Wsp system to reveal both predictable and surprising evolutionary innovations to modulate c-di-GMP production. For the sake of an  Table S1. Isolate 1.0 was constructed from the ancestral M isolate to lack the wspC gene, which encodes a methyltransferase, known to activate the Wsp system. Each isolate was sampled from naturally emerging fans from their respective parent colony with the opposite phenotype. This phylogeny depicts the bidirectional transitions between mucoid and dry states, which we show are c-di-GMP dependent. (C) C-di-GMP levels of evolved isolates. C-di-GMP was quantified through LC-MS/MS and is reported as nM per OD 600 . Samples are categorized by phenotype, and we report the triplicate mean 6 standard deviation (SD) for each phenotypic classification: mucoid n = 14 and dry n = 21 (***, Student's t test P , 0.001).
Cyclic-di-GMP Modulation in Pseudomonas fluorescens mSystems accessible narrative, we will primarily refer to five simplified lineages (Fig. 1B), with corresponding c-di-GMP measurements (Fig. 1C) and colony morphologies (Fig. S2). Diverse mutations dynamically rewire signal flow to produce a spectrum of functional states of the Wsp signal transduction system. The Wsp signal transduction system was first reported 2 decades ago in P. fluorescens (30) and continues to serve as the primary model system of c-di-GMP regulation in pseudomonads. The Wsp system comprises seven proteins encoded in a monocistronic operon under the control of an unidentified transcriptional regulator (31), and the current functional model is based on a combination of empirical studies and predictions based on homologies to the chemotaxis system in Escherichia coli (19). The main difference with the latter is that the signal cascade is initiated by the extracellular sensory domain of WspA detecting surface contact (32,33) to ultimately drive c-di-GMP production by the DGC WspR (Fig. 3). Given that a mutation in any one of the seven Wsp protein-coding genes could potentially increase or decrease c-di-GMP production through WspR, we expected wsp genes to be a common target in our parallel evolution experiments. Indeed, the vast majority of the first transitions from the ancestral M to D occurred through either a missense or an in-frame deletion mutation in wsp genes (Table S1).
Across our subset of five lineages (Fig. 4), initial mutations appear to artificially activate WspA, followed by surprisingly similar patterns of sequential mutations in various wsp genes to oscillate c-di-GMP production. Lineages 2 and 3 began with the first respective D isolate possessing a unique in-frame deletion that partially spans the predicted methylation helices of WspA. Lineage 5 began with a missense mutation in the phosphoreceiver domain of Only the missense and in-frame mutations identified in this study are shown, and they are mapped to the consensus sequence of each Wsp protein (19). Previously annotated mutations (19) are shown in the gray box below each protein. Mutations indicated in the gray box are annotated to represent specific residue positions that were affected in the respective protein sequence of a given organism, but the gray horizontal bars indicate their actual positions within the consensus sequence. Therefore, the horizontal bar locations of the mutations identified in this study are directly comparable to previously identified mutations through the consensus sequence. Mutations that stimulate c-di-GMP production are indicated by orange horizontal bars, and those that reduce c-di-GMP production are indicated by blue horizontal bars. WspF (P105S), which likely disrupts its methyl-esterase activity to keep WspA in an activated state. This activation is undone in the following M isolate by an in-frame deletion in wspA that removes the amino acid residue S329, and then c-di-GMP production is restored again in the next D isolate through an in-frame deletion of residues 283 to 296 in WspA. This deletion falls within the same deleted region as those observed in lineages 2 and 3 ( Fig. 2, Fig. 4).
Although methylation of WspA remains to be empirically demonstrated, its predicted methylation helices clearly have a functional role in Wsp signaling, as they were commonly mutated to activate c-di-GMP production (Fig. 2). Based on our mutation data and homologies to Tsrmethylation in the chemotaxis system of E. coli (34), we predict that the residues Q285, E292, E299, Q502, and E512 are methylated and the deletion of these sites mimics a methylated state. We also set up an independent lineage to commence with an engineered DwspC M isolate ( Fig. 4, isolate 1.0), under the expectation that the deletion of wspC would steer downstream mutations away from the Wsp system. However, we detected a missense mutation within the signaling domain of WspA (A381V) in the first D progeny isolate 1.1, and the same A381V mutation was previously observed in an independently evolved D isolate with an unaltered wspC gene (29). We detected additional missense mutations in WspA's signaling domain to independently stimulate c-di-GMP production ( Fig. 2), which likely bypass the methylation requirement to activate WspE. Mutations that stimulate WspA were typically followed by a missense mutation in WspB, WspD, or WspE (Table S1). A missense mutation in WspD (W90L) in lineage 2 reduced c-di-GMP production like the WspB (L136P) mutation in lineage 1 (Fig. 4). The functional roles of WspB and WspD in signal relay are loosely modeled after homology to CheW of E. coli, which remains the least characterized protein in the chemotaxis system. Mutations in CheW that successfully disrupt signal transduction are largely limited to five unique amino acid clusters that situate on the protein surface and are essential for the docking of CheW to Tsr (WspA homolog) (35). WspB and WspD also possess six and seven highly conserved amino acid clusters, respectively (19). Our WspB L135P and WspD W90L mutations fall within these conserved clusters that likely situate on the surface of WspB and WspD and natively act as docking sites for WspA and WspE to permit signal relay. We also observed additional missense and in-frame deletion mutations in WspB and WspD that independently impact c-di-GMP production (Fig. 2, Table S1), confirming that both WspB and WspD are necessary for Wsp signal transduction through domains and interactions that remain to be characterized. Homology with the chemotaxis system predicts that once WspA is stimulated by surface contact, it undergoes a conformational change to expose cytosolic methylation helices that are necessary for signal transduction (34). The methyl-transferase WspC is then predicted to dock at the 59 tail of WspA and methylate the WspA trimer-of-dimers complex (58). The signal is then relayed through the coupling proteins WspB and WspD, and subsequently to WspE, which functions as a histidine kinase (59). WspE phosphorylates the DGC WspR to stimulate c-di-GMP production (60). Simultaneously, WspE phosphorylates the methyl-esterase WspF, which demethylates WspA and terminates the signaling cascade (31). Significant work has been carried out to characterize the biochemical functions of WspE and WspR in particular, and WspE is a common mutational target in clinical settings (19). More recently, studies have begun to explore the functional role of WspA during signal activation and subcellular localization, revealing that both properties are dependent on WspB and WspD proteins (32,61). These studies also suggest that WspB and WspD likely possess unique functions, but they remain unresolved. Loss-of-function mutations (i.e., nonsense and frameshift mutations) in wspE were exclusively associated with D to M transitions throughout all lineages to reduce c-di-GMP production (Table S1). These observations are not particularly surprising, since WspR cannot be phosphorylated in the absence of a functional WspE (Fig. 3). In contrast, multiple missense and in-frame deletions were observed in the phosphoreceiver domain of WspE to exclusively stimulate c-di-GMP production ( Fig. 2). Accordingly, mutations within the conserved phosphoreceiver domain of CheY, a WspE homolog in E. coli, can activate the protein without phosphorylation (36,37). In addition, the DGC activity of WspR in P. aeruginosa could be modulated in a phosphorylation-independent manner (18,38). Artificial activation of WspE would deregulate WspA-dependent activation and was often the evolutionary solution to recover c-di-GMP production when the upstream signaling components of the Wsp system became functionally crippled (Table S1). The M isolate 1.4 of lineage 1 carries a frameshift mutation in wspE, which should terminally shut down the Wsp system, and the c-di-GMP level was reduced as predicted (Fig. 4). However, c-di-GMP production was restored in this lineage once again in the following D isolate 1.5 through a missense mutation in DgcH (V459G), a bioinformatically predicted DGC that impacts biofilm formation in P. aeruginosa PAO1 (39). DgcH is capable not only of producing c-di-GMP, but of producing it at one of the highest levels detected in our study (Fig. 4, Table S1). In contrast to the Wsp system, it is difficult to predict functional consequences of mutations in other DGCs, since so little is known beyond their conserved catalytic domain.
Many DGCs are functionally capable of causing the M to D shift independently from WspR. Sequential patterns of mutations and c-di-GMP quantification across all 33 lineages indicate that the DGC DgcH is also a dominant contributor to the intracellular c-di-GMP pool under our experimental conditions. Other known or predicted DGCs were typically targeted to modulate c-di-GMP levels after the Wsp system and DgcH became genetically dysfunctional. One exception is MorA, which represents the first mutation target in two lineages (Table S1). The first D isolate 4.1 of lineage 4 carries an in-frame deletion mutation in dgcH that removes amino acid residues 524 to 537 ( Fig. 4 and 5). Interestingly, the  Table S1, and mutations reflect the affected amino acid residues, except that frameshift (FS) and promoter region (PR) mutations refer to the specific nucleotides.
GGDEF diguanylate cyclase domain of DgcH spans residues 537 to 693, and the removal of the immediately adjacent residues appears to activate c-di-GMP production. We observed a reduction in c-di-GMP in the following M isolate 4.2, which has a missense mutation in DgcH (L634P). A search against NCBI's conserved domain database identifies that residues 619 to 643 contain the inhibitory I-site (24,40), with residue N637 as the predicted interaction site with c-di-GMP. Our sequence analysis supports this annotation and the likelihood that DgcH is under allosteric control by c-di-GMP; the L634P mutation in DgcH could mimic a c-di-GMP bound state to cause the observed reduction in c-di-GMP. We cannot discount the possibility of additional regulatory mechanisms, since numerous missense and in-frame deletion mutations were detected across the coding sequence of DgcH (Fig. 5). As already discussed above, the V459G mutation in DgcH greatly increased c-di-GMP production in lineage 1 (Fig. 4), and in lineage 27, the Q454P mutation caused a shift to the D phenotype and then reverted back to the M phenotype through another missense mutation (Q148P) ( Table S1). Similar to DgcH, specific missense mutations in other DGCs caused bidirectional shifts, and most activating mutations were found outside their conserved catalytic domain (Fig. 5). We also detected the I-site in other DGCs, with MorA being the lone exception, but no associated mutations were observed (Fig. 5). To the best of our knowledge, no functional studies have been carried out to characterize the bioinformatically predicted domains of these DGCs beyond their catalytic domain. We have clearly established that mutations Cyclic-di-GMP Modulation in Pseudomonas fluorescens mSystems observed in this study impact c-di-GMP production, which should stimulate future studies to characterize their catalytic and/or regulatory consequences.
Mutations in genes with previously unknown associations with c-di-GMP production integrate with Wsp or other DGC mutations to modulate c-di-GMP levels. We also observed mutations in bioinformatically annotated genes with no known associations to c-di-GMP modulation but with clear associations with Wsp components and other DGCs. This pattern is best represented in lineage 5 (Fig. 4). Like most lineages, mutations initially accumulate within wsp genes and deregulate the Wsp system. Subsequent isolates D 5.5 and M 5.6 carry mutations that activate and then deactivate the DGC SadC, respectively. Surprisingly, isolate 5.7 returns this lineage to the D phenotype through a missense mutation in IlvH (G150S), which is predicted to function in the biosynthesis of branched-chain amino acids (BCAA). We observed three independent missense mutations in IlvH across our study that exclusively cause a phenotypic shift from M to D and reduced motility (Table S1). Accordingly, quantification of c-di-GMP shows that each of these three mutations significantly increases c-di-GMP respective to their immediate M parent (Fig. 6). These results indicate that a previously unknown relationship exists between BCAA biosynthesis and c-di-GMP production, and we explore this relationship in detail below.
Lineage 5 returns to the M phenotype in isolate 5.8 through a missense mutation in the hypothetical protein Pfl01_0730 (D115N) (Fig. 4). This protein, here ParA (PAP2 associated c-di-GMP regulator), shares homology with the PAP2 superfamily of proteins, suggesting that it may be involved in phospholipid turnover and impact the bioavailability of phosphate within the cell (41). Although we found the relative decrease in c-di-GMP in isolate 5.8 to be not statistically significant from its parent likely due to technical noise (Fig. 6), the observed mutation clearly drives the phenotypic transition from D to M. Furthermore, if this mutation does not impact c-di-GMP levels, then we would not expect a further transition through elevated c-di-GMP production. Indeed, we captured two independent transitions back to the D phenotype, both with substantially increased c-di-GMP levels (Fig. 4). Isolates 5.9a and 5.9b are sister isolates (i.e., they share a parent) that have developed unique chimeric Wsp proteins. Both chimeric proteins (WspA::D and WspB::D) contain identical partial segments of WspD commencing at residue 158, situated immediately adjacent to the in-frame deletion detected in their common ancestral isolate 5.4, which previously deactivated the Wsp system (Fig. 4). To determine whether these chimeric proteins alone are sufficient to drive the phenotypic shift from M to D, we replaced the corresponding native Wsp proteins in the ancestral M isolate with each chimeric protein. In parallel, we also replaced the corresponding naturally modified Wsp proteins in the immediate parental M isolate 5.8 with each chimeric protein. Both engineered isolates in the ancestral M background retained the M phenotype, but both engineered isolates in the M isolate 5.8 background recapitulated the D transition Cyclic-di-GMP Modulation in Pseudomonas fluorescens mSystems that was observed in isolates 5.9a and 5.9b (Table 1). These results indicate that the two chimeric Wsp proteins require additional mutation(s) observed within this lineage to produce c-di-GMP. Both the WspA::D and WspB::D chimeras have lost WspC, which is not surprising since the methylation sites of WspA were removed previously in the D isolate 5.3. Similarly, WspB is completely removed in the WspA::D chimera, and nearly 75% of WspB has been removed in the WspB::D chimera. At this stage, it remains unclear whether the observed mutations in ilvH and/or parA are required for the function of the two Wsp chimeras. However, we present additional evidence below that strengthens the relationship between IlvH and the Wsp system. The WspB::D chimera isolate 5.9b was next used to evolve M progeny isolates. Sister isolates 5.10a and 5.10b were each found to carry frameshift deletions in ilvI that decreased c-di-GMP levels (Fig. 4). The two ilvI mutants carry a deletion of 17 or 23 nucleotides that start at the same position, which should abolish IlvI's function by causing a frameshift. IlvI and IlvH form a complex that catalyzes the first step of BCAA biosynthesis, where IlvH serves as the regulatory subunit and IlvI serves as the catalytic subunit (42). IlvH in E. coli has been bioinformatically predicted to respond to c-di-GMP levels due to the presence of the inhibitory I-site (40), which we also detected in most DGCs (Fig. 5). We identified that IlvH of P. fluorescens Pf0-1 also contains the I-site RXXD motif at residues 103 to 106. None of the independently identified IlvH mutations (A36E, S149R, and G150S) in our study fall within this motif, but all three mutations significantly increase c-di-GMP levels (Fig. 6). All identified mutations in IlvH fell outside the core RXXD motif. However, it does not necessarily mean these are not the key residues contributing to the binding affinity of c-di-GMP. Beyond the potential allosteric regulation of IlvH by c-di-GMP, there is no known relationship between BCAA biosynthesis and c-di-GMP production to the best of our knowledge. However, we subsequently observed frameshift mutations in the catalytic subunit (IlvI) that return this lineage back to the M phenotype with reduced c-di-GMP levels.
The frameshift mutations in ilvI should lead to BCAA auxotrophy, which is likely masked by the nutrient-rich medium utilized in our experimental evolution system. To test this prediction, we grew isolates 5.10a and 5.10b in a chemically defined minimal medium without BCAA. We did not observe any growth unless BCAA was exogenously supplemented (Fig. S3), confirming that our ilvI mutants are unable to synthesize BCAA. We next tested if exogenous supplementation of BCAA alone could restore c-di-GMP production in the ilvI mutant (M isolate 5.10a). Since the ilvI mutant cannot grow in the minimal medium without BCAA supplementation, we also included our M and D ancestors (Fig. 1A) as relative controls in addition to the ilvH mutant (D isolate 5.7) and the parA mutant (M isolate 5.8) from the same lineage. Exogenous  (Fig. 1). c Each mutation was reconstructed in the direct parent and the corresponding ancestor isolate. WGS confirmed the reconstruction of each mutation.
supplementation of BCAA in the ilvI mutant indeed restored c-di-GMP production at comparable levels to both D and the ilvH mutant (D isolate 5.7) (Fig. 7). In contrast, c-di-GMP levels remained the same in both M and the parA mutant (M isolate 5.8) whether or not BCAA was added (Fig. 7). These results show that exogenous supplementation of BCAA does not correlate with c-di-GMP production, at least across the concentration range we have tested (50 to 150 mM). We chose this range for our experiment since higher concentrations did not impact the growth of either the ancestral M or the ilvI mutants in the minimal medium (Fig. S3).
Given that there is no evidence of exogenous BCAA dose dependence on c-di-GMP production across all tested strains of either the M or D phenotype (Fig. 7), it is unlikely that higher concentrations of BCAA will amplify c-di-GMP production. Although BCAA supplementation restores c-di-GMP production in the ilvI mutant in the minimal medium, the reduction of c-di-GMP levels observed in the same ilvI mutant under our nutrient-rich experimental evolution conditions (Fig. 6) is not likely associated with losing the ability to synthesize BCAA. Accordingly, reconstructing the two frameshift ilvI mutations independently in the ancestral D isolate did not cause a shift to the M phenotype, but the initially observed shift to the M phenotype was recapitulated when these mutations were introduced into their direct parental D isolates (Table 1). Since IlvI forms a complex with IlvH, the observed frameshift mutations in ilvI likely decouple the complex and consequently impact IlvH's function to modulate c-di-GMP production. Lineage 5 returns to the D phenotype in four sister branches following the ilvI mutations through mutations in wspE or dgcH (Fig. 4). As already discussed above, we commonly observed activating wspE mutations when the upstream signaling Wsp components are disrupted, and the in-frame deletion in dgcH is also identical to the deletion reported in the D isolate 4.1. These four independent transition steps into the D phenotype clearly indicate that the catalytic activity of IlvI is not required for c-di-GMP production. Lineage 27 was commenced with an engineered DwspCD M isolate (Table S1, isolate 27.0), where we also observed the shift to the D phenotype through a missense mutation in IlvH (A36E) by restoring c-di-GMP production (Fig. 6). In contrast to the M-shifting ParA (D115N) mutation that followed the D-shifting IlvH (G150S) mutation in lineage 5, the transition from IlvH (A36E) to the M phenotype was associated with a mutation which converts the stop codon of wspB with a cysteine residue to extend its open reading frame, but out of frame with wspE (both wspC and wspD were already deleted). A common trend exists  (Table S1) in that multiple mutations first occur in Wsp proteins, which are proceeded by mutations in a non-Wsp DGC (SadC or DgcH respectively), and then unique missense mutations occur in IlvH to stimulate c-di-GMP production. Reconstructing the three IlvH missense mutations in the ancestral M isolate does not lead to a shift to the D phenotype, but reconstructing these mutations in their direct parental M isolate recapitulates the initially observed shift to the D phenotype (Table 1). The common targets of mutations between the two lineages that converge on IlvH mutations are Wsp proteins, which implies that the c-di-GMP-stimulating mutations in IlvH specifically associate with a deregulated Wsp system. Furthermore, the undoing of the IlvH (A36E)induced c-di-GMP production by the stop codon mutation in wspB implies potential IlvH-Wsp interactions that are unlikely to be purely metabolic in nature.
A macro assessment of sequential mutation patterns across all 33 parallel lineages captures potential associations within the c-di-GMP regulome. A long-standing question in the field is why microbes possess so many independent DGC and/or PDE machineries. A systematic analysis of genes that are bioinformatically predicted to be DGCs in P. fluorescens demonstrated that most have little impact on the intracellular c-di-GMP pool (9), and a significant number of DGCs/PDEs in E. coli are present in an inactive form (8). In addition, independent DGCs could regulate a subcellularly localized role (11) that could also contribute to the total intracellular pool to regulate diverse physiological pathways across a gradient (10). Our study indicates that many known or predicted DGCs could be independently stimulated to trigger a common phenotypic shift, which appears to operate across a broad range of c-di-GMP levels (Fig. 1C). We also captured mutations in genes with unknown associations with c-di-GMP modulation that significantly impact the c-di-GMP pool (Table S1). Repetitive patterns of sequential mutations across independently evolving lineages could reveal previously unknown connectivity between discrete regulatory and catalytic elements.
To visualize conserved patterns of mutations across our complex data set, we compiled all sequences of mutation events across the 33 parallel lineages into a single node-edge plot (43) (Fig. 8A). The nodes represent the mutated genes in our experiment. The edges (i.e., lines) of the plot report if the mutation causes an M (blue) or D (orange) phenotype, and the Each node represents a protein, and each edge represents a relationship where one mutation precedes or proceeds another. Orange edges represent M to D shifts, and blue edges represent D to M shifts. The relative size of a node and thickness of an edge is proportional to frequencies within the data set. Proteins that collectively form the ring represent commonly mutated targets in this study and likely have a substantial role in c-di-GMP regulation. Proteins with limited connectivity fall outside the main ring. (B) LC-MS/MS validates that proteins with significant network connectivity, but not previously associated with c-di-GMP, have a significant impact on intracellular c-di-GMP levels. C-di-GMP measurements are reported as the triplicate mean 6 SD. C-di-GMP values of each mutant and its respective parent are summarized in Table S1 with the corresponding ID tag. D isolates are shown in orange, and M isolates are shown in blue. (Student's t test: *, P , 0.05; **, P , 0.01). thickness of the edge reflects the frequency of the observed pattern (Table S2). Proteins that are integrated into this plot-meaning that many edges connect to a specific node-collectively form the main ring. Nodes associated with a single edge fall out of this ring, and mutations in these proteins typically occur at the end of a lineage or a branched sublineage that was not sequenced. Since we did not sequence any progeny isolates from these 10 mutants that fall outside the main ring, we will not discuss them further.
The largest nodes in Fig. 8A are represented by Wsp proteins, DgcH, and other known DGCs such as YfiN, GcbC, MorA, and SadC. The majority of DGCs have blue recursive edges that cause an M phenotype, which indicates that mutations occur in the same protein to shut down its activity, immediately following a mutation that initially stimulated its activity. The high frequency of blue recursive edges associated with each DGC and yellow edges between distinct DGCs collectively indicates that each DGC simply replaces another to restore c-di-GMP production. Blue edges between DGCs are extremely rare, which suggests that individual DGCs do not directly interact with each other. Notably, there are several proteins with both yellow and blue edges that are integrated in this plot and have not been shown to have a role in c-di-GMP regulation or production. We next infer potential functions of these proteins with respect to c-di-GMP modulation.
NarA (N-ethylammeline chlorohydrolase-associated c-di-GMP regulator). Pfl01_3147 is annotated as a hypothetical protein, but we found homology to an N-ethylammeline chlorohydrolase. In isolate 2.3, we observe that the A383S mutation in NarA increases c-di-GMP levels (Fig. 8B), which confirms that this hypothetical protein is indeed biologically active. This NarA mutation follows the M-shifted W90L mutation in WspD, which suggests that NarA is either able to restore Wsp signaling or activate other DGCs to increase c-di-GMP levels. We have no reason to suspect that NarA directly catalyzes c-di-GMP, since it lacks any c-di-GMP-associated domains or motifs. The progeny of this NarA mutant carries a nonsense mutation in wspD, reverting this isolate back to an M phenotype (Table S1). Thus, the A383S mutation in NarA associates specifically with wspD mutations to both restore and shut down c-di-GMP production. N-ethylammeline is a derivative of ammeline, and a recent study theorized that ammeline is likely an essential building block required for the synthesis of cyclic complexes (44). Although the metabolic connections are unknown, there is a strong indication that NarA interacts with WspD at some capacity to modulate c-di-GMP levels.
CdrB. Pfl01_0652 encodes a protein that shares homology with the transporter protein CdrB in P. aeruginosa (PA4624) (45)(46)(47). In our lineage 6, the missense mutation YfiN (G157C) increases c-di-GMP production, which is immediately followed by the CdrB A290T missense mutation to decrease c-di-GMP levels (Fig. 8B). This suggests that an unaltered CdrB or its native transport activity is necessary for YfiN's DGC activity. Mutations that follow the CdrB mutation occur in morA or dgcH, each increasing c-di-GMP as indicated by the D-phenotype (Table S1). Therefore, the negative impact of our CdrB mutation on YfiN could be recovered through the genetic activation of other DGCs, which could also reflect interactions with CdrB. Expression of the cdrB gene increases proportionally with intracellular levels of c-di-GMP during biofilm formation (47), but it has no identified role in c-di-GMP modulation.
CalM (c-di-GMP-associated LysM). The hypothetical gene Pfl01_1895, here calM, is a commonly mutated gene in our study (Table S1). We observed a duplication of five nucleotides within calM that results in a significant decrease in c-di-GMP (Fig. 8B). Two repetitive clusters of CAGCG shift to three clusters to cause a frameshift. This mutation follows a dgcH inframe deletion that dramatically increases c-di-GMP production, indicating that CalM impacts DgcH's activity. Remarkably, three independently evolved progeny isolates in this lineage carry mutations that remove the initial five-nucleotide insertion, returning each isolate back to the D phenotype, likely by restoring DgcH's activity (Table S1). CalM contains domains with homology to the LysM peptidoglycan-binding superfamily and the FimV C-terminal transmembrane domain that is implicated in modulating motility (48,49). The presence of such domains suggests a potential link to peptidoglycan stress, and peptidoglycan stress has been associated with surface-sensing through the Wsp system in P. aeruginosa (33).
RndA (RND transporter). Pfl01_2749 is bioinformatically predicted to encode an RND transporter, which functions as an efflux pump to remove toxins, antibiotics, quorum sensing molecules, and other metabolites (50). We observed that the nonsense mutation Q51* in rndA significantly decreases c-di-GMP levels (isolate 11.2, Fig. 8B). The rndA mutant then shifts to the D phenotype through an in-frame deletion that removes amino acid residues 524 to 537 in DgcH. This specific mutation was observed in multiple lineages to exclusively stimulate c-di-GMP production. Interestingly, a recent study speculates that quorum sensing molecules could contribute to a feedback loop to inhibit c-di-GMP production (51).
DgcX and DgcY (bioinformatically predicted DGCs). Pfl01_3550 and Pfl01_3508 are both bioinformatically predicted to encode DGCs, but their enzymatic function has not yet been demonstrated. Here, we establish that they are both functionally active in c-di-GMP catalysis and annotate them as dgcX and dgcY, respectively, since genes dgcA through dgcW already exist. Pfl01_3550, here dgcX, is a commonly mutated gene in our data set which increases c-di-GMP (Table S1). We observed a blue recursive loop in the node-edge plot for DgcX, which is consistent with other known DGCs (Fig. 8A). We isolated two independently evolved D mutants with the same M350I mutation and an M mutant with an inframe deletion of L289. Phenotypic and motility data support that the observed mutations in dgcX influence the intracellular c-di-GMP pool. The genomic location of dgcY is in a region of hypothetical genes that appear to form an operon. Interestingly, our dgcY mutants carry large deletions of the hypothetical genes within this genomic region (Fig. 5). The deletion of Pfl01_3505, Pfl01_3506, and the 59 untranslated region (UTR) region of dgcY results in a significant increase in c-di-GMP (isolate 6.5a, Fig. 8B), which likely replaces the native promoter region of dgcY with that of Pfl01_3505. C-di-GMP is significantly decreased in the progeny isolate 6.6a (Fig. 8B) after the deletion of Pfl01_3500-Pfl01_3510, which includes dgcY. It is unclear whether the encoded proteins of the hypothetical genes that flank dgcY are functionally related to DgcY. We had initially suspected that these hypothetical genes could encode a complex signal transduction system like Wsp, but we did not detect any homologous Wsp domains or other c-di-GMP-associated motifs.
To confirm whether the observed mutations in these atypical c-di-GMP-associated genes alone could modulate c-di-GMP levels, we reconstructed representative mutations in their respective ancestral backgrounds. None of the engineered isolates produced a phenotypic shift, but reconstructing the same mutations in their direct parental isolates recapitulated the initially observed phenotypic shifts (Table 1). These results clearly indicate that specific mutations in the associated DGCs are first required to impact c-di-GMP levels. Although functional studies are required to prescribe specific mechanisms, diverse cellular processes appear to be integrated with numerous catalytic machineries to modulate the c-di-GMP pool in P. fluorescens, and the missense and in-frame deletion mutations observed in Wsp proteins and other DGCs (Fig. 5) could reflect potential interaction sites.
Conclusion. Although we have extensive knowledge of how specific enzymes make or break c-di-GMP and how this molecule interacts with a large number of regulatory elements to modulate diverse physiological processes, relatively little is known about how C-di-GMP production itself is regulated. This study takes advantage of a repeatedly evolving social behavior in P. fluorescens to demonstrate how numerous proteins could independently and collectively produce remarkable sequences of genetic innovation to regulate c-di-GMP levels. Our extensive library of functionally critical residues reinforces the current model of the Wsp system, which heavily draws functional analogies to the E. coli chemotaxis system and particularly lacks domain resolution in WspB and WspD. Beyond the Wsp system, we have confirmed that several bioinformatically predicted DGCs are indeed enzymatically active and identified unique missense and in-frame deletion mutations in multiple DGCs that fall outside their conserved catalytic domain. These mutations likely reflect regulatory domains that are critical for modulating enzymatic activity. We have also uncovered a set of unexpected or hypothetical proteins that clearly have the capacity to influence c-di-GMP turnover. Three missense mutations in ilvH were identified to independently stimulate c-di-GMP production, which appears to be independent of IlvH's primary role in BCAA biosynthesis but, rather, operates through an unknown interaction with the Wsp system. We have identified the same I-site in IlvH as in most DGCs, which suggests that IlvH physically interacts with c-di-GMP. Similarly, CdrB and RndA homologs in other organisms have been implicated to respond to c-di-GMP levels, but they have never been associated with modulating c-di-GMP production. We have confirmed that the mutations observed in these proteins alone do not impact c-di-GMP levels but do so specifically in the presence of genetically modified DGCS or the Wsp signal transduction system.
An emerging concept here is that the intracellular pool of c-di-GMP, whether it be subcellularly localized or in sum, appears to be modulated through the catalytic activities of many independent DGCs that are in tune with diverse proteins. Such intricacies are masked by the activities of DGCs with high output, but they are clearly hardwired in the genome of P. fluorescens and reflect a complex and underappreciated interconnectivity. Similar relationships could also manifest in other Pseudomonas and related species that share the same DGCs and unexpected proteins described in this study. In particular, very little is known about the functional domains of most DGCs beyond their conserved catalytic domain. The vast majority of the missense mutations or in-frame deletions captured in this study occur outside the catalytic domain to clearly activate or shut down c-di-GMP production. Our extensive library of functionally important mutations presented here should facilitate mechanistic studies of phylogenetically conserved DGCs that remain heavily underexplored.
Parallel experimental evolution. Pure cultures of the ancestral M isolate or engineered M isolates were grown overnight in 5 mL LB broth, and 25-mL volumes were spotted on PAF plates and then incubated at room temperature for 3 to 8 days to allow spreading fans to emerge. To ensure consistency across experiments, PAF plates were always prepared and left to dry at room temperature for 2 days before being inoculated. One fan was randomly sampled from each colony and streaked on LB plates in duplicate. Sampling a fan always produced two isolates with either a mucoid (M) or dry (D) colony phenotype. A single isolated colony exhibiting the D phenotype was randomly chosen to start a new liquid culture and was also frozen down. This process was repeated to sample six independent fans (i.e., six lineages). Three of the six independent D isolates were randomly chosen and cultured as described above to start the next round of experimental evolution. We repeated the same sampling process as above except that new M phenotypes were isolated instead. We bidirectionally selected M and D phenotypes for as many as 13 rounds for a single lineage, but the number of sampled fans at each step was reduced over time in a varying manner across independent lineages. In total, approximately 600 independently evolved isolates were frozen down across 33 parallel lineages. Every isolate was subjected to a motility assay in an LB plate with 0.25% agar, and the M and D ancestors were used as controls (29). In addition, 25 mL of pure culture of each isolate, the ancestral isolate with the opposite phenotype, and a 1:1 mixture of both were spotted onto PAF plates in triplicate and assessed for cooperative radial spreading (Fig. 1A) to filter out hypermotile and noncooperative mutants (29). Colony diameters were measured and recorded after 4 days of incubation at room temperature and considered cooperative if the measured diameter of the 1:1 colony exceeded the diameters of both the mutant and the respective ancestral isolate. Isolates that failed this test were removed from the study.
Reconstruction of naturally derived mutations in ancestral and direct parental backgrounds. Experimentally evolved isolates were whole-genome sequenced to identify the causal mutations as described below. We reconstructed select naturally acquired mutations in the respective ancestral background and also in the direct parental background using previously described methods (29). PCR primers were developed from the sequence data to amplify mutations of interest and are reported in Table S3. Isolates were grown overnight and underwent DNA extractions (Qiagen). Dream Taq DNA polymerase (Thermo Fisher) was used in PCRs with purified DNA to generate insert fragments carrying the mutations, which ranged in size from 456 bp to 934 bp. PCR was carried out using the recommended conditions of the Dream Taq amplification kit for amplicons below 2,000 bp. Inserts were cloned into the pGEM-T Easy vector system (Promega) and then subcloned into the NotI site of the suicide plasmid pSR47s (54). Suicide plasmids were heat-shocked into E. coli S17.1lpir, mated with the target isolates reported in Table S3 on LB plates at 30°C for 16 h, and then plated onto PMM plates supplemented with 50 mg/mL kanamycin. We observed that ilvI mutants could not grow on PMM and therefore replaced the PMM selection plate with an LB plate supplemented with 50 mg/mL kanamycin and 100 mg/mL ampicillin for ilvI mutant constructs. Isolated colonies were grown overnight, serially diluted, and plated out on LB plates supplemented with 15% sucrose (wt/vol) for counterselection. Isolated colonies were placed in cryo-storage as described above. Crude DNA extracts were made for each isolate, where 100 mL of overnight culture was boiled at 100°C and centrifuged at 15,000 relative centrifugal force (RCF) for 30 s, and 10 mL of supernatant was used for PCR using the primers reported in Table S3. PCR Cyclic-di-GMP Modulation in Pseudomonas fluorescens mSystems products underwent gel extraction and subsequent Sanger sequencing (primers in Table S3) to confirm that the mutation was successfully engineered into the target isolate. Upon PCR confirmation, pure cultures of cells from cryo-storage underwent DNA extraction (Qiagen) and were submitted for whole-genome sequencing (WGS), which confirmed that no additional mutations had accumulated in the target isolates during mutant construction. DNA extraction, whole-genome sequencing, and variant calling. Pure cultures were grown overnight in 5 mL LB broth at 30°C. Then, 1 mL was spun down to form a pellet for DNA extraction. Pelleted cells were mechanically lysed with the Qiagen TissueLyser LT at 50 oscillations/s for 10 min. The Qiagen DNeasy UltraClean microbial DNA extraction kit was used following the standard operating procedures. Samples were then eluted into 30 mL H 2 O after a 5-min incubation at 4°C with the elution medium and stored at 220°C prior to shipment. Genomic DNA library preparations were carried out by the Microbial Genome Sequencing Center (Pittsburgh, PA) following standard operating procedures. Pooled and indexed samples were sequenced on the NextSeq 2000 platform with 150-bp paired-end reads and were demultiplexed prior to data delivery. The Trimmomatic algorithm was used to group and quality filter paired reads at a Phred score of 20 using a 4-bp sliding window. The BreSeq algorithm was used to align the paired reads to the P. fluorescens Pf0-1 reference genome (CP000094.2) and evaluate for single nucleotide polymorphisms (SNPs), insertions, deletions, or rearrangements with a set minimum of 20Â coverage. Output sequence data from BreSeq was manually evaluated for each isolate to determine the order of mutations for a given lineage, where each M mutation was inferred by a common mutation found between its two sister progeny D isolates.
Alignment of mutation data to annotated domains. Domain analyses of Wsp proteins and DGCs were conducted and rendered as previously described (19). Peptide sequences of the predicted or known DGCs were obtained from the Pseudomonas Genome Database (55). NCBI's Conserved Domain Database (CDD) was used to evaluate the functional domains (56). Mutation data of Wsp proteins from this study were manually mapped to the rendered consensus sequence of the respective Wsp protein (19).
C-di-GMP extraction and LC-MS/MS. A colony of each isolate was transferred to PF broth and grown overnight. Overnight cultures were diluted to an optical density at 600 nm (OD 600 ) of 0.04 in triplicate and incubated for 2 to 4 h until the samples reached an OD 600 of 0.5 and then promptly processed for c-di-GMP extraction. We verified by microscopy that there was no clumping of cells at the time of c-di-GMP extraction. C-di-GMP extraction and quantification procedures followed the protocols established by the Michigan State University Research Technology Support Facility (MSU-RTSF): MSU_MSMC_009 protocol for di-nucleotide extractions and MSU_MSMC_009a for LC-MS/MS, with the following modifications. All the steps of c-di-GMP extraction were carried out at 4°C unless noted otherwise. A sample of 1 mL of 0.5 OD 600 culture was centrifuged at 15,000 RCF for 30 s, and cell pellets were resuspended in 100 mL of ice-cold acetonitrile/methanol/ water (40:40:20 [vol/vol/vol]) extraction buffer supplemented with a final concentration of 0.01% formic acid and 25 nM c-di-GMP-fluorinated internal standard (InvivoGen catalog [cat. no.] 1334145-18-4). Pelleted cells were mechanically disturbed with the Qiagen TissueLyser LT at 50 oscillations/s for 2 min. Resuspended slurries were incubated at 220°C for 20 min and pelleted at 15,000 RCF for 15 min at 4°C. The supernatant was transferred to a prechilled tube and then supplemented with 4 mL of 15% wt/vol ammonium bicarbonate (Sigma-Aldrich) buffer for stable cryo-storage. Extracts were stored at 280°C for less than 2 weeks prior to LC-MS/MS analysis at MSU-RTSF. Sample degradation was observed with repeated freeze-thaw cycles, so extracts were never thawed until LC-MS/MS analysis. All samples were evaporated under vacuum (SpeedVac, no heat) and redissolved in 100 mL mobile phase (10 mM tryptose blood agar [TBA]/15 mM acetic acid in water/methanol, 97:3 vol/vol, pH 5). LC-MS/MS quantification was carried out with the Waters Xevo TQ-S instrument with the following settings: c-di-GMP-F at m/z transition of 693 to 346 with cone voltage of 108 and collision voltage of 33; c-di-GMP at m/z transition of 689 to 344 with cone voltage of 83 and collision voltage of 31. C-di-GMP data were normalized to 25 nM c-di-GMP-F by MSU-RTSF to account for sample matrix effects and sample loss during preparation. OD 600 measurements of the initial samples were used to report the quantified c-di-GMP as nM/OD and visualized with GraphPad Prism 9.
BCAA supplementation assay and c-di-GMP quantification. Pure cultures were grown overnight in 5 mL LB broth. Then, 5 mL of each culture was transferred into 5 mL PMM. Cultures were supplemented with the BCAA valine, leucine, and isoleucine as previously described (57). Assessment of growth with BCAA supplementation contained either 0 mM (loading control), 180 mM, 200 mM, 220 mM, 250 mM, or 300mM in triplicate. Samples were incubated at 30°C for 24 h with 250 rpm shaking. Cultures were manually pipetted to resuspend cells, and 300 mL of culture was transferred to a 96-well black-wall clear-bottom plate for OD 600 quantification. OD was quantified using the SpectraMax plate reader, and background (determined from blanks) was removed from each quantified sample. Data analysis and visualization of growth during BCAA supplementation were conducted in Excel. To evaluate the relationship between BCAA and c-di-GMP, pure cultures were grown for 24 h at 30°C with 250 rpm shaking in PMM liquid medium containing either 0 mM, 50 mM, 100 mM, or 150 mM of supplemented BCAAs. OD 600 measurements were taken and samples were diluted to an OD 600 of 0.04 in PMM at the same BCAA concentrations. Samples were incubated for 24 h under the same conditions stated and then sonicated (time 0.00:20, pulse 01 01, amplitude 20%) to disrupt cell aggregates. OD 600 measurements occurred every hour until the OD 600 was approximately equal to 0.25. C-di-GMP analyte extraction was carried out as stated above with the modification that 2-mL cultures were used for cell lysis. Quantified values of c-di-GMP are reported as nM/OD 600 and were visualized with GraphPad Prism 9.
Node-edge analysis and visualization. The data reported in Table S1 were used to generate the "to-from" table, which assesses mutation patterns at the gene/protein level (43). Counts were determined by the frequency of the to-from patterns calculated from Table S1, where a mutation in a given gene/protein (to) occurred after a mutation in another gene/protein (from). This table was uploaded to Cytoscape to visualize the nodes (to genes/proteins) and edges (counts). Default renderings of the plot reported node centrality (commonly mutated targets) through coordinate placement on the graph. This metric was converted to node size, with larger nodes representing more commonly mutated targets. Nodes were manually aligned for visual clarity, with nodes connected by two or more edges being placed in the main ring and nodes with a single edge being placed outside the ring.
Data availability. Raw sequence data of the evolved isolates from this study are accessible at NCBI's Sequence Read Archive under BioProject accession number PRJNA877395. The raw data from c-di-GMP measurements, BCAA-supplemented growth measurements, and BreSeq summaries of the evolved isolates are available at https://github.com/wook-kim-lab/cyclic-di-gmp-modulation.

SUPPLEMENTAL MATERIAL
Supplemental material is available online only.