The mechanisms of gene regulatory networks constrain evolution: A lesson from synthetic circuits

Phenotypic variation is the raw material of adaptive Darwinian evolution. The phenotypic variation found in organismal development is biased towards certain phenotypes, but the molecular mechanisms behind such restrictions are still poorly understood. Gene regulatory networks have been proposed as one cause of constrained phenotypic variation. However, most of the evidence for this is theoretical rather than experimental. Here, we study evolutionary biases in two synthetic gene regulatory circuits expressed in E. coli that produce a gene expression stripe - a pivotal pattern in embryonic development. The two parental circuits produce the same phenotype, but create it through different regulatory mechanisms. We show that mutations cause distinct novel phenotypes in the two networks and use a combination of experimental measurements, mathematical modelling and DNA sequencing to understand why mutations bring forth only some but not other novel gene expression phenotypes. Our results reveal that the regulatory mechanisms of networks restrict the possible phenotypic variation upon mutation. Consequently, seemingly equivalent networks can indeed be distinct in how they constrain the outcome of further evolution.


37
The ability of biological systems to bring forth novel and beneficial phenotypes as a 38 consequence of genetic mutations is essential for evolutionary adaptation and innovation. This 39 ability is encapsulated in the concept of evolvability (Kirschner and Gerhart, 1998;Wagner, 40 2005b). Evolvability can be limited by evolutionary constraints, which are biases or limitations 41 in the production of novel phenotypes (Smith et al., 1985). An example of such constraints 42 comes from laboratory selection experiments with butterfly populations for enhanced wing 43 eyespot colours (Allen et al., 2008). Selection was able to increase the amount of black or gold 44 colouring in the two eyespots simultaneously, but was unable to do so for the two different 45 colours independently in the two eyespots. Constrained variation can have multiple genetic and 46 developmental causes that can be difficult to disentangle in a complex developing organism 47 (Arnold, 1992;Wagner, 2011). Therefore, few experimental demonstrations of evolutionary 48 constraints exist. What is more, thirty years after this concept rose to prominence (Smith et al.,49 1985), we still do not understand the mechanistic causes of evolutionary constraints. The instructions for an organism's development are encoded in gene regulatory networks 52 (GRNs)networks of interacting transcription factors that control gene expression in both time 53 and space (Davidson, 2006). Mutations in the cis-regulatory regions of GRNs play an important 54 part in evolutionary adaptation and innovation (Payne and Wagner, 2014; Prud'homme et al., 55 Briscoe, 2017; Stanojevic et al., 1991;Wolpert, 1996). The question of which regulatory 89 mechanisms can produce stripes is therefore itself crucial for developmental genetics (Cotterell 90 and Sharpe, 2010; Francois et al., 2007). Here we go beyond this question to ask whether 91 different GRNs that have the same phenotype (a "stripe" of gene expression) can produce 92 different novel (i.e., "non-stripe") gene expression phenotypes in response to mutations, and if 93 so, why. Specifically, we use here two synthetic circuits that employ different regulatory mechanisms to 96 produce a striped gene expression pattern. Both of these circuits are hosted by Escherichia coli 97 bacteria. When these bacteria are grown as a lawn in the presence of a concentration gradient 98 of the morphogen analogue, they display a spatially striped gene expression pattern (Fig. 1 c). 99 We introduced random mutations into the regulatory regions of these circuits, and analysed the 100 resulting phenotypes. The two circuits indeed produce a different spectrum of novel gene 101 expression phenotypes. That is, the gene expression variation they produce is constrained. To 102 identify the mechanistic causes of these constraints, we combined experimental DNA sequence 103 and phenotypic data with a mathematical model of gene expression dynamics. networks (see Box 1). The expression level of the "green" gene is the phenotypic "output" of 144 the network. Circles: Bacterial lawns display green fluorescent rings as a function of arabinose 145 gradients from central paper disks (white). Figure adapted from (Schaerli et al., 2014).
Box 1. Two starting circuits producing stripes through two different mechanisms 149 Opposing gradients mechanism (Incoherent FFM type 2): The "red" gene (with the open 150 reading frames (ORFs) for LacI and TetR encoded on the same transcript) is activated by the 151 "morphogen" arabinose (vertical arrow). Its products thus form a gradient of increasing 152 concentration with increasing arabinose concentration. The "blue" gene (LacI) and the "green" 153 gene (GFP) are expressed from constitutive promoters. However, the "blue" gene is also 154 repressed by the "red" gene product (TetR). Thus, the "blue" gene product forms an opposing 155 gradient with respect to the gradient of the "red" gene product. Both the "blue" (LacI) and "red" 156 (LacI) gene products repress the "green" gene. The GFP thus reaches a high expression only at 157 medium morphogen concentration where the repression from the "red" and "blue" genes are 158 low.

