Adaptive evolution of feed-forward loops versus diamonds to filter out short spurious signals

Transcriptional regulatory networks (TRNs) are enriched for certain network motifs. This could either be the result of natural selection for particular hypothesized functions of those motifs, or it could be a byproduct of mutation (e.g. of the prevalence of gene duplication) and of less specific forms of selection. We have developed a powerful new method for distinguishing between adaptive vs. non-adaptive causes, by simulating TRN evolution under different conditions. We simulate mutations to transcription factor binding sites in enough mechanistic detail to capture the high prevalence of weak-affinity binding sites, which can complicate the scoring of motifs. Our simulation of gene expression is also highly mechanistic, capturing stochasticity and delays in gene expression that distort external signals and intrinsically generate noise. We use the model to study a well-known motif, the type 1 coherent feed-forward loop (C1-FFL), which is hypothesized to filter out short spurious signals. We found that functional C1-FFLs evolve readily in TRNs under selection for this function, but not in a variety of negative controls. Interestingly, a new “diamond” motif also emerged as a short spurious signal filter. Like the C1-FFL, the diamond integrates information from a fast pathway and a slow pathway, but their speeds are based on gene expression dynamics rather than topology. When there is no external spurious signal to filter out, but only internally generated noise, only the diamond and not the C1-FFL evolves. Author Summary Frequently occurring motifs are thought to be fundamental building blocks of biological networks, conducting specific functions. However, we still lack definitive evidence that these motifs have evolved “adaptively” (to perform the particular function proposed for them), rather than “non-adaptively” (as byproducts of some other function, or as an artifact of patterns of mutations). Here we develop a powerful null model that captures important non-adaptive factors that can shape the evolution of transcriptional regulatory networks, and use it to provide the missing piece of evidence of adaptive origin in the case of the most studied motif, a feed-forward loop that is hypothesized to filter out short spurious signals. We also find evidence for an alternative solution to this problem, where the functionality of the feed-forward loop is encoded not in network topology, but in the dynamics of gene expression. Our model is suitable for studying whether other network features have evolved adaptively vs. non-adaptively.


21
Transcriptional regulatory networks (TRNs) are enriched for certain network motifs. This could 22 either be the result of natural selection for particular hypothesized functions of those motifs, or 23 it could be a byproduct of mutation (e.g. of the prevalence of gene duplication) and of less 24 specific forms of selection. We have developed a powerful new method for distinguishing 25 between adaptive vs. non-adaptive causes, by simulating TRN evolution under different 26 conditions. We simulate mutations to transcription factor binding sites in enough mechanistic 27 detail to capture the high prevalence of weak-affinity binding sites, which can complicate the 28 scoring of motifs. Our simulation of gene expression is also highly mechanistic, capturing 29 stochasticity and delays in gene expression that distort external signals and intrinsically generate 30 noise. We use the model to study a well-known motif, the type 1 coherent feed-forward loop 31 (C1-FFL), which is hypothesized to filter out short spurious signals. We found that functional C1-32 FFLs evolve readily in TRNs under selection for this function, but not in a variety of negative 33 controls. Interestingly, a new "diamond" motif also emerged as a short spurious signal filter. Like 34 Introduction internal fluctuations. Stochasticity in gene expression also shapes how external spurious signals 100 are propagated. Stochasticity is a constraint on what TRNs can achieve, but can be adaptively 101 co-opted in evolution [26]; either way, it might underlie the evolution of certain motifs. Most 102 computational models of TRN evolution that consider gene expression as the major phenotype 103 do not simulate stochasticity in gene expression (see [27][28][29] for three notable exceptions). The 104 genotype to phenotype map we develop here does include intrinsic stochasticity in gene 105 expression. 106

