Predicting mutational routes to new adaptive phenotypes

Predicting evolutionary change poses numerous challenges. Here we take advantage of the model bacterium Pseudomonas fluorescens in which the genotype-to-phenotype map determining evolution of the adaptive ‘wrinkly spreader’ (WS) type is known. We present mathematical descriptions of three necessary regulatory pathways and use these to predict both the rate at which each mutational route is used and the expected mutational targets. To test predictions, mutation rates and targets were determined for each pathway. Unanticipated mutational hotspots caused experimental observations to depart from predictions but additional data led to refined models. A mismatch was observed between the spectra of WS-causing mutations obtained with and without selection due to low fitness of previously undetected WS-causing mutations. Our findings contribute toward the development of mechanistic models for forecasting evolution, highlight current limitations, and draw attention to challenges in predicting locus-specific mutational biases and fitness effects.


Introduction
Adaptation requires the realization of beneficial mutations. As self-evident as this may be, predicting the occurrence of beneficial mutations and their trajectories to improved fitness is fraught with challenges (Lässig et al., 2017). Nonetheless progress has been made for phenotypically diverse asexual populations subject to strong selection. Effective approaches have drawn upon densely sampled sequence data and equilibrium models of molecular evolution to predict amino acid preferences at specific loci (Luksza and Lässig, 2014). Predictive strategies have also been developed based on selection inferred from the shape of coalescent trees (Neher et al., 2014). In both instances the models are coarse-grained and sidestep specific molecular and mutational details.
There is reason to by-pass molecular details: mutation, being a stochastic process, means that for the most part details are likely to be idiosyncratic and unpredictable. But an increasing number of investigations give reason to think otherwise -that adaptive molecular evolution might follow rules (Yampolsky and Stoltzfus, 2001;Pigliucci, 2010;Stern, 2013;Laland et al., 2015). This is particularly apparent from studies of parallel molecular evolution (Colosimo et al., 2005;Woods et al., 2006;Ostrowski et al., 2008;Flowers et al., 2009;Meyer et al., 2012;Tenaillon et al., 2012;Zhen et al., 2012;Herron and Doebeli, 2013;Galen et al., 2015;Bailey et al., 2017;Kram et al., 2017;Stoltzfus and McCandlish, 2017), and particularly from studies that show parallel evolution to be attributable -at least in part -to factors other than selection (McDonald et al., 2009;Lind et al., 2015;Bailey et al., 2017).
A standard starting position for predicting adaptive evolution recognises the importance of population genetic parameters including mutation rate, generation time, population size, selection and more recently information on the distribution of beneficial fitness effects. An often used approach appeals to 'origin-fixation' models that estimate the probability that selection will realise particular mutational trajectories based on the assumption that the rate of evolution is a function of the rate at which mutations originate, multiplied by their fixation probability (for review, see (McCandlish and Stoltzfus, 2014), for application, see for example (Weinreich et al., 2006)). Leaving aside a number of simplifying and restrictive assumptions, population genetic approaches such as those encompassed by origin-fixation models require pre-existing knowledge of (or assumptions concerning) mutation rate to a particular phenotype and selection coefficients (obtained in the absence of frequency dependent effects) and thus have predictive power only in specific situations and with the benefit of hindsight (McCandlish and Stoltzfus, 2014).
Looking to the future, there is hope that one day it might be possible to predict the course of evolutionary change in response to specific environmental challenges in much the same way as gene function, regulation and interactions can be predicted today based on knowledge of nucleotide sequence data. A central issue is to define the necessary information. Mechanistic understanding of the connection between genotype and phenotype combined with knowledge of the genomic bases of mutational bias offer opportunities for progress.
Mutations arise randomly with respect to utility, but genetic architecture can influence the translation of mutation into phenotypic variation: the likelihood that a given mutation generates phenotypic effects depends on the genotype-to-phenotype map (Alberch, 1991;Gompel and Prud'homme, 2009;Stern and Orgogozo, 2009;Rainey et al., 2017). Thus, the functions of gene products and their regulatory interactions provide information on likely mutational targets underpinning particular phenotypes. For example, consider a hypothetical structural determinant subject to both positive and negative regulation and whose over-expression generates a given adaptive phenotype. Assuming a uniform distribution of mutational events, mutations in the negative regulator (and not the positive activator) will be the primary cause of the adaptive phenotype. This follows from the fact that loss-of-function mutations are more common than gain-of-function mutations. Indeed, an emerging rule indicates that phenotypes determined by genetic pathways that are themselves subject to negative regulation are most likely to arise by loss-of-function mutations in negative regulatory components (McDonald et al., 2009;Tenaillon et al., 2012;Lind et al., 2015;Fraebel et al., 2017). eLife digest Predicting evolution might sound like an impossible task. The immense complexity of biological systems and their interactions with the environment has meant that many biologists have abandoned the idea as a lost cause. But despite this, evolution often repeats itself. This repeatability offers hope for being able to spot in advance how evolution will happen. To make general predictions, it is necessary to understand the mechanisms underlying evolutionary pathways, and studying microbes in the laboratory allows for real-time experiments in evolution.
One of the best studied microbes for experimental evolution is Pseudomonas fluorescens, which repeatedly evolves flattened wrinkled colonies instead of round smooth ones when there is limited oxygen. The underlying molecular pathways that lead to this change have been studied in detail. Lind et al. developed mathematical models to predict how often the three most common pathways would be used and which genes were most likely to be mutated. After controlling for the effects of natural selection and refining the models to take into account mutation hotspots, Lind et al. were able to accurately predict the genes that would be targeted by mutations.
The findings suggest that biologists need not lose hope when it comes to the goal of predicting evolution. A deep understanding of the molecular mechanisms of evolutionary changes are essential to predicting the mutations that lead to adaptive change. The results are an important first step towards forecasting organisms' responses to changing conditions in the future. In the short term, this is important for medical issues, including antibiotic resistance, cancer and immune receptors. In the long term, predicting the course of evolution could be essential for survival of life on the planet.
Mutation is not equally likely at each nucleotide of a given genome (Lind and Andersson, 2008;Lynch, 2010;Seier et al., 2011;Foster et al., 2015;Reijns et al., 2015;Sankar et al., 2016;Stoltzfus and McCandlish, 2017). Numerous instances of mutational bias have been reported. Prime examples are simple sequence repeats such as homopolymeric nucleotide tracts or di-, triand tetrameric repeats that mutate at high frequency via slipped strand mispairing (Levinson and Gutman, 1987). These readily identifiable sequences define contingency loci in obligate human pathogens and commensals (Moxon et al., 1994) and are widespread in eukaryotic genomes (Tautz and Renz, 1984). The behaviour of contingency loci can be further modulated by defects in components of methyl-directed mismatch repair systems (Richardson and Stojiljkovic, 2001;Martin et al., 2004;Hammerschmidt et al., 2014;Heilbron et al., 2014).
Certain palindromic structures also lead to mutational bias (Viswanathan et al., 2000;Lovett, 2004) and promote amplification events that can increase target size for additional mutations (Roth et al., 1996;Kugelberg et al., 2010;Reams and Roth, 2015) Transition-transversion bias (Stoltzfus and McCandlish, 2017) and elevated mutation rates at CpG sites (Galen et al., 2015) can further skew the distributions of mutational effects. Additional bias arises from the chromosomal neighbourhood of genes under selection (Steinrueck and Guet, 2017), the location of genes with regard to interactions with DNA replication/transcription machineries (Sankar et al., 2016), and environmental factors that affect not only mutation rate but also the spectra of mutational events (Krasˇovec et al., 2017;Maharjan and Ferenci, 2017;Shewaramani et al., 2017).
Beyond the genotype-to-phenotype map and mutational biases, predicting adaptive evolution requires ability to know a priori the fitness effects of specific mutations. At the present time there is much theoretical and empirical interest in the distribution of fitness effects (DFE) (Eyre-Walker and Keightley, 2007) -and particularly the DFE of beneficial mutations (Orr, 2005) -because of implications for predicting the rate of adaption and likelihood of parallel evolution (de Visser and Krug, 2014), but knowledge of the shape of the distribution is insufficient to connect specific mutations to their specific effects, or to their likelihood of occurrence. Such connections require a means of knowing the connection between mutations and their environment-specific fitness effects. This is a tall order. A starting point is to understand the relationship between all possible mutational routes to a particular phenotype and the set that are realised by selection.
Here we take a bacterial system in which the genetic pathways underpinning evolution of the adaptive 'wrinkly spreader' (WS) type are known and use this to explore the current limits on evolutionary forecasting. Pseudomonas fluorescens SBW25 growing in static broth microcosms rapidly depletes available oxygen establishing selective conditions that favour mutants able to colonise the air-liquid interface. The most successful mutant-class encompasses the WS types (Ferguson et al., 2013;Lind et al., 2017b). These types arise from mutational activation of diguanylate cyclases (DGCs) that cause over-production of the second messenger c-di-GMP (Goymer et al., 2006;McDonald et al., 2009), over-production of an acetylated cellulose polymer (Spiers et al., 2002;Spiers et al., 2003) and ultimately formation of a self-supporting microbial mat ( Figure 1A). McDonald et al. (McDonald et al., 2009) showed that each time the tape of WS evolution is rerun mutations generating the adaptive type arise in one of three DGC-encoding pathways (Wsp, Aws, or Mws) ( Figure 1A). Subsequent work revealed that when these three pathways are eliminated from the ancestral type evolution proceeds along multiple new pathways (Lind et al., 2015) resulting in WS types whose fitnesses are indistinguishable from those arising via mutations in Wsp, Aws, or Mws. Parallel evolution of WS involving preferential usage of the Wsp, Aws and Mws pathways is therefore not explained by selection: repeated use of Wsp, Aws and Mws stems from the fact that these pathways are subject to negative regulation and thus, relative to pathways subject to positive regulation, or requiring promoter-activating mutations, gene fusion events, or other rare mutations, present a large mutational target (Lind et al., 2015).
Given repeatability of WS evolution, knowledge of the Wsp/Aws/Mws pathways, plus genetic tools for mechanistic investigation -including capacity to obtain WS mutants in the absence of selection -the WS system offers a rare opportunity to explore the extent to which knowledge of the genotype-to-phenotype map can lead mechanistic models for evolutionary forecasting. Our findings show that short-term mechanistic-level predictions of mutational pathways are possible, but also draw attention to challenges that stem from current inability to a priori predict locus-specific mutational biases and environment-specific fitness effects.

Results
Obtaining an unbiased measure of pathway-specific mutation rates to WS Knowledge of the rate at which mutation generates WS types via each of the Wsp, Aws and Mws pathways -unbiased by the effects of selection -provides a benchmark against which the predictive power of null models can be appraised. To achieve such measures we firstly constructed a set of genotypes containing just one of the three focal pathways: PBR721 carries the Wsp pathway but is devoid of Aws and Mws, PBR713 carries the Aws pathway but is devoid of Wsp and Mws, while PBR712 harbours the Mws pathway but is devoid of Wsp and Aws. Into each of these genotypes a promoterless kanamycin resistance gene was incorporated immediately downstream of the promoter of the cellulose-encoding wss operon and transcriptionally fused to an otherwise unaffected wss operon ( Figure 1B).
In the ancestral SM genotype the cellulose promoter is inactive in shaken King's Medium B (KB) broth (Spiers et al., 2002) and thus the strain is sensitive to kanamycin. When a WS-causing mutation occurs, the wss promoter becomes active resulting in a kanamycin-resistant WS type (Fukami et al., 2007;McDonald et al., 2011). Individual growth of this set of three genotypes in shaken KB, combined with plating to detect kanamycin-resistant mutants, makes possible a fluctuation assay (Luria and Delbrück, 1943;Hall et al., 2009) from which a direct measure of the rate at which WS mutants arise can be obtained. Importantly, because WS types are maladapted in shaken broth culture, the screen for kanamycin-resistant clones allows the pathway-specific mutation rate to WS to be obtained without the biasing effects of selection for growth at the air-liquid interface ( Figure 1B). The results are shown in The mutation rate was highest for the Aws pathway (6.5 Â 10 À9 ); approximately double that of Wsp (3.7 Â 10 À9 ) and an order of magnitude higher than that of the Mws pathway (0.74 Â 10 À9 ) ( Figure 2). The rate at which WS mutants arose from the ancestral genotype in which the three pathways are intact (11.2 Â 10 À9 ) was approximately the sum of the rates for the three pathways (11.0 Â 10 À9 ) confirming that the Wsp, Aws and Mws pathways are the primary routes by which WS types evolve (Lind et al., 2015). That the Aws pathway has the greatest capacity to generate WS is surprising given the smaller target size (three genes and 2.3 kb compared to seven genes (8.4 kb) in the Wsp pathway).
The genotype-to-phenotype map underpinning WS evolution Much is known about the function and interactions among components of each of the three focal pathways. This knowledge allows development of models that capture the dynamic nature of each pathway and thus allow predictions as to the likelihood that evolution will precede via each of the three mutational routes. An unresolved issue is the extent to which these models match experimental findings. Following a brief description of each pathway we describe the models.
The 8.4 kb Wsp pathway is a chemotaxis-like system (Goymer et al., 2006;Güvener and Harwood, 2007;Rö mling et al., 2013;Micali and Endres, 2016) composed of seven genes with the first six genes (wspA-wspF) being transcribed as a single unit and the last (wspR from its own promoter (Bantinaki et al., 2007). WspA (PFLU1219) is a methyl-accepting chemotaxis (MCP) protein that forms a complex with the CheW-like scaffolding proteins WspB (PFLU1220) and WspD (PFLU1222). WspA senses environmental stimuli and transmits the information via conformational changes in the WspA/WspB/WspD complex to effect activity of WspE (PFLU1223), a CheA/Y hybrid histidine kinase response regulator. WspE activates both the WspR (PFLU1225) diguanylate cyclase (DGC) and the CheB-like methylesterase WspF (PFLU1224) following transference of an active phosphoryl group. The activity of WspA is modulated by methylation by the constitutively active CheRlike methyltransferase WspC (PFLU1221) that transfers methyl groups to conserved glutamine residues on WspA. The demethylase WspF serves to remove these groups when in the phosphorylated active form. WS mutants are known to arise by mutations in the WspF negative regulator and also in the WspE kinase (McDonald et al., 2009). In vitro manipulations of WspR that abolish repression of the DGC domain by the response regulator domain are known, but have never been observed to occur in experimental populations (Goymer et al., 2006). The 2.3 kb aws operon contains three genes transcribed from a single promoter (awsXRO). Homologous genes in Pseudomonas aeruginosa (yfiRNB, PA1121-1119) have been characterised in detail (Malone et al., 2010;Malone et al., 2012;Xu et al., 2016). The outer membrane lipoprotein AwsO (PFLU5209) has an OmpA domain, a signal peptide and binds to peptidoglycan. AwsO is thought to be the sensor whose activity is modulated in response to envelope stress (Malone et al., 2012). AwsO sequesters the periplasmic protein AwsX (PFLU5211) at the outer membrane. AwsX functions as a negative regulator of the DGC AwsR (PFLU5210) in the inner membrane. Both increased binding of AwsX to AwsO or loss of negative regulation by inactivation of the interaction between AwsX and AwsR can lead to WS (McDonald et al., 2009;Malone et al., 2010;Malone et al., 2012).
The 3.9 kb mwsR gene (PFLU5329) is known as morA (PA4601) in Pseudomonas aeruginosa, and encodes a predicted membrane protein with both a DGC domain that produces c-di-GMP and a phosphodiesterase (PDE) domain that degrades c-di-GMP. Little is known of the molecular details determining its function, but both catalytic domains appear to be active (Phippen et al., 2014). Deletion of the PDE domain results in a WS phenotype with activity being dependent on a functional DGC domain (McDonald et al., 2009).

Development of mechanistic models predicting mutational routes to WS
If the specific effects of changing each nucleotide (and sets of nucleotides) were known then models for each pathway would not be required. Here, we show how knowledge of genetic architecture can be used to build models that predict the likelihood that mutations generating WS types will arise in a given pathway -and even in specific genes. In the following section we present four null models that incorporate increasing levels of information concerning the genotype-to-phenotype map. The goal of these models is two-fold: firstly to demonstrate that incorporation of knowledge of genetic architecture allows development of models with explanatory value and secondly, to define minimal necessary information for reliable forecasting. The results are summarised in Figure 3, which displays the experimental data from Figure 2 ( Figure 3A), along with predictions from each of the null models.
Null Model I is intentionally naïve. It uses sequence length as a proxy for mutational target size, but ignores genetic organisation, function of predicted proteins and interaction among proteins. The model assumes that mutational target size is proportional to the number of nucleotides at a given locus and thus the probability that a given pathway is used to generate WS relative to another is simply the ratio of the probability of generating WS for the two focal pathways. For any one pathway the probability that a mutation generates WS is given by 1-(1-p) n , where p is the probability of a mutation at a nucleotide and n is the number of nucleotides in the pathway. If the mutation probability is low such that the expected number of mutations in a pathway is below 1, that is, np <<1 then the binomial approximation can be used: 1-(1-p) n = np. Thus the probability that evolution follows the Wsp pathway over the Aws pathway is: 8400 p 1 / 2300 p 2 = 3.65 p 1 /p 2 , where p 1 and p 2 are the mutation rates for each pathway. Assuming equivalency of mutation rate, p 1 = p 2 , evolution is predicted to proceed via the Wsp pathway 3.65 times more often than via Aws, with evolution predicted to proceed via Mws 1.65 times more often than Aws. Comparison with experimental data shows a departure both in terms of the priority of pathways used by evolution and the frequency of pathway usage ( Figure 3A versus 3B).
Null Model II builds on Model I but only in a marginal sense. It recognises that nucleotides defining loci of interest are organised into genes, and therefore adopts gene number as a proxy for mutational target size. As above, the probability that a pathway is used is 1-(1-p) n = np, but in this instance n is the number of genes. Model II predicts that mutations in Wsp generate WS types 2.33 times more often than mutations in Aws, with mutations in Mws being 3-fold less likely to generate WS compared to mutations in Aws ( Figure 3C). This marginal adjustment makes little difference to the fit between experimental data and predictions.
Past work has shown the explanatory value of information that comes from knowledge of gene (protein) function and interactions (McDonald et al., 2009;Lind et al., 2015). These form the basis of Null Model III. The relevant functions and interactions are depicted in Figure 4 as reaction diagrams that reflect how changes in different interactions affect the production of WS types. Organizing interactions within each pathway according to reactions has the advantage that it allows for a standard mathematical description of the biochemical dynamics using differential equations (see ). An additional advantage is that such an approach allows ready incorporation of new experimental data including interaction type and interaction strength.
One immediate consequence of this approach is that production of WS types is described entirely through biochemical reactions. Knowing whether a WS-type is generated by a particular mutation amounts to determining which reaction rates are altered by the mutation and whether those changes affect the likelihood of producing WS types. This leads to an approach where mutations in components can be classified according to their effects on reaction rates: enabling mutations increase reaction rates, whereas disabling mutations decrease reaction rates.
Based on the mathematical models for the pathways, each reaction rate can be altered (either increased or decreased) to generate higher levels of the key WS effector (the active DGC). For example, in the Wsp pathway, if reactions 2 or 6 experience a disabling mutation, or alternatively, if any of the reactions 1, 3, 4, or 5 experience an enabling mutation, then the level of activated (phosphorylated) WspR is increased. For now, it is assumed that increased activity of the relevant DGC  generates a WS type. Reactions whose enabling/disabling mutations increase the amount of WS effector are termed enabling/disabling reactions. By comparing the number of disabling and enabling reactions in different pathways it is possible to calculate the relative likelihood that evolution uses each mutational pathway. The Wsp pathway encompasses two disabling and four enabling reactions, whereas the Aws pathway is defined by one disabling and three enabling reactions. Mws encompasses a single disabling reaction and two enabling reactions. If disabling mutations are as likely as enabling mutations, then the likelihood that a pathway will be used by evolution is simply the ratio of the total number of enabling and disabling reactions: 6:4:3 (Wsp:Aws:Mws). If instead, disabling mutations are much more likely than enabling mutations, then enabling reactions can be ignored and the likelihood that a pathway will be used by evolution becomes simply the ratio of the number of disabling reactions: 2:1:1 (Wsp:Aws:Mws) ( Figure 3D left hand panel). If the reverse is true, then the likelihood of use is a ratio of the number of enabling reactions: 4:3:2 (Wsp:Aws:Mws) ( Figure 3D right hand panel). In all cases, mutations in the Wsp pathway are predicted to be 1.3-2 times more likely to generate WS than mutations in the Aws pathway, with the Aws pathway being the target of evolution 1-1.5 times more often than the Mws pathway. Although still inflating the importance of the Wsp pathway relative to Aws, Mws is more prominent than in null Model II. Null Model IV takes a further step toward mechanistic accuracy by endowing interactions with pleiotropic and continuous effects. This removes two simplifying assumptions of null Model III likely to limit predictive power. In Model III mutations affect one reaction at a time, but mutations in certain components, for example wspA, can affect more than a single reaction (WspA appears in reactions 1-4). Additionally, in Model III, changes in reaction rate are assumed to be binary, however reactions may have a range of effects on the evolution of WS types.
To accommodate pleiotropic effects, null Model IV systematically considers all combinations of enabling and disabling changes to reaction rates and determines the likelihood that a WS type is generated. An example of one set of the possible mutations (m i ) in Wsp is 1,-1, 0, 0, 0, 0 (an increase in r 1 , a decrease in r 2 , but no change in r 3 , r 4 , r 5 , or r 6 ( Figure 3A)). Since the Wsp pathway has six reaction rates this amounts to 3 6 or 729 total combinations. However, note that reaction 3 does not share any reactants or products with reactions 5 or 6. Thus, mutations such as 0, 0, 1, 0, 1, 0 or 0, 0, 1, 0, 0, 1 are not considered because they require mutations in two separate genetic components.
To accommodate a range of effects null Model IV simulates enabling/disabling changes of different magnitudes and determines the resulting effect on the respective effector DGC (see Materials and methods). Briefly, the approach addresses the lack of information concerning biochemical reaction rates and molecular concentrations in the mathematical models describing WS-producing pathway dynamics. By repeatedly sampling from the space of all possible reaction rates, initial concentrations, and magnitudes of effects, this approach computes the probability that a particular set of mutations (m i ), for example 1,-1, 0, 0, 0, 0, results in a wrinkly spreader. This probability is represented as the conditional probability P (WS |m i 2 Wsp), which motivates a Bayesian formulation to compare the relative probability that the different pathways produce WS. To this end, the probability that a particular pathway will be used is decomposed into two terms: the probability that a particular set of mutations (m i ) occurs in Wsp (or Aws, or Mws) represented as P (m i 2 Wsp) and the probability that those mutations give rise to a wrinkly spreader represented as P (WS |m i 2 Wsp) (or Aws, or Mws).
To estimate P (m i 2 Wsp) we assume fixed probabilities of enabling and disabling changes and compute the product. Thus, the probability of m i = 1,-1, 0, 0, 0, 0 is p e p d (1 À p e À p d ) 4 , where p e is the probability of a mutation with an enabling effect and p d is the probability of a mutation with a disabling effect. Recognising the value of accommodating the possibility of localised mutational bias we note that p e and p d can be adjusted for the affected reactants. The second term, P (WS |m i 2 Wsp), relies on our sampling methodology and describes the probability that a set of disabling/ enabling changes of different magnitudes will yield a WS type (see Materials and methods).
Despite the mechanistic advances incorporated into null Model IV the Wsp pathway is still predicted to be the pathway most commonly used by evolution. The extreme cases in which disabling mutations are more probable than enabling mutations (and vice versa) are shown in Figure 3E. The results of simulations for the full range of p d and p e values are shown in Figure 5. Figure 5A shows that the Wsp pathway is predicted to be the target of mutation 1.5-1.9 times more often than the Aws pathway while Figure 5B shows that the Mws pathway is predicted to be the target of mutation 0.4-0.5 times less often that the Aws pathway. While these results agree with the experimental data showing Mws to be least likely pathway to be followed, the predictions are at odds with the data showing WS types to be twice as likely to arise from mutation in Aws, versus Wsp. The causes of this discrepancy are accounted for in the next section 'analysis of mutants reveals sources of mutational bias'.
The Bayesian approach of null Model IV has additional utility: it allows prediction of genes likely to be affected by mutation. Predictions as to specific genetic targets come from appraisal of the relative importance of each reaction in terms of generating WS types (Figure 5-figure supplements 1 and 2). While it is recognised that a majority of WS mutations arise from defects in negative regulators of DGCs, such as WspF and AwsX (McDonald et al., 2009;Lind et al., 2015), further predictions are possible based on impacts of alterations in gene function on reaction rates. For example, with reference to the Wsp pathway (Figure 5-figure supplement 1), there are two reactions (2 and 6) that are affected by WspF function: r 2 describes the rate of removal of methyl groups from the WspA signalling complex and r 6 the rate at which WspF is activated by transfer of active phosphoryl groups from the WspE kinase. Loss-of-function (disabling) mutations being much more common than gain-of-function (enabling) mutations means that WspF, WspA and WspE are all likely targets. The null model predicts that in the area of parameter space in which Wsp is most likely compared to Aws, 45% of the time WS will be generated when the second reaction, r 2 , is altered ( Figure 5-figure supplement 1). The same is true for reaction r 6 . Thus the presence of a negative regulator is predicted to extend the mutational target size well beyond the gene itself. This is also true for Aws, where r 3 is the main contributor to the WS phenotype in the case where disabling change is more common than enabling change. Here mutations are predicted not only in the negative regulator AwsX, but also in the interacting region of the DGC AwsR ( Figure 5-figure supplement 2).
Loss-of-function mutations in negative regulators and their interacting partners are not the only predicted targets. For Wsp r 1 , r 3 , r 4 , and r 5 are altered approximately 5% of the time in the parameter region where disabling mutations are more common than enabling mutations, but contribute more when the rate of enabling mutations is increased ( Figure 5-figure supplement 1). Enabling mutations based on null Model IV are likely to be found in WspC (increasing r 1 ), WspABD (increasing r 3 ), WspABD/WspE (increasing r 4 ) and WspR (increasing r 5 ) ( Figure 4A). For Aws, enabling mutations are predicted to increase reaction r 1 by mutations causing constitutive activation of AwsO, r 2 A B C 10 -7 10 -6 10 -5 10 -4 10 -3 10 -2 Probability of enabling change In summary, high rates of WS mutations are predicted for wspF, wspE, wspA, awsX and awsR with lower rates for wspC, wspR and awsO. Several of these predictions sit in accord with previous experimental observations, however, notable are predictions that evolution might also target wspA and wspR, two genes that have not previously been identified as mutational causes of WS types (McDonald et al., 2009).

Analysis of mutants reveals sources of mutational bias
There are several reasons why predictions from the models might be out of kilter with experimental data on mutation rates. We firstly looked to the distribution of WS generating mutations among the 109 mutants collected during the course of the fluctuation assays. Of the 109 mutants, 105 harboured a mutation in wsp (46 mutants), aws (41 mutants) or mws (18 mutants) ( Figure 6A, Figure 6source data 1). The remaining four had mutations in previously described rare pathways (PFLU0085, PFLU0183), again confirming that these non-focal pathways produce just a fraction of the total set of WS mutants (Lind et al., 2015).
The distribution of mutations for each of the three pathways is indicative of bias. As shown in Figure 6B, almost 29% of all WS-causing mutations (adjusted for differences in mutation rates between the three pathways) were due to an identical 33 base pair in-frame deletion in awsX (Dt229-g261, DY77-Q87), while a further 13% were due to an identical mutation (79 a->c, T27P) in awsR. At least 41 different mutations in Aws lead to WS: if mutation rates were equal for these sites the probability of observing 20 identical mutations would be extremely small. In fact 10 million random samplings from the observed distribution of mutations failed to recover this bias. While the Wsp pathway also contains sites that were mutated more than once (six positions were mutated twice, one site three times and one five times), sources of mutational bias in Wsp were less evident than in Aws ( Figure 6B).

Incorporating mutational bias into null models
The mathematical models presented above assumed no mutational bias, but the null models can be modified to incorporate such bias. With focus on the AwsX hotspot, we show this for models II-IV.

Null Model II
The probability of a mutation at the hotspot is kp where k > 1 and p is the probability of a mutation in non-hotspot locations. The probability that a pathway is used is: kp + (1-kp)(1-(1-p) (n-1) ) where kp is the probability of a mutation at the hotspot and (1-kp)(1-(1-p )(n-1) ) is the probability there is no mutation at the hotspot, but that there is a mutation elsewhere in the pathway. Using the binomial approximation, the probability a given pathway is used is: kp + (1-kp)(n-1)p, which for Aws is kp+(1kp)2p. Figure 7A shows that if the mutation rate at the hotspot increases 10-fold (k = 10.3) then Aws is used 1.8 times more often than Wsp. While this matches experimental data, the same 10-fold increase applied to Mws results in mismatch. Only with a reduced mutation rate at the same locus (k = 6.8) do predicted and observed data agree ( Figure 7B).

Null Model III
The mutational hotspot in the Aws pathway occurs in AwsX, which is involved in both an enabling and a disabling reaction. If the hotspot increases the rate of a disabling reaction by a factor k where k > 1, then the probability the pathway is used is: which is the sum of probabilities of the events the hotspot is used, the hotspot is not used but a disabling change occurs, and the hotspot is not used but only an enabling change occurs. Figure 7C shows that if k is between 4 and 5.3, then there is a good fit to experimental data. The actual value depends on the probabilities of enabling/disabling change and which reaction has the hotspot. For Mws the mutational hotspot must have a factor between 10 and 15 ( Figure 7D).

Null Model IV
The mutational hotspot in awsX affects reaction rates r 2 and r 3 in the Aws system of differential equations ( Figure 4B, Figure 4-figure supplement 2). If the mutation factor is between 3.3-5.6 for the Wsp/Aws comparison ( Figure 7E) or 4.8-7.6 for the Mws/Aws comparison ( Figure 7F), then predictions match observed data. Thus Model IV is the only one where the same five-fold change in the probability of enabling/disabling change predicts observed data for all three mutational pathways ( Figure 7G). This model, updated with hotspot data, predicts that the Aws pathway is more likely to generate WS types than Wsp for all probability values (see Figure 7H). If disabling changes are more probable than enabling changes, then r 3 dominates the predicted rate, whereas when enabling changes are more common r 2 is the major contributor ( Figure 7-figure supplement 1). Knowledge of the functional interaction between AwsX and AwsR allows interpretation of the likely effect of the hotspot mutations in AwsX and AwsR as disabling change in r 3 ( Table 1).

Analysis of mutants reveals mutational targets and effects
Null Model IV of the Wsp, Aws, and Mws pathways allow estimation of the relative probability that a WS is produced by enabling or disabling changes in a certain reaction rate. This means that in order to evaluate the success of our predictions, the 105 WS mutations must be connected to their likely effects on reaction rates. For some mutations this is straightforward. If a mutation completely disables the negative regulator WspF, this will reduce r 2 and r 6 , thereby producing a wrinkly spreader. In other cases interpretations of likely effects are more difficult and requires knowledge of the molecular functions of the proteins and their interactions, which needs to be obtained from literature, bioinformatics and structural modelling. For example, mutations in WspA can potentially produce WS by disruption of the interaction with WspF (reducing r 2 ), increase methylation by WspC (increasing r 1 ), increase signalling rate (r 3 ) or increase phosphorylation of WspE (r 4 ). A summary of the interpretation of the molecular effects of mutations and how they are connected to reaction rates is available in Table 1 with a more detailed description in Appendix 1. Mutations were found in five genes in the Wsp pathway. The majority of the mutations were found in the negative regulator WspF or its interacting partners WspE and WspA. These mutations likely reduced reaction rates of r 2 and r 6 in accord with predictions of null Model IV. The remaining mutations in WspC and WspR were less common, in accord with predictions and likely due to increases in reaction rates r 1 and r 5 , respectively. Results for the Aws pathway were also in line with predictions with the majority of the mutations arising in the negative regulator AwsX or the interacting part of AwsR (both decreasing r 3 ), with rare mutations in AwsO and other parts of AwsR. For MwsR, loss-of-function mutations were expected primarily in the phosphodiesterase (PDE) domain, but mutations were found in both diguanylate cyclase (DGC) and PDE domains. A structural analysis revealed that mutations were clustered in the interface between the domains and unlikely to disrupt PDE function. These mutations most likely change the dynamics between the DGC and EAL domains in a way that increases production of c-di-GMP.

Differences of mutational spectra with and without selection
The null models -and especially null Model IV -successfully predicted the mutational targets and explained mutation rates to WS when knowledge of mutational hot spots was included. To what and r 3 in the Aws differential equation system. Middle pie chart use p d = 0.001, p e = 0.0001 and right pie chart p d = 0.0001, p e = 0.001 to allow comparison to Figure 3F and G. (H) Null Model IV predicted probability of Wsp relative to Aws with a hotspot in AwsX that increases the mutation rate five-fold for p e and p d for r 2 and r 3 in the Aws system of differential equation (I) Null Model IV predicted probability of Aws relative to Mws with a hotspot in AwsX that increase the mutation rate five-fold for p e and p d for r 2 and r 3 in the Aws system of differential equation. Figure 7-figure supplement 1 shows the relative contribution of each reaction rate in the Aws network to the production of WS for the revised null model IV with a 5 times increase for p e and p d for r 2 and r 3 . DOI: https://doi.org/10.7554/eLife.38822.016 The following figure supplement is available for figure 7: Figure supplement 1. Relative contributions of reaction rates for Aws for null Model IV with a five-fold increase in mutation factor for r 2 and r 3 . DOI: https://doi.org/10.7554/eLife.38822.017 degree can such knowledge allow prediction of the outcome of the original experimental evolution under selection (McDonald et al., 2009)? A comparison between the frequencies of mutations isolated here (without selection) with those isolated under selection, reveals several notable differences (Figure 8). The most obvious difference is in use of the Wsp pathway, which is most commonly used (15/24) under selection and yet produces WS types at a lower rate than the Aws pathway. Differences are also apparent in the spectrum of wsp mutations, with no wspA mutations being found under selection despite being the most commonly mutated gene without selection (15/46), and the previous failure to detect wspR mutants in a screen of 53 WS mutants (Goymer et al., 2006).

Fitness of WS types
The most obvious explanation for the differences in mutational spectra between WS isolated with or without selection (Figure 8) is that certain mutants have a lower fitness and thus their relative frequencies will be lower in the original evolution experiment (McDonald et al., 2009). We measured the fitness of representative WS types with mutations in each of the mutated genes (wspA, wspC/D, wspE, wspF, wspR, awsX, awsR, awsO, mwsR) in 1:1 competitions against a reference WspF DT226-G275 deletion mutant marked with GFP ( Figure 9). This type of fitness data should be interpreted with caution because the fitness of WS mutants are frequency-dependent (Rainey and Travisano, 1998) and some WS mutants are superior in early phase attachment as opposed to growth at the air-liquid interface (Lind et al., 2015). Nevertheless, the competition experiments provide an estimate of fitness when several different WS mutants compete at the air-liquid interface (a likely situation given a~10 À8 mutation rate to WS and a final population size of >10 10 ). The fitness data account for the over-or under-representation of some WS mutants when grown under selection (McDonald et al., 2009) compared to those uncovered without selection (as reported here). The three wspF mutants, the wspC-wspD fusion, and the wspE mutants have similar fitness (p>0.38, except for WspE K734E that has slightly higher fitness (p=0.027), two-tailed t-tests). In contrast, both wspA mutants are slightly less fit (p<0.0214, two-tailed t-test) and both wspR mutants are severely impaired (p<0.00007, two-tailed t-test) (Figure 9). This sits in accord with previous work in which mutations generating WS obtained with selection have been detected in wspF and wspE, but not wspA or wspR (Goymer et al., 2006;McDonald et al., 2009). All awsXRO mutants have similar lower fitness (p<10 À6 , two-tailed t-test) compared to the wspF reference strain (Figure 9), which explain why under selection these are found at lower frequencies compared to mutations in the wsp pathway (McDonald et al., 2009) despite a roughly two-fold higher mutation rate to WS.

Discussion
The issue of evolutionary predictability and the relative importance of stochastic events compared to deterministic processes have a long history in evolutionary biology (Darwin, 1872;Simpson, 1949;Jacob, 1977;Gould, 1989;Conway Morris, 2003;Orgogozo, 2015). Recent interest has been sparked by an increasing number of observations that both phenotypic and molecular evolution, under certain circumstances, can be remarkably repeatable (Colosimo et al., 2005;Shindo et al., 2005;Jost et al., 2008;Barrick et al., 2009;Lee and Marx, 2012;Meyer et al., 2012;Zhen et al., 2012;Herron and Doebeli, 2013), but whether these cases are representative for evolutionary processes in general remains to be determined. A related question, with greater potential for practical applications, is whether it is possible to forecast short-term evolutionary events and if so, then the challenge is to stipulate the data necessary to make successful predictions. Our uniquely detailed knowledge of the WS experimental evolution system has provided a rare opportunity to determine the contributions of mutational bias and genetic architecture to the generation of new adaptive phenotypes, and consequently explore the limits of evolutionary forecasting. A thorough understanding of the function of the molecular species and their interactions allowed development of null models (see especially Model IV) that capture essential features of the genotype-to-phenotype map sufficient to predict the relative likelihood that evolution will follow each of the three principle pathways, along with specific mutational (genetic) targets. Even the most sophisticated of these models failed to forecast outcomes that matched the experimental data shown in Figure 2, however, the reason became apparent upon characterisation of the set of WS mutants obtained without selection: the presence of mutational hotspots. Armed with knowledge of sources of mutational bias it was a simple matter to refine Model IV leading to predictions matching those observed by experiment.
Problematic at the current time is inability to a priori detect all causes of mutational bias, however, it is likely that this will improve as understanding of the biochemical causes of bias improves and algorithms trained to recognise and detect nucleotide patterns indicative of bias are implemented. One specific deletion (DY77-Q87) in awsX was found to account for nearly half (20/41) of the mutations in the Aws pathway. Thus, despite the existence of hundreds of possible mutations In frame deletion Fusion Amino acid substitution WspC/D M1-A324 fused to D18-S232 Figure 9. Fitness of different WS mutants. Competitive fitness against a WspF DT226-G275 reference strain was measured for representative mutations in the Wsp, Aws, Mws pathways. Pairwise competitions were performed in quadruplicates and error bars represent ±one standard deviation. Full data, including statistical tests are available in Figure 9-source data 1. DOI: https://doi.org/10.7554/eLife.38822.020 The following source data is available for figure 9: Source data 1. Data from fitness assay. DOI: https://doi.org/10.7554/eLife.38822.021 leading to WS (this work and (McDonald et al., 2009;McDonald et al., 2011;Lind et al., 2015)) one single mutation accounts for more than one quarter of all WS mutations. While the six base pair direct repeat flanking the deletion provides a convincing explanation for its increased rate, it is not clear why this deletion would be ten times more common than the DP34-A46 deletion in the same gene that is flanked by ten base pair repeats and contains five base pairs identical to those from the DY77-Q87 deletion ( Figure 6-source data 1). There are also instances where single base pair substitutions are overrepresented: the AwsR T27P mutation is found in nine cases, while eight other single pair substitutions in Aws were found only once. Consider further the fact that wspE (a gene of~2.3 kb), where changes to only four specific amino acids repeatedly cause WS, and wspF (a gene of~1 kb) where any mutation that disrupts function results in WS ( Figure 6A) contribute equally to the rates of at which new WS types arise. It is evident from these findings and from related studies (Pollock and Larkin, 2004) that there is need for detailed experimental measurement of local mutation rates in specific systems. Such investigations stand to contribute to understanding of the causes of mutational bias and the extent to which biases might be conserved among related or even unrelated organisms. If local nucleotide sequence is the major determinant, an estimate of mutation rate will apply strictly to very closely related species, but if the dynamics of molecular processes, such as transcription and replication (Sankar et al., 2016), are major influences then estimates might be applicable to a wider range of species.
Forecasting mutational routes to new phenotypes is one component of a comprehensive forecasting strategy. The second requirement is ability to a priori predict fitness effects. Solving this problem is the Holy Grail for predicting evolution, but at the current time this is not possible. That it matters is made clear by the observation that a subset of all possible WS mutants was found after experimental evolution (Figure 8) due to inability of some to successfully compete with high fitness WS types (Figure 9).
Direct measurement of the fitness effects of large numbers of mutations is difficult, time-consuming, typically only possible for microbial species and against the spirit of identifying a priori predictors of fitness effects. Future success likely rests on ability to infer fitness from parameters such as estimated effects of mutations on thermodynamic stability (Capriotti et al., 2005;Dean and Thornton, 2007;Rodrigues et al., 2016), molecular networks, evolutionary conservation of amino acid residues (Glaser et al., 2003;Ng and Henikoff, 2003;Landau et al., 2005;Kumar et al., 2009;Celniker et al., 2013;Choi and Chan, 2015) or machine-learning methods combining several types of information (Bromberg and Rost, 2007;Li et al., 2009;Capriotti et al., 2013;Yates et al., 2014;Hecht et al., 2015). Possibilities likely also exist to extrapolate findings from a small number of mutations that are either directly constructed and assayed in the laboratory or through fitness estimates of polymorphism data from natural populations. Recent work on the prediction of the fitness effects of random mutations in several genes suggests that in many cases large effect mutations can be predicted using methods based on evolutionary conservation (Lind et al., 2017a).
On a less ambitious scale it may be possible to take advantage of the fact that the distribution of fitness effects associated with mutations in single genes is often bimodal -a consequence of many mutations causing complete loss-of-function rather than intermediate deleterious; effects (Sanjuan, 2010;Hietpas et al., 2011;Jacquier et al., 2013;Firnberg et al., 2014;Lind et al., 2017a;Sarkisyan et al., 2016;Lundin et al., 2017). As such it may be sufficient to know gene function plus consequences arising from loss of function mutations and thus treat mutations within single genes as having equivalent fitness (Sanjuan, 2010;Jacquier et al., 2013;Sarkisyan et al., 2016;Lind et al., 2017a). Interestingly WS mutations occurring in the same gene typically show similar fitness effects (Figure 9). While still requiring experimental data, this simplification may fuel exploration of the relationship between large numbers of mutations and their fitness effects and thus understanding of the extent to which fitness effects are transferable between strains with different genetic backgrounds or closely related species (Ungerer et al., 2003;Pearson et al., 2012;Wang et al., 2014).
Taken at face value, our findings give reason to question the value of aspiring to forecast evolution from first principles based on mechanistic understanding. But we argue against such pessimism and point firstly to the value that stems from a clear understanding of current limitations on forecasting: defining what is known, what is not known, and what needs to be known. That it is possible to take knowledge of the genotype-to-phenotype map and forecast with good accuracy targets and rates is an important advance. There is good reason to suspect that the principles outlined here and previously (Lind et al., 2015) are transferable to other systems and even generalizable (Lind, 2018).
Evolutionary forecasting is in its infancy. In the short term it is likely to be most successful for biological systems where there are experimental data on a large number of independent evolutionary events, such as is the case for influenza, HIV and cancer (Kouyos et al., 2012;Fraser et al., 2014;Lawrence et al., 2014;Luksza and Lässig, 2014;Neher et al., 2014;Eirew et al., 2015). Evolution might appear idiosyncratic indicating that every specific system requires detailed investigation, but our hope is that deeper knowledge of the genotype-to-phenotype map, distribution of fitness effects and mutational biases will allow short term forecasts to be produced using modelling without the need for large-scale experimental studies. A major boost to further refinement of evolutionary forecasting is likely to come from combining coarse (top down) and fine-grained (bottom-up) approaches. Our demonstration that simple null models of functional networks can produce quantitative predictions is a step forward allowing predictions to be directly tested in other experimental systems (Lind, 2018).

Strains and media
The strains used in the study are all Pseudomonas fluorescens SBW25 (Silby et al., 2009) or derivatives thereof. The reporter construct (pMSC), used for isolation of WS mutants before selection, fused the P wss promoter to a kanamycin resistance marker (nptII) (Fukami et al., 2007;McDonald et al., 2011). P. fluorescens strains with deletions of the wsp (PFLU1219-1225), aws (PFLU5209-5211) and mws (PFLU5329) operons were previously constructed as described by McDonald et al. (McDonald et al., 2011). All experiments used King's medium B (KB) (King et al., 1954), solidified with 1.5% agar and incubation was at 28˚C. All strains were stored in glycerol saline KB at À80˚C.

Fluctuation tests and isolation of WS mutants before selection
Strains with the pMSC reporter construct and either wild type genetic background or double deletions of aws/mws, wsp/mws or wsp/aws were used to estimate mutation rates to WS before selection. Overnight cultures were diluted to approximately 10 3 cfu/ml and 60 independent 110 ml cultures were grown for 16-19 hr (OD600 = 0.9-1.0) with shaking (200 rpm) in 96-well plates before plating on KB plates with 30 mg/l kanamycin. Viable counts were estimated by serial dilution and plating on KB agar. One randomly chosen colony per independent culture with WS colony morphology was restreaked once on KB agar. The assay was repeated at least four times for the double deletion mutants and twice for the wild type strain in order to obtain enough mutants to allow estimation of mutation rates. Mutations rates and confidence intervals were estimated using the Ma-Sandri-Sarkar Maximum Likelihood Estimator ) available at www.keshavsingh.org/ protocols/FALCOR.html. The mutation rates between the different strains were statistically evaluated using a t-test as previously described (Rosche and Foster, 2000), but this method has only been shown to be valid in cases where total population size is not significantly different for the strains used. In our assay this was not the case, as determined by ANOVA, and therefore the results of the statistical analysis should be interpreted with caution. As the estimated number of mutants per well was <0.5 for all strains, the biasing effect of differences in fitness between WS mutants is minimal.

Sequencing
Mutations causing the WS phenotype were identified by Sanger sequencing of candidate genes in the remaining common pathway to WS, for example the wsp operon for the aws/mws deletion strain. In a few cases where no mutations were identified in the previously established WS target genes, we used genome sequencing (Illumina HiSeq, performed by Macrogen Korea).

Fitness assays
Competition assays were performed as previously described (Lind et al., 2015) by mixing the WS mutant 1:1 with a reference strain labelled with green fluorescent protein and measuring the ratio of each strain before and after static growth for 24 hr using flow cytometry (BD FACS Canto). We used a WspF DT226-G275 deletion mutant as the reference strain because WspF mutants are the most commonly found WS type when grown under selective conditions (McDonald et al., 2009) and the in frame deletion of 50 amino acids most likely represents a complete loss-of-function mutation with minimal polar effects on the downstream wspR. Selection coefficients per generation were calculated as s = [ln(R(t)/R(0))]/[t], as previously described (Dykhuizen, 1990) where R is the ratio of alternative WS mutant to WspF DT226-G275 GFP and t the number of generations. Viable counts on KB plates of initial and final populations were performed to calculate the number of generations. Stability of colony morphologies was confirmed and data from microcosms with >5% smooth colonies were excluded (two cases). Control competition experiments with isogenic WspF DT226-G275 reference strains with and without GFP were used to correct for the cost of the GFP marker. Control competitions were also used to determine the cost of the double deletions and the reporter construct relative to a wild type genetic background, for example an AwsX DY77-Q87 mutant in Dwsp/Dmws background with pMSC was competed with a GFP labeled AwsX DY77-Q87 mutant in wild type background. Competitions were performed in independently inoculated quadruplicates for each strain with the number of replicates based on previous work (Lind et al., 2015).

Probability estimation in the null model IV
The differential equation models describe the interactions between proteins in each of the three WS pathways. In order to solve the differential equations, two pieces of information are required: (i) the initial concentrations of the molecular species and (ii) the reaction rates. Although this information is unavailable a random-sampling approach was used to generate different random sets of initial concentrations and reaction rates. Each random set was used to establish a baseline of potential WS expression making it possible to evaluate whether a set of mutations results in a WS type. Effectively, this approach allows sampling of the probability distribution P (WS |m i 2 Wsp) used in our Bayesian model. We randomly sample 1000 different sets of reaction rates and initial concentrations from uniform priors: reaction rates were sampled randomly from a uniform distribution on log space (i.e. 10 U[À2,2] ) and initial concentrations of reactants were sampled from a uniform distribution U[0,10]. For each set, the appropriate differential equation model was integrated and the steady state concentration of the compounds that correspond to a wrinkly spreader (RR in Aws, R* in Wsp and D* for Mws) computed. This served as a baseline for the non-WS phenotype that was used for comparison to determine whether combinations of mutations result in increased WS expression. After obtaining the baseline, we implemented particular combinations of enabling/disabling mutations (a m i ). Ideally, a distribution linking enabling/disabling mutations to a fold change in reaction rates would be used, but this information is unavailable. In order to progress the effect sizes for enabling and disabling mutations were sampled from 10 U[0,2] and 10 U[À2,0] , respectively, and then multiplied by the reaction rates. The differential equations were then solved for the same time that it took the baseline to reach steady state. The final concentration of R* ( Figure 4A), RR ( Figure 4B) and D* ( Figure 4C) was then compared to the baseline and the number of times out of 1000 that the WS-inducing compound increased served as an estimate of P (WS|m i 2 Wsp). The probability distribution stabilized by 500 random samples and additional sampling did not produce significant changes (data not shown).
The absence of empirical data on reaction rates, initial concentrations, and expected mutation effect size meant using a random sampling approach requiring estimates for parameter ranges. Parameter ranges were chosen to be broad enough to capture differences spanning several orders of magnitudes while allowing numerical computations for solving the differential equations. To assess the effect of these ranges on the results, the sampling procedure was repeated for WSP for three different parameter regimes (i) an expanded range for initial concentrations [0-50], (ii) an expanded range for reaction rates 10 [-3,3] , (iii) a compressed range for mutational effect size 10^ [ -1,1] . This analysis shows that qualitative results are robust to these changes (see Figure 5-figure supplement 3).
hotspot produced in-frame deletions likely mediated by 6 bp direct repeats ( Figure 6-source data 1). The deletions are consistent with a loss of function and a decrease in r 3 ( Figure 4B) that would leave the partially overlapping open reading frame of the downstream gene (awsR) unaffected.
The DGC AwsR, was mutated in 14 cases with an apparent mutational hot spot at T27P (nine mutants) in a predicted transmembrane helix (amino acids 19-41). The remaining mutations were amino acid substitutions in the HAMP linker and in the PAS-like periplasmic domain between the two transmembrane helices. These amino acid substitutions are distant to the output DGC domain ( Figure 6) and their effects are difficult to interpret, but they could cause changes in dimerization (Malone et al., 2012) or the packing of HAMP domains, which could, in turn, alter transmission of conformational changes in the periplasmic PAS-like domain to the DGC domain causing constitutive activation (Parkinson, 2010). Such effects would increase r 4 in Figure 4B. Mutations in the N-terminal part of the protein are easier to interpret based on the existing functional model (Malone et al., 2012) and most likely disrupt interactions with the periplasmic negative regulator AwsX resulting in a decrease in r 3 in Figure 4B.
Two mutations were found in the outer membrane lipoprotein protein AwsO between the signal peptide and the OmpA domain ( Figure 6). Both mutations were glutamine to proline substitutions (Q34P, Q40P), which together with a previously reported G42V mutation (McDonald et al., 2009) suggest that multiple changes in this small region can cause a WS phenotype. This is also supported by data from Pseudomonas aeruginosa in which mutations in nine different positions in this region lead to a small colony variant phenotype similar to WS (Malone et al., 2012). A functional model based on the YfiBNR in P. aeruginosa (Malone et al., 2012;Xu et al., 2016), suggest that AwsO sequesters AwsX at the outer membrane and that mutations in the N-terminal part of the protein lead to constitutive activation and increased binding of AwsX. This would correspond to an increase in r 2 in Figure 4B, which would relieve negative regulation of AwsR.
Mws pathway: The MwsR pathway (comprising just a single gene) harbours mutations in both DGC and phosphodiesterase (PDE) domains. Mainly mutations in the C-terminal phosphodiesterase (PDE) domain were predicted ( Figure 4C). Eleven of 18 mutations were identical in-frame deletions (DR1024-E1026) in the PDE domain, mediated by 8 bp direct repeats ( Figure 6, Figure 6-source data 1). It has been shown previously that deletion of the entire PDE domain generates the WS phenotype (McDonald et al., 2009), suggesting a negative regulatory role that causes a decrease of r 2 in the model in Figure 4C. One additional mutation was found in the PDE domain (E1083K) located close to R1024 in a structural homology model made with Phyre2 (Kelley et al., 2015), but distant to the active site residues (E1059-L1061). Previously reported mutations (A1018T, ins1089DV) (McDonald et al., 2009) are also distant from the active site and cluster in the same region in a structural homology model. This suggests that loss of phosphodiesterase activity may not be the mechanism leading to the WS phenotype. This is also supported by the high solvent accessibility of the mutated residues, which indicates that major stability-disrupting mutations are unlikely and changes in interactions between domains or dimerization are more probable. Thus, it is likely that the WS phenotype resulting from a deletion in the PDE domain is caused by disruption of domain interactions or dimerization rather than loss of phosphodiesterase activity.
The remaining mutations within mwsR are amino acid substitutions in the GGDEF domain, close to the DGC active site (G927-F931) with the exception of a duplication of I978-G985. While it is possible that these mutations directly increase the catalytic activity of the DGC, increasing r 1 in Figure 4C, such enabling mutations are considered to be rare. An alternative hypothesis is that these mutations either interfere with c-di-GMP feedback regulation or produce larger conformational changes that change inter-domain or inter-dimer interactions, similar to the mutations in the PDE domain. Based on these data we reject the current model of Mws function, which predicted mutations decreasing r 2 ( Figure 4C) through mutations inactivating the PDE domain. We instead suggest that the mutations are likely to disrupt the conformational dynamics between the domains and could be seen either as activating mutations causing constitutive activation or disabling mutations with much reduced mutational target size that must specifically disrupt the interaction surface between the domains. In both cases the previous model leads to an overestimation of the rate to WS for the Mws pathway.