159
Concurring gradients mechanism (Incoherent FFM type 3): The "red" gene (with the ORFs for 160 SP6 RNA polymerase (RNAP) and LacI encoded on the same transcript) is activated by the 161 "morphogen" arabinose, just as in the previous circuit. Its expression thus also mimics the 162 arabinose gradient. However, in this circuit the "red" gene product SP6 RNAP activates the 163 "blue" gene, which thus forms a concurring gradient with respect to the gradient of the "red" 164 gene product. The "green" gene is activated by the "blue" gene (T7 RNAP) and repressed by 165 LacI of the "red" gene. Its maximum expression occurs at medium arabinose concentration 166 where there is already activation from the "blue" gene, but not yet a high level of repression of 167 the "red" gene. We introduced mutations into the regulatory regions of these two networks by replacing the 170 wild-type regulatory sequence with semi-randomised weighted oligonucleotides (Isalan, 2006).

171
Resulting average mutation rates per regulatory regions ranged from 2.6 to 3.5 mutations 172 (mainly point mutations and <5% of insertions and deletions) per regulatory region with 173 individual mutants carrying 1 to 9 mutations (SI sequences, SI Table 4). For each of our two 174 networks, we first generated three libraries of mutant networks in which mutations were 175 restricted to regulatory regions of the "red", "blue" or "green" gene ( Fig. 1). After plating cells 176 from a population whose members harboured a synthetic network variant, we randomly picked 177 colonies, grew them in liquid culture, and measured their GFP expression at low (0%), middle 178 (0.0002%) and high (0.2%) arabinose concentrations (SI expression). We classified the 179 observed fluorescence phenotypes into six categories ( Fig. 2a) (see Methods for exact 180 definitions): "stripe", "increase", "decrease", "flat" and "broken" (all expression values below 181 a threshold) and "other" (phenotypes that do not fall in any of the previous categories).  Fig. 2b summarises the spectrum of phenotypes we observed after mutagenesis. We first note 184 that both networks are to some extent robust to mutations; that is, a considerable fraction of 185 mutations do not change the "stripe" phenotype (black sectors in Fig. 2b). What is more, the 186 two types of networks we study differ in their robustness. Averaged across the three genes, 45.5 187 % of analysed mutants preserve the stripe phenotype in the concurring gradients network, 188 whereas only 32.9% do so in the opposing gradients network. The concurring gradients network 189 is thus significantly more robust to mutations [Chi-square goodness of fit test, X2 (1, N = 215) 190 = 13.67, p =0.0002]. Next, we note that within any one of the two networks the novel 191 phenotypes do not occur at the same frequency, providing evidence for the biased production 192 of novel phenotypes, where certain types of phenotypes are more common than others.  at medium arabinose concentration is compared to the GFP expression levels at low (x axis) 203 and high arabinose (y axis) concentrations. The numbers written close to each phenotype group 204 are the average mutation rates for that group. We omitted the "broken" phenotype from this 205 analysis, as the networks with this phenotype do not show any significant GFP expression. d 206 Experimentally observed phenotype distributions as displayed in b and c, grouped according to 207 the mutated gene. We also observed differences in the types of novel phenotypes between the two networks. For  Fig. 2b). In contrast, mutations in the concurring gradients network did not produce a single 214 such phenotype. In addition, mutations in the opposing gradients network are more likely to 215 create a "decrease" phenotype (purple, 29.8% of all novel phenotypes) rather than an "increase" 216 phenotype (orange, 15.4%). For the concurring gradients network, the opposite is true:

219
Next, we analysed the GFP expression levels of the measured phenotypes quantitatively (Fig.   220 2c). To this end we compared the GFP expression at medium arabinose concentration to those 221 at high (y axis) and at low arabinose concentrations (x axis). We note that the previously 222 classified phenotypes (Fig. 2a) form well-separated clusters in this analysis. For example, 223 networks in the bottom-right quadrant correspond to "stripe" phenotypes, because their pattern 224 is described as an increase (positive x-axis) followed by a decrease (negative y-axis) in 225 expression. Consequently, "decrease" and "increase" phenotypes occupy the upper-right and 226 bottom-left quadrants, respectively. We also sequenced the mutated regulatory regions of all 227 analysed networks, and find a weak association between the number of mutations a network 228 carries, and the extent to which its observed phenotype differs from the starting stripe phenotype 229 (as quantified through the Euclidean distance) (SI Figure 2).

231
Subsequently, we analysed the differences in novel phenotypes created by mutations in specific 232 regulatory regions (i.e. of the "red", "blue" or "green" gene). Within any one of the two network 233 types, regulatory mutations in the "red" gene most often create "increase" phenotypes ( Fig. 2d, 234 pie charts left to the "red" genes), whereas those in the "blue" gene most often create "decrease" 235 phenotypes (Fig. 2d, pie charts at the bottom of the "blue" genes), and those in the "green" gene 236 preferably create "broken" phenotypes ( Fig. 2d, pie charts to the right of the "green" genes). As 237 a consequence, not all phenotypes can be reached by introducing mutations in the regulatory 238 region of any of the three genes. For example, in the opposing gradient network, the "increase" 239 phenotype is only reachable by introducing mutations into the "red" gene, but not in the "blue" 240 and "green" genes.

241
The two networks differ in the spectrum of novel phenotypes that mutations in individual genes 242 create, which is especially obvious for mutations in the "green" gene: Unless regulatory 243 mutations in this gene lead to a complete loss of expression ("broken"), the opposing gradients 244 network is >5 times more likely to create a "flat" phenotype (23.2 %) than a "decrease" 245 phenotype (4.1 %). In contrast, the concurrent gradients network does not produce any "flat" 246 phenotype at all, but readily produces "increase" phenotypes (4.5 %). In sum, mutations in 247 networks which start with the same phenotype (single stripe formation), but which have 248 alternative topologies and regulatory mechanisms, create different kinds of novel phenotypes.

249
Hence, phenotypic variation is subject to constraints, and these constraints differ between 250 regulatory regions and networks. Differences in constrained variation can be explained by differences in the regulatory 253 mechanisms behind stripe formation 254 We next asked whether the regulatory mechanisms contributing to stripe formation can help 255 explain these phenotypic constraints. In doing so, we focused on novel phenotypes produced 256 by regulatory mutations in the "green" gene, because such mutations produced the most distinct 257 spectrum of novel phenotypes (Fig. 2d). Also, the regulation of this gene is most complex,  Table 1 and SI Tables 1 -2 for details). The unmutated ("wild-type") model for each circuit 265 used parameter values determined in our previous study (Schaerli et al., 2014). Into these 266 models we now introduced quantitative changes in the parameters relating to the promoter    the "stripe" phenotype. Its area is therefore a measure for a network's robustness to parameter 289 changes. The value of those diagrams is that they provide information on which parameters 290 must be mutated, and by how much, in order to access a given phenotype.