107
Here we use this model to ask whether AND-gated C1-FFLs evolve as a response to selection for 108 filtering out short and spurious external signals, compared to conditions that control for both 109 mutational biases and for less specific forms of selection. We find that they evolve far more 110 often under these specific selection conditions than under control conditions, providing long-111 awaited support for the adaptive hypothesis. We also ask whether there are alternative motifs 112 that evolve to solve the same selective challenge. We find that a "diamond" [30] is such a motif, 113 filtering out short spurious signals by requiring them to arrive not through both a long and a 114 short path, but through both a fast and a slow path of equal topological lengths. We also 115 compare motifs that evolve to filter out external spurious signals to those that evolve in 116 response to intrinsic stochastic noise in gene expression. We find that while both diamonds and 117 C1-FFLs evolve in response to the former, only diamonds evolve in response to the latter. regulation is allowed, but not shown. Following Milo et al. [1], we exclude the case in which A 123 and B regulate one another, rather than treating this case as two overlapping FFLs. 124 Models 125 Overview of the model 126 We simulate the dynamics of TRNs as the TFs activate and repress one another's transcription. 127 For each moment in developmental time (i.e. on the timescale of one cell responding to stimuli), 128 we simulate the numbers of nuclear and cytoplasmic mRNAs in a cell, the protein 129 concentrations, and the chromatin state of each transcription start site. Transitions between 130 three possible chromatin states --Repressed, Intermediate, and Active --are a stochastic 131 function of TF binding, and transcription initiation from the Active state is also stochastic. An 132 overview of the model is shown in Fig 2. The pattern of TF binding affects chromatin, which 133 affects transcription rates, eventually affecting the concentration of TFs and so completing 134 regulatory feedback loops. The genotype is specified by a set of cis-regulatory sequences that 135 contain TFBSs to which TFs may bind (which, as nucleotide sequences, are subject to realistic 136 mutational parameters), by which consensus sequence each TF recognizes and with what 137 affinity, and by 5 gene-specific parameters that control gene expression as a function of TF 138 binding: mean duration of transcriptional bursts, mRNA degradation, protein production, and 139 protein degradation rates, and gene length which affects delays in transcription and translation. 140 An external signal is treated like another TF, and the concentration of an effector gene in 141 response is a primary determinant of fitness, combined with a cost associated with gene 142 expression (Fig 2). Mutants replace resident genotypes as a function of the difference in 143 estimated fitness. Parameter values, taken as far as possible from Saccharomyces cerevisiae, are 144 summarized in Table 1 Thus, about 90% of total TFs are bound non-specifically, leaving about 10% free. The relatively 200 small number of specific TFBSs is not enough to significantly perturb the proportion of free TFs, 201 and so for the specific TFBSs with m<3 that are of interest in our model, we simply use Kd*(m) = 202 10Kd(m) to account for the reduction in the amount of available TF due to non-specific binding. 203 We also rescale Kd* from moles/liter to the more convenient number of molecules per cell by 204 multiplying by 3×10 -15 liter × 6.02×10 23 molecules/mole = 1.8×10 9 molecules cell -1 M -1 , for a 205 total multiplication factor of 1. where Ni is the per-cell number of molecules of TF i; note that we assume all TF molecules are 211 located in the nucleus. 212