292
For each of the novel phenotypes observed experimentally when mutating the regulatory region 293 of the "green" gene ( Fig. 2b), some "mutated" model parameter values exist which reproduce 294 the phenotype (Fig. 3a). However, the phenotype diagrams of the two networks ( Fig. 3a) are 295 visually very distinct, indicating that the two networks have different potentials to access 296 specific phenotypes. Specifically, for the opposing gradients network we find regions 297 corresponding to the "broken" (grey), "decrease" (purple), "flat" (yellow) and "other" (beige) 298 phenotypes, whereas for the concurring gradients network we find regions for the "broken" and 299 "increase" (orange) phenotypescorresponding to the phenotypes observed experimentally 300 when mutating the "green" gene ( Fig. 2d). Especially instructive are mutants with strongly 301 decreased repressor binding (i.e. reduced operator activity, arrows in Fig. 3a). Such mutants 302 produce a "flat" phenotype (yellow region) in the opposing gradient network, but an "increase" 303 phenotype (orange region) in the concurring gradient network (Fig. 3a). In other words, even 304 though both networks contain the same operator (LacO) in the "green" gene, the model predicts  categorisation reveals that mutations in the operator produce a "flat" phenotype in the opposing 339 gradient network, but an "increase" phenotype in the concurring gradient network (Fig. 4a,   340 smaller diagrams), thus validating model predictions (Fig. 3). In addition, the model predicts 341 that mutations in the operator of the opposing gradient network are able to produce a "decrease" 342 phenotype (Fig. 3a). Even though we did not find a circuit with a "decrease" phenotype that has 343 only operator mutations, all circuits with this phenotype carry at least one mutation in the 344 operator (and additional mutations in the promoter, SI sequences). So far, we demonstrated that each of the two analysed networks yields a biased spectrum of 381 novel phenotypes after mutation, and that two networks with different regulatory mechanisms 382 yield different spectra of novel phenotypes. However, these spectra may not be influenced only 383 by a network's regulatory mechanisms. They may also differ among networks with the same 384 topology and the same regulatory mechanism, but with quantitative differences in the 385 biochemical parameters that determine a networks gene expression pattern. To find out whether 386 this is the case, we performed the following experiments: We took two mutant stripe-forming 387 networks of the concurring gradient mechanism with mutations in all three genes (mutants A 388 and B) and introduced further mutations into their "green" regulatory regions. Fig. 5a shows 389 the resulting phenotype distributions and compares them to the initial ("wild-type", WT) 390 network. As in the WT network, we observe "stripe", "broken" and "increase" phenotypes in 391 the mutants. However, the figure also shows that the proportions of these phenotypes differ 392 among the networks. In addition, 3% (mutant A) and 1.3% (mutant B) of the two concurring 393 gradient network variants displayed a "decrease" phenotype. This suggests that by making 394 "neutral" or "silent" genetic changes in a regulatory network that do not affect its ("stripe") of their wild-type parameter value), because a mutation is likely to affect these specific 420 parameters in a similar way (SI Table 3). The changed parameters were drawn from a uniform 421 distribution, and for each parameter we aimed to identify upper and lower bounds for this 422 distribution that give the best possible agreement between the experimental and model data (SI 423  also has an upper bound higher than 100% of the unmutated (WT) value (SI Table 3). This was 440 unexpected and led us to discover a context dependent effect in the plasmid we used to express

Opposing gradients
Concurring gradients red blue green red blue green 23 regulatory regions (SI sequences). Fig. 6b shows the resulting distribution of phenotypes.