213
The transition rates between chromatin states (see section below) are a function of the 214 numbers of activators A and repressors R bound to a cis-regulatory region. Note that in our 215 model, each TF is either always an activator, or always a repressor, independently of binding 216 context. The joint probability distribution of A and R is derived in S1 Text section 1. TFs bound is for them to be different TFs (Fig 3B). All other genes are AND-gate-incapable, 237 meaning that their activation requires only one TFBS to be occupied by an activator. denotes 238 the probability of having at least one activator bound for an AND-gate-incapable gene, or two 239 for an AND-gate-capable gene. denotes the probability of having at least one repressor 240 bound. 241 242 Noise in yeast gene expression is well described by a two step process of transcriptional 243 activation [40,41], e.g. nucleosome disassembly followed by transcription machinery assembly. 244 We denote the three possible states of the transcription start site as Repressed, Intermediate, 245 and Active (Fig 2). Transitions between the states depend on the numbers of activator and 246 repressor TFs bound (e.g. via recruitment of histone-modifying enzymes [42,43] Brown et al. [40] reported that the rate of transcription machinery assembly is 263 3.3 min -1 for a constitutively active PHO5 promoter, and 0.025 min -1 when the Pho4 activator of 264 the PHO5 promoter is knocked out. We use this range to set 265 where PA_no_R is the probability of having no repressors and either one (for an AND-gate-269 incapable gene) or two (for an AND-gate-capable gene) activators bound, and _ _ is the 270 probability of having no TFs bound (for AND-gate-incapable genes) or having no repressors and 271 not more than one activator bound (for AND-gate-capable genes). 272

273
The promoter sequence not only determines which specific TFBSs are present, but also 274 influences non-specific components of the transcriptional machinery [47,48]. We capture this 275 via gene-specific but TF-binding-independent rates rAct_to_Int with which the machinery disassembles and a burst of transcription ends. In other words, we let TF binding regulate the 277 frequency of "bursts" of transcription, while other properties of the cis-regulatory region 278 regulate their duration. E.g., yeast transcription factor Pho4 regulates the frequency but not 279 duration of bursts of PHO5 expression, by regulating the rates of nucleosome removal and of 280 transition to but not from a transcriptionally active state [40]. We estimate the distribution of 281 rAct_to_Int from the observed rates of mRNA production of 255 yeast genes [49] that are likely to 282 have similarly low nucleosome occupancy [50] and thus are constitutively open to expression 283 (see S1 Text section 2 for details and also for the bounds of rAct_to_Int). For modeling simplicity, we 284 assume that the core promoter sequence responsible for the value of rAct_to_Int is distinct from 285 the 150-bp sequences in which our TFBSs are found. 286 287 mRNA and protein dynamics 288 Once in the Active state, a gene initiates new transcripts stochastically at rate rmax_transc_init = 6.75 289 mRNA/min [40]. There is a delay before transcription is completed, of duration 1 + L / 600 290 minutes, where L is the length of the ORF in codons (see S1 Text section 3). 291

292
We model a second delay between the completion of a transcript and the production of the first 293 protein from it. The delay comes from a combination of translation initiation and elongation; it 294 ends when the mRNA is fully loaded with ribosomes all the way through to the stop codon and 295 the first protein is produced. We ignore the time required for mRNA splicing; introns are rare in 296 yeast [51]. mRNA transportation from nucleus to cytosol, which is likely diffusion-limited [52, 297 53], is fast even in mammalian cells [54] let alone much smaller yeast cells, and the time it takes 298 is also ignored. The median time in yeast for initiating translation is 0.5 minute [ Table 1 in 55], 299 and the genomic average peptide elongation rate is 330 codon/min [55]. After an mRNA is produced, we therefore wait for 0.5 + L / 330 minutes, and then model protein production as 301 continuous at a gene-specific rate rprotein_syn (see S1 Text section 4 for details of rprotein_syn). 302 303 Protein transport into the nucleus is rapid [56] and is approximated as instantaneous and 304 complete, so that the newly produced protein molecules immediately increase the probability of 305 TF binding. Each gene has its own mRNA and protein decay rates, initialized from distributions 306 taken from data (see S1 Text section 5). 307 308 All the rates regarding transcription and translation are listed in times are simulated between transcription initiation and completion, and between transcript 317 completion and the production of the first protein. Protein production and degradation are 318 described deterministically with ODEs, and updated frequently in order to recalculate TF 319 concentrations and hence chromatin transition rates. We initialize developmental simulations 320 with no mRNA or protein (except for the signal), and all genes in the Repressed state. Details of 321 our simulation algorithm are given in the S1 Text section 6. 322 323

Selection conditions
Filtering out short spurious signals is a special case of signal recognition more generally. In 325 environment 1, expressing the effector is beneficial, and in environment 2 it is deleterious. We 326 select for TRNs that take information from the signal and correctly decide whether to express 327 the effector. In our control condition, the signal is "on" at a constant level when the effector is 328 beneficial in environment 1, and off in environment 2. Fitness is a weighted average across 329 these two environments. In our test condition (Fig 4), the signal is constantly on in environment 330 1 and briefly on (for the first 10 minutes) in environment 2 -selection is to ignore this short 331 spurious signal. The signal is treated as though it were an activating TF whose concentration is 332 controlled externally, with an "off" concentration of zero and an "on" concentration of 1,000 333 molecules per cell, which is the typical per-cell number of a yeast TF [58]. 334

335
We make fitness quantitative in terms of a "benefit" ( ) as a function of the amount of 336 effector protein Ne(t) at developmental time t. Our motivation is the scenario in which the 337 effector protein directs resources from metabolic program I to II. When program II produces 338 benefits, 339 where bmax is the maximum benefit if all resources were redirected to program II, and Ne_sat is 343 the minimum of amount of effector protein to achieve this. Similarly, when program I is 344 beneficial, 345 ( 3)  347   348 We set Ne_sat to 10,000 molecules, which is about the average molecule number of a 349 metabolism-associated protein per cell in yeast [58]. Without loss of generality given that fitness 350 is relative, we set bmax to 1. 351 352 A second contribution to fitness comes from the cost of gene expression C(t) (Fig 2, bottom  353 center). We make this cost proportional to the total protein production rate. We estimate a 354 fitness cost of gene expression of 2×10 -6 per protein molecule translated per minute, based on 355 the cost of expressing a non-toxic protein in yeast [59; see S1 Text section 7 for details]. 356

357
We simulate gene expression for 90 minutes of developmental time (Fig 4), and calculate 358 "cellular fitness" in a given environment as the average instantaneous fitness (B(t)-C(t)) over 359 these 90 minutes. We consider environment 2 to be twice as common as environment 1 (a 360 "signal" should be for an uncommon event rather than the default), and take the appropriate 361 weighted average. This differs from Kimura's [60] equation for fixation probability, but captures the same flavor; 382 due to stochasticity in ̂, fixation probability is a monotonic function of the true difference in 383 fitness. Note that it is possible, especially at the beginning of an evolutionary simulation, for 384 relative fitness to be paradoxically negative. In this rare case, for simplicity, we use the absolute 385 value of ̂ on the denominator. 386

387
If 2000 successive mutants are all rejected, the simulation is terminated; upon inspection, we 388 found that these resident genotypes had evolved to not express the effector in either 389 environment. We refer to each change in resident genotype as an evolutionary step. We stop 390 the simulation after 50,000 evolutionary steps; at this time, most replicate simulations seem to 391 have reached a fitness plateau (S2 Fig); we use all replicates except those terminated early. To 392 reduce the frequency of early termination in the case where the signal was not allowed to 393 directly regulate the effector, we used a burn-in phase selecting on a more accessible 394 intermediate phenotype (see S1 Text section 9). In this case, burn-in occurred for 1000 395 evolutionary steps, followed by the usual 50,000 evolutionary steps with selection for the 396 phenotype of interest (S2 Fig). 397 398 Genotype Initialization 399 We initialize genotypes with 3 activator genes, 3 repressor genes, and 1 effector gene. Cis-400 regulatory sequences and consensus binding sequences contain As, Cs, Gs, and Ts sampled with 401 equal probability. Rate constants associated with the expression of each gene, are sampled from 402 the distributions described above and summarized in Table 1. 403 404 A genotype is subjected to 5 broad classes of mutation, at rates summarized in Table 2 and 406 justified in S1 Text section 8. First are single nucleotide substitutions in the cis-regulatory 407 sequence; the resident nucleotide mutates into one of the other three types of nucleotides with 408 equal probability. Second are single nucleotide changes to the consensus binding sequence of a 409 TF, with the resident nucleotide mutated into one of the other three types at equal probability. 410 Both of these can affect the number and strength of TFBSs. 411 412  (0) 3.5×10 -9 per gene k = 0.5, µ = -5 [2] , σ = 0.776 Mutation to L 1.2×10 -11 per codon Mutation to rprotein_syn 9.5×10 -12 per codon k = 0.5, µ = 0.021 [2] , σ = 0.760 Mutation to rprotein_deg 9.5×10 -12 per codon k = 0.5, µ = -1.88, σ = 0.739 Mutation to rAct_to_Int 9.5×10 -12 per codon k = 0.5, µ = 1.57 [2] , σ = 0.773 Mutation to rmRNA_deg 9.5×10 -12 per codon k = 0.5, µ = -1.19, σ = 0.396 Mutation to these quantitative rates takes the form log 10 ′ = log 10 + Normal( ( − 414 log 10 ), ), where x is the original value of a rate and x' is the value after mutation. See S1 Text 415 section 8 for details. 416 2 The value of this parameter is different during burn-in. See S1 Text section 8 for details. 417 Third are gene duplications or deletions. Because computational cost scales steeply (and non-418 linearly) with network size, we do not allow effector genes to duplicate once there are 5 copies, 419 nor TF genes to duplicate once the total number of TF gene copies is 19. We also do not allow 420 the signal, the last effector gene, nor the last TF gene to be deleted. 421 422 Fourth are mutations to gene-specific expression parameters. Most of these (L, rAct_to_Int, 423 rprotein_syn, rmRNA_deg, and rprotein_deg) apply to both TFs and effector genes, while mutations to the 424 gene-specific values of Kd(0) apply only to TFs. Each mutation to L increases or decreases it by 1 425 codon, with equal probability unless L is at the upper or lower bound. Effect sizes of mutations 426 to the other five parameters are modeled in such a way that mutation would maintain log-normal stationary distributions for these values, in the absence of selection or arbitrary bounds 428 (see S1 Text section 8 for details). Upper and lower bounds (S1 Text section 8) are used to 429 ensure that selection never drives these parameters to unrealistic values. 430 431 Fifth is conversion of a TF from being an activator to being a repressor, and vice versa. The signal 432 is always an activator, and does not evolve. 433 434 Importantly, this scheme allows for divergence following gene duplication. When duplicates 435 differ due only to mutations of class 4, i.e. protein function is unchanged, we refer to them as 436 "copies" of the same gene, encoding "protein variants". Mutations in classes 2 and 5 can create

441
Functional AND-gated C1-FFLs evolve readily under selection for filtering out a short 442 spurious signal 443 We begin by simulating the easiest case we can devise to allow the evolution of C1-FFLs for their 444 purported function of filtering out short spurious signals. The signal is allowed to act directly on 445 the AND-gate-capable effector, so all that needs to evolve is a single activating TF between the 446 two, as well as AND-logic for the effector. We score motifs at the end of a set number of 447 generations (see Methods). Evolved C1-FFLs are scored and classified into subtypes based on 448 the presence of non-overlapping TFBSs (Fig 3B). The important subtype comparison for our 449 purposes being the AND-gated C1-FFL vs. the next three non-AND-gated C1-FFL types combined vanishingly rare. The adaptive hypothesis predicts the evolution of the subtype with AND-452 regulatory logic, which requires both the effector to be stimulated both by the signal and by the 453 slow TF. While all replicates show large increases in fitness, a multimodal distribution of final 454 fitness states is observed, indicating whether or not the replicate was successful at evolving the 455 phenotype of interest rather than becoming stuck at an alternative locally optimal phenotype 456 ( Fig 5A). AND-gated C1-FFLs frequently evolve in the high fitness outcomes, but not the low 457 fitness outcomes (Fig 5B). 458

459
We also see C1-FFLs that, contrary to expectations, are not AND-gated; while found primarily in 460 the low fitness replicates, some are also in the high fitness genotypes (Fig 5B). However, this is 461 based on scoring motifs and their logic gates on the basis of all TFBSs, even those with two 462 mismatches and hence low binding affinity. Unless these weak TFBSs are deleterious, they will 463 appear quite often by chance alone. A random 8-bp sequence has probability ( 8 2 ) × 0.25 6 × 464 0.75 2 = 0.0038 of being a two-mismatch binding site for a given TF. In our model, a TF has the 465 potential to recognize 137 different sites in a 150-bp cis-regulatory sequence (taking into 466 account steric hindrance at the edges), each with 2 orientations. Thus, by chance alone a given 467 TF will have 0.0038 × 137 × 2 ≈ 1 two-mismatch binding sites in a given cis-regulatory 468 sequence (ignoring palindromes for simplicity), compared to only ~0.1 one-mismatch TFBSs. 469 Excluding two-mismatch TFBSs when scoring motifs significantly reduces the non-AND-gated C1-470 FFLs, while only modestly reducing the observed frequency of adaptively evolved AND-gated C1-471 FFLs in the high fitness mode (Fig 5C).  Because one TRN can contain multiple C1-FFLs of different subtypes, "Any logic" will generally 481 be less than the sum of the occurrences of all seven subtypes. See S1 Text section 10 for details 482 on the calculation of the y-axis. (C) The over-representation of AND-gated C1-FFLs becomes 483 even more pronounced relative to alternative logic-gating when weak (two-mismatch) TFBSs are 484 excluded while scoring motifs. Data are shown as mean±SE of the occurrence over replicate 485 evolution simulations. n = 23 for high-fitness group, and n = 24 for low-fitness group. 486 487 To confirm the functionality of these AND-gated C1-FFLs, we mutated the evolved genotype in 488 two different ways (Fig 6A) to remove the AND regulatory logic. As expected, this lowers fitness 489 in the presence of the short spurious signal but increases fitness in the presence of constant 490 signal, with a net reduction in fitness (Fig 6B). This is consistent with AND-gated C1-FFLs 491 representing a tradeoff, by which a more rapid response to a true signal is sacrificed in favor of 492 the greater reliability of filtering out short spurious signals. replicate with AND-gated C1-FFLs and lacking other potentially confounding motifs (see S1 Text 500 section 11 for details). (B) Destroying the AND-logic slightly increases the ability to respond to 501 the signal, but leads to a larger loss of fitness when short spurious signals are responded to. 502 Data are shown as mean±SE over replicate evolutionary simulations. 503 504 To test the extent to which C1-FFLs can evolve non-adaptively, we simulated evolution under 505 three negative control conditions: 1) neutrality, i.e. all mutations are accepted to become the 506 new resident genotype; 2) no spurious signal, i.e. the effector should be expressed under a 507 constant "ON" signal and not under a constant "OFF" signal; 3) harmless spurious signal, i.e. the 508 effector should be expressed under a constant "ON" environment whereas effector expression 509 in the "OFF" environment with short spurious signals is neither punished nor rewarded beyond 510 the cost of unnecessary gene expression. AND-gated C1-FFLs evolve much less often under all 511 three negative control conditions (Fig 7). Non-AND-gated C1-FFLs do evolve under the negative 512 control conditions (Fig 7A), but disappear when weak TFBSs are excluded during motif scoring 513 (Fig 7B). To enforce indirect regulation, we ran simulations in which the signal was not allowed to bind to 533 the cis-regulatory sequence of effector genes. The fitness distribution of the evolutionary 534 replicates has only one mode (S4 Fig), so we compared the highest fitness, lowest fitness, and 535 median fitness replicates. In agreement with results when direct regulation is allowed, 536 genotypes of low and medium fitness contain few AND-gated C1-FFLs, while high fitness 537 genotypes contain many (Fig 8A, left). While visually examining the network context of these C1-FFLs, we discovered that many were 551 embedded within AND-gated "diamonds" to form "FFL-in-diamonds" (Fig 8A right). This led us 552 to discover that AND-gated diamonds also occurred frequently without AND-gated C1-FFLs to 553 form "isolated diamonds" (Fig 8A middle). Note that it is in theory possible, but in practice 554 uncommon, for diamonds to be part of more complex conjugates. Systematically scoring the 555 AND-gated isolated diamond motif confirmed its high occurrence (Fig 8B, middle). 556 557 An AND-gated C1-FFL integrates information from a short/fast regulatory pathway with 558 information from a long/slow pathway, in order to filter out short spurious signals. A diamond 559 achieves the same end of integrating fast and slow transmitted information via differences in 560 the gene expression dynamics of the two regulatory pathways, rather than via topological length 561 (Fig 9). 562 563 Note that a simple transcriptional cascade, signal -> TF -> effector, has also been found 564 experimentally to filter out short spurious signals, e.g. when the intermediate TF is rapidly 565 degraded, dampening the effect of a brief signal [62]. Two such transcriptional cascades 566 involving different intermediate TFs form a diamond, so the utility of a single cascade is a 567 potential explanation for the high prevalence of double-cascade diamonds. However, in this case we would have no reason to expect marked differences in expression dynamics between 569 the two TFs, as illustrated in Fig 9. We will also see below that AND-gates evolve between the 570 two cascades. Each TFs is a different protein, and each is encoded by 3 gene copies, shown separately in 577 colors, with the total in thick black. The expression of one TF plateaus faster than that of the 578 other; this is characteristic of the AND-gated diamond motif, and leads to the same functionality 579 as the AND-gated C1-FFL. The two TFs are indistinguishable topologically, but can be easily and 580 reliably assigned identities as "fast" and "slow" by using the fact that the fast TF degrades faster 581 (has higher rprotein_deg). We use the geometric mean rprotein_deg over gene copies of a TF in order to 582 differentiate the two TFs for analysis in Fig 9 and elsewhere. 583 584 Weak TFBSs make motif scoring more difficult 585 Results depend on whether we include weak TFBSs when scoring motifs. Weak TFBSs can either 586 be in the effector's cis-regulatory region, affecting how the regulatory logic is scored, or 587 upstream, affecting only the presence or absence of motifs. When a motif is scored as AND-588 gated only when two-mismatch TFBSs in the effector are excluded, we call it a "near-AND-589 gated" motif. Recall from Fig 3B that effector expression requires two TFs to be bound, with 590 only one TFBS of each type creating an AND-gate. When a second, two-mismatch TFBS of the 591 same type is present, we have a near-AND-gate. TFs may bind so rarely to this weak affinity TFBS 592 that its presence changes little, making the regulatory logic still effectively AND-gated. A near-593 AND-gated motif may therefore evolve for the same adaptive reasons as an AND-gated one. Fig  594   8B and C shows that both AND-gated and near-AND-gated motifs are enriched in the high fitness 595 genotypes. 596 597 When we exclude upstream weak TFBSs while scoring motifs, FFL-in-diamonds are no longer 598 found, while the occurrence of isolated C1-FFLs and diamonds increases (Fig 8C). This makes 599 sense, because adding one weak TFBS, which can easily happen by chance alone, can convert an 600 isolated diamond or C1-FFL into a FFL-in-diamond (added between intermediate TFs, or from 601 signal to slow TF, respectively). 602 603 AND-gated isolated C1-FFLs appear mainly in the highest fitness outcomes, while AND-gated 604 isolated diamonds appear in all fitness groups (Fig 8C), suggesting that diamonds are easier to 605 evolve. 18 out of 30 high-fitness evolutionary replicates are scored as having a putatively 606 adaptive AND-gated or near-AND-gated motif in at least 50% of their evolutionary steps when 607 upstream weak TFBSs are ignored (close to addition of bars in Fig 8C, because these two AND-608 gated motifs rarely coexist in a high-fitness genotype). The remaining 12 have more complex 609 arrangements of weak TFBSs that mimic a single strong one. If we add a two-mismatch TFBS instead, this converts an AND-gated motif to a near-AND-gated 619 motif. This lowers fitness only when the extra link is from the slow TF to the effector, and not 620 when the extra link is from the fast TF to the effector (Fig 10B.ii, 10C.ii). Indeed, these extra 621 links are tolerated during evolution too: if we take the 7 high-fitness replicates that contain a 622 near-AND-gated C1-FFL in at least 5% of the evolutionary steps, in all 7 cases this motif is near-623 AND-gated rather than AND-gated because of an extra weak TFBS for the fast TF, while this is 624 never due to a weak TFBS for the slow TF in C1-FFLs. Similarly, out of the 20 high-fitness 625 replicates that contain a near-AND-gated diamond, 11 cases are primarily because of an extra 626 weak TFBS of the fast TF, 9 cases (all of them OR-gated) are because of weak TFBSs for both TFs, and no cases are primarily due to an extra TFBS for the slow TF. By chance alone, fast and slow 628 TF should be equally likely to contribute the weak TFBS that makes a motif near-AND-gated 629 rather than AND-gated. This non-random occurrence of weak TFBSs creating near-AND-gates 630 illustrates how even weak TFBSs can be shaped by selection against some (but not all) motif-631 breaking links. 632