475
Similar to the one-gene mutants (Fig. 2b, repeated in Fig. 6a), some phenotypes occur more 476 frequently than others, and the opposing and concurring gradient networks produce any one 477 phenotype at different frequencies. Regardless of the regulatory mechanism, the frequencies of novel phenotypes differed 490 significantly between networks with mutations in multiple versus single genes ( Fig. 6a and 6b; to multiple mutations, but the concurring gradient networks did not produce this phenotype in 494 response to single-gene mutations (compare yellow sectors in Fig. 6a and 6b).  In SI Fig. 6 we show an experimental example of how mutations in the "green" and "blue" 502 genes can interact to produce a "flat" phenotype in the opposing gradient network: The network 503 with the mutated "green" gene maintains the "stripe" phenotype (SI Fig. 6a) and the network 504 with the mutated "blue" gene leads to a "decrease" phenotype (SI Fig. 6b). When these two 505 mutations are combined, the resulting phenotype is "flat" (SI Fig. 5c). Importantly, this new 506 phenotype cannot just be explained as an additive superposition of the two individual 507 phenotypes.

508
To understand the phenotype distributions of the multiple-gene mutants and in particular, the 509 non-additive interactions, we turned again to our model. Analogous to the experiments, we now 510 changed parameters of multiple genes simultaneously, within the exact same ranges as used to 511 model single-gene mutants (SI table 3). The resulting phenotype distributions (Fig. 6c, SI Table   512 5) predict the experimentally observed distributions (Fig. 6b) (Fig. 1) and analysed the resulting distributions of novel gene expression phenotypes (Fig. 2). to understand the differences between accessible novel phenotypes for the two networks ( Fig.   532 2 and 3). The model predictions are also supported by DNA sequencing data (Fig. 4). We thus 533 provide for the first time empirical evidence that GRNs with different regulatory mechanisms 534 can cause different constrained variation, as was recently proposed (Jimenez et al., 2015). We 535 also provide experimental evidence that the mechanism by which a network produces a stripe 536 constrains the origin of novel expression phenotypes more than quantitative parameters driving 537 gene expression dynamics (Fig. 5).

538
Comparisons of GRNs in related species indicate that they indeed solve the problem of 539 producing a specific adaptive phenotype in many different ways, and that these solutions  (Carroll, 2008). In 547 an attempt to understand the rules that govern this selection, Savageau formulated its "demand 548 rule" (Savageau, 1977). He observed that activators and repressors can achieve the same 549 regulatory goals, but that frequently expressed genes tend to be regulated by activators (positive 550 mode), whereas rarely expressed genes tend to be regulated by repressors (negative mode).

551
These differences can be explained by the fact that the negative and positive regulatory modes 552 can lead to different phenotypes upon mutation and the error-load caused by a mutation is 553 different, favouring one or the other mode of regulation (Savageau, 1977(Savageau, , 1983(Savageau, , 1998a. does not capture the complexity of a developing animal, this reduced complexity also allowed 588 us to study the potential of GRNs to bias phenotype production without confounding effects. Oligonucleotides covering the regulatory regions of these circuits were synthesised with 648 weighted base mixtures (Isalan, 2006). That is, four custom-weighted phosphoramidite 649 mixtures were prepared, with the WT base pair constituting 95% and each of the other bases 650 constituting 1.67% of any one mixture. These mixtures were used to randomise the regulatory 651 regions (Table 3) during oligonucleotide synthesis (Microsynth). These semi-randomised 652 weighted single-stranded oligonucleotides (2 µM) were annealed to a reverse primer (Table 4,   653 2.4 µM) to render them double-stranded by primer extension (2 min 95º C, cooling down to 72º 654 C in 7 min, 10 min 72º C). The resulting library of double-stranded oligonucleotides was then 655 purified with the QIAquick nucleotide removal kit (QIAGEN). Next, these double-stranded 656 oligonucleotides encoding the mutated regulatory regions were cloned into the plasmid 657 encoding the gene whose regulatory region was to be mutated. For the "blue" and "green" The cloning of the "red" regulatory regions by Gibson assembly was performed as follows: The 669 plasmid into which the region was to be inserted was first amplified using the primers 670 pBAD_Gibson_f and pBAD_Gibson_r (Table 4), and was then assembled with the double-671 stranded oligonucleotide library using the Gibson assembly master mix (NEB). Due to a 672 mistake in primer design, Gibson assembly produced a two-nucleotide deletion (CT) 673 downstream of the pBAD promoter. Therefore, the library contains sequences with and without 674 this deletion. We confirmed that this change did not affect the initial "stripe" phenotype.

675
Ligation ("blue" and "green" regulatory regions) and Gibson assembly reaction products ("red" 676 regulatory region) were transformed into electrocompetent MK01 cells (Kogenaru and Tans,677 2014) that carried already the two other plasmids necessary to complete the synthetic circuit.

678
Transformants were plated out on SM-agar plates.    "Phenotype classification, metabolic load" for details).

717
The pipetting steps for this part of the experiment were carried out with a manual pipetting  All expression data is listed in the file "SI expression". Flat: The average nF at the lowest, medium and highest arabinose concentration was calculated.

778
If all three nFs differed by less than 10% from this average, the circuit was assigned the "flat" 779 phenotype.

781
Decrease: If the following three statements were true, the mutant circuit was assigned the 782 "decrease" phenotype:

783
1) The nF at the lowest arabinose concentration was higher than 90% of the nF at the 784 medium arabinose concentration.

785
2) The nF at the medium arabinose concentration was higher than 90% of the nF at the 786 highest arabinose concentration.

787
3) The nF at the lowest arabinose concentration was higher than 120% of the nF at the 788 highest arabinose concentration.

790
Increase: If the following three statements were true, the mutant was assigned the "increase" 791 phenotype:

792
1) The nF at the highest arabinose concentration was higher than 90% of the nF at the 793 medium arabinose concentration.