Fig 10. Perturbation analysis shows that AND-gated C1-FFLs (A) and diamonds (B) filter out 633
short spurious signals. We add a strong TFBS (i) or a two-mismatch TFBS (ii) or (iii); the latter 634 creates near-AND-gated motifs. Allowing the effector to respond to the slow TF alone slightly 635 increases the ability to respond to the signal, but leads to a larger loss of fitness when effector 636 expression is undesirable. Allowing the effector to respond to the fast TF alone does so only 637 when the conversion uses a strong TFBS not a two-mismatch TFBS. short spurious external signals (Fig 8B), (iii) is based on 18 of the 31 replicates from selection for 641 signal recognition in the absence of an external spurious signal (Fig 11B). The 26 and 31 642 replicates were the ones with AND-gated diamond. Replicate exclusion was based on the co-643 occurrence of other motifs with the potential to confound results (see S1 Text section 11 for 644 details). Data are shown as mean±SE of the averaged fitness over replicate evolutionary 645 simulations. 646 647 AND-gated isolated diamonds also evolve in the absence of external spurious signals 648 We simulated evolution under the same three control conditions as before, this time without 649 allowing the signal to directly regulate the effector. In the "no spurious signal" and "harmless 650 spurious signal" control conditions, motif frequencies are similar between low and high fitness 651 genotypes (S5 Fig, S6 Fig), and so our analysis includes all evolutionary replicates. When weak 652 (two-mismatch) TFBSs are excluded, AND-gated isolated C1-FFLs are seen only after selection 653 for filtering out a spurious signal, and not under other selection conditions (Fig 11A). However, 654 AND-gated isolated diamonds also evolve in the absence of spurious signals, indeed at even 655 higher frequency (Fig 11B). Results including weak TFBSs are similar (S7 Fig). Perturbing the AND-gate logic in these isolated diamonds reduces fitness via effects in the 657 environment where expressing the effector is deleterious (Fig 10B.iii). Even in the absence of 658 external short spurious signals, the stochastic expression of intermediate TFs might effectively 659 create short spurious signals when the external signal is set to "OFF". It seems that AND-gated 660 diamonds evolve to mitigate this risk, but that AND-gated C1-FFLs do not. The duration of 661 internally generated spurious signals has an exponential distribution, which means that the 662 optimal filter would be one that does not delay gene expression [63]. The two TFs in an AND-663 gated diamond can be activated simultaneously, but they must be activated sequentially in an 664 AND-gated C1-FFL; the shorter delays possible with AND-gated diamonds might explain why 665 only diamonds and not FFLs evolve to filter out intrinsic noise in gene expression. "Neutral", n = 50 for "No spurious signal", and n = 60 for "Harmless spurious signal". We reused 675 data from Fig 8 for Knabe et al. [23] suggested that including a cost for gene expression may suppress 695 unnecessary links and promote motifs. However, we found AND-gated C1-FFLs still evolve in the 696 high-fitness genotypes under selection for filtering out a spurious signal, even when there is no 697 cost of gene expression (S8 Fig). 699 AND-gated C1-FFLs express an effector after a noise-filtering delay when the signal is turned on, 700 but shut down expression immediately when the signal is turned off, giving rise to a "sign-701 sensitive delay" [3,7]. Rapidly switching off has been hypothesized to be part of their selective 702 advantage, above and beyond the function of filtering out short spurious signals [63]. We 703 selected only for filtering out a short spurious signal, and not for fast turn-off, and found that 704 this was sufficient for the adaptive evolution of AND-gated C1-FFLs. 705 706 Most previous research on C1-FFLs has used an idealized implementation (e.g. a square wave) of 707 what a short spurious signal entails [4,63,69]. In real networks, noise arises intrinsically in a 708 greater diversity of forms, which our model does more to capture. Even when a "clean" form of 709 noise enters a TRN, it subsequently gets distorted with the addition of intrinsic noise [70]. 710 Intrinsic noise is ubiquitous and dealing with it is an omnipresent challenge for selection. 711 Indeed, we see adaptive diamonds evolve to suppress intrinsic noise, even when we select in 712 the complete absence of extrinsic spurious signals. 713 714 Our model, while complex for a model and hence capable of capturing intrinsic noise, is 715 inevitably less complex than the biological reality. However, we hope to have captured key 716 phenomena, albeit in simplified form. E.g., a key phenomenon is that TFBSs are not simply 717 present vs. absent but can be strong or weak, i.e. the TRN is not just a directed graph, but its 718 connections vary in strength. Our model, like that of Burda et al. [68] in the context of circadian 719 rhythms, captures this fact by basing TF binding affinity on the number of mismatch deviations 720 from a consensus TFBS sequence. While in reality, the strength of TF binding is determined by 721 additional factors, such as broader nucleic context and cooperative behavior between TFs 722 (reviewed in Inukai et al. [71]), these complications are unlikely to change the basic dynamics of 723 frequent appearance of weak TFBSs and enhanced mutational accessibility of strong TFBSs from 724 weak ones. Similarly, AND-gating can be quantitative rather than qualitative [72], a 725 phenomenon that weak TFBSs in our model provide a simplified version of. Note that our 726 model, while powerful in some ways, is computationally limited to small TRNs. 727 Core links in adaptive motifs involve strong not weak TFBSs. However, weak (two-mismatch) 728 TFBSs can create additional links that prevent an adaptive motif from being scored as such. 729 Some potential additional links are neutral while others are deleterious; the observed links are 730 thus shaped by this selective filter, without being adaptive. Note that there have been 731 experimental reports that even weak TFBSs can be functionally important [73,74]; these might, 732 however, better correspond to 1-mismatch TFBSs in our model than two-mismatch TFBSs. 733 Ramos et al. [74] and Crocker et al. [73] identified their "weak" TFBSs in comparison to the 734 strongest possible TFBS, not in comparison to the weakest still showing affinity above baseline. 735 736 A striking and unexpected finding of our study was that AND-gated diamonds evolved as an 737 alternative motif for filtering out short spurious external signals, and that these, unlike FFLs, 738 were also effective at filtering out intrinsic noise. Diamonds are not overrepresented in the TRNs 739 of bacteria [2] or yeast [75], but are overrepresented in signaling networks (in which post-740 translational modification plays a larger role) [76], and in neuron networks [1]. In our model, we 741 treated the external signal as though it were a transcription factor, simply as a matter of 742 modeling convenience. In reality, signals external to a TRN are by definition not TFs (although 743 they might be modifiers of TFs). This means that our indirect regulation case, in which the signal 744 is not allowed to directly turn on the effector, is the most appropriate one to analyze if our 745 interest is in TRN motifs that mediate contact between the two. Note that if we were to score the signal as not itself a TF, we would observe adaptive C1-FFLs but not diamonds in this case, in 747 agreement with the TRN data. However, this TRN data might miss functional diamond motifs 748 that spanned levels of regulatory organization, i.e. that included both transcriptional and other 749 forms of regulation. The greatest chance of finding diamonds within TRNs alone come from 750 complex and multi-layered developmental cascades, rather than bacterial or yeast [77]. Multiple 751 interwoven diamonds are hypothesized to be embedded with multi-layer perceptrons that are 752 adaptations for complex computation in signaling networks [30]. 753

754
The function of a motif relies ultimately on its dynamic behavior, with topology merely a means 755 to that end. The C1-FFL motif is based on two pathways between signal and effector, one much 756 faster than the other, which is achieved by making them different lengths. This same function 757 was achieved non-topologically in our adaptively evolved diamond motifs. Multiple motifs have 758 previously been found capable of generating the same steady state expression pattern [21]; 759 here we find multiple motifs for a much more complex function. 760 761 It is difficult to distinguish adaptations from "spandrels" [8]. Standard procedure is to look for 762 motifs that are more frequent than expected from some randomized version of a TRN [2,78]. 763 For this method to work, this randomization must control for all confounding factors that are 764 non-adaptive with respect to the function in question, from patterns of mutation to a general 765 tendency to hierarchy -a near-impossible task. Our approach to a null model is not to 766 randomize, but to evolve with and without selection for the specific function of interest. This 767 meets the standard of evolutionary biology for inferring the adaptive nature of a motif [13][14][15][16][17][18][19][20][21][22][23]. 768 769 Acknowledgements 771 Work was supported by the University of Arizona and by a Pew Scholarship to JM, John 772 Templeton Foundation grant 39667 to JM and KX, and by National Institutes of Health grants 773 R35GM118170 to MLS and R01GM076041 to JM and AKL. We thank Hinrich Boeger for helpful 774 discussions and careful reading of the manuscript, Jasmin Uribe for early work on this project, 775 and the high-performance computing center at the University of Arizona for generous 776 allocations. 777