794
2) The nF at the medium arabinose concentration was higher than 90% of the nF at the 795 lowest arabinose concentration.

796
3) The nF at the highest arabinose concentration was higher than 120% of the nF at the 797 lowest arabinose concentration.

799
Stripe: If the nF at medium arabinose concentration was higher than 120% of the nF at the 800 lowest and highest arabinose concentrations, the mutant was assigned the "stripe" phenotype.

802
Other: Any phenotype that did not fall into one of the previous categories. The mutagenised region(s) of all circuits whose phenotype fell into one of our six main 806 categories, and did so consistently in three independent measurements were sent for Sanger 807 sequencing (High-throughput service, Microsynth, see Table 5 for primers). Experimental confirmation of epistasis (SI Fig. 6) 828 The plasmids from mutant 4_4_f (see SI sequences) were isolated. For SI Fig. 5a the mutated 829 plasmid coding for the "green" gene was transformed together with the WT plasmids for the 830 "blue" and "red" genes into electrocompetent MK01 cells (Kogenaru and Tans, 2014). For SI 831 Fig. 5b the mutated plasmid coding for the "blue" gene was transformed together with the WT 832 plasmids for the "green" and "red" genes into electrocompetent MK01 cells (Kogenaru and 833 Tans, 2014). For SI Fig. 5c the initial 4_4_f mutant was assayed (the "red" gene is not mutated).

834
The fluorescent phenotypes were measured in a 96-well plate assay as described (   In order to fit quantitatively the distributions of novel phenotypes for the single-gene mutants 862 and predict the distributions of novel phenotypes for the multiple-gene mutants a custom-made 863 Python scripts was used. We first discuss the single-gene mutants:

864
Multiple iterations of a procedure were performed that consisted of the following three steps:

865
First, a series of simulated mutants was created. For each gene, a random binary vector whose 866 length corresponds to the number of nucleotides in the regulatory sequence (SI Table 4) was 867 generated. In this binary vector 0's and 1's indicate whether a nucleotide is not mutated (0) or 868 mutated (1), and the probability to obtain either 0 or 1 is given by the average mutation rate 869 extracted from the experimental sequencing data (SI Table 4). For every network, it was then 870 assessed which genes were mutated, according to the entries of this vector. For the single-gene 871 mutants, mutants that had only one gene mutated were selected. This process was repeated until 872 1000 single-gene mutants had been obtained. For a particular mutant, all parameters related to 873 the mutated sequence were varied (see SI Tables 1-3). A given single-gene mutant had either 874 only its promoter mutated, only its operator mutated, or both. If a mutation affected the 875 promoter (operator), all parameters determining promoter (operator) activity were changed (SI 876   Table 3). A subset of the changed parameters was varied jointly and to the same extent (i.e. all 877 of them were changed to same percentage of their wild-type parameter value), because a 878 mutation is likely to affect these parameters in a similar way (SI Table 3). New (mutant) 879 parameters were chosen according to a standard uniform distribution between upper and lower 880 ranges which were kept constant for a given model iteration.

881
Second, for each mutant the phenotype of the model (SI Tables 1 and 2) was evaluated at 3 882 arabinose concentrations (0%, 0.000195%, 0.2% for opposing gradients network and 0%, 883 0.000195%, 0.1% for concurring gradients network). In order to allocate the resulting GFP 884 expression pattern to a phenotype category the same rules as described above for the 885 experimental data ("Phenotypic categories") were used.

41
Third, the obtained phenotype distribution of the 1000 assessed mutants was compared to the 887 experimentally observed phenotype distribution.

888
After each iteration of these three steps, the upper and lower ranges of each parameter were 889 manually adjusted to best fit the results of the model to the phenotype distributions observed in 890 the experimental data (see SI "Discussion of intervals"). Finally, the best upper and lower 891 ranges for each mutant parameter distribution were kept (SI Table 3) and used to produce Fig.   892 2c.

893
To predict the distributions of novel phenotypes for the multiple-gene mutants (Fig. 2e) the 894 same procedure was used, with the following modifications: 1) only mutants containing more 895 than one mutated gene were kept for analysis.
2) The upper and lower ranges of the parameter 896 distributions were not adjusted, but the intervals derived from the single-gene mutants were 897 used (SI Table 3). 3) Only one iteration of the three steps above was performed.