The adaptive landscapes of three global Escherichia coli transcriptional regulators

eLife Assessment

This study maps the genotype-phenotype landscapes of three E. coli transcription factors and the topographical features of these landscapes. It shows that ruggedness and epistasis do not hinder the evolution of strong transcription factor binding sites. These convincing findings contribute important insights into fitness landscape theories and highlight the role of chance, contingency, and evolutionary biases in gene regulation. The authors then study the topographical features of these landscapes, especially the number and distribution of local maxima, as well as the statistical properties of evolutionary paths on these landscapes.

https://doi.org/10.7554/eLife.103774.3.sa0

Significance of the findings:

Important: Findings that have theoretical or practical implications beyond a single subfield

Landmark
Fundamental
Important
Valuable
Useful

Strength of evidence:

Convincing: Appropriate and validated methodology in line with current state-of-the-art

Exceptional
Compelling
Convincing
Solid
Incomplete
Inadequate

During the peer-review process the editor and reviewers write an eLife Assessment that summarises the significance of the findings reported in the article (on a scale ranging from landmark to useful) and the strength of the evidence (on a scale ranging from exceptional to inadequate). Learn more about eLife Assessments

Abstract
eLife digest
Introduction
Results
Discussion
Materials and methods
Appendix 1
Data availability
References
Article and author information
Metrics

Abstract

The evolution of gene regulation is a major source of evolutionary adaptation and innovation, particularly when organisms encounter new or changing environments. Central to this process is the emergence of new transcription factor binding sites (TFBSs). Adaptive landscapes provide a powerful framework to study such emergence by linking regulatory DNA sequences to their transcriptional outputs. Although several landscapes have been characterized for DNA, RNA, and proteins, large-scale in vivo adaptive landscapes for bacterial TFBSs remain scarce. Here, we address this gap by experimentally mapping the first comprehensive in vivo regulatory landscapes for three global transcription factors in Escherichia coli: cAMP receptor protein, Fis, and IHF. Using a massively parallel reporter assay, we quantify the regulation strength of more than 30,000 TFBS variants for each factor, and reconstruct their adaptive landscapes. All three landscapes are highly rugged and exhibit pervasive epistasis, with thousands of local peaks distributed broadly across sequence space. This ruggedness contrasts sharply with the much smoother TFBS landscapes of eukaryotes. It suggests greater constraints on the evolution of prokaryotic gene regulation. Nonetheless, evolutionary simulations show that ~10% of evolving populations can reach a peak of strong regulation, a proportion that is significantly greater than in comparable random landscapes. Adaptive evolution starting from the same DNA sequence can attain different high peaks, and some peaks are reached more frequently than others. Together, our results show that de novo adaptive evolution of new gene regulation in bacteria is feasible, but subject to a blend of chance, historical contingency, and evolutionary biases.

eLife digest

Organisms regulate the activity of their genes by turning them up or down at the right time and place. To do so, they use proteins called transcription factors, which bind to short stretches of DNA near a gene known as binding sites.

How tightly a transcription factor binds to its binding site determines how strongly it regulates the gene's activity. As genomes evolve, mutations can create new binding sites, establishing entirely new connections within a cell's regulatory network. This is an important way in which organisms adapt to new environments during evolution. Yet the evolutionary steps that transform a random DNA sequence into a strong, functional binding site remain poorly understood, especially in bacteria. Mapping all possible binding sites for a transcription factor can reveal how readily evolution can generate new regulatory interactions.

To test how evolution builds a strong binding site from scratch, Westmann, Goldbach and Wagner studied three major transcription factors in the bacterium Escherichia coli. For each transcription factor, they measured how strongly 30,000 different DNA sequences regulated a gene. To do this, they linked each sequence to a reporter gene that emitted light when activated and measured light emission in millions of cells.

Using these data, the researchers constructed genetic landscapes, in which each DNA sequence occupied a position in the landscape and its height (or elevation, analogous to the height of a hill) represented the strength of gene regulation. Higher peaks, therefore, corresponded to greater gene activation, likely reflecting stronger transcription factor binding.

All three landscapes contained thousands of distinct peaks, far more than comparable landscapes reported for animals and plants. Most of these peaks corresponded to only weak or moderate gene regulation, while peaks of strong regulation were rare. Despite this complexity, the analyses showed that evolution could still reach these rare, strong peaks through a series of small beneficial mutations more often than expected by chance. In other words, even highly complex landscapes can contain accessible routes to strong new binding sites.

The work of Westmann et al. sheds light on how evolution creates and refines the simplest elements of gene regulation. These findings improve our understanding of how bacteria adapt to new environments, including the evolution of antibiotic resistance. In addition, the experimental and analytical approaches developed in this study provide valuable tools for researchers investigating the evolution of gene regulation and for those engineering synthetic genetic circuits in the laboratory.

Introduction

Transcriptional regulation controls the ability of RNA polymerase to initiate the transcription of a gene into mRNA (Browning and Busby, 2016). It is crucial in the life of all organisms, orchestrating the expression of genes in response to environmental cues and cellular states (Cases et al., 2003; Seshasayee et al., 2006). Transcriptional regulation is mediated by transcription factors (TFs), proteins that bind DNA near a gene at short DNA words known as transcription factor binding sites (TFBSs). The binding of a TF to its binding site on DNA can either hinder or facilitate transcription initiation by RNA polymerase (Browning and Busby, 2016; Barnard et al., 2004; Browning and Busby, 2004). The stronger the binding of a TF to its TFBS is, the more strongly the TF can regulate a nearby gene (Aguilar-Rodríguez et al., 2017; Sharon et al., 2012; Shultzaberger et al., 2010; Haldane et al., 2014). In bacteria, TFs operate within a hierarchically organized gene regulatory network. At the bottom of this hierarchy are local TFs that regulate the expression of one or few genes and modulate specific biological processes (Shen-Orr et al., 2002; Martínez-Antonio and Collado-Vides, 2003; Lozada-Chávez et al., 2008). At the top are global TFs that may regulate hundreds of genes and control many cellular functions (Martínez-Antonio and Collado-Vides, 2003; Kurafeiski et al., 2019; Browning et al., 2019; Visweswariah and Busby, 2015). Consistent with their broader role, global regulators are typically more highly expressed and bind to a broader range of TFBSs than local regulators (Lozada-Chávez et al., 2008).

A genotype–phenotype (GP) map is a conceptual analog of a physical landscape, where each location corresponds to one genotype in sequence space, and is associated with a quantitative phenotype. If the phenotype is related to gene expression, the map is also called a regulatory landscape (Aguilar-Rodríguez et al., 2017; Vaishnav et al., 2022; Schweizer and Wagner, 2021; Westmann et al., 2024b; Westmann et al., 2024a). Another special case of such a map is a fitness landscape or adaptive landscape, in which the phenotype is fitness and is interpreted as an elevation (Wright, 1932; Wright, 1931). Over the last decade, GP maps and fitness landscapes have become central tools for understanding how molecular systems evolve under mutation and selection (Blanco et al., 2019; Yi and Dean, 2019; Bank, 2022; Fragata et al., 2019). Such maps and landscapes have been experimentally studied for DNA (Aguilar-Rodríguez et al., 2017; Shultzaberger et al., 2010; Westmann et al., 2024b; Westmann et al., 2024a; Otwinowski and Nemenman, 2013; Chattopadhyay et al., 2025), protein (Herrera-Álvarez et al., 2025; Weeks and Ostermeier, 2023; Steinberg and Ostermeier, 2016; Papkou et al., 2023a; Sarkisyan et al., 2016), and RNA (Schuster, 2002; García-Galindo et al., 2023; Bendixsen et al., 2019) molecules, revealing key topographical properties that shape evolutionary outcomes, including epistasis (Bank, 2022; de Visser et al., 2011) – the non-additive effects of multiple mutations on phenotype – landscape ruggedness, reflected in the number and distribution of fitness peaks, and constraints on adaptive evolution. For example, one large-scale study adopted CRISPR-Cas9 technology to measure the fitness of more than 200,000 Escherichia coli genotypes that encode variants of the bacterial antibiotic resistance gene dihydrofolate reductase (DHFR). It showed that this landscape is highly rugged (multi-peaked; Papkou et al., 2023a).

For TFBSs, most pertinent large-scale studies are based on in vitro binding assays, such as protein-binding microarrays (PBMs), and they focus predominantly on eukaryotic TFs (Aguilar-Rodríguez et al., 2017). While these studies have been instrumental in characterizing TF binding preferences, they typically do not measure regulatory output in a native cellular context. In contrast, comprehensive in vivo data for bacterial TFBSs remain extremely rare. To our knowledge, only two high-resolution in vivo landscapes have been previously mapped for bacterial regulators, those of the local regulators TetR (Westmann et al., 2024b) and LacI (Chattopadhyay et al., 2025). As a result, it remains unclear whether principles inferred from protein landscapes, eukaryotic TFBSs, or in vitro binding assays generalize to transcriptional regulation in bacteria, particularly for global regulators (Martínez-Antonio and Collado-Vides, 2003) that integrate multiple physiological signals.

Both TFs and their TFBSs evolve, but TFBSs evolve more rapidly. The reason is that a mutation in a TF can affect the expression of many genes, whereas a mutation in a TFBS may affect only one gene and is thus less likely to be deleterious (Kurafeiski et al., 2019; Wray, 2007; Signor and Nuzhdin, 2018; Wittkopp et al., 2004). A special case of TFBS evolution is the evolution of a strong TFBS from a DNA word with weak or no regulatory activity. Such de novo evolution of TFBSs can create new regulatory interactions and change the structure of gene regulatory networks (Yona et al., 2018; Fuqua and Wagner, 2023; McAdams et al., 2004; Dorman et al., 2018). Unfortunately, we know little about the ability of Darwinian evolution to create TFBSs de novo. A strong TFBS may have to arise from a non-binding site through an evolutionary path of multiple mutational steps. Darwinian evolution can favor this process only if strong binding is adaptive and if each mutational step in a path increases binding strength, that is if the path is accessible to Darwinian evolution. To find whether such paths exist may require the analysis of multiple paths. This is challenging because sequence space contains an astronomical number of potential TFBSs for any one TF. The number of evolutionary paths to strong transcriptional regulation is thus also astronomical. For each such path, the strength of each TFBS along the path has to be measured experimentally (Kinney and McCandlish, 2019; de Visser et al., 2018; Draghi and Ogbunugafor, 2023; Louis, 2016).

In principle, one could attempt to construct such landscapes in silico using commonly employed models of TF–DNA interactions, such as position weight matrices (PWMs; D’haeseleer, 2006; Stormo, 2000; Stormo and Zhao, 2010). However, PWMs are derived from a limited set of naturally occurring binding sites and primarily reflect sequence conservation rather than quantitative regulatory output. Moreover, PWMs assume independent and additive contributions of individual nucleotide positions to DNA binding (D’haeseleer, 2006; Stormo, 2000; Stormo and Zhao, 2010). They therefore cannot capture the influences of epistatic interactions between positions, which can dramatically affect landscape topography and evolutionary accessibility (Bank, 2022). Lastly, PWMs do not account for important biological effects that modulate gene regulation such as DNA shape (Zhou et al., 2015; Mathelier et al., 2016), cooperative interactions (Ibarra et al., 2020; Bintu et al., 2005), and chromosomal context (Jones et al., 2014; Scholz et al., 2019). Thus, experiments are necessary not only to obtain quantitative measurements of gene regulatory activity, but also to refine and inform PWM-based models using thousands of experimentally characterized sequences.

Building on our previous work on a local TF (Westmann et al., 2024b), we address this challenge for three global regulators in E. coli by performing three independent and massively parallel experiments (Kinney and McCandlish, 2019) for each TF. The experiments use a synthetic biology platform to quantify how strongly each of more than 30,000 binding sites for a TF can regulate the expression of a nearby reporter gene.

The first TF we study is the cAMP receptor protein (CRP). In the absence of glucose, CRP modulates the expression of genes mostly involved in carbon metabolism. It allows E. coli to efficiently switch between sources of carbon and energy (Kim et al., 2018; Borirak et al., 2015; Pal et al., 2022; Khankal et al., 2009). The second TF is the factor for inversion stimulation (Fis), which helps to regulate the expression of genes involved in growth phase transitions. It also modulates the supercoiling of DNA (Nowak-Lovato et al., 2013; Kahramanoglou et al., 2011; Gawade et al., 2020) and influences DNA replication, recombination, and repair (Nowak-Lovato et al., 2013; Kahramanoglou et al., 2011; Gawade et al., 2020). The third factor is the integration host factor (IHF). IHF regulates genes involved in stress responses and stationary phase survival (Prieto et al., 2012; Arfin et al., 2000). Similar to Fis, it is involved in DNA compaction, replication, and recombination, but also in the assembly of complex nucleoprotein structures (Prieto et al., 2012; Arfin et al., 2000). We chose these factors because they are the most global regulators in E. coli, and they are diverse, belonging to different protein families.

We use our experimental data for each TF’s binding sites to map the relationship between binding site genotype and gene expression. Our first aim is to characterize the resulting regulatory landscapes for global bacterial regulators, and to find out whether these landscapes are different or similar. When strong regulation of a gene is adaptive, a regulatory landscape becomes a fitness landscape (Wright, 1932; Wright, 1931; Fragata et al., 2019; Taylor and Provine, 1987). Our second aim is to understand how populations would evolve on each of our three landscapes when they are viewed as fitness landscapes. Specifically, we study how evolving populations would traverse each landscape through individual mutational steps that change a TFBS's ability to regulate a gene via its cognate TF. A peak in such a landscape is a TFBS conveying stronger regulation than all neighboring TFBSs in sequence space. If such a landscape is rugged (has multiple peaks), natural selection alone may not enable a population to discover the highest peaks, that is the strongest TFBSs. The reason is that the peaks may be separated by valleys of low regulation strength that cannot be traversed by natural selection alone (Fragata et al., 2019; Lynch et al., 2016; Iwasa et al., 2004; Taylor and Higgs, 2000). In other words, high peaks may not be easily accessible through Darwinian evolution – they may be reachable by few or no evolutionary paths of single DNA mutations in which each mutational step increases TFBS strength (Aguilar-Rodríguez et al., 2017). One factor that can reduce peak accessibility is epistasis, which can reduce the predictability of evolutionary trajectories towards a peak (Bank, 2022; de Visser et al., 2011).

To accomplish both aims, we first studied the topography of the three landscapes, and then simulated the dynamics of evolving populations on them. All three landscapes are highly rugged. They contain more than 2000 peaks that are scattered through genotype space and are rife with epistatic interactions, in striking contrast to the comparatively smooth TFBS landscapes described for eukaryotic systems (Aguilar-Rodríguez et al., 2017). Despite these features, evolving populations can reach the strongest TFBSs in all three landscapes. Which of several high peaks is reached is contingent on chance events during adaptive evolution.

Results

Landscape mapping

We constructed a modular plasmid system based on the common backbone plasmid pCAW-Sort-Seq-V2 (Westmann et al., 2024b; Supplementary file 1; Figure 1—figure supplements 1–2; Appendix 1, The design of plasmid pCAW-Sort-Seq-V2, Construction of the plasmid pCAW-Sort-Seq-V2 and its variants). This backbone contains all shared regulatory and reporter elements required for fluorescence-based measurements, but it encodes neither a TF nor its binding site(s) for any of the regulators studied here. From this backbone, we generated three TF-specific plasmid derivatives. Each of these plasmids encodes one of the global TFs CRP, Fis, or IHF under inducible control (hereafter pCAW-Sort-Seq-V2-CRP, pCAW-Sort-Seq-V2-Fis, and pCAW-Sort-Seq-V2-IHF; Figure 1—figure supplement 2). In each of the TF-specific plasmids, a TFBS insertion site is positioned immediately upstream of the gfp reporter gene, such that TF binding represses the transcription of gfp. Consequently, stronger TF–DNA binding results in lower GFP expression, which enables a quantitative readout of a binding site’s regulation strength.

For each TF-specific plasmid, we constructed three kinds of variants. The first carries a wild-type (WT) TFBS for the corresponding TF upstream of gfp. It serves as a reference conferring wild-type regulation of the reporter gene. The second contains a complete TFBS library, in which we randomized the eight most information-rich base-pair positions of the respective binding site (Appendix 1, Library design, synthesis, and cloning; Supplementary file 3), as determined from alignments of experimentally characterized binding sites curated in RegulonDB (Salgado et al., 2024). Each library comprised 4⁸=65,536 unique TFBS sequences, and we constructed three independent biological replicate plasmid libraries from them per TF (Appendix 1, Library design, synthesis, and cloning Analysing and sorting cells). The third variant lacks a promoter upstream of gfp and serves as a negative control that allows us to quantify cellular autofluorescence during fluorescence measurements (Supplementary file 1).

We introduced each TF-specific plasmid into an E. coli host strain in which the corresponding chromosomal TF gene had been deleted (Δcrp, Δfis, or Δihf; Supplementary file 2). This genetic background ensures that the focal TF is expressed exclusively from the plasmid. Although the mutant strains grow more slowly than the WT, they reached similar cell densities during late exponential or early stationary phase, the growth phase at which we performed all measurements (Figure 1—figure supplement 3). TF expression is controlled by a tetracycline-inducible promoter and can be precisely tuned using anhydrotetracycline (aTc), allowing us to regulate TF abundance independently of growth conditions and to isolate the effects of TF–DNA binding on transcriptional regulation (Figure 1a).

Figure 1 with 9 supplements see all

Download asset Open asset

Mapping transcription factor binding site (TFBS) regulatory landscapes.

(**a–c**) Sort-Seq procedure. We utilized a plasmid-based fluorescence reporter system followed by sort-seq to map TFBS regulatory landscapes. (a) Library generation. For generating our *E. coli* libraries, we cloned TFBS sequence variants (4⁸=65,536), represented as black circles, into our plasmid, between a σ (Salgado et al., 2024) constitutive promoter (shown as a black right-facing arrow) and a *gfp* gene (shown as a green right-facing arrow), to measure transcriptional repression through fluorescence intensity. When a transcription factor (TF) binds to the TFBS, it blocks RNA polymerase activity through steric hindrance, reducing *gfp* transcription. Thus, lower fluorescence levels (light green) represent stronger TF-TFBS binding. The TFBS variants produce a range of regulation strengths resulting in variable green fluorescence intensities in bacterial cells (green-colored rounded rectangles). (b) Sorting procedure. We sorted libraries into expression bins based on fluorescence intensity using a fluorescence-activated cell sorting (FACS, Figure 1—figure supplements 4–6). (c) Sequencing and phenotyping. We sequenced TFBS variants from each fluorescence bin and used these data to calculate a continuous regulation strength for each genotype. Regulation strength (S) was computed as a weighted average of fluorescence across bins, based on the distribution of sequencing reads (see Methods). Regulation strength is visualized using a color gradient from dark green (low-affinity TFBSs, high fluorescence) to light green (high-affinity TFBSs, low fluorescence). (d) Genotype-phenotype mapping. To construct a regulatory landscape, we connected TFBS genotypes (colored circles) that differed by a single nucleotide via edges, thereby establishing an interconnected genotype network. Each genotype conveys a regulatory phenotype (strength of regulation, heatmap colors), which can be viewed as the elevation dimension (z-axis) in a landscape.

We then mapped these TFBS genotypes to their respective regulatory phenotypes using a well-established technique known as sort-seq (Vaishnav et al., 2022; Kinney and McCandlish, 2019; Peterman and Levine, 2016; de Boer et al., 2020b; Feng et al., 2023; Figure 1b), which combines fluorescence-activated cell sorting (FACS) with high-throughput sequencing (Figure 1c). In sort-seq, one first sorts cells into multiple fluorescence ‘bins’ depending on the level of GFP expression. A cell’s GFP fluorescence serves as a proxy for GFP expression levels (Garcia et al., 2011) and quantifies how strongly the TFBS library member in this cell can regulate GFP expression (Methods, Appendix 1, Analysing and sorting cells). For each TFBS genotype, we quantified regulation strength (S) as a weighted average of its sequencing counts across the different fluorescence bins, yielding a single continuous measure of regulatory activity (see Methods).

The results of our sort-seq experiments are three maps – one for each TF – from each of more than 30,000 genotypes (TFBS variants) to regulatory phenotypes (regulation strength). One can view each map as a regulatory landscape (Shultzaberger et al., 2010; Haldane et al., 2014; Vaishnav et al., 2022; Westmann et al., 2024b; Friedlander et al., 2017; Le et al., 2018; Maerkl and Quake, 2007; Poelwijk et al., 2007; Nghe et al., 2020; Figure 1d) that becomes a fitness landscape whenever strong regulation entails high fitness (Aguilar-Rodríguez et al., 2017; Fragata et al., 2019). For the purpose of analyzing the landscape quantitatively, we represent it as a network of genotypes (TFBS variants). Each node in this network corresponds to a genotype. Edges link neighboring genotypes, which differ in a single nucleotide (Figure 1d).

Landscapes exhibit diverse regulation strengths and distribution breadths

To evaluate the ability of our library to regulate gene expression, we first measured the distribution of GFP fluorescence intensities across the bins in two conditions, i.e., in the presence or absence of the TF (with or without the atc inducer, Figure 1—figure supplements 4–6). In the presence of the TF, GFP expression was lower on average and showed a broader distribution than in the absence of the TF (Figure 1—figure supplements 4–6). This indicates that the TFBSs in each library can indeed downregulate GFP expression, but some library variants convey stronger regulation than others, hence the broader fluorescence distribution.

We then pooled barcoded DNA sequences extracted from cells in each bin and biological replicate. We sequenced at least 250 unique TFBS genotypes from each bin, with an average of 3992.3±4293.9 unique sequences per bin for the three TFs (as detailed in Methods, Appendix 1, Data analysis, and Figure 1—figure supplements 4–6). The resulting sequences covered 95%, 90%, and 93% of the total library sizes (N=65,536 genotypes) for CRP, Fis, and IHF, respectively. To ensure the reliability of our data, we excluded sequences with low read coverage and sequences that were not present in all triplicates (Methods, Appendix 1, Combining data from triplicates). This quality filtering step resulted in library sizes of 31,975 genotypes for CRP (49% of all 4⁸ genotypes), 43,222 genotypes for Fis (66%), and 41,325 genotypes for IHF (63%). The correlation in read counts for each variant across replicates was very high for this quality-filtered data (Pearson’s R=0.98–0.99; Figure 1—figure supplements 7–9). We used this data for all further analyses.

The fluorescence and sequence data from each bin allowed us to quantify the variant’s ability to regulate gene expression for each TFBS library variant. We refer to the resulting metric as the regulation strength S conveyed by the variant (Methods, Appendix 1, Calculating regulation strengths, Figure 1d). It ranges between S=0 (no regulation) to S=1 (strongest regulation among all variants in the library). Although our experiments directly quantify expression regulation, and not the affinity or binding strength of a TF to a TFBS variant, we also use binding strength as a proxy for regulation strength, because TF-DNA binding is necessary for regulation (Browning and Busby, 2016; Barnard et al., 2004; Martínez-Antonio and Collado-Vides, 2003; Bintu et al., 2005; Struhl, 1999; Barnes et al., 2019). To each of our three TFs, we also assigned a reference TFBS that is naturally occurring and conveys strong regulation by the TF, as proven by previous experimental work (Khankal et al., 2009; Kolb et al., 1983; Shao et al., 2008; Huo et al., 2009). We refer to this TFBS as the wild-type (WT_CRP Kolb et al., 1983; Gunasekera et al., 1992_, WT_Fis Shao et al., 2008; Aiyar et al., 2002 and WT_IHF Huo et al., 2009; Ho, 2013, Supplementary file 4). We refer to TFBSs with regulation strengths below and above the wild-type (WT_CRP: S_WT = 0.71, WT_Fis: S_WT = 0.97, WT_IHF: S_WT = 0.95) as weak and strong, respectively.

We observed a broad range of regulation strengths S for each TF landscape, with varying dispersions. The CRP landscape exhibited the lowest average regulation strength and a narrower distribution of S compared to Fis and IHF, which had similar distributions (mean ± SD, 0.37±0.1 for CRP, 0.57±0.13 for Fis, and 0.52±0.14 for IHF; see Figure 2a, Figure 2—figure supplement 1). Next, we analyzed the strongly regulating TFBSs to identify nucleotides that may be particularly frequent, and thus potentially important for strong regulation (Figure 2b). We discovered a moderate association between the most frequent nucleotides (highlighted in yellow in Figure 2b) and the most informative nucleotides from the available PWMs for these TFs (Salgado et al., 2024). (Pearson correlation coefficients: R=0.51, R=0.43, R=0.47 for Fis, CRP, and IHF, respectively, with p-values smaller than 10^–16, rejecting the null hypothesis of an absence of association). This observation suggests that our logos capture different information compared to available PWMs. This is expected because our approach allows us to filter sequences by regulation strength thresholds to construct our matrices, unlike the traditional method of aligning genomic TFBSs without considering their regulation strengths (D’haeseleer, 2006; Stormo, 2000).

Figure 2 with 12 supplements see all

Download asset Open asset

Regulation strength and sequence features of the CRP, Fis, and IHF binding site landscapes.

(a) Genotypes in each landscape vary broadly in their regulation strength. The violin plots show the distribution of regulation strengths S (vertical axis) for each transcription factor (TF) landscape (horizontal axis). The width of a plot at a given value of S represents the frequency of TFBSs at this value. The vertical length of the box in each violin plot covers the range between the first and third quartiles (IQR). The horizontal line within the box represents the median value, and whiskers span 1.5 times the IQR. The white circle shows the regulation strength of the wild type for each landscape (cAMP receptor protein [CRP]: 0.71, Fis: 0.97, integration host factor [IHF]: 0.95). The landscape sizes are as follows: CRP: N=31,975 transcription factor binding site (TFBS) variants; Fis: N=43,222; IHF: N=41,325. (b) Peak genotypes are stronger regulators than non-peak genotypes. Dual violin plots show the distribution of regulation strength (vertical axis) for the three TF landscapes (horizontal axis), stratified by non-peak genotypes (gray, CRP, N=29,821, Fis, N=40,910, IHF, N=38,872), and peak genotypes (red, CRP, N=2154, Fis, N=2312, IHF, N=2453). The black tick-mark in each plot indicates the mean regulation strength of both non-peak and peak genotypes taken together (mean ± SD, CRP: 0.37±0.1, Fis: 0.57±0.13, IHF: 0.52±0.14). The white circle on each plot marks the regulation strength of the wild-type (CRP: 0.71, Fis: 0.97, IHF: 0.95). (c) Sequence logos and nucleotide frequency matrices for strong CRP, Fis, and IHF binding sites. Each sequence logo (D’haeseleer, 2006; Stormo, 2000) is based on an alignment of TFBSs with greater regulation strength than the wild-type for each of the three TFs IHF, Fis, and CRP. Each logo also shows the non-varying position (gray) of the TFBS genotype from which our libraries were created. In each logo, the height of each letter at each TFBS position indicates the information content at that nucleotide position – the taller the letter, the more frequent the nucleotide is in strongly regulating TFBSs (D’haeseleer, 2006; Stormo, 2000). Similar information is conveyed by the frequency matrices displayed as heat maps below each logo. They represent the variability of each nucleotide at each position (horizontal axis) through a color gradient (see color legend). Tall letters in the sequence logo and yellow letters in the frequency matrix indicate frequent, and thus likely important nucleotides for strong regulation.

To further validate the data from our sort-seq experiments, we isolated cells harboring 10 different TFBS variants from each of 13 bins of each library (i.e. 10×13 = 130 variants per library), and determined their regulation strength more directly by quantifying GFP expression with a microplate reader (Methods and Figure 2—figure supplement 2). This comparison validates the sort seq approach by revealing a strong association of the two independent quantifications of regulation strength (Pearson’s R=−0.81 for CRP, R=−0.74 for Fis, and R=−0.73 for IHF). (Methods and Figure 2—figure supplement 2).

All three regulatory landscapes are highly rugged

The study of our landscapes in a network framework (Figure 1d) can help to quantify different aspects of landscape topography (Supplementary file 5). One of them is the ruggedness of each landscape. It can be quantified by the number of peaks (Papkou et al., 2023a; Song and Zhang, 2021; Obolski et al., 2018; Kauffman and Levin, 1987). In the network framework, a peak is a TFBS whose neighbors are all weaker regulators (with lower S) than the TFBS itself. We find that all three landscapes are highly rugged, with 2154, 2312, and 2453 peaks for the CRP, Fis, and IHF landscapes, respectively (Figure 2c, Figure 2—figure supplement 3, Supplementary file 5). Not surprisingly, peak genotypes generally are stronger regulators than non-peak genotypes (Figure 2c). Only a small fraction of peak genotypes regulate expression more strongly than the wild-type (Figure 2c; 61 for CRP, 172 for Fis, and 199 for IHF). We refer to such peaks as high peaks and distinguish them from low peaks (S<S_WT).

The prevalence of sign epistasis (Supplementary file 5) supports the notion that our landscapes are indeed rugged (see Appendix 1, Creation of genotype networks and determining network metrics for further details on epistasis and its evolutionary consequences). Independent evidence for landscape ruggedness comes from comparing the ruggedness of our landscapes with that of a well-established theoretical model of uncorrelated random landscapes, in which each sequence is assigned a fitness at random from the same fitness distribution, and neighbors have uncorrelated fitness values (Kauffman and Levin, 1987; Weinberger, 1990; Orr, 2003). Such landscapes are maximally rugged (Kauffman and Levin, 1987; Weinberger, 1990; Orr, 2003). We created 10³ uncorrelated random landscapes for each of our three TF landscapes by randomly shuffling the measured regulation strengths among all genotypes (Appendix 1, Generating randomly shuffled landscapes). The number of peaks in our TF landscapes lies within 93%, 96%, and 98% of that of a maximally rugged random landscape, which, on average (mean ± SD), has 2308±133, 2405±89, and 2373±97 peaks for CRP, Fis, and IHF, respectively. This analysis underlines that our landscapes are indeed highly rugged.

Because we use only quality-filtered genotype data, our landscapes lack regulatory information for about 40% of the 4⁸ genotypes. While we cannot exclude that this undersampling of genotypes has led to systematic biases in our determination of landscape ruggedness, we note that the sampling of landscape genotypes by our experiments was not strongly biased with respect to regulation strength. Specifically, when we analyzed the relative connectivity of genotypes in our landscapes – the fraction of each genotype’s 24 possible neighbors for which our experiments yielded regulatory data – we found that it is only weakly correlated with regulation strength (R=−0.1,–0.1, 0.01 for the CRP, Fis, and IHF landscapes, Figure 2—figure supplements 4–6). Similarly, the relative connectivity of peak genotypes is only weakly correlated with their regulation strength (R=−0.05, –0.04, 0.06 for the CRP, Fis, and IHF landscapes).

Landscape peaks are widely scattered in genotype space

In a rugged adaptive landscape, reaching high fitness peaks can be challenging, because such peaks are separated from other genotypes by valleys of low fitness that cannot be traversed by natural selection alone (Fragata et al., 2019). If natural selection favors strong gene regulation, the ruggedness of our regulatory landscapes may thus present a challenge for adaptive evolution. To better understand the potential magnitude of this challenge, we next analyzed our landscapes’ topography in greater detail. We began by studying the distribution of peaks in genotype space. If a landscape’s highest peaks are widely scattered through genotype space, then they may be accessible via fitness-increasing evolutionary paths from diverse non-peak genotypes. This may facilitate adaptive evolution compared to a landscape where peaks are clustered in a small region of genotype space.

For each of our three landscapes, we determined the genetic distance between peaks, that is the minimum number of mutations needed to transition from one peak to another, regardless of their effect on regulation strength. We compared the distribution of these distances to the distribution of genetic distances for an equal number of randomly selected non-peak variants. In all three landscapes, the distances between peaks are almost indistinguishable from those of random genotypes, differing on average by fewer than 0.1 mutational steps. In other words, peaks are about as widely dispersed in each landscape as random genotypes (Figure 2—figure supplement 7). This is the case for both low and high peaks (Figure 2—figure supplements 8–9). A principal component analysis further underscores this dispersion (Appendix 1, Principal component analysis, Figure 2—figure supplements 10–12).

Accessible paths to a peak are often not the shortest possible paths

Next, we focused on mutational paths to peaks that are evolutionarily accessible, that is a series of mutational steps where each step increases regulation strength. Specifically, we enumerated all accessible paths that exist from each non-peak genotype to each peak genotype for all three landscapes. We found that these paths are generally longer than the shortest genetic distance between a non-peak and the peak genotypes. Specifically, the mean length of accessible paths exceeded the shortest distance by 1.5, 1.8, and 1.8 steps for the CRP, Fis, and IHF landscapes (Figure 3—figure supplement 1). In other words, accessible mutational paths are often not the most direct paths.

The existence of indirect paths implies that some mutational steps are evolutionarily prohibited because they decrease regulation strength. The reason is closely linked to non-additive (epistatic) effects of two or more mutations on regulation strength. (see Appendix 1, Creation of genotype networks and determining network metrics for further details on epistasis). More specifically, such inaccessible steps are a result of sign epistasis. In this kind of epistasis, a double mutant of a TFBS regulates expression more strongly than the TFBS itself, even though one or both constituent single mutants regulate expression more weakly than the TFBS (Bank, 2022; Saona et al., 2022; Poelwijk et al., 2011; Kvitek and Sherlock, 2011). Indeed, epistasis is prevalent in all three landscapes (Supplementary file 5). Specifically, we observe sign epistasis in 62%, 66%, and 65% of interactions between single mutant pairs in the CRP, Fis, and IHF landscapes.

The existence of accessible paths alone does not tell us how easily a high peak can be found through Darwinian evolution. The reason is that only a very small fraction of all paths to that peak may be accessible, and an evolving population may not find any one of those paths. To quantify the fraction of accessible paths, we determined the total number of paths from each non-peak genotype to each high-fitness peak, and computed the fraction of those paths that are accessible. Figure 3a shows how this fraction depends on the number of mutational steps in a path. It behaves similarly for all three landscapes (Figure 3a). The majority of two-step paths (a fraction greater than 80%) are accessible in all three landscapes (Figure 3a), but the fraction of accessible paths dwindles rapidly with an increasing number of mutational steps. It reaches a minimum below 1% for all three landscapes at eight mutational steps (Figure 3a).

Figure 3 with 4 supplements see all

Download asset Open asset

Peaks and their basins of attractions.

(a) The fraction of accessible paths declines with increasing path length. The figure illustrates how the fraction of accessible paths to high peaks (vertical axis, logarithmic scale) decreases with the length of the shortest accessible path (horizontal axis). Accessible paths are defined as paths where each step increases regulation, starting from a specified initial genotype. Each colored line corresponds to data from a different TF landscape: cAMP receptor protein (CRP; green), Fis (orange), and integration host factor (IHF; blue). Circles indicate the mean fraction of accessible paths for a given path length. (b) Higher peaks have larger basins of attraction. The split-violin plots display the distribution of basin sizes (vertical axis) for high (red) and low (blue) peaks in the three TF landscapes CRP, Fis, and IHF (horizontal axis). High peaks have significantly larger basins of attraction (Welch Two-Sample t-tests; CRP: t-value=9.0898, df = 63.509, and p-value = 4.22 × 10^–13, with mean basin sizes of x_high=9528.426 and x_low = 5655.035. Fis: t-value=7.617, df = 172.35, and p-value = 1.645 × 10^–12, with mean basin sizes of x_high=14206.70 and x_low=10931.94. IHF: t-value=6.1777, df = 201.99, and p-value = 3.521 × 10^–9, with mean basin sizes of x_high = 10965.872 and x_low = 8664.953.). The violin plots show the distribution of basin sizes for each of the two kinds of peaks. Their width represents the frequency of a given basin size. The vertical length of the box in each plot covers the range between the first and third quartiles (IQR). The horizontal line within the box represents the median value, and whiskers span 1.5 times the IQR. (c) Peak genotypes with larger basins of attraction regulate expression more strongly. The regulation strength of peak genotypes (horizontal axis) is plotted against the size of their basins of attraction (vertical axis), shown as a percentage of all (non-peak genotypes) for peaks in all three landscapes (color legend). Color-coded values of R at the top of the graph indicate Pearson correlation coefficients for each landscape (color legend). Curved lines are based on a linear regression analysis (in linear space), with gray shading indicating 95% confidence intervals (CRP: R²=0.28, N=2154; Fis: R²=0.24, N=2312; IHF: R²=0.23, N=2434). (d) Basins of attractions share many TFBS genotypes. The violin plots with embedded boxplots illustrate the fraction of genotypes shared between basins of attractions (vertical axis) for all pairs of high peaks in the CRP (green), Fis (orange), and IHF (purple) TF landscapes (horizontal axis). Specifically, we computed basin overlap as one minus the Jaccard coefficient (Appendix 1, Creation of genotype networks and determining network metrics). A value of one (zero) indicates that the basin of attraction of two peaks share all (no) genotypes (CRP: 61 high peaks, N=1830 comparisons, Fis: 172 high peaks, N=14,706 comparisons, IHF: 199 high peaks, N=19,701 comparisons). Their width represents the frequency of a given basin overlap. The vertical length of the box in each box covers the range between the first and third quartiles (IQR). The horizontal line within the box represents the median value, and whiskers span 1.5 times the IQR.

High peaks have large basins of attraction that share many TFBS genotypes

In the following analysis, we computed the basin of attraction of each peak, defined as the set of non-peak genotypes from which accessible paths to the peak exist (Papkou et al., 2023a; Conrad, 1990; Krug and Oros, 2023). In other words, a peak’s basin of attraction comprises all genotypes from which adaptive evolution can access the peak. The peaks in our landscapes have widely different basin sizes. They include peaks accessible from a mere fraction of genotypes to those accessible by a majority of genotypes (Figure 3b). Notably, high peaks generally had larger basins of attraction than low peaks (Figure 3c, Figure 3—figure supplement 2). In addition, we found that many genotypes are part of multiple basins of attraction. Adaptive evolution may reach different high peaks starting from any such genotype. For all pairs of high peaks, we computed the pairwise overlap (intersection) of the basins of attraction, that is the fraction of genotypes that are part of both basins of attraction. (Figure 3d, Figure 3—figure supplements 3 and 4). The number of genotypes in this intersection varies widely, but the basins of attraction of high peaks generally share a substantial fraction of genotypes (mean shared genotypes: 48%, 42%, and 36% among all pairs of high peaks for the CRP, Fis, and IHF landscapes, Figure 3d).

Genetic drift facilitates and clonal interference reduces the attainment of high peaks

Taken together, our analyses thus far indicate that high peaks in the CRP, Fis, and IHF landscapes are more accessible than low peaks (Figure 3b), and their basins of attraction also share more genotypes (Figure 3d). However, the accessible evolutionary paths to high peaks are often indirect. In addition, the fraction of accessible paths to any one high peak dwindles rapidly with the distance from that peak (Figure 3a).

To understand how these topographical features jointly affect adaptive evolution, we simulated the evolutionary dynamics on our three landscapes under the assumption that natural selection would favor strong regulation and that fitness is proportional to regulation strength.

Because high peaks constitute only a small fraction of all peaks (2.8%, 0.4%, and 0.5% in the CRP, Fis, and IHF landscapes), we hypothesized that most evolving populations would reach only low peaks. To test this hypothesis, we took advantage of the fact that E. coli has a low genomic mutation rate of 2.2×10⁻¹⁰ per base pair per generation (Lynch et al., 2016), and that we study evolution only within an eight-nucleotide segment of a TFBS. In addition, E. coli has a large effective population size (>10⁸ individuals Lynch et al., 2016), which means that genetic drift is weak and even minor differences in fitness are visible to natural selection (Lynch et al., 2016; Kimura, 1962). These conditions imply that a population on our landscape would evolve in the well-studied strong-selection weak-mutation regime (SSWM; Vaishnav et al., 2022; Bank et al., 2016; Orr, 2002; Gillespie, 1984). In this regime, a population is monomorphic most of the time, until a beneficial mutation arises, which usually sweeps rapidly to fixation. In other words, adaptive evolution can be modeled as an adaptive walk, in which each mutational step is taken with a fixation probability that has been derived by Kimura, 1962. We thus call the resulting adaptive walks Kimura walks (Appendix 1, Simulated adaptive walks). We simulated 10³ such adaptive walks starting from each of 15,000 randomly and uniformly selected starting (non-peak) TFBS genotypes from each landscape. Each random walk comprised up to 25 mutational fixation steps, unless it reached a fitness peak earlier. Each adaptive walk also accounted for known E. coli mutational biases (Lind and Andersson, 2008; Horton and Taylor, 2023; Lee et al., 2012; Appendix 1, Simulated adaptive walks).

As we hypothesized, most adaptive walks indeed reach only low peaks (90%, 85%, and 87% for the CRP, Fis, and IHF landscapes). However, the percentage of walks that terminate at high peaks is several-fold greater than the proportion of high peaks. Specifically, 10% of walks terminate at high peaks in the CRP landscape, which is 3.6-fold higher than the percentage (2.8%) of high peaks in this landscape. In the Fis and IHF landscape, adaptive walks are 2- and 1.6-fold more likely to terminate at high peaks than expected from the proportion of these peaks among all peaks (7.4% and 8.1%) (Figure 4a).

Figure 4 with 10 supplements see all

Download asset Open asset

Peak accessibility, contingency, and biases in adaptive walks.

(a) More than 10% of adaptive walks lead to high peaks. The graph shows the cumulative distribution function (CDF, vertical axis) of regulation strength (horizontal axis) attained at the end of 15 million adaptive walks (15,000 starting genotypes ×1000 adaptive walks each) for each transcription factor (TF) landscape (color legend: cAMP receptor protein (CRP), green; Fis, orange; integration host factor [IHF], purple). Dashed lines intersect the CDF at points equivalent to the wild-type regulation strength for each TF (color legend). They help to infer the proportion y_TF of walks that terminate at high peaks (S>S_WT), which is indicated by the numerical values at the graph’s top for each of the landscapes (color legend). (b) Paths to high peaks are not much longer than genetic distances from the starting genotype. Colored boxplots summarize, for each of the three landscapes (horizontal axis), the distribution of shortest genetic distances (blue) and actual path lengths (green) for adaptive walks terminating at high peaks. Path lengths are typically less than one mutational step longer than genetic distances. Each box covers the range between the first and third quartiles (IQR). The horizontal line within each box represents the median value, and whiskers span 1.5 times the IQR. Numbers atop each plot indicate means ± 1 SD. (c) Different numbers of high peaks are reached from different starting variants. Violin plots integrated with boxplots show the distribution in the number of high peaks reached (vertical axis) from each starting variant that attained any high peak during 1000 adaptive walks. (d) Some high peaks are reached more often than others. We randomly and uniformly sampled 10 starting genotypes from the CRP landscape, started 103 adaptive walks from each, and recorded the number and frequency of distinct high peaks attained in these random walks. Results for each starting genotype are symbolized by a vertical bar. The number of stacks within each bar (delineated by horizontal lines, also indicated by an integer above each bar) indicates the number of high peaks reached by the 103 adaptive walks. Starting variants are ordered in ascending order based on this number of attained peaks. Stack height indicates the fraction of walks that reached the same peak, and is indicated in red, orange, and yellow for the three most frequently attained peaks. See Figure 4—figure supplement 8 for the Fis and IHF landscapes.

To determine whether high regulatory peaks are more accessible due to chance alone, we compared the empirical landscapes to randomized ‘null’ landscapes. Specifically, we generated these randomized landscapes by permuting regulation strengths across genotypes while preserving the sampled genotype network (Appendix 1, Generating randomly shuffled landscapes). On each randomized landscape, we then performed adaptive-walks (15,000 walks for each of 10³ starting genotypes) with the same parameters as for the biological landscape. For all three landscapes, the fraction of adaptive walks reaching high regulatory peaks in the empirical landscapes exceeds that for the randomized landscapes by almost threefold, a difference that is statistically significant (Figure 4—figure supplement 5, Appendix 1, Randomized landscape null model for peak accessibility). In sum, rugged regulatory landscapes strongly constrain evolutionary trajectories, yet do not render the evolution of strong regulation vanishingly rare. Instead, strong regulatory phenotypes remain evolutionarily attainable at levels that exceed null expectations, even though they are reached by only a minority of evolving populations.

The adaptive walks that reached a high peak were short (Figure 4b). On average, they were also no more than half a mutational step longer than the shortest genetic distance between each starting genotype and peak (Figure 4b, CRP landscape: path length 3.2±1.6 [mean ± SD] vs. genetic distance 2.8±1.2; Fis landscape: path length 3.1±1.5 vs. genetic distance 2.7±1.2; IHF landscape: path length 3.1±1.5 vs. genetic distance 2.7±1.2).

These short paths can be explained by our previous observation that peaks are widely distributed across the genotype spaces of the three landscapes. This distribution makes it easier for populations starting from any (non-peak) genotype to find a local peak through few mutations. Additionally, path lengths realized during adaptive walks tend to be shorter than the average lengths of accessible paths. (Figure 4b).

Because genetic drift can cause evolving populations to escape a low fitness peak and attain a higher fitness peak, we also asked how small population sizes affect the likelihood that a population attains a high fitness peak (Appendix 1, Simulated adaptive walks, Figure 4—figure supplements 1–3). When simulating adaptive evolution in populations of N=10² individuals, we found that the likelihood that a population reaches a high peak increased by at least 10% (to 18%, 20%, and 21% for the CRP, Fis, and IHF landscapes, Figure 4—figure supplements 1–3).

We also studied how deviations from the strong selection weak mutation regime would affect evolutionary dynamics. In large populations or at high mutation rates, populations tend to be polymorphic and are subject to clonal interference, where the most advantageous of several co-occurring mutations typically dominates and becomes fixed (Melissa et al., 2022; Li et al., 2024; Stolyarova et al., 2022). To approximate this dynamic, we conducted simulations using ‘greedy’ adaptive walks (Appendix 1, Simulated adaptive walks; Fragata et al., 2019; Park et al., 2016). In such an adaptive walk, it is always a genotype’s mutational neighbor with the highest increase in regulation strength that becomes fixed (Orr, 2002; Park et al., 2016; Kauffman and Levin, 1987; Orr, 2003). In other words, a greedy walk starting from a given genotype is deterministic. We, therefore, simulated only one walk for each of the 15,000 randomly chosen (non-peak) starting genotypes. We found that the fraction of greedy walks reaching high regulatory peaks is slightly lower than in the SSWM regime, with 1%, 2%, and 5% fewer walks achieving such peaks in the CRP, Fis, and IHF landscapes, respectively (Figure 4—figure supplement 4). This outcome is expected, because greedy walks prioritize immediate fitness gains at the expense of the ability to discover higher fitness peaks (Park et al., 2016).

The attainment of peaks is highly contingent on chance events

Because different basins of attraction share many TFBS genotypes (Figure 3d, Figure 3—figure supplements 3 and 4), it is possible that adaptive evolution starting from any one genotype can reach different peaks, depending on which mutational paths it takes. In other words, the structure of the landscapes we study may lead to evolutionary contingency (Blount et al., 2018; Nonoyama and Chiba, 2019) – the dependence of an outcome of evolution on unpredictable chance events (Louis, 2016; Blount et al., 2018; Xie et al., 2021). Indeed, non-peak genotypes can access on average around 30% of the total number of high peaks in each landscape (Figure 3—figure supplement 4). To assess the prevalence of evolutionary contingency, we determined how many different peaks are attained by 10³ adaptive walks starting from the same randomly chosen genotype. We performed this analysis on 15,000 starting genotypes, focusing on the subset of those starting genotypes from which at least one high peak is reached during the 10³ walks. Specifically, 28.6%, 23%, and 24.5% of starting genotypes for the CRP, Fis, and IHF landscapes reached at least one high peak. Notably, these percentages represent the fraction of starting genotypes that can access at least one high peak, not the total fraction of adaptive walks that lead to a high peak, which is lower (as shown in Figure 4a). For instance, while adaptive walks starting from 28.6% of genotypes reached at least one high peak in the CRP landscape, only 10% of all adaptive walks in this landscape end at a high peak (Figure 4a). Importantly, for any one starting genotype where random walks reached at least one high peak, they typically reached more than one high peak (Figure 4c, median ±IQR, 4±4, 4±8, 6±6 for the CRP, Fis, and IHF landscapes, respectively). The number of high peaks attained from any one starting genotype was as large as 21, 30, and 33 for the three respective landscapes. Thus, evolutionary contingency is pervasive in all three landscapes. In addition, not only the end point of evolution but also the evolutionary trajectories leading to this end point show contingency. That is, adaptive walks starting from a specific genotype and ending at a specific high peak often take multiple paths to this peak (Figure 4—figure supplement 6).

Among the different peaks that can be reached from a starting genotype, some may be reached preferentially. Such biases towards specific evolutionary outcomes have been documented in both empirical and computational studies (Westmann et al., 2024b; Papkou et al., 2023a; Louis, 2016; Blount et al., 2018; Spor et al., 2014; Lenski, 2017; Lind et al., 2015). They also exist in our landscapes (Figure 4d, Figure 4—figure supplement 6). Figure 4d illustrates both contingency and bias for 10³ adaptive walks starting from each of 10 different genotypes in the CRP landscape (thus a total of 10⁴ adaptive walks). Depending on the starting variant, the adaptive walks reached between 1 and 16 high peaks. Whenever different peaks are reached, they are reached by markedly different proportions of adaptive walks. The same holds in the Fis and IHF landscapes (Figure 4—figure supplement 8). Additionally, while multiple routes to each peak can be traversed (Figure 4—figure supplement 6), we observed in all three landscapes biases towards specific paths that are taken more frequently than others (Figure 4—figure supplement 9).

The number of these variants is shown as an integer on top of the panel, and its percentage among all 15,000 starting variants is shown in parentheses (color legend). Underneath it, the panel shows the mean ±1 SD of the number of attained high peaks. a. d. Some high peaks are reached more often than others. We randomly and uniformly sampled 10 starting genotypes from the CRP landscape, started 10³ adaptive walks from each, and recorded the number and frequency of distinct high peaks attained in these random walks. Results for each starting genotype are symbolized by a vertical bar. The number of stacks within each bar (delineated by horizontal lines, also indicated by an integer above each bar) indicates the number of high peaks reached by the 10³ adaptive walks. Starting variants are ordered in ascending order based on this number of attained peaks. Stack height indicates the fraction of walks that reached the same peak and is indicated in red, orange, and yellow for the three most frequently attained peaks. See Figure 4—figure supplement 8 for the Fis and IHF landscapes.

Lastly, because sort-seq measurements are subject to experimental uncertainty (Peterman and Levine, 2016; Trippe et al., 2022; Gilliot and Gorochowski, 2023), an important question is whether such noise affects the inferred structure of a regulatory landscape and the resulting evolutionary dynamics. To address this issue, we explicitly incorporated empirically estimated measurement uncertainty into our fitness comparisons and repeated pertinent analyses under two conditions: a noise-free landscape ( $G^{S}$ ) and an ‘uncertainty-aware’ landscape $G^{s, τ}$ using genotype-specific noise estimates (τ, see Figure 4—figure supplement 10, Appendix 1, Incorporating experimental uncertainty into adaptive walks). We examined whether noise-induced changes in landscape structure influence evolutionary trajectories. To do so, we compared adaptive walk dynamics between the two kinds of landscapes using identical population genetic parameters. Despite differences in peak counts and local topology, the overall pattern of genotype visitation during adaptive walks was highly similar between the two kinds of landscapes. Specifically, the visitation frequency profiles of genotypes were strongly correlated between $G^{S}$ and $G^{S, τ}$ landscapes (Spearman’s ρ<0.001 shown in Figure 4—figure supplement 10). This means that genotypes frequently accessed in the noise-free landscape remain frequently accessed in the noisy landscape.

Discussion

We evaluated the ability of three E. coli global regulators to control gene expression through each of more than 30,000 TFBSs for each regulator. To this end, we utilized a synthetic plasmid-based system that facilitates high-throughput fluorescence measurements. This system allowed us to quantify the regulation strength of individual TFBSs by measuring gene expression through GFP fluorescence (Garcia et al., 2011). Additionally, the system insulates the control of transcription from any direct effects library sequences might have on transcription, such as the transcriptional and translational impacts of 5’ untranslated regions on gene expression (Evfratov et al., 2017; Cuperus et al., 2017). Recent studies have investigated large empirical datasets of cis-regulatory genotypes, examining the impact of sequence variation on gene regulation in both eukaryotes and prokaryotes (Aguilar-Rodríguez et al., 2017; Vaishnav et al., 2022; de Boer et al., 2020b; Barnes et al., 2019; Belliveau et al., 2018; Kinney et al., 2010; Lagator et al., 2022; Urtecho et al., 2019). However, few studies have focused on these interactions from an evolutionary perspective and studied adaptive landscapes of gene regulation, as we do here (Aguilar-Rodríguez et al., 2017; Vaishnav et al., 2022; Lagator et al., 2022).

We showed that the regulatory landscapes of all three TFs are highly rugged and have multiple peaks. The ruggedness of all three landscapes is also supported by the prevalence of epistasis between pairs of TFBS mutations (Supplementary file 5). A particularly important form of epistasis is sign epistasis (Bank, 2022; Saona et al., 2022; Poelwijk et al., 2011), because it can lead to multiple adaptive peaks (Bank, 2022; Saona et al., 2022; Poelwijk et al., 2011) (see Appendix 1, Creation of genotype networks and determining network metrics). Our landscapes contain up to 65% of mutation pairs with sign epistasis, a value that is especially high compared to the almost exclusively additive interactions of mutations in eukaryotic TFs (Aguilar-Rodríguez et al., 2017; Aguilar-Rodríguez and Payne, 2021). However, the TFs we study are not exceptions among other prokaryotic TFs. Recent smaller-scale studies have shown that epistatic interactions are common in binding sites of local prokaryotic regulators, such as AraC (sign epistasis in more than 50% of 20 mutants Lagator et al., 2016) and Cl (sign epistasis in 85% of 113 mutants Lagator et al., 2017). A possible reason for this greater incidence of epistasis lies in the nature of prokaryotic TFBSs. Specifically, prokaryotic TFBSs are at approximately 20 bps twice as long as eukaryotic TFBSs (Struhl, 1999; Stewart et al., 2012) and exhibit symmetries that reflect the dimeric state of their cognate TFs (Madan Babu and Teichmann, 2003; Huffman and Brennan, 2002; Perez-Rueda et al., 2018). These factors may increase the likelihood of intramolecular epistasis. Our observations raise important questions for future work, such as why the landscapes of prokaryotic TFBSs differ so dramatically from those of eukaryotic ones. And what do these differences imply for the evolutionary dynamics of gene regulation?

Despite the high ruggedness and pervasive epistasis of these landscapes, we found that a modest fraction of evolving populations can still access the highest regulatory peaks. Specifically, our evolutionary simulations show that 10% of populations with a size typical of E. coli reach one of the highest peaks. This percentage is significantly higher than in randomized landscapes (Appendix 1, Randomized landscape null model for peak accessibility; Figure 4—figure supplement 5), which shows that the structure of our regulatory landscapes facilitates access to stronger regulation than expected by chance. We speculate that this property reflects the biological role of global regulators, which coordinate the expression of many target genes (Martínez-Antonio and Collado-Vides, 2003) and thus must operate across a wide spectrum of regulation strengths, while still permitting the evolution of strong regulation when necessary.

The clonal interference that may occur in even larger populations or in populations with high mutation rates reduces the accessibility of high peaks from 10% to 5% of evolving populations (Figure 4—figure supplement 4). Conversely, in small populations where genetic drift can help a population escape from a low peak, this percentage increases to 18%. Numbers like these render the de novo evolution of strong TFBSs plausible for our three global regulators. Once such binding sites have originated in one population, they can also spread to others through horizontal gene transfer (Oren et al., 2014; Price et al., 2008).

In addition to increasing the ruggedness of a landscape, epistasis can also influence the evolution of populations on the landscape (Bank, 2022; Poelwijk et al., 2011; Weinreich et al., 2005). The regulation strength of a TFBS is partially determined by the combined interactions of the nucleotides within the TFBS (Stormo and Zhao, 2010; Gerland et al., 2002; Zhao et al., 2012; O’Flanagan et al., 2005), which in turn affects TF-TFBS binding affinities. Different combinations of mutations can result in similar regulation strengths and create multiple evolutionary pathways to achieve optimal or near-optimal regulation (Aguilar-Rodríguez et al., 2017; Payne and Wagner, 2014; Figure 4—figure supplement 6). This can lead to overlapping basins of attraction among different peak TFBSs conveying strong regulation, such that even evolution starting from the same genotypes can take different paths and reach different peaks with similar regulation strengths. Indeed, we observed such overlapping basins in our landscapes (Figure 3d and Figure 3—figure supplement 3). Moreover, epistasis can create plateaus of regulation strength where various combinations of nucleotides yield TFBSs with intermediate regulation strengths. These plateaus further contribute to the overlapping basins of attraction we observe, and may serve as common evolutionary intermediates that multiple starting sequences can traverse on their way to different peaks.

In a landscape with substantial epistasis, the sequence of mutations that occurs in an evolving population can also render adaptive evolution highly contingent on this sequence (Blount et al., 2018; Nonoyama and Chiba, 2019). For example, some beneficial mutations within a TFBS might only be accessible after other specific mutations have occurred. Different starting sequences or early mutations can also change the spectrum of accessible mutations, leading to different peaks with similarly strong regulation (Blount et al., 2018; Nonoyama and Chiba, 2019). Indeed, we observed in all three landscapes that different evolving populations starting from the same genotypes in the landscape attain different peaks (Figures 3d and 4d and Figure 4—figure supplement 8; Blount et al., 2018; Xie et al., 2021; Blount et al., 2008; Palmer et al., 2015). Such contingency reduces the predictability of evolution (Louis, 2016; Blount et al., 2018; Vermeij, 2006).

Despite the prevalence of contingency in all three landscapes, evolving populations starting from the same genotype more often attain some peaks than others (Figure 4—figure supplement 8). Moreover, for each attained peak with multiple possible evolutionary routes, we also observed that some paths are more frequently transversed than their alternatives (Figure 4—figure supplement 9). That is, evolution is biased towards traversing some paths and attaining some peaks more often than others (Lässig et al., 2017; de Visser and Krug, 2014). These observations emphasize the complex interplay between chance, contingency, and evolutionary biases in shaping the outcomes of adaptive evolution (Louis, 2016; Blount et al., 2018; Xie et al., 2021; Lenski, 2017).

The three TFs we studied here belong to different protein families, yet they have regulatory landscapes with similar topography. All three landscapes are highly rugged, highly epistatic, and harbor multiple peaks that are widely scattered, and this holds for peaks of low, intermediate, and high regulation strength. In addition, peak accessibility in all three landscapes increases with peak height (Figure 4—figure supplement 7). On the one hand, these commonalities may be caused by common biological or biochemical properties of the three TFs. For example, CRP has been suggested to possess nucleoid-associated protein properties similar to Fis and IHF, due to its ability to bend and loop DNA (Heyde et al., 2021). On the other hand, the commonalities may reflect general characteristics of global gene regulation. One of them is that global TFs often bind unspecifically to multiple TFBSs (Martínez-Antonio and Collado-Vides, 2003; Kurafeiski et al., 2019). Also, they may bind DNA with a broad range of different affinities centered around intermediate affinity, rather than bind few sites but very strongly (Kurafeiski et al., 2019). This possibility is supported by comprehensive analyses of in vitro eukaryotic TF binding affinities for thousands of TFBS variants, showing essentially continuous binding affinity distributions for hundreds of sites (Aguilar-Rodríguez et al., 2017).

Despite broad similarities among the three landscapes we study, we also found some differences. Most notable is the narrower distribution of regulation strengths for the binding sites of CRP compared to those of Fis and IHF (Figure 2a). This is not unexpected, given that CRP binding sites are known for being quasi-symmetric and less degenerate than those of Fis and IHF (Kim et al., 2018; Khankal et al., 2009). Previous studies have also shown that Fis and IHF binding sites are less conserved in their nucleotide sequence and more biased toward AT-rich content (Dorman et al., 2018; Monteiro et al., 2020; Dillon and Dorman, 2010). Differences in the distribution of regulation strengths may also result from differences in functions among the three proteins. CRP’s function is mostly restricted to gene regulation. As a result, it may have evolved to finely tune its regulatory output, resulting in a narrower distribution of regulation strengths. In contrast, Fis and IHF are involved in functions beyond gene regulation, such as DNA replication and genomic structural maintenance (Dorman, 2013; Dorman et al., 2020), which might impose additional constraints or require a broader range of DNA binding strengths.

One limitation of our work stems from our rigorous quality filtering of sort-seq data. As a result, we lack regulatory data for some 40% of the TFBSs in each landscape. This limited diversity of reliable data is a common feature in mutational library studies. It has several technical causes, such as biases in library synthesis (Chen et al., 2020), PCR amplification (Kebschull and Zador, 2015), cloning, and loss of sequence diversity after cell sorting (Kinney et al., 2010). In future work, this limitation could be overcome by a combination of strategies, such as subsampling complete genotype spaces or combining different molecular methods to overcome biases and diversity loss during PCR amplification (Wong et al., 2006; Aird et al., 2011), sorting (Peterman and Levine, 2016; Trippe et al., 2022; Gilliot and Gorochowski, 2023), and high-throughput sequencing (Lagator et al., 2022). Importantly, although undersampling of genotype space is a limitation, standard approaches such as random subsampling or predictive modeling are not straightforward remedies. Several of our core analyses – including peak identification, quantification of epistasis, and assessment of evolutionary accessibility – rely on combinatorially complete local neighborhoods in genotype space. Random subsampling of genotypes would remove mutational neighbors, and thereby confound pairwise comparisons and the interpretation of landscape topology. Predictive modeling could be used to infer missing genotypes and reconstruct more complete landscapes (Wagner, 2022), but it requires additional assumptions that introduce their own limitations. In addition, developing, validating, and benchmarking such models would be beyond the scope of this study, which is focused on empirical landscape mapping, and can serve as a starting point for future modeling work.

A second limitation comes from our use of the sort-seq method, which is best-suited for our work, because it allows the high-throughput measurement and sorting of millions of individual cells in a standardized and straightforward manner. However, the method’s accuracy depends on the binning procedure. Other studies have used between two and several dozen bins, with or without unbiased sampling (Peterman and Levine, 2016; Feng et al., 2023; Barnes et al., 2019; Trippe et al., 2022; Belliveau et al., 2018; Kinney et al., 2010; Boer et al., 2018). Recommendations emerging from this work include sorting of cells into at least four logarithmically (log₂) equally-spaced bins, each covering approximately 12.5–15% of the fluorescence distribution (Barnes et al., 2019; Kinney et al., 2010; de Boer et al., 2020a). We followed these recommendations, using 13 bins. In addition, we computed a weighted average to calculate expression values from sequences appearing in multiple bins, a straightforward method validated by robust studies in transcriptional regulation analysis (Vaishnav et al., 2022; de Boer et al., 2020b). In addition, we validated individual regulation strengths with an independent method to demonstrate its reliability.

A third limitation of our study is that we only examined variation at eight positions for each TF. We selected these positions based on their importance for DNA binding and regulation by our TFs, as evidenced by their high information content. We cannot exclude the possibility that selecting a different set of positions could yield different landscape topographies. However, we speculate that less information-rich positions would not reduce but rather increase the potential for the de novo evolution of strong TFBSs. For example, they may facilitate landscape navigability through extradimensional bypasses (Conrad, 1990; Wu et al., 2016) or provide small-incremental changes in regulation strengths. These small changes help populations reach peaks via diminishing returns effects, meaning that as a population gets closer to a peak, each subsequent mutation contributes progressively smaller improvements to regulation strength (Diaz-Colunga et al., 2023; MacLean et al., 2010). Investigating a wider array of positions and larger regulatory landscapes remains an important task for future work.

Fourth, we use a simplified empirical system for the study of gene regulation – a promoter followed by a single TFBS. Although this regulatory architecture exists for some genes (Rydenfelt et al., 2014), global regulators often form part of more complex, combinatorial architectures (Barnard et al., 2004; Rydenfelt et al., 2014; Ishihama, 2010) that are influenced by environmental factors (Cases et al., 2003; Cases and de Lorenzo, 2005; Shepherd et al., 2023), concentrations of active TFs (González Pérez et al., 2009; Ali Azam et al., 1999; Weinert et al., 2014), chromosome structure (Lagomarsino et al., 2015; Dorman and Dorman, 2016; El Houdaigui et al., 2019; Sobetzko et al., 2012; Meyer et al., 2018), methylation states (Sánchez-Romero et al., 2015), and more. Studies on more complex promoter architectures may reveal regulatory landscapes with different topographies.

Lastly, we also analyzed the adaptive evolution of TFBSs in a simplified manner, performing data-driven evolutionary simulations rather than experimental evolution. Although many studies assume that regulation strength and fitness are correlated, this is not always the case. Low binding affinities can be adaptive during development (Crocker et al., 2016), and many genes exhibit a nonlinear fitness-expression function (Srivastava and Payne, 2022) with a plateau of maximal fitness across a wide range of expression levels (Rest et al., 2013; Duveau et al., 2017; Bergen et al., 2016). For instance, most mutations and polymorphisms in the promoter of the yeast gene TDH3 do not significantly affect fitness in a glucose-rich medium (Duveau et al., 2017). For these reasons, we focused our observations and interpretations on regulatory phenotypes (regulation strength). Our simplified approach allowed us to model population evolution in a large genotype space and avoid monitoring thousands of evolving regulatory sequences simultaneously in vivo (Yi and Dean, 2019; Ornelas et al., 2023). However, it cannot replace experimental evolution of TFBSs, which also remains an important challenge for future work.

TFBSs are among the simplest units of biological organization. Our work provides the first large-scale analysis of the regulatory landscapes formed by such sites for global transcriptional regulators. It shows that strong binding sites can readily evolve de novo, even though prokaryotic TFBSs are much larger than their eukaryotic counterparts and have more rugged regulatory landscapes. In addition, the evolution of these simple sequences also displays phenomena that have been characterized in much more complex systems, such as evolutionary contingency and evolutionary biases.

Materials and methods

Key resources table

Reagent type (species) or resource	Designation	Source or reference	Identifiers	Additional information
Strain, strain background (Escherichia coli)	SIG10-MAX	Sigma-Aldrich	Cat# CMC0004	Cloning strain (high transformation efficiency)
Strain, strain background (E. coli)	DH5α	Kelly et al., 2009	–	Cloning strain (genotype in Supplementary file 2)
Strain, strain background (E. coli)	BW25113	KEIO collection (Baba et al., 2006)	–	Wild-type parent of KEIO deletion mutants
Genetic reagent (E. coli)	JW5702-4 (Δcrp)	KEIO collection (Baba et al., 2006)	Δcrp-765::kan	Sort-seq host for CRP
Genetic reagent (E. coli)	JW1702-1 (ΔihfA)	KEIO collection (Baba et al., 2006)	ΔihfA786::kan	Sort-seq host for IHF
Genetic reagent (E. coli)	JW3229-1 (Δfis)	KEIO collection (Baba et al., 2006)	Δfis-779::kan	Sort-seq host for Fis
Recombinant DNA reagent	pCAW-Sort-Seq (plasmid)	Westmann et al., 2024b	–	Parental sort-seq plasmid
Recombinant DNA reagent	pCAW-Sort-Seq-V2 and TF-specific derivatives (plasmids)	This study; https://doi.org/10.5281/zenodo.13838265	–	Full list in Supplementary file 1
Sequence-based reagent	TFBS libraries (CRP, Fis, IHF)	This study; IDT	–	Ultramer ssDNA; see Supplementary file 3
Sequence-based reagent	Reference (wild-type) TFBS	This study	–	See Supplementary file 4
Sequence-based reagent	PCR and barcoding primers	This study; IDT; Eurofins	–	See Supplementary file 6
Commercial assay or kit	Q5 High-Fidelity DNA Polymerase	New England Biolabs	Cat# M0491L	PCR amplification
Commercial assay or kit	NEBuilder HiFi DNA Assembly Master Mix	New England Biolabs	Cat# E2621L	Gibson assembly
Commercial assay or kit	Monarch DNA Gel/PCR Extraction Kit	New England Biolabs	Cat# T1020L	DNA purification
Commercial assay or kit	QIAprep Spin Miniprep Kit	Qiagen	–	Plasmid extraction
Commercial assay or kit	HindIII-HF	New England Biolabs	Cat# R3104	Restriction digestion
Commercial assay or kit	BamHI	New England Biolabs	Cat# R3136	Restriction digestion
Commercial assay or kit	T4 DNA Ligase	New England Biolabs	Cat# M0202L	Ligation
Commercial assay or kit	Quick CIP	New England Biolabs	–	Dephosphorylation
Commercial assay or kit	DpnI	New England Biolabs	Cat# R0176L	Template removal before sequencing
Commercial assay or kit	Exonuclease I	New England Biolabs	Cat# M0293L	ssDNA removal
Chemical compound, drug	Anhydrotetracycline (aTc)	Cayman Chemicals	Cat# 10009542	Induction of TF expression
Chemical compound, drug	Dulbecco’s PBS	Sigma-Aldrich	Cat# D8537	FACS buffer
Chemical compound, drug	D-glucose	Sigma	Cat# G8270	Growth medium component
Chemical compound, drug	Magnesium sulfate	Sigma	Cat# 230391	Growth medium component
Chemical compound, drug	LB medium	Sigma-Aldrich	Cat# L3522	Growth medium
Software, algorithm	Cutadapt	Martin, 2011	RRID:SCR_011841	Adapter trimming
Software, algorithm	FLASH	Magoč and Salzberg, 2011	RRID:SCR_005531	Paired-end read merging
Software, algorithm	FASTX-Toolkit	Hannon, 2010	RRID:SCR_005534	Quality filtering
Software, algorithm	R	R Development Core Team, 2026	RRID:SCR_001905	Statistical analysis and plotting
Software, algorithm	igraph (Python/R)	Csárdi and Nepusz, 2006	RRID:SCR_019225	Genotype-network analysis
Software, algorithm	ggplot2	Wickham, 2016	RRID:SCR_014601	Plotting
Software, algorithm	SnapGene	GSL Biotech; https://snapgene.com	RRID:SCR_015052	Cloning design
Software, algorithm	BD FACSDiva v9.0	BD Biosciences	RRID:SCR_001456	Cell sorting control
Software, algorithm	Custom analysis code	This study; https://doi.org/10.5281/zenodo.13838265	–
Other	FACSAria III cell sorter	BD Biosciences	RRID:SCR_016695	Fluorescence-activated cell sorting

Appendix 1 contains extended details of experimental procedures and data analysis.

Strains and plasmids

Request a detailed protocol

Bacterial strains and plasmids used in this work are listed in Supplementary files 1 and 2. We obtained electrocompetent E. coli cells of strain SIG10-MAX from Sigma Aldrich (CMC0004). We used this strain for molecular cloning and library generation due to its high transformation efficiency. The genotype of this strain (Supplementary file 2) is similar to DH5α (Sigma Aldrich commercial information, see Supplementary file 2). The strain is resistant to the antibiotic streptomycin.

We amplified plasmid libraries in SIG10-MAX, extracted their DNA, and transformed them into mutants derived from E. coli K-12 strain BW25113 that harbor chromosomal deletions of the crp, fis, or ihfa gene. We obtained these mutant strains from the KEIO collection (Baba et al., 2006) and used them for sort-seq experiments. The design, genetic parts, and assembly of the plasmid vectors we used in this study are available in Appendix 1, as are all primers, TFBS sequences/libraries, strains, and plasmids.

Sort-Seq procedure

Request a detailed protocol

To explore the regulatory effects of each TF on binding sites in the corresponding library, we constructed three plasmids, each of which enables the inducible expression of one of our three TFs. These are plasmids pCAW-Sort-Seq-V2-CRP, pCAW-Sort-Seq-V2-Fis, and pCAW-Sort-Seq-V2-IHF (Appendix 1, Construction of the plasmid pCAW-Sort-Seq-V2 and its variants and Supplementary file 1). We cloned the TFBS libraries (Supplementary file 3) into their respective plasmids and then transformed them into mutant strains lacking the corresponding TFs (Δcrp, Δfis, and Δihf, as listed in Supplementary file 2). We induced TF expression using anhydrotetracycline (Atc) and, after overnight growth, performed cell sorting for cell populations. During sorting, we distributed cells into 13 equally spaced logarithmic bins based on their fluorescence levels. We replicated each sort-seq experiment three times from three separate library transformations for each TF. To mitigate the impact of extrinsic noise (gene expression variation among cells Elowitz et al., 2002), we adhered to standard protocols by normalizing our GFP fluorescence measurements against mScarlet-I fluorescence values obtained from flow-cytometry assays (Sharon et al., 2012; Peterman and Levine, 2016; Sharon et al., 2014; Rudge et al., 2016). We subsequently recovered the sorted cells from each bin in 50 mL Falcon tubes containing 10 mL of LB medium supplemented with chloramphenicol and incubated them overnight at 37 °C with shaking at 220 rpm. After this growth period, we re-sorted cell cultures from each bin to eliminate potential contaminants and ensure that the cell populations had preserved their fluorescence distributions. Following re-sorting, we extracted plasmids from the cell population of each bin. We amplified and barcoded the TFBS region from each population through a polymerase chain reaction (PCR). Lastly, we sequenced barcoded amplicons containing TFBS sequences, and used the sequencing results to calculate the regulation strengths of each TF to its TFBSs in the corresponding library. More details are provided in the Supplementary Material.

Regulation strengths

Request a detailed protocol

Due to gene expression and measurement noise, individual TFBS variants in a sort-seq experiment usually appear in more than a single bin, and their read count (frequency) varies among bins (Peterman and Levine, 2016; de Boer et al., 2020b; Belliveau et al., 2018; Gilliot and Gorochowski, 2023). Following established practice (Westmann et al., 2024b; de Boer et al., 2020b; Lagator et al., 2022), we used a weighted average of these frequencies for each variant to represent the mean expression level caused by the variant. To facilitate the interpretation of this quantity, we converted this expression level into a regulation strength relative to the highest observed regulation strength for a given TF, to which we assigned a value of one. From each library, we selected a single naturally occurring binding site for each TF that was previously characterized in the literature as a strong binder. We called this TFBS the WT sequence and used it as a baseline to separate weakly (higher GFP expression) from strongly (lower GFP expression) regulating TFBSs.

Validating regulation strengths with plate reader measurements

Request a detailed protocol

To further validate our regulation strength data, we chose 10 DNA binding sites from each bin (10 variants ×13 bins=130 variants in total per TF library, Supplementary file 2), covering a wide range of measured regulation strengths. We cloned these sequences into the appropriate vector pCAW-Sort-Seq-V2-TF (TF: CRP, Fis, or IHF), and transformed them into the appropriate mutant strain. We picked individual colonies and grew them overnight (16 hr, 37 °C, 220 rpm) in liquid LB supplemented with 50 μg/mL of chloramphenicol and anhydrotetracycline. We diluted the cultures to 1:10 (v/v) in cold Dulbecco’s PBS (Sigma-Aldrich #D8537) to a final volume of 1 mL. We transferred 200 μl of the diluted cultures into individual wells in 96-well plates and measured GFP fluorescence (emission: 485 nm/excitation: 510 nm, bandpass: 20 nm, gain: 50), as well as the optical density at 600 nm (OD₆₀₀) as an indicator of cell density. We then normalized fluorescence by the measured OD₆₀₀ value to account for differences in cell density among cultures and compared the obtained ratios to the previously inferred regulation strengths for the 130 selected variants. We performed all such measurements in biological and technical triplicates (three colonies per sample, and three wells per colony, respectively).

Code availability and data analysis

Request a detailed protocol

All code used for processing data and plotting, as well as the final processed data, plasmid sequences, and primer sequences are available in the Zenodo public repository and is accessible via the following DOI: https://doi.org/10.5281/zenodo.13838265.

Supplementary information is linked to the online version of the paper.

Appendix 1

General procedures

Although all general procedures have been previously described (Westmann et al., 2024b), we describe them again below for completeness.

Media and reagents

To prepare SOB medium, we mixed 25.5 g of solid medium stock (VWR J906) with 960 ml of purified water and subjected the resulting suspension to autoclaving. To prepare the SOC medium, we dissolved 20 ml of 1 M D-glucose (Sigma G8270) and 20 ml of 1 M magnesium sulfate (Sigma 230391) in 960 ml of pre-prepared SOB solution. For the LB medium, we combined 25 g of solid medium stock (Sigma-Aldrich L3522) with 1 liter of purified water and then autoclaved it. To prepare the M9 minimal medium, we diluted M9 minimal salt sourced from Sigma (M6030) in distilled water as per the manufacturer’s guidelines, autoclaved the solution, and added 0.4% glucose (Sigma G8270), 0.2% casamino acids (Merk Millipore, 2240), 2 mM magnesium sulfate (Sigma 230391), and 0.1 mM calcium chloride (Sigma C7902). Where required, we supplemented growth media with chloramphenicol (50 µg/mL), anhydrotetracycline (100 ng/mL, Cayman Chemicals #10009542), and/or glucose (0.4% w/v final concentration). We prepared anhydrotetracycline by diluting the dried chemical in absolute ethanol, from a stock concentration (1000 X) of 100 µg/mL to a working concentration of 100 ng/mL.

Overnight incubation of cultures in liquid and solid medium

We cultivated bacteria in liquid LB medium (using either 15 mL or 50 mL Falcon tubes), enriched with chloramphenicol at a concentration of 50 µg/mL. We incubated these cultures for a period of 16 hr at a temperature of 37 °C, with a shaking speed of 200 rpm and 50 mm orbital motion, using an Infors HT Multitron Incubator Shaker. Similarly, for cultures in solid medium, we grew bacterial colonies on LB-agar plates (using sterile plastic Petri dishes of 90 mm × 15 mm dimensions), also supplemented with chloramphenicol at a concentration of 50 µg/mL, and incubated them for the same time and at the same temperature.

PCR reactions

Except where specifically mentioned, we amplified DNA fragments through polymerase chain reactions (PCR), employing Q5 high-fidelity polymerase (NEB #M0491L) to minimize mutation introduction into the amplicons. We followed the protocol recommended by NEB, aiming for a final reaction mixture of 50 µL. We conducted each PCR reaction twice and combined the products from these duplicates after the completion of the reaction. We determined the primer melting temperatures (Tm) using the NEB Tm calculator (accessible at https://tmcalculator.neb.com/#!/main), with primers at a concentration of 500 nM.

Verifying PCR products through gel electrophoresis

Unless indicated otherwise, we verified successful PCR amplification via gel electrophoresis to ensure the presence of singular-band amplicons and the absence of non-specific bands. This process involved separating PCR products in a 0.8% agarose Tris-EDTA (TAE) gel. We conducted electrophoresis for 45 min at 120 V or until the bands had progressed beyond the halfway point of the gel’s total length.

DNA purification with commercial kits

Once we had confirmed a successful PCR through gel electrophoresis, we purified the PCR products utilizing the Monarch DNA PCR/Gel Extraction Kit (NEB #T1020L), adhering to the manufacturer’s protocol. In situations requiring further purification (such as the occurrence of non-specific bands post-PCR), we performed gel purification. This involved mixing 10 μL of 6 x NEB DNA dye with each 50 μL PCR product, then loading all 60 μL onto a 1% agarose gel. We carried out electrophoresis for 45 min at 120 V or until the bands had moved more than halfway through the gel. We then excised the DNA band corresponding to the amplified sequence using a scalpel. For extracting DNA from the gel, we used the Monarch DNA Gel Extraction Kit (NEB #T1020L).

Gibson assembly (Gibson et al., 2009)

We assembled PCR-amplified fragments using the NEBuilder-HiFi DNA Assembly Master Mix kit (NEB #E2621L). We determined the molarity required for assembly based on the guidelines outlined by the Barrick Lab (details available at https://barricklab.org/twiki/bin/view/Lab/ProtocolsGibsonCloning). We incubated the assembly mix for 1 hr at 50 °C in a dry bath incubator, followed by cooling on ice for subsequent steps.

Preparation of electrocompetent cells

We utilized glycerol/mannitol step centrifugation to prepare electrocompetent cells (Warren, 2011). To this end, we first cultured the appropriate E. coli strain (Supplementary file 2) in 5 mL SOB medium at 37 °C with a shaking speed of 250 rpm overnight. The next day, we transferred 3 mL of the culture to 300 mL SOB medium and incubated under the same conditions until the OD600 reached a value between 0.4 and 0.6 (measured at an optical path length of 1 cm), which took approximately 2–4 hr. After cooling the culture on ice for 15 min, we centrifuged the cells at 4 °C and 1500 × g for 15 min. We then resuspended the cells in 60 mL of ice-cold distilled water (dH₂O) and divided them into three 50 mL tubes. Gradually, we added 10 mL of an ice-cold glycerol/mannitol solution (consisting of 20% glycerol (w/v) and 1.5% mannitol (w/v)) to each tube using a 10 mL pipette. We centrifuged the tubes at 1500 × g and 4 °C for 15 min in an Eppendorf 5810/5810 R centrifuge with acceleration/deceleration set to zero. After discarding the supernatant, we resuspended cells in 3.0 mL of the same glycerol/mannitol solution. We transferred these suspensions to pre-cooled 1.5 mL tubes and incubated them in a dry ice-ethanol bath for about 1 min. Finally, we stored the suspensions at –80 °C for future transformation experiments.

Electroporation

In all transformation experiments described in this study, we used 100 µL of electrocompetent cells for electroporation, employing 0.2 cm cuvettes (EP202, Cell Projects, UK) and a Micropulser electroporator (Bio-Rad) set to the EC3 setting (15 k V/cm). Post-electroporation, we recovered cells in 1 mL of SOC media, warmed in advance in 15 mL Falcon tubes. This recovery step lasted for 1.5 hr at 37 °C with a shaking speed of 220 rpm. Unless specified otherwise, we spread 300 µL of each recovered culture on an LB agar plate that contained 50 μg/mL chloramphenicol. We then incubated this plate overnight for 16 hr at 37 °C. Subsequently, we confirmed the identity of the clones on the plate via Sanger sequencing.

The design of plasmid pCAW-Sort-Seq-V2

The plasmid pCAW-Sort-Seq-V2 is a derivative of the plasmid pCAW-Sort-Seq18. Briefly, plasmid pCAW-Sort-Seq harbors a pBBR1 replication origin, which ensures a broad host range and a low copy number—typically between 5 and 10 copies per cell (Jahn et al., 2016). It also harbors a chloramphenicol resistance gene, a TetR repression system (Lutz and Bujard, 1997), and a TFBS measuring module. The latter consists of an interchangeable TFBS that lies between a constitutive promoter and a superfolder GFP (sfgfp Pédelacq et al., 2006) reporter gene. In this system, TF-TFBS interactions can decrease GFP production by physically obstructing the bacterial RNA polymerase, a phenomenon known as steric hindrance. The reporter gene sfgfp is insulated by a transcriptional insulator named RiboJ, a synthetic ribozyme that removes 5’ UTR interferences from variable TFBS sequences in the mRNA by self-cleavage (Lou et al., 2012). More details and features have been previously described in Westmann et al., 2024b. Here, we have integrated a bicistronic expression cassette into pCAW-Sort-Seq to enable the regulated expression of a TF of choice. This cassette harbors the gene encoding one of our focal TFs situated upstream of a mscarlet-I (Bindels et al., 2017) reporter gene to monitor the expression of this bicistronic operon via fluorescence. Regulation is achieved through the pLtetO-1 (Lutz and Bujard, 1997) promoter, a synthetic promoter tightly repressed by TetR. Expression is initiated by the addition of anhydrotetracycline (Cayman Chemicals, catalog #10009542), which relieves TetR-mediated repression and activates the expression of the entire cassette.

Construction of the plasmid pCAW-Sort-Seq-V2 and its variants

We initiated the construction of pCAW-Sort-Seq-V2 by synthesizing (Twist Biosciences, California, USA) and cloning a codon-optimized mscarlet-I gene expression cassette into the pCAW-Sort-Seq plasmid. This gene is oriented antiparallel to the tetracycline resistance gene (tetr). This cloning step resulted in the formation of the intermediate pCAW-Sort-Seq-V2t plasmid (Figure 1—figure supplement 1). Subsequently, we amplified for each of our three TFs the encoding gene (crp, fis, and ihf) from the genome of the E. coli strain BW25113 and inserted it upstream of the mscarlet-I gene in a bicistronic operon configuration using Gibson assembly. This process produced TF-specific expression plasmid variants (Figure 1—figure supplement 2). The final phase involved removing the core promoter to create negative controls for each TF and cloning libraries into each plasmid variant upstream of the sfgfp gene. This yielded the plasmids used in our experiments (Supplementary file 1, Figure 1—figure supplement 2). The specific sets of primers used for each PCR reaction are listed in Supplementary file 6.

Library design, synthesis, and cloning

We based the design of the TFBS library for each TF on consensus sequences available in the RegulonDB database (Salgado et al., 2024), the most comprehensive database for E. coli transcriptional studies. We designed each library by randomizing the eight most important TFBS positions, such that each of the four nucleotides had an equal probability (0.25) to occur at each position. The expected library size for this approach is (4⁸)=65,536 sequences. Library compositions can be found in Supplementary file 3. We designed libraries in the Snapgene software (https://snapgene.com), and had them synthesized by IDT (Coralville, USA) as single-stranded DNA Ultramers of 140 bp (4 nmol). We resuspended each library in nuclease-free distilled water and serially diluted it to a concentration of 50 ng/µL. We used Ultramers as templates in a PCR reaction for the formation of dsDNA fragments and library amplification. We amplified the Ultramer DNA molecules with the following PCR program: 98 °C/30 s; 25 cycles of 98 °C/10 s, 60 °C/15 s and 72 °C/80 s; and 1 cycle of 72 °C/5 min. We opted for a maximum of 20 cycles to reduce amplification biases. We analyzed amplification products through gel electrophoresis to confirm that only a single product band with the expected size of around 140 bp was present for each PCR reaction. After confirming the presence of single bands, we purified the products using the Monarch DNA gel extraction kit (NEB #T1020L). Whenever unspecific bands appeared during electrophoresis, we gel-purified the PCR products using the same kit.

For the construction of each plasmid-based library, we digested 1 μg of the purified library with HindIII-HF (NEB #R3104) and BamHI (NEB #R3136) restriction enzymes in a 100 μL reaction, followed by overnight incubation at 37 °C. We isolated the appropriate cloning plasmid from its host strain using the QIAprep spin miniprep kit (Qiagen, Germany), and digested it with the same enzymes. Post-digestion, we added 3 μL of Quick CIP (calf intestinal alkaline phosphatase) to prevent self-ligation by dephosphorylating DNA ends. We purified the ligated DNA using the Monarch DNA gel extraction kit (NEB #T1020L).

We performed the ligation with a 10:1 molar ratio of insert-to-vector, using 100 ng of vector, 10 units of T4 DNA ligase (NEB #M0202L), and 2 µL of 10 X ligation buffer in a 20 µL reaction. We incubated the mixture at 20–22°C for approximately 16 hr, followed by a 15-min inactivation of the ligase at 65 °C. We purified the ligation product, resulting in 10 µL of DNA resuspended in dH₂O. We transformed E. coli SIG10-MAX cells by electroporation with this purified product.

Post-transformation, we plated 50 µL of each recovered culture on LB agar for colony-forming unit counting (cfu) and transformation efficiency estimation. From the agar plates, we selected 30 colonies for colony PCR and Sanger sequencing to assess library diversity (NightSeq service, Microsynth, Switzerland). Our mean transformation efficiency was 10⁶ cells per transformation. We diluted the remaining 950 µL of transformants in 9 mL of LB medium with chloramphenicol and cultured it overnight. We then aliquoted the overnight culture into 1 mL cryotubes with 20% glycerol and stored it at –80 °C.

After transforming plasmid libraries into the SIG10-MAX strains for reasons of transformation efficiency and plasmid maintenance, we extracted the plasmids using a QIAprep spin miniprep kit (Qiagen, Germany), and transformed them into the appropriate host mutant host strains Δcrp, Δfis, and Δihf (see the strain genotypes in Supplementary file 2). We cultured each library-transformed strain overnight, aliquoted it in 1 mL cryotubes with 20% glycerol, and stored it at –80 °C for further experimentation.

Analysing and sorting cells

In preparation for cell sorting, we cultivated cells harboring each library in liquid LB medium enriched with chloramphenicol. Specifically, we cultured 1 mL of transformed cell aliquots and a streak of cells with a control plasmid (a promoterless pCAW-Sort-Seq-V2 plasmid lacking sfGFP expression) in 50 ml Falcon tubes with 9 mL LB medium (containing 50 µg/mL chloramphenicol) overnight. Subsequently, we diluted these overnight cultures at a 1:100 ratio (v/v) in LB medium with chloramphenicol and divided them into two aliquots. We supplemented one aliquot with the anhydrotetracycline (Atc) inducer to express the plasmid-encoded TF. The other aliquot remained unchanged. We incubated both cultures for 5 hr until late-exponential/early-stationary phase (200 RPM, 37 °C). Subsequently, we diluted 20 µL of the cultures in 1 mL of cold filtered Dulbecco’s PBS (Sigma-Aldrich #D8537) in 15 mL FACS tubes.

We performed FACS-sorting on a FACS Aria III flow cytometer (BD Biosciences, San Jose, CA) using a 70 µm nozzle. We utilized a 488 nm laser for detecting forward scatter (FSC) and side scatter (SSC) with a 488 nm/10 nm band-pass filter. We set the flow rate to 1.0, adjusting sample dilution as necessary to achieve no more than ≈10,000 events/second. Considering the small size of bacterial cells, we reduced the particle detection threshold to the lowest feasible setting (200 arbitrary units on FSC and SSC channels), increasing it to a maximum of 500 units if background noise was excessive. We then adjusted FSC-H and SSC-H for cells with the negative control plasmid to center the bacterial population in the cytometer’s software (BD FACSDiva Software v9.0). We based the sorting and binning of cells on both mScarlet-I (PE-Texas-Red channel, excitation laser: YellowGreen 561 nm, LP filter: 600 nm, BP filter: 610/20 nm) and sfGFP fluorescence (FITC channel, excitation laser: 488 nm, LP filter: 502 nm, BP filter: 530 nm/30 nm), setting both the PE-Texas-Red and FITC channels voltages so that the median fluorescence of the negative control was between 0 and 100 (arbitrary units) on the FITC-H and PE-Texas-Red-H axes.

Initially, the sorting process involved establishing a gate for identifying cells that are red-fluorescence-positive, that is reporter and thus TF-expressing. To establish this gate, we began by measuring the autofluorescence of our negative control culture on the PE-Texas-Red-H axis. Thereafter, we examined the fluorescence of a positive control that had been induced with Atc and was expressing the mScarlet-I protein. We then configured a gate around the positive mScarlet-I-expressing population, ensuring that all cells sorted through this gate exhibited a level of reporter expression higher than the negative control.

Next, we established the green-fluorescent sorting gates. For setting sorting gates on the FITC-H axis, we first recorded the autofluorescence of the negative control culture. This median autofluorescence defined the upper boundary of the lowest bin (B1) for the experimental population. We then analyzed the fluorescence of 10⁶ cells expressing sfGFP and containing the library, without sorting, in order to establish the boundaries of the next binning gates. We took the lower bound of the highest bin (B13) for the experimental population to correspond to the 95th percentile of the fluorescence distribution of this population. We chose boundaries between the remaining 11 intermediate bins with equidistant spacing on a binary logarithmic (log₂) scale. After we had determined the gates in this way, we calculated the fractions of the previously 10⁶ recorded cells that fell inside each bin.

To maintain statistical robustness and minimize sampling error in our downstream analyses, we wanted to ensure that each of the 65,536 unique sequences in our library was represented by at least 30 cells after sorting. In addition, it was crucial to account for an anticipated diversity loss of up to 70% of the cell population during sorting, due to factors such as cell death and the dilution of low-frequency and lower-fitness genotypes during the post-sorting recovery phase (Kinney et al., 2010). To establish the number of cells required to be sorted initially, we thus determined a multiplicative factor $f$ based on the minimally needed number of 30 cells per sequence and the expected cell retention rate post-sorting of 30%, that is, $f = \frac{30}{0.30} = 100$ cells/sequence. We thus multiplied our initial library size by 100 for sorting purposes, that is, we sorted a total of 6,553,600 cells. This calculation aims to ensure that despite a substantial reduction in cell numbers post-sorting, each unique sequence remains adequately represented in the surviving cell population. We sorted cells into 1.5 mL Eppendorf tubes, each containing 500 µL of LB medium, and kept at 4 °C to prevent growth during sorting and sample processing. We replicated the entire sorting procedure three times, each based on independent library transformations.

After sorting, we added 1 mL of LB without antibiotics to each tube and removed 20 µL of the resulting volume for serial dilutions. We transferred the remaining liquid culture (980 µL) to 50 mL falcon tubes and allowed each culture to recover for 2 hr (37 °C, 220 rpm). After recovery, we added 9 mL of LB supplemented with chloramphenicol, and grew the cultures overnight for freezing part of them in glycerol stock aliquots, for validating our binning procedure (see below), and for extracting plasmid DNA for subsequent PCR and sequencing steps. We used the 20 µL initially removed from each 1 mL culture for preparing two serial dilutions (10^–4 and 10^–6) that we plated on LB-Cm agar plates (200 µL per plate) to estimate the post-sorting viability through colony forming unit (cfu) counting. By knowing how many cells were sorted into each bin, we can estimate how many cells would be expected in our dilutions and compare this number with the cfu counts we observed. In this way, we estimated that on average (across bins) 77% of cells remained viable (standard deviation: 18.15%). We also estimated the genetic diversity of the library through Sanger sequencing of DNA from individual colonies (NightSeq service, Microsynth, Switzerland).

To validate our binning procedure, we re-grew binned cultures from either overnight recovered cultures or frozen aliquot stocks. We then measured their expression distributions by flow cytometry, which reproduced the original expression measurements. This approach allowed us to compare the fluorescence distributions of the re-grown cultures with the pre-sorting distributions. Specifically, we assessed whether the geometric mean of the fluorescence distribution for each sorting bin matched those recorded during the initial sorting procedure. Next, we re-sorted cells derived from each bin in order to eliminate cell cross-contaminations and other factors that could alter the bin distributions. The geometric mean is often preferred over the arithmetic mean in this type of analysis, because it is less affected by the presence of outliers in the data (Beal, 2017; Beal et al., 2016). It is also a more accurate representation of the central tendency of data that is log-normally distributed, which is often the case with flow cytometry fluorescence measurements (Beal, 2017; Beal et al., 2016).

DNA extraction and sequencing

We diluted 500 uL of individual glycerol stocks of each replicate subpopulation of sorted cells (i.e., cells from each bin of fluorescence intensity) in 5 mL of LB supplemented with chloramphenicol in 15 mL Falcon tubes, and grew the resulting cell culture overnight (16 hr, 37 °C, 220 rpm). On the next day, we isolated plasmids from each culture using a QIAprep spin miniprep kit (Qiagen, Germany). In order to allow the sequencing of multiple pooled samples (multiplexing), we barcoded our regulatory region through PCR with specific HPLC-purified primers (Supplementary file 6) provided by Eurofins (Konstanz, Germany). We added barcodes to the 5' region of the amplicon through a PCR reaction. We performed this PCR reaction with the Q5 high-fidelity polymerase, and did so in triplicate for each miniprep-isolated plasmid library. To calculate primer melting temperatures (Tm), we used the NEB Tm calculator (https://tmcalculator.neb.com/#!/main) for a primer concentration of 500 nM. We performed the PCR with the following program: 98 °C/30 s; 25 cycles of 98 °C/10 s, 64 °C/30 s and 72 °C/30 s; and 1 cycle of 72 °C/2 min.

After PCR amplification, we digested the reaction products with the restriction enzymes DpnI (NEB #R0176L) and Exonuclease I (NEB #M0293L) in order to remove traces of genomic DNA, plasmids, and single-stranded DNA that could interfere with sequencing. The Master Mix we used for a single digestion harbored 1 µL of 10 x CutSmart Buffer (NEB #B6004S), 1 µL of Exonuclease I (NEB #M0293L), 1 µL of DpnI restriction enzyme (NEB #R0176L), and 7 µL of distilled nuclease-free water. For each PCR product, we added 10 µL of the Master Mix. We incubated the reaction for 1 hr at 37 °C, following 15 min at 80 °C for deactivation of the enzymes.

We purified the digestion products using the Monarch DNA PCR/Gel Extraction Kit (NEB #T1020L) and fractionated them through gel electrophoresis to confirm that only a single band with a size ≈150 bp was present. After this confirmation, we pooled the purified PCR products from the different bins of each replicate sorting equimolarly to a total mass of 2,600 ng and a volume of 100 µL (26 ng/µL of DNA) in 1.5 mL Eppendorf tubes. We then sent the pooled purified amplicons for adapter ligation and sequencing at Eurofins (NGSelect Amplicons on Illumina HiSeq), obtaining 15 million paired-end reads (2×150 bp) for all samples.

We quantified the total number of reads (t) required for sequencing based on several factors, including the number of bins (b), biological replicates (r), strains (s), genotypes (g), and the reads per genotype (p). Using these variables, we first calculated the total number $n_{b c}$ of needed barcodes as follows:

n_{b c} = b \times r \times s

For each of our three TFs, the values of these parameters are b=13 bins, r=3 biological replicates, and s=1 strains, yielding n_bc = 39. The total number of required sequence reads can be estimated through the equation:

t = r \times g \times p

where we aimed for P=30 paired-reads per genotype. The number g of expected genotypes per library is the same (65,536) for each TF. These values lead to a total of t=3 × 65,536×30 = 7,667,712 required paired-end reads for each library.

Data analysis

Filtering and preparing sequencing reads

We processed the sequencing data with a blend of custom python and awk scripts, complemented by established bioinformatics utilities. Firstly, we trimmed sequences by computationally excising Illumina adapters, followed by merging paired-end reads. We then organized the paired-end reads into distinct files, each tagged with a sequence barcode identifying the sequence bin from which the corresponding reads originated.

For the removal of Illumina adapter sequences from the paired-end reads, we employed Cutadapt (Martin, 2011) with the following parameters:

cutadapt -j 8 -e 0.1 --no-indels --overlap=8 --discard-untrimmed \
-a "^\$FWD...\$REV_RC;max_error_rate = 0.2;min_overlap = 6" \
-A "^\$REV...\$FWD_RC;max_error_rate = 0.2;min_overlap = 6" --pair-filter=any \
-o 'Sample_${sample}_trimmed_1.fastq.gz' \
-p 'Sample_${sample}_trimmed_2.fastq.gz' \
--max-ee=2 -l=114 \
$reads

Here, ‘ADAPTER_FWD’ and ‘ADAPTER_REV_RC’ are placeholders for the actual adapter sequences used, which are 5’-AGATCGGAAGAGCACACGTCTGAACTCCAGTCA-3’ for read 1, and 5’-AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGT-3’ for read 2.

Because our amplicons were short (153 bp) and each sequencing sample was divided into two files with overlapping single-end reads, merging the two files (single-end reads) was necessary, for which we used the FLASH software (Magoč and Salzberg, 2011). After the filtering and merging of reads, we demultiplexed paired-end reads with the following FLASH parameters:

flash -t $ncores $reads -O -m 60 M 140 -z -o 'Sample_${sample}_merged'

Next, we utilized the FastX toolkit (http://hannonlab.cshl.edu/fastx_toolkit/) (Hannon, 2010) to discern and preserve only those sequences that exceeded a high-quality threshold of Q=33. Subsequently, we refined our dataset by removing sequences that contained undesirable mutations or insertions/deletions (indels) within the region extending from the start of the core promoter to the end of the variable TFBS library. We performed this filtering step using a custom R (R Development Core Team, 2026) script. Thereafter, we used custom awk and R (R Development Core Team, 2026) scripts to convert the data into a table. Each row of this table contains data from one TFBS in the library, and each column contains the number of reads for this TFBS from one of the 13 bins into which we had sorted cells. We estimated the average regulation strength of each TFBS from this data.

Calculating regulation strengths

We calculated regulation strengths as previously described18. Briefly, in sort-seq experiments, sequences often appear in multiple fluorescence bins due to random mis-sorting (Peterman and Levine, 2016; Meyer et al., 2018) and the stochastic nature of gene expression (Elowitz et al., 2002). Following methods established in previous research (Vaishnav et al., 2022; de Boer et al., 2020b; Lagator et al., 2022), we determined the reporter expression level driven from each TFBS variant by calculating a weighted average of bins in which the variant occurred. This involved multiplying the frequency of each sequence (x_i) in a given bin i by a numerical value representing that bin (w_i = 1, 2, 3,…, 13), and then averaging these products over the total sequence count. Mathematically, this weighted average calculates as

e = \sum_{i = 1}^{n} \frac{(x_{i} * w_{i})}{\sum_{i = 1}^{n} x_{i}},

where e is the expression level driven by the sequence.

This approach yields a continuous spectrum of expression levels e within the range of 1–13. High expression values indicate robust GFP expression and, thus, weak binding of a TF to a TFBS variant. To facilitate interpretation, we inverted this scale so that higher values indicate stronger repression and lower GFP expression, that is we define the regulation strength b (strength of regulation) as follows:

b = e_{m a x} + 1 - e

where e_max = 13 is the maximal expression level (maximal bin value). We then normalized b by the maximal regulation strength among all TFBSs we studied for a given TF. The resulting normalized regulation strength scores b vary from 0 to 1, where low scores denote TFBS variants with weak TF binding (weak reporter repression), and high scores indicate variants with strong TF binding (strong reporter repression). A score of b=1 signifies the highest regulation strength (repression) calculated among all TFBSs for a given TF.

Combining data from triplicates

As a quality-filtering step, we eliminated TFBS variants that did not appear in all three replicates or that were represented by less than 30 reads in total (summing their read counts over the thirteen bins). Given that we have 13 bins and most sequences are typically found within an average of 3.4 bins, we wanted a minimum of 10 reads per bin to prevent misinterpretation of repression levels and also guide our threshold’s choice. This procedure reduced the fraction of TFBS variants for which we had sort-seq-based sequence data from 95%, 90%, and 93% of the total library size (65,536 sequences) for CRP, Fis, and IHF, respectively, to 49%, 66%, and 63%.

Subsequently, we calculated the regulation strengths b of the remaining TFBS variants (Appendix 1, Calculating regulation strengths). Then, we determined the coefficient of variation of regulation strengths across replicates for each TFBS variant, which estimate the consistency of measured regulation strengths among replicates. Notably, most TFBS variants exhibited a coefficient of variation (CV) below 0.5, which was the threshold we set for filtering sequences. Sequences with a CV above 0.5 were excluded from the analysis. This threshold is commonly used in transcriptional studies employing fluorescent reporters to account for transcriptional noise and measurement variability (Elowitz et al., 2002). Finally, we averaged regulation strengths for each TFBS variant across all replicates and normalized the resulting averages by the maximal observed regulation strength among all TFBSs of a given TF.

Frequency matrices and sequence logos

We generated frequency matrices of nucleotides that occur in TFBSs from each fluorescence bin by counting the frequency of each nucleotide at each variable position of the TFBS library. From this data, we generated heatmaps and DNA sequence logos representing the frequency matrices graphically. A sequence logo consists of a stack of the letters A, C, G, and T, at each position of a DNA sequence, where the relative size of each letter indicates its frequency in the sequence. The total height of the stack corresponds to the information content of that position, in bits (Stormo, 2000; Schneider and Stephens, 1990). A sequence logo is a graphical representation of the informational properties of a TFBS. When mutated, nucleotides with high information content are more likely to lead to a loss of binding (repression) than nucleotides with low information content.

Creation of genotype networks and determining network metrics

We used in-house Python script and the Python package igraph (Csárdi and Nepusz, 2006) to generate directed genotype networks. These are graphs in which TFBS variants are nodes (vertices, genotypes), and variants that differ in a single nucleotide are connected by an edge. Each node of this network is associated with the corresponding DNA sequence and the associated regulation strength. Each edge is directed, that is it corresponds to a binding-score-increasing mutation, and points from a TFBS variant with lower regulation strength to a neighbor with higher regulation strength. We extracted the largest weakly connected subgraph (the ‘giant component’ Aguilar-Rodríguez et al., 2017) of the network, and used it for all further analyses. This giant component comprises the vast majority (99%) of sequenced genotypes for all of our TF landscapes. We used in-house Python and R scripts for all network analyses described below.

Epistasis

Epistasis refers to non-additive interactions between two or more mutations. It can impose severe constraints on molecular evolution, because the mutations that are beneficial in one genetic background may be deleterious in another (Bank, 2022). Epistasis between two mutations can be classified as magnitude, simple sign, or reciprocal sign epistasis, depending on the sign (i.e. positive or negative) of the fitness effect of individual mutations and their combinations (Bank, 2022). In magnitude epistasis, the effect of a mutation on regulation strength varies depending on the genetic background but the sign of this effect (increasing or decreasing regulation strength) does not. Simple sign epistasis occurs if one single mutant has a lower regulation strength than both the wild type and the double mutant, while the other single mutant has a regulation strength that is intermediate to the wild type and double mutant. Reciprocal sign epistasis occurs when both mutations independently decrease regulation strength, but their combination increases regulation strength. The presence of reciprocal sign epistasis is a necessary condition for the existence of multiple peaks in an adaptive landscape (Saona et al., 2022; Poelwijk et al., 2011).

To determine the incidence of epistasis in our landscapes, we employed a method that involves identifying all ‘squares’ in a genotype network with the motifs function from the igraph library in R. Each square consists of a ‘wild-type’ sequence, a double-nucleotide mutant, and the corresponding two single mutants. We assessed epistasis for each square along a single axis by selecting the highest-regulation strength sequence as the double mutant. Our analysis placed each square into one of three categories: no sign epistasis, simple sign epistasis, and reciprocal sign epistasis. The no sign epistasis category included both magnitude epistasis and additivity (no epistasis) without differentiating between them, because neither affects the accessibility of landscape peaks through only binding-increasing mutations (Aguilar-Rodríguez et al., 2017). We determined the proportion of all squares that fell into each category (see Supplementary file 5).

Peaks

A peak is a genotype (TF binding site variant) whose neighbors all convey lower regulation strength than itself. Two peaks are connected if they are neighbors and convey the same regulation strength. We refer to the genotype with the highest regulation strength as the summit or global peak (Aguilar-Rodríguez et al., 2017; Khalid et al., 2016).

Genotype connectivity

To evaluate the sparsity of our networks and the likelihood of peak misclassification due to the incompleteness of our landscapes, we used a metric called ‘genotype connectivity’. This represents the number of immediate (1-mutant) neighbors of a given genotype for which our experiment produced repression strength data. A perfect experiment creating a complete landscape data set would provide repression data for each of our 4⁸ genotypes and all of its 24 immediate neighbors. In practice, however, data for some genotypes and some of their neighbors is unavailable. We quantified the ‘relative connectivity’ of any one genotype as the ratio of the actual number of adjacent genotypes for which we have repression data to the theoretical maximum of 24 neighbors. This measure helps us to address questions about the sparsity of our landscape and its implications. Specifically, it allows us to determine if our sample is biased, with some regions or genotypes being more connected than others. This is particularly important for assessing whether peaks and high peaks are less connected than non-peaks, which could potentially lead to peak misassignments.

Accessible paths

In the context of our genotype networks, we call a mutational path accessible if the regulation strength increases with each mutational step (Aguilar-Rodríguez et al., 2017; Khalid et al., 2016). We systematically enumerated all shortest accessible paths, distinguishing them by their path length, that is by the number of mutational steps in a path. Notably, there can be multiple alternative accessible paths with the same number of mutational steps from a single starting genotype to a peak.

Basins of attraction

The basin of attraction of a peak comprises all TFBS variants from which accessible paths to the peak exist. We refer to the basin’s size as the number of variants in the basin. We determined basin sizes by exhaustive enumeration.

Overlap between basins

The basins of attraction of different peaks may comprise overlapping sets of variants. To determine the overlap between two basins B₁ and B₂, we used the Jaccard index J (Jaccard, 1912; Papkou et al., 2023b), which is equal to the size of the intersection between two sets of variants divided by the size of their union:

J = \frac{B 1 \cap B 2}{B 1 \cup B 2}

Generating randomly shuffled landscapes

To generate uncorrelated random landscapes for our three TF landscapes, we performed a random shuffling of regulation strength values from our genotypes, followed by landscape analysis to identify peaks. To this end, we employed an ad-hoc Python script to randomly permute experimentally measured regulation strength values among all genotypes in each of our landscapes. This procedure preserves the distribution of regulation strengths in each landscape, but also creates an uncorrelated random landscape (Weinberger, 1990; Krug and Oros, 2023). Following random shuffling, we analyzed each shuffled landscape’s topography with the same methods we had applied to the original (unshuffled) landscape to identify its peaks. To ensure statistical robustness, we repeated the random shuffling and peak identification 10³ times. We then compared the number of peaks between the 10³ shuffled landscapes and the original (unshuffled) landscape.

Principal component analysis

In order to investigate which nucleotides in a TFBS are most important for changes in regulation strengths, we performed a principal component analysis for each landscape (Figure 2—figure supplements 10–12). To this end, we first one-hot-encoded our data. This encoding represents each categorical value (nucleotide at each position of a DNA string) as a binary vector of length 4. It thus converts an entire DNA string of length L into a 4xL binary matrix. We used this binary matrix to perform PCA with the R base function prcomp.

Simulated adaptive walks

We simulated the adaptive evolution of a population on our landscapes by performing two different types of random walks, ‘greedy’ adaptive random walks, and random walks based on Kimura’s model of fixation probabilities.

For all simulations, we assumed that only point mutations occur and that the time it takes for a point mutation to become fixed in a population is much shorter than the time it takes for a new mutation to appear that will eventually go to fixation. This scenario is also called the strong selection weak mutation (SSWM) scenario (Vaishnav et al., 2022; Bank et al., 2016; Orr, 2002; Gillespie, 1984). Under this scenario, evolving populations are monomorphic most of the time, that is all individuals have the same genotype. This scenario is realistic when the product of effective population size N and mutation rate μ is small (Nμ<1), which is the case for E. coli (N=1.8 × 10⁸, μ=2 × 10^–10; Lynch et al., 2016). It allows us to model adaptive evolution as an adaptive random walk in our landscapes. Our simulations also account for mutation bias, which means that different types of mutation have a different probability of occurring in a population. We use mutation biases that were experimentally determined for E. coli (Lee et al., 2012).

In a greedy adaptive random walk, starting from any one genotype, only the mutational neighbor that conveys the largest fitness advantage is fixed in a population. We initiated one greedy random walk from each non-peak genotype and terminated the walk once a fitness peak was reached, that is once every neighbor of the current genotype had lower fitness than the genotype itself. Because every greedy random walk is deterministic in the absence of neutral mutations, it is sufficient to initiate one greedy random walk per starting genotype.

Another well-established model for random walks allows for genetic drift. It uses fixation probabilities computed by Kimura, 1962; Kimura, 1983; Crow and Kimura, 2009 that is f_ij = (1 – e^-2s) / (1 – e^-2Ns), where f_ij is the probability of fixing mutation j in the background of genotype i, N is the effective population size, and s is the selection coefficient, that is the difference in fitness between genotypes i and j (Kimura, 1962; Kimura, 1983). For a given pair of genotypes, the only parameter of this model is the effective population size N, for which we explored values of N=10⁸, N=10⁵, and N=10². Generally, the smaller the population size is, the larger is the probability that neutral or deleterious mutations become fixed in a population.

For these ‘Kimura’ adaptive walks, we chose 15,000 random starting genotypes for each of the three values of N and simulated 1000 random walks for each of them. At each step, we randomly picked a mutation j, generated a random number in the interval [0, 1], and considered the mutation to become fixed if the random number fell within the interval [0, f_ij]. Computationally, we accelerated this process by precomputing all fixation probabilities for all genetic backgrounds, and then used the Python function numpy.choice to generate a random sample from the multinomial distribution of fixation probabilities at each step (Harris et al., 2020). Because of genetic drift, Kimura’s random walks do not necessarily terminate when they reach a fitness peak. For this reason, we simulated each such random walk for a maximum of 25 mutational steps.

Randomized landscape null model for peak accessibility

To evaluate whether the accessibility of strong regulatory peaks observed in the empirical landscapes exceeds random expectations, we constructed randomized ‘null’ landscapes and repeated evolutionary simulations on them under identical conditions as for the empirical landscapes.

Specifically, for each transcription factor landscape, we generated 10³ independently shuffled landscapes (see Appendix 1, Generating randomly shuffled landscapes) by permuting repression-strength values across genotypes while preserving the network of genotypes and its mutational neighborhood structure. This procedure maintains the number of genotypes, their connectivity, and the distribution of repression-strength values while removing correlations between genotype and phenotype. All other aspects of the analysis were kept identical to those used for the empirical landscapes.

For each randomized landscape, we simulated adaptive evolution using the same Kimura random-walk framework as for the empirical data. We assumed a population size of 10⁸ individuals, initiated adaptive walks from the same number of starting genotypes, and used identical stopping criteria, fixation probabilities, and mutational bias parameters as described in Appendix 1, Simulated adaptive walks (‘Simulated adaptive walks’).

For each randomized landscape, we computed the fraction of adaptive walks that reached a high regulatory peak. This yielded a null distribution of peak-accessibility values for each transcription factor. We tested the null hypothesis that the observed fraction from the corresponding empirical landscape is greater than the mean fraction derived from this null distribution using a one-sided Monte Carlo permutation test ( $N = 10^{3}$ randomized landscapes). For this test, we used Z-scores computed as the difference between the observed accessibility and the mean of the shuffled-landscape distribution, normalized by the standard deviation of that distribution.

Incorporating experimental uncertainty into adaptive walks

Fluorescence-based sort-seq measurements are subject to experimental uncertainty arising from finite fluorescence bin resolution, variability in fluorescence measurements, and sequencing noise. To account for this uncertainty when defining fitness relationships between genotypes and identifying peaks in the landscape, we incorporated empirically estimated noise into pairwise repression strength comparisons.

For each genotype, we computed repression strength ( $S$ ) as a continuous quantity derived from the distribution of sequencing reads across fluorescence bins (see Methods). We quantified experimental uncertainty in $S$ by a genotype-specific noise parameter ( $τ_{G}$ ), calculated as the standard deviation of repression strength across three biological replicates. This parameter defines an uncertainty interval within which repression strength values cannot be reliably distinguished.

We compared the repression strengths of pairs of genotypes while explicitly accounting for their associated experimental uncertainty. Let $S_{A}$ and $S_{B}$ denote the repression strengths of genotypes $A$ and $B$ , with corresponding uncertainty estimates $τ_{A}$ and $τ_{B}$ . We considered genotype $A$ to exhibit significantly stronger repression than genotype $B$ if

S_{A} - τ_{A} > S_{B} + τ_{B} .

In this case, we assigned a directed edge from $A$ to $B$ in the genotype network. Conversely, genotype $A$ was considered to exhibit significantly weaker repression than genotype $B$ if

S_{B} - τ_{B} > S_{A} + τ_{A},

in which case we assigned a directed edge from $B$ to $A$ . If neither condition was satisfied – that is if the uncertainty intervals of the two genotypes overlapped – we considered their repression strengths to be indistinguishable within experimental uncertainty. Such genotype pairs were treated as neutral and connected by two directed edges, one in each direction.

Applying these rules to all genotype pairs yielded the uncertainty-aware genotype network, denoted $G^{S, τ}$ . For comparison, we also constructed a noise-free genotype network, denoted $G^{S}$ , in which experimental uncertainty was ignored ( $τ_{G} = 0$ ) and pairwise comparisons were based solely on point estimates of repression strength $S$ . In both networks, peaks were defined as genotypes without outgoing edges.

Data availability

Sequencing data has been deposited in the NCBI database under the BioProject accession code: PRJNA1162449. The data generated in this study have been deposited in the Zenodo public repository and are accessible via the following DOI: https://doi.org/10.5281/zenodo.13838265.

The following data sets were generated

1. Westmann C
(2024) NCBI BioProject
ID PRJNA1162449. Sort-seq data from CRP, Fis and IHF.

https://www.ncbi.nlm.nih.gov/bioproject/?term=PRJNA1162449
1. Leander G
2. Andreas W
(2024) Zenodo
Scripts and Datasets for Manuscript: The adaptive landscapes of three global Escherichia coli transcriptional regulators.

https://doi.org/10.5281/zenodo.13838265

References

(2017) A thousand empirical adaptive landscapes and their navigability
Nature Ecology & Evolution 1:45.

https://doi.org/10.1038/s41559-016-0045
- PubMed
- Google Scholar
Book
1. Aguilar-Rodríguez J
2. Payne JL
(2021) Robustness and evolvability in transcriptional regulation
In: Crombach A, editors. Evolutionary Systems Biology. Springer. pp. 197–219.

https://doi.org/10.1007/978-3-030-71737-7_9
- Google Scholar
1. Aird D
2. Ross MG
3. Chen W-S
4. Danielsson M
5. Fennell T
6. Russ C
7. Jaffe DB
8. Nusbaum C
9. Gnirke A
(2011) Analyzing and minimizing PCR amplification bias in Illumina sequencing libraries
Genome Biology 12:1–14.

https://doi.org/10.1186/gb-2011-12-2-r18
- PubMed
- Google Scholar
1. Aiyar SE
2. McLeod SM
3. Ross W
4. Hirvonen CA
5. Thomas MS
6. Johnson RC
7. Gourse RL
(2002) Architecture of Fis-activated transcription complexes at the Escherichia coli rrnB P1 and rrnE P1 promoters
Journal of Molecular Biology 316:501–516.

https://doi.org/10.1006/jmbi.2001.5390
- PubMed
- Google Scholar
1. Ali Azam T
2. Iwata A
3. Nishimura A
4. Ueda S
5. Ishihama A
(1999) Growth phase-dependent variation in protein composition of the Escherichia coli nucleoid
Journal of Bacteriology 181:6361–6370.

https://doi.org/10.1128/JB.181.20.6361-6370.1999
- PubMed
- Google Scholar
1. Arfin SM
2. Long AD
3. Ito ET
4. Tolleri L
5. Riehle MM
6. Paegle ES
7. Hatfield GW
(2000) Global gene expression profiling in Escherichia coli K12. The effects of integration host factor
The Journal of Biological Chemistry 275:29672–29684.

https://doi.org/10.1074/jbc.M002247200
- PubMed
- Google Scholar
1. Baba T
2. Ara T
3. Hasegawa M
4. Takai Y
5. Okumura Y
6. Baba M
7. Datsenko KA
8. Tomita M
9. Wanner BL
10. Mori H
(2006) Construction of Escherichia coli K-12 in-frame, single-gene knockout mutants: the Keio collection
Molecular Systems Biology 2:2006.0008.

https://doi.org/10.1038/msb4100050
- PubMed
- Google Scholar
(2016) On the (un)predictability of a large intragenic fitness landscape
PNAS 113:14085–14090.

https://doi.org/10.1073/pnas.1612676113
- PubMed
- Google Scholar
1. Bank C
(2022) Epistasis and adaptation on fitness landscapes
Annual Review of Ecology, Evolution, and Systematics 53:457–479.

https://doi.org/10.1146/annurev-ecolsys-102320-112153
- Google Scholar
(2004) Regulation at complex bacterial promoters: how bacteria use different promoter organizations to produce different regulatory outcomes
Current Opinion in Microbiology 7:102–108.

https://doi.org/10.1016/j.mib.2004.02.011
- PubMed
- Google Scholar
(2019) Mapping DNA sequence to transcription factor binding energy in vivo
PLOS Computational Biology 15:e1006226.

https://doi.org/10.1371/journal.pcbi.1006226
- PubMed
- Google Scholar
(2016) Reproducibility of fluorescent expression from engineered biological constructs in E. coli
PLOS ONE 11:e0150182.

https://doi.org/10.1371/journal.pone.0150182
- Google Scholar
1. Beal J
(2017) Biochemical complexity drives log‐normal variation in genetic expression
Engineering Biology 1:55–60.

https://doi.org/10.1049/enb.2017.0004
- Google Scholar
(2018) A systematic and scalable approach for dissecting the molecular mechanisms of transcriptional regulation in bacteria
Biophysical Journal 114:151a.

https://doi.org/10.1016/j.bpj.2017.11.849
- Google Scholar
(2019) Genotype network intersections promote evolutionary innovation
PLOS Biology 17:e3000300.

https://doi.org/10.1371/journal.pbio.3000300
- PubMed
- Google Scholar
(2016) Divergent MLS1 promoters lie on a fitness plateau for gene expression
Molecular Biology and Evolution 33:1270–1279.

https://doi.org/10.1093/molbev/msw010
- PubMed
- Google Scholar
1. Bindels DS
2. Haarbosch L
3. van Weeren L
4. Postma M
5. Wiese KE
6. Mastop M
7. Aumonier S
8. Gotthard G
9. Royant A
10. Hink MA
11. Gadella TJ
(2017) mScarlet: a bright monomeric red fluorescent protein for cellular imaging
Nature Methods 14:53–56.

https://doi.org/10.1038/nmeth.4074
- PubMed
- Google Scholar
1. Bintu L
2. Buchler NE
3. Garcia HG
4. Gerland U
5. Hwa T
6. Kondev J
7. Phillips R
(2005) Transcriptional regulation by the numbers: models
Current Opinion in Genetics & Development 15:116–124.

https://doi.org/10.1016/j.gde.2005.02.007
- PubMed
- Google Scholar
1. Blanco C
2. Janzen E
3. Pressman A
4. Saha R
5. Chen IA
(2019) Molecular fitness landscapes from high-coverage sequence profiling
Annual Review of Biophysics 48:1–18.

https://doi.org/10.1146/annurev-biophys-052118-115333
- PubMed
- Google Scholar
(2008) Historical contingency and the evolution of a key innovation in an experimental population of Escherichia coli
PNAS 105:7899–7906.

https://doi.org/10.1073/pnas.0803151105
- PubMed
- Google Scholar
(2018) Contingency and determinism in evolution: Replaying life’s tape
Science 362:eaam5979.

https://doi.org/10.1126/science.aam5979
- PubMed
- Google Scholar
Preprint
1. Boer C
2. Sadeh R
3. Friedman N
4. Regev A
(2018) Deciphering cis-regulatory logic with 100 million random promoters
bioRxiv.

https://doi.org/10.1101/224907
- Google Scholar
1. Borirak O
2. Rolfe MD
3. de Koning LJ
4. Hoefsloot HCJ
5. Bekker M
6. Dekker HL
7. Roseboom W
8. Green J
9. de Koster CG
10. Hellingwerf KJ
(2015) Time-series analysis of the transcriptome and proteome of Escherichia coli upon glucose repression
Biochimica et Biophysica Acta 1854:1269–1279.

https://doi.org/10.1016/j.bbapap.2015.05.017
- PubMed
- Google Scholar
1. Browning DFDF
2. Busby SJWS
(2004) The regulation of bacterial transcription initiation
Nature Reviews. Microbiology 2:57–65.

https://doi.org/10.1038/nrmicro787
- PubMed
- Google Scholar
1. Browning DF
2. Busby SJW
(2016) Local and global regulation of transcription initiation in bacteria
Nature Reviews. Microbiology 14:638–650.

https://doi.org/10.1038/nrmicro.2016.103
- PubMed
- Google Scholar
(2019) Bacterial transcription factors: regulation by pick “N” mix
Journal of Molecular Biology 431:4067–4077.

https://doi.org/10.1016/j.jmb.2019.04.011
- Google Scholar
(2017) Reducing DNA context dependence in bacterial promoters
PLOS ONE 12:e0176013.

https://doi.org/10.1371/journal.pone.0176013
- PubMed
- Google Scholar
(2003) Transcription regulation and environmental adaptation in bacteria
Trends in Microbiology 11:248–253.

https://doi.org/10.1016/s0966-842x(03)00103-3
- PubMed
- Google Scholar
1. Cases I
2. de Lorenzo V
(2005) Promoters in the environment: transcriptional regulation in its natural context
Nature Reviews. Microbiology 3:105–118.

https://doi.org/10.1038/nrmicro1084
- PubMed
- Google Scholar
Preprint
(2025) The fitness landscape of the E. coli Lac operator is highly rugged in two different environments
bioRxiv.

https://doi.org/10.1101/2025.07.23.666252
- Google Scholar
1. Chen YJ
2. Liu P
3. Nielsen AAK
4. Brophy JAN
5. Clancy K
6. Peterson T
7. Voigt CA
(2013) Characterization of 582 natural and synthetic terminators and quantification of their design constraints
Nature Methods 10:659–664.

https://doi.org/10.1038/nmeth.2515
- PubMed
- Google Scholar
1. Chen YJ
2. Takahashi CN
3. Organick L
4. Bee C
5. Ang SD
6. Weiss P
7. Peck B
8. Seelig G
9. Ceze L
10. Strauss K
(2020) Quantifying molecular bias in DNA data storage
Nature Communications 11:3264.

https://doi.org/10.1038/s41467-020-16958-3
- PubMed
- Google Scholar
1. Conrad M
(1990) The geometry of evolution
Biosystems 24:61–81.

https://doi.org/10.1016/0303-2647(90)90030-5
- Google Scholar
(2016) The soft touch: low-affinity transcription factor binding sites in development and evolution
Current Topics in Developmental Biology 117:455–469.

https://doi.org/10.1016/bs.ctdb.2015.11.018
- PubMed
- Google Scholar
Book
1. Crow J
2. Kimura M
(2009)
An Introduction to Population Genetics Theory

The Blackburn Press.
- Google Scholar
1. Csárdi G
2. Nepusz T
(2006)
The igraph software package for complex network research

InterJournal, Complex Systems 1695:1–9.
- Google Scholar
1. Cuperus J
2. Groves B
3. Kuchina A
4. Rosenberg AB
5. Jojic N
6. Fields S
7. Seelig G
(2017) Deep learning of the regulatory grammar of yeast 5’ untranslated regions from 500,000 random sequences
Genome Res 163:1–10.

https://doi.org/10.1101/137547
- Google Scholar
1. de Boer CG
2. Ray JP
3. Hacohen N
4. Regev A
(2020a) MAUDE: inferring expression changes in sorting-based CRISPR screens
Genome Biology 21:134.

https://doi.org/10.1186/s13059-020-02046-8
- PubMed
- Google Scholar
1. de Boer CG
2. Vaishnav ED
3. Sadeh R
4. Abeyta EL
5. Friedman N
6. Regev A
(2020b) Deciphering eukaryotic gene-regulatory logic with 100 million random promoters
Nature Biotechnology 38:56–65.

https://doi.org/10.1038/s41587-019-0315-8
- PubMed
- Google Scholar
(2011) The causes of epistasis
Proceedings. Biological Sciences 278:3617–3624.

https://doi.org/10.1098/rspb.2011.1537
- PubMed
- Google Scholar
1. de Visser J
2. Krug J
(2014) Empirical fitness landscapes and the predictability of evolution
Nature Reviews. Genetics 15:480–490.

https://doi.org/10.1038/nrg3744
- PubMed
- Google Scholar
(2018) The utility of fitness landscapes and big data for predicting evolution
Heredity 121:401–405.

https://doi.org/10.1038/s41437-018-0128-4
- PubMed
- Google Scholar
1. D’haeseleer P
(2006) What are DNA sequence motifs?
Nature Biotechnology 24:423–425.

https://doi.org/10.1038/nbt0406-423
- PubMed
- Google Scholar
(2023) Global epistasis on fitness landscapes
Philosophical Transactions of the Royal Society B 378:20220053.

https://doi.org/10.1098/rstb.2022.0053
- PubMed
- Google Scholar
1. Dillon SC
2. Dorman CJ
(2010) Bacterial nucleoid-associated proteins, nucleoid structure and gene expression
Nature Reviews. Microbiology 8:185–195.

https://doi.org/10.1038/nrmicro2261
- PubMed
- Google Scholar
1. Dorman CJ
(2013) Genome architecture and global gene regulation in bacteria: making progress towards a unified model?
Nature Reviews. Microbiology 11:349–355.

https://doi.org/10.1038/nrmicro3007
- PubMed
- Google Scholar
1. Dorman CJ
2. Dorman MJ
(2016) DNA supercoiling is a fundamental regulatory principle in the control of bacterial gene expression
Biophysical Reviews 8:89–100.

https://doi.org/10.1007/s12551-016-0238-2
- PubMed
- Google Scholar
Book
(2018) The evolution of gene regulatory mechanisms in bacteria
In: Dorman CJ, editors. Molecular Mechanisms of Microbial Evolution. Springer. pp. 125–152.

https://doi.org/10.1007/978-3-319-69078-0_6
- Google Scholar
(2020) When is a transcription factor a NAP?
Current Opinion in Microbiology 55:26–33.

https://doi.org/10.1016/j.mib.2020.01.019
- PubMed
- Google Scholar
1. Draghi JA
2. Ogbunugafor CB
(2023) Exploring the expanse between theoretical questions and experimental approaches in the modern study of evolvability
Journal of Experimental Zoology. Part B, Molecular and Developmental Evolution 340:8–17.

https://doi.org/10.1002/jez.b.23134
- PubMed
- Google Scholar
(2017) Fitness effects of cis-regulatory variants in the Saccharomyces cerevisiae TDH3 promoter
Molecular Biology and Evolution 34:2908–2912.

https://doi.org/10.1093/molbev/msx224
- PubMed
- Google Scholar
(2019) Bacterial genome architecture shapes global transcriptional regulation by DNA supercoiling
Nucleic Acids Research 47:5648–5657.

https://doi.org/10.1093/nar/gkz300
- PubMed
- Google Scholar
(2002) Stochastic gene expression in a single cell
Science 297:1183–1186.

https://doi.org/10.1126/science.1070919
- PubMed
- Google Scholar
(2017) Application of sorting and next generation sequencing to study 5΄-UTR influence on translation efficiency in Escherichia coli
Nucleic Acids Research 45:3487–3502.

https://doi.org/10.1093/nar/gkw1141
- Google Scholar
1. Feng H
2. Li F
3. Wang T
4. Xing X-H
5. Zeng A-P
6. Zhang C
(2023) Deep-learning-assisted Sort-Seq enables high-throughput profiling of gene expression characteristics with high precision
Science Advances 9:eadg5296.

https://doi.org/10.1126/sciadv.adg5296
- PubMed
- Google Scholar
(2019) Evolution in the light of fitness landscape theory
Trends in Ecology & Evolution 34:69–82.

https://doi.org/10.1016/j.tree.2018.10.009
- PubMed
- Google Scholar
(2017) Evolution of new regulatory functions on biophysically realistic fitness landscapes
Nature Communications 8:216.

https://doi.org/10.1038/s41467-017-00238-8
- PubMed
- Google Scholar
Preprint
1. Fuqua T
2. Wagner A
(2023) Mobile DNA Is Replete with Hotspots for the de Novo Emergence of Gene Regulation
bioRxiv.

https://doi.org/10.1101/2023.10.22.563463
- Google Scholar
(2011) Comparison and calibration of different reporters for quantitative analysis of gene expression
Biophysical Journal 101:535–544.

https://doi.org/10.1016/j.bpj.2011.06.026
- PubMed
- Google Scholar
(2023) The non-deterministic genotype-phenotype map of RNA secondary structure
Journal of the Royal Society, Interface 20:20230132.

https://doi.org/10.1098/rsif.2023.0132
- PubMed
- Google Scholar
1. Gawade P
2. Gunjal G
3. Sharma A
4. Ghosh P
(2020) Reconstruction of transcriptional regulatory networks of Fis and H-NS in Escherichia coli from genome-wide data analysis
Genomics 112:1264–1272.

https://doi.org/10.1016/j.ygeno.2019.07.013
- PubMed
- Google Scholar
1. Gerland U
2. Moroz JD
3. Hwa T
(2002) Physical constraints and functional characteristics of transcription factor–DNA interaction
PNAS 99:12015–12020.

https://doi.org/10.1073/pnas.192693599
- PubMed
- Google Scholar
1. Gibson DG
2. Young L
3. Chuang RY
4. Venter JC
5. Hutchison CA
6. Smith HO
(2009) Enzymatic assembly of DNA molecules up to several hundred kilobases
Nature Methods 6:343–345.

https://doi.org/10.1038/nmeth.1318
- PubMed
- Google Scholar
1. Gillespie JH
(1984) Molecular evolution over the mutational landscape
Evolution; International Journal of Organic Evolution 38:1116–1129.

https://doi.org/10.1111/j.1558-5646.1984.tb00380.x
- PubMed
- Google Scholar
1. Gilliot PA
2. Gorochowski TE
(2023) Effective design and inference for cell sorting and sequencing based massively parallel reporter assays
Bioinformatics 39:btad277.

https://doi.org/10.1093/bioinformatics/btad277
- PubMed
- Google Scholar
(2009) From sequence to dynamics: the effects of transcription factor and polymerase concentration changes on activated and repressed promoters
BMC Molecular Biology 10:92.

https://doi.org/10.1186/1471-2199-10-92
- Google Scholar
(1992)
DNA sequence determinants for binding of the Escherichia coli catabolite gene activator protein

The Journal of Biological Chemistry 267:14713–14720.
- PubMed
- Google Scholar
(2014) Biophysical fitness landscapes for transcription factor binding sites
PLOS Computational Biology 10:e1003683.

https://doi.org/10.1371/journal.pcbi.1003683
- PubMed
- Google Scholar
Software
1. Hannon GJ
(2010) FASTX-toolkit
Cold Spring Harbor Laboratory.

http://hannonlab.cshl.edu/fastx_toolkit
1. Harris CR
2. Millman KJ
3. van der Walt SJ
4. Gommers R
5. Virtanen P
6. Cournapeau D
7. Wieser E
8. Taylor J
9. Berg S
10. Smith NJ
11. Kern R
12. Picus M
13. Hoyer S
14. van Kerkwijk MH
15. Brett M
16. Haldane A
17. Del Río JF
18. Wiebe M
19. Peterson P
20. Gérard-Marchant P
21. Sheppard K
22. Reddy T
23. Weckesser W
24. Abbasi H
25. Gohlke C
26. Oliphant TE
(2020) Array programming with NumPy
Nature 585:357–362.

https://doi.org/10.1038/s41586-020-2649-2
- PubMed
- Google Scholar
(2025) The structure of an ancient genotype–phenotype map shaped the functional evolution of a protein family
Nature Ecology & Evolution 9:1656–1669.

https://doi.org/10.1038/s41559-025-02777-6
- PubMed
- Google Scholar
(2021) Restoring global gene regulation through experimental evolution uncovers a NAP (Nucleoid-Associated Protein)-like behavior of Crp/Cap
mBio 12:e0202821.

https://doi.org/10.1128/mBio.02028-21
- PubMed
- Google Scholar
Book
1. Ho S
(2013) Characteristics of IHF Binding to Holliday Junction DNA
Middletown, CT: Wesleyan University.

https://doi.org/10.14418/wes01.1.932
- Google Scholar
1. Horton JS
2. Taylor TB
(2023) Mutation bias and adaptation in bacteria
Microbiology 169:001404.

https://doi.org/10.1099/mic.0.001404
- PubMed
- Google Scholar
1. Huffman JL
2. Brennan RG
(2002) Prokaryotic transcription regulators: more than just the helix-turn-helix motif
Current Opinion in Structural Biology 12:98–106.

https://doi.org/10.1016/s0959-440x(02)00295-6
- PubMed
- Google Scholar
1. Huo YX
2. Zhang YT
3. Xiao Y
4. Zhang X
5. Buck M
6. Kolb A
7. Wang YP
(2009) IHF-binding sites inhibit DNA loop formation and transcription initiation
Nucleic Acids Research 37:3878–3886.

https://doi.org/10.1093/nar/gkp258
- PubMed
- Google Scholar
1. Ibarra IL
2. Hollmann NM
3. Klaus B
4. Augsten S
5. Velten B
6. Hennig J
7. Zaugg JB
(2020) Mechanistic insights into transcription factor cooperativity and its impact on protein-phenotype interactions
Nature Communications 11:124.

https://doi.org/10.1038/s41467-019-13888-7
- PubMed
- Google Scholar
1. Ishihama A
(2010) Prokaryotic genome regulation: multifactor promoters, multitarget regulators and hierarchic networks
FEMS Microbiology Reviews 34:628–645.

https://doi.org/10.1111/j.1574-6976.2010.00227.x
- PubMed
- Google Scholar
(2004) Stochastic tunnels in evolutionary dynamics
Genetics 166:1571–1579.

https://doi.org/10.1534/genetics.166.3.1571
- PubMed
- Google Scholar
1. Jaccard P
(1912) The distribution of the flora of the alpine zone
New Phytologist 11:37–50.

https://doi.org/10.1111/j.1469-8137.1912.tb05611.x
- Google Scholar
1. Jahn M
2. Vorpahl C
3. Hübschmann T
4. Harms H
5. Müller S
(2016) Copy number variability of expression plasmids determined by cell sorting and Droplet Digital PCR
Microbial Cell Factories 15:211.

https://doi.org/10.1186/s12934-016-0610-8
- PubMed
- Google Scholar
(2014) Promoter architecture dictates cell-to-cell variability in gene expression
Science 346:1533–1536.

https://doi.org/10.1126/science.1255301
- PubMed
- Google Scholar
(2011) Direct and indirect effects of H-NS and Fis on global gene expression control in Escherichia coli
Nucleic Acids Research 39:2073–2091.

https://doi.org/10.1093/nar/gkq934
- PubMed
- Google Scholar
1. Kauffman S
2. Levin S
(1987) Towards a general theory of adaptive walks on rugged landscapes
Journal of Theoretical Biology 128:11–45.

https://doi.org/10.1016/s0022-5193(87)80029-2
- PubMed
- Google Scholar
1. Kebschull JM
2. Zador AM
(2015) Sources of PCR-induced distortions in high-throughput sequencing data sets
Nucleic Acids Research 43:e143.

https://doi.org/10.1093/nar/gkv717
- PubMed
- Google Scholar
1. Kelly JR
2. Rubin AJ
3. Davis JH
4. Ajo-Franklin CM
5. Cumbers J
6. Czar MJ
7. de Mora K
8. Glieberman AL
9. Monie DD
10. Endy D
(2009) Measuring the activity of BioBrick promoters using an in vivo reference standard
Journal of Biological Engineering 3:4.

https://doi.org/10.1186/1754-1611-3-4
- PubMed
- Google Scholar
(2016) Genonets server-a web server for the construction, analysis and visualization of genotype networks
Nucleic Acids Research 44:W70–W76.

https://doi.org/10.1093/nar/gkw313
- PubMed
- Google Scholar
1. Khankal R
2. Chin JW
3. Ghosh D
4. Cirino PC
(2009) Transcriptional effects of CRP* expression in Escherichia coli
Journal of Biological Engineering 3:1–14.

https://doi.org/10.1186/1754-1611-3-13
- PubMed
- Google Scholar
1. Kim D
2. Seo SW
3. Gao Y
4. Nam H
5. Guzman GI
6. Cho BK
7. Palsson BO
(2018) Systems assessment of transcriptional regulation on central carbon metabolism by Cra and CRP
Nucleic Acids Research 46:2901–2917.

https://doi.org/10.1093/nar/gky069
- PubMed
- Google Scholar
1. Kimura M
(1962) On the probability of fixation of mutant genes in a population
Genetics 47:713–719.

https://doi.org/10.1093/genetics/47.6.713
- PubMed
- Google Scholar
Book
1. Kimura M
(1983) The Neutral Theory of Molecular Evolution
Cambridge University Press.

https://doi.org/10.1017/CBO9780511623486
- Google Scholar
1. Kinney JB
2. Murugan A
3. Callan CG
4. Cox EC
(2010) Using deep sequencing to characterize the biophysical mechanism of a transcriptional regulatory sequence
PNAS 107:9158–9163.

https://doi.org/10.1073/pnas.1004290107
- PubMed
- Google Scholar
1. Kinney JB
2. McCandlish DM
(2019) Massively parallel assays and quantitative sequence-function relationships
Annual Review of Genomics and Human Genetics 20:99–127.

https://doi.org/10.1146/annurev-genom-083118-014845
- PubMed
- Google Scholar
1. Kolb A
2. Spassky A
3. Chapon C
4. Blazy B
5. Buc H
(1983) On the different binding affinities of CRP at the lac, gal and malT promoter regions
Nucleic Acids Research 11:7833–7852.

https://doi.org/10.1093/nar/11.22.7833
- PubMed
- Google Scholar
Preprint
1. Krug J
2. Oros D
(2023) Evolutionary accessibility of random and structured fitness landscapes
arXiv.

https://arxiv.org/abs/2311.17432v2
- Google Scholar
(2019) Evolutionary potential of cis-regulatory mutations to cause rapid changes in transcription factor binding
Genome Biology and Evolution 11:406–414.

https://doi.org/10.1093/gbe/evy269
- PubMed
- Google Scholar
1. Kvitek DJ
2. Sherlock G
(2011) Reciprocal sign epistasis between frequently experimentally evolved adaptive mutations causes a rugged fitness landscape
PLOS Genetics 7:e1002056.

https://doi.org/10.1371/journal.pgen.1002056
- PubMed
- Google Scholar
1. Lagator M
2. Igler C
3. Moreno AB
4. Guet CC
5. Bollback JP
(2016) Epistatic interactions in the arabinose cis-regulatory element
Molecular Biology and Evolution 33:761–769.

https://doi.org/10.1093/molbev/msv269
- PubMed
- Google Scholar
1. Lagator M
2. Paixão T
3. Barton NH
4. Bollback JP
5. Guet CC
(2017) On the mechanistic nature of epistasis in a canonical cis-regulatory element
eLife 6:e25192.

https://doi.org/10.7554/eLife.25192
- PubMed
- Google Scholar
(2022) Predicting bacterial promoter function and evolution from random sequences
eLife 11:e64543.

https://doi.org/10.7554/eLife.64543
- PubMed
- Google Scholar
(2015) From structure to function of bacterial chromosomes: Evolutionary perspectives and ideas for new experiments
FEBS Lett 589:2996–3004.

https://doi.org/10.1016/j.febslet.2015.07.002
- PubMed
- Google Scholar
(2017) Predicting evolution
Nature Ecology & Evolution 1:0077.

https://doi.org/10.1038/s41559-017-0077
- Google Scholar
1. Le DD
2. Shimko TC
3. Aditham AK
4. Keys AM
5. Longwell SA
6. Orenstein Y
7. Fordyce PM
(2018) Comprehensive, high-resolution binding energy landscapes reveal context dependencies of transcription factor binding
PNAS 115:E3702–E3711.

https://doi.org/10.1073/pnas.1715888115
- PubMed
- Google Scholar
1. Lee H
2. Popodi E
3. Tang H
4. Foster PL
(2012) Rate and molecular spectrum of spontaneous mutations in the bacterium Escherichia coli as determined by whole-genome sequencing
PNAS 109:E2774–E2783.

https://doi.org/10.1073/pnas.1210309109
- PubMed
- Google Scholar
1. Lenski RE
(2017) Convergence and divergence in a long-term experiment with bacteria
The American Naturalist 190:S57–S68.

https://doi.org/10.1086/691209
- PubMed
- Google Scholar
1. Li J
2. Amado A
3. Bank C
(2024) Rapid adaptation of recombining populations on tunable fitness landscapes
Molecular Ecology 33:e16900.

https://doi.org/10.1111/mec.16900
- PubMed
- Google Scholar
1. Lind PA
2. Andersson DI
(2008) Whole-genome mutational biases in bacteria
PNAS 105:17878–17883.

https://doi.org/10.1073/pnas.0804445105
- PubMed
- Google Scholar
(2015) Experimental evolution reveals hidden diversity in evolutionary pathways
eLife 4:e07074.

https://doi.org/10.7554/eLife.07074
- PubMed
- Google Scholar
1. Lou C
2. Stanton B
3. Chen YJ
4. Munsky B
5. Voigt CA
(2012) Ribozyme-based insulator parts buffer synthetic circuits from genetic context
Nature Biotechnology 30:1137–1142.

https://doi.org/10.1038/nbt.2401
- PubMed
- Google Scholar
1. Louis AA
(2016) Contingency, convergence and hyper-astronomical numbers in biological evolution
Studies in History and Philosophy of Biological and Biomedical Sciences 58:107–116.

https://doi.org/10.1016/j.shpsc.2015.12.014
- PubMed
- Google Scholar
(2008) The role of DNA-binding specificity in the evolution of bacterial regulatory networks
Journal of Molecular Biology 379:627–643.

https://doi.org/10.1016/j.jmb.2008.04.008
- PubMed
- Google Scholar
1. Lutz R
2. Bujard H
(1997) Independent and tight regulation of transcriptional units in Escherichia coli via the LacR/O, the TetR/O and AraC/I1-I2 regulatory elements
Nucleic Acids Research 25:1203–1210.

https://doi.org/10.1093/nar/25.6.1203
- PubMed
- Google Scholar
1. Lynch M
2. Ackerman MS
3. Gout JF
4. Long H
5. Sung W
6. Thomas WK
7. Foster PL
(2016) Genetic drift, selection and the evolution of the mutation rate
Nature Reviews. Genetics 17:704–714.

https://doi.org/10.1038/nrg.2016.104
- PubMed
- Google Scholar
(2010) Diminishing returns from beneficial mutations and pervasive epistasis shape the fitness landscape for rifampicin resistance in Pseudomonas aeruginosa
Genetics 186:1345–1354.

https://doi.org/10.1534/genetics.110.123083
- PubMed
- Google Scholar
1. Madan Babu M
2. Teichmann SA
(2003) Functional determinants of transcription factors in Escherichia coli: protein families and binding sites
Trends in Genetics 19:75–79.

https://doi.org/10.1016/S0168-9525(02)00039-2
- PubMed
- Google Scholar
1. Maerkl S
2. Quake SR
(2007) A systems approach to measuring the binding energy landscapes of transcription factors
Science 315:233–237.

https://doi.org/10.1126/science.1131007
- PubMed
- Google Scholar
1. Magoč T
2. Salzberg SL
(2011) FLASH: fast length adjustment of short reads to improve genome assemblies
Bioinformatics 27:2957–2963.

https://doi.org/10.1093/bioinformatics/btr507
- PubMed
- Google Scholar
1. Martin M
(2011) Cutadapt removes adapter sequences from high-throughput sequencing reads
EMBnet.Journal 17:10.

https://doi.org/10.14806/ej.17.1.200
- Google Scholar
1. Martínez-Antonio A
2. Collado-Vides J
(2003) Identifying global regulators in transcriptional regulatory networks in bacteria
Current Opinion in Microbiology 6:482–489.

https://doi.org/10.1016/j.mib.2003.09.002
- PubMed
- Google Scholar
1. Mathelier A
2. Xin B
3. Chiu TP
4. Yang L
5. Rohs R
6. Wasserman WW
(2016) DNA shape features improve transcription factor binding site predictions in vivo
Cell Systems 3:278–286.

https://doi.org/10.1016/j.cels.2016.07.001
- PubMed
- Google Scholar
(2004) The evolution of genetic regulatory systems in bacteria
Nature Reviews. Genetics 5:169–178.

https://doi.org/10.1038/nrg1292
- PubMed
- Google Scholar
1. Melissa MJ
2. Good BH
3. Fisher DS
4. Desai MM
(2022) Population genetics of polymorphism and divergence in rapidly evolving populations
Genetics 221:iyac053.

https://doi.org/10.1093/genetics/iyac053
- PubMed
- Google Scholar
(2018) Chromosomal organization of transcription: in a nutshell
Current Genetics 64:555–565.

https://doi.org/10.1007/s00294-017-0785-5
- PubMed
- Google Scholar
(2019) Escherichia coli “Marionette” strains with 12 highly optimized small-molecule sensors
Nature Chemical Biology 15:196–204.

https://doi.org/10.1038/s41589-018-0168-3
- PubMed
- Google Scholar
(2020) Unraveling the complex interplay of Fis and IHF through synthetic promoter engineering
Frontiers in Bioengineering and Biotechnology 8:510.

https://doi.org/10.3389/fbioe.2020.00510
- PubMed
- Google Scholar
1. Nghe P
2. de Vos MGJ
3. Kingma E
4. Kogenaru M
5. Poelwijk FJ
6. Laan L
7. Tans SJ
(2020) Predicting evolution using regulatory architecture
Annual Review of Biophysics 49:181–197.

https://doi.org/10.1146/annurev-biophys-070317-032939
- PubMed
- Google Scholar
1. Nonoyama T
2. Chiba S
(2019) Phenotypic determinism and contingency in the evolution of hypothetical tree-like organisms
PLOS ONE 14:e0211671.

https://doi.org/10.1371/journal.pone.0211671
- PubMed
- Google Scholar
(2013) Binding of nucleoid-associated protein fis to DNA is regulated by DNA breathing dynamics
PLOS Computational Biology 9:e1002881.

https://doi.org/10.1371/journal.pcbi.1002881
- PubMed
- Google Scholar
1. Obolski U
2. Ram Y
3. Hadany L
(2018) Key issues review: evolution on rugged adaptive landscapes
Reports on Progress in Physics 81:012602.

https://doi.org/10.1088/1361-6633/aa94d4
- PubMed
- Google Scholar
(2005) Non-additivity in protein-DNA binding
Bioinformatics 21:2254–2263.

https://doi.org/10.1093/bioinformatics/bti361
- PubMed
- Google Scholar
1. Oren Y
2. Smith MB
3. Johns NI
4. Kaplan Zeevi M
5. Biran D
6. Ron EZ
7. Corander J
8. Wang HH
9. Alm EJ
10. Pupko T
(2014) Transfer of noncoding DNA drives regulatory rewiring in bacteria
PNAS 111:16112–16117.

https://doi.org/10.1073/pnas.1413272111
- PubMed
- Google Scholar
(2023) Evolution and synthetic biology
Current Opinion in Microbiology 76:102394.

https://doi.org/10.1016/j.mib.2023.102394
- PubMed
- Google Scholar
1. Orr HA
(2002) The population genetics of adaptation: the adaptation of dna sequences
Evolution; International Journal of Organic Evolution 56:1317–1330.

https://doi.org/10.1111/j.0014-3820.2002.tb01446.x
- PubMed
- Google Scholar
1. Orr HA
(2003) A minimum on the mean number of steps taken in adaptive walks
Journal of Theoretical Biology 220:241–247.

https://doi.org/10.1006/jtbi.2003.3161
- PubMed
- Google Scholar
1. Otwinowski J
2. Nemenman I
(2013) Genotype to phenotype mapping and the fitness landscape of the E. coli lac promoter
PLOS ONE 8:e61570.

https://doi.org/10.1371/journal.pone.0061570
- PubMed
- Google Scholar
(2022) Global pleiotropic effects in adaptively evolved Escherichia coli lacking CRP reveal molecular mechanisms that define the growth physiology
Open Biology 12:210206.

https://doi.org/10.1098/rsob.210206
- PubMed
- Google Scholar
1. Palmer AC
2. Toprak E
3. Baym M
4. Kim S
5. Veres A
6. Bershtein S
7. Kishony R
(2015) Delayed commitment to evolutionary fate in antibiotic resistance fitness landscapes
Nature Communications 6:1–8.

https://doi.org/10.1038/ncomms8385
- PubMed
- Google Scholar
(2023a) A rugged yet easily navigable fitness landscape
Science 382:eadh3860.

https://doi.org/10.1126/science.adh3860
- PubMed
- Google Scholar
Preprint
(2023b) A rugged yet easily navigable fitness landscape of antibiotic resistance
bioRxiv.

https://doi.org/10.1101/2023.02.27.530293
- Google Scholar
(2016) Greedy adaptive walks on a correlated fitness landscape
Journal of Theoretical Biology 397:89–102.

https://doi.org/10.1016/j.jtbi.2016.02.035
- PubMed
- Google Scholar
1. Payne JL
2. Wagner A
(2014) The robustness and evolvability of transcription factor binding sites
Science 343:875–877.

https://doi.org/10.1126/science.1249046
- PubMed
- Google Scholar
(2006) Engineering and characterization of a superfolder green fluorescent protein
Nature Biotechnology 24:79–88.

https://doi.org/10.1038/nbt1172
- PubMed
- Google Scholar
(2018) Abundance, diversity and domain architecture variability in prokaryotic DNA-binding transcription factors
PLOS ONE 13:e0195332.

https://doi.org/10.1371/journal.pone.0195332
- PubMed
- Google Scholar
1. Peterman N
2. Levine E
(2016) Sort-seq under the hood: implications of design choices on large-scale characterization of sequence-function relations
BMC Genomics 17:206.

https://doi.org/10.1186/s12864-016-2533-5
- PubMed
- Google Scholar
(2007) Empirical fitness landscapes reveal accessible evolutionary paths
Nature 445:383–386.

https://doi.org/10.1038/nature05451
- PubMed
- Google Scholar
(2011) Reciprocal sign epistasis is a necessary condition for multi-peaked fitness landscapes
Journal of Theoretical Biology 272:141–144.

https://doi.org/10.1016/j.jtbi.2010.12.015
- PubMed
- Google Scholar
(2008) Horizontal gene transfer and the evolution of transcriptional regulation in Escherichia coli
Genome Biology 9:R4.

https://doi.org/10.1186/gb-2008-9-1-r4
- PubMed
- Google Scholar
(2012) Genomic analysis of DNA binding and gene regulation by homologous nucleoid-associated proteins IHF and HU in Escherichia coli K12
Nucleic Acids Research 40:3524–3537.

https://doi.org/10.1093/nar/gkr1236
- PubMed
- Google Scholar
Software
1. R Development Core Team
(2026) R: a language and environment for statistical computing
R Foundation for Statistical Computing, Vienna, Austria.

https://www.R-project.org/
1. Rest JS
2. Morales CM
3. Waldron JB
4. Opulente DA
5. Fisher J
6. Moon S
7. Bullaughey K
8. Carey LB
9. Dedousis D
(2013) Nonlinear fitness consequences of variation in expression level of a eukaryotic gene
Molecular Biology and Evolution 30:448–456.

https://doi.org/10.1093/molbev/mss248
- PubMed
- Google Scholar
1. Rudge TJ
2. Brown JR
3. Federici F
4. Dalchau N
5. Phillips A
6. Ajioka JW
7. Haseloff J
(2016) Characterization of intrinsic properties of promoters
ACS Synthetic Biology 5:89–98.

https://doi.org/10.1021/acssynbio.5b00116
- PubMed
- Google Scholar
(2014) The influence of promoter architectures and regulatory motifs on gene expression in Escherichia coli
PLOS ONE 9:e114347.

https://doi.org/10.1371/journal.pone.0114347
- PubMed
- Google Scholar
(2024) RegulonDB v12.0: a comprehensive resource of transcriptional regulation in E. coli K-12
Nucleic Acids Research 52:D255–D264.

https://doi.org/10.1093/nar/gkad1072
- PubMed
- Google Scholar
(2015) DNA methylation in bacteria: from the methyl group to the methylome
Current Opinion in Microbiology 25:9–16.

https://doi.org/10.1016/j.mib.2015.03.004
- PubMed
- Google Scholar
(2022) Relation between the number of peaks and the number of reciprocal sign epistatic interactions
Bulletin of Mathematical Biology 84:74.

https://doi.org/10.1007/s11538-022-01029-z
- PubMed
- Google Scholar
1. Sarkisyan KS
2. Bolotin DA
3. Meer MV
4. Usmanova DR
5. Mishin AS
6. Sharonov GV
7. Ivankov DN
8. Bozhanova NG
9. Baranov MS
10. Soylemez O
11. Bogatyreva NS
12. Vlasov PK
13. Egorov ES
14. Logacheva MD
15. Kondrashov AS
16. Chudakov DM
17. Putintseva EV
18. Mamedov IZ
19. Tawfik DS
20. Lukyanov KA
21. Kondrashov FA
(2016) Local fitness landscape of the green fluorescent protein
Nature 533:397–401.

https://doi.org/10.1038/nature17995
- PubMed
- Google Scholar
1. Schneider TD
2. Stephens RM
(1990) Sequence logos: a new way to display consensus sequences
Nucleic Acids Research 18:6097–6100.

https://doi.org/10.1093/nar/18.20.6097
- PubMed
- Google Scholar
1. Scholz SA
2. Diao R
3. Wolfe MB
4. Fivenson EM
5. Lin XN
6. Freddolino PL
(2019) High-resolution mapping of the Escherichia coli chromosome reveals positions of high and low transcription
Cell Systems 8:212–225.

https://doi.org/10.1016/j.cels.2019.02.004
- PubMed
- Google Scholar
Book
1. Schuster P
(2002) A testable genotype-phenotype map: modeling evolution of RNA molecules
In: Lässig M, Valleriani A, editors. Biological Evolution and Statistical Physics. Springer. pp. 55–81.

https://doi.org/10.1007/3-540-45692-9_4
- Google Scholar
1. Schweizer G
2. Wagner A
(2021) Both binding strength and evolutionary accessibility affect the population frequency of transcription factor binding sequences in Arabidopsis thaliana
Genome Biology and Evolution 13:evab273.

https://doi.org/10.1093/gbe/evab273
- PubMed
- Google Scholar
(2006) Transcriptional regulatory networks in bacteria: from input signals to output responses
Current Opinion in Microbiology 9:511–519.

https://doi.org/10.1016/j.mib.2006.08.007
- PubMed
- Google Scholar
(2008) Functional characterization of the Escherichia coli Fis-DNA binding sequence
Journal of Molecular Biology 376:771–785.

https://doi.org/10.1016/j.jmb.2007.11.101
- PubMed
- Google Scholar
1. Sharon E
2. Kalma Y
3. Sharp A
4. Raveh-Sadka T
5. Levo M
6. Zeevi D
7. Keren L
8. Yakhini Z
9. Weinberger A
10. Segal E
(2012) Inferring gene regulatory logic from high-throughput measurements of thousands of systematically designed promoters
Nature Biotechnology 30:521–530.

https://doi.org/10.1038/nbt.2205
- PubMed
- Google Scholar
1. Sharon E
2. van Dijk D
3. Kalma Y
4. Keren L
5. Manor O
6. Yakhini Z
7. Segal E
(2014) Probing the effect of promoters on noise in gene expression using thousands of designed sequences
Genome Research 24:1698–1706.

https://doi.org/10.1101/gr.168773.113
- PubMed
- Google Scholar
1. Shen-Orr SS
2. Milo R
3. Mangan S
4. Alon U
(2002) Network motifs in the transcriptional regulation network of Escherichia coli
Nature Genetics 31:64–68.

https://doi.org/10.1038/ng881
- PubMed
- Google Scholar
(2023) Transcription factor expression levels and environmental signals constrain transcription factor innovation
Microbiology 169:001378.

https://doi.org/10.1099/mic.0.001378
- PubMed
- Google Scholar
(2010) The fitness landscapes of cis-acting binding sites in different promoter and environmental contexts
PLOS Genetics 6:e1001042.

https://doi.org/10.1371/journal.pgen.1001042
- PubMed
- Google Scholar
1. Signor SA
2. Nuzhdin SV
(2018) The evolution of gene expression in cis and trans
Trends in Genetics 34:532–544.

https://doi.org/10.1016/j.tig.2018.03.007
- PubMed
- Google Scholar
(2012) Gene order and chromosome dynamics coordinate spatiotemporal gene expression during the bacterial growth cycle
PNAS 109:E42–E50.

https://doi.org/10.1073/pnas.1108229109
- PubMed
- Google Scholar
1. Song S
2. Zhang J
(2021) Unbiased inference of the fitness landscape ruggedness from imprecise fitness estimates
Evolution; International Journal of Organic Evolution 75:2658–2671.

https://doi.org/10.1111/evo.14363
- PubMed
- Google Scholar
1. Spor A
2. Kvitek DJ
3. Nidelet T
4. Martin J
5. Legrand J
6. Dillmann C
7. Bourgais A
8. de Vienne D
9. Sherlock G
10. Sicard D
(2014) Phenotypic and genotypic convergences are influenced by historical contingency and environment in yeast
Evolution; International Journal of Organic Evolution 68:772–790.

https://doi.org/10.1111/evo.12302
- PubMed
- Google Scholar
1. Srivastava M
2. Payne JL
(2022) On the incongruence of genotype-phenotype and fitness landscapes
PLOS Computational Biology 18:e1010524.

https://doi.org/10.1371/journal.pcbi.1010524
- PubMed
- Google Scholar
1. Steinberg B
2. Ostermeier M
(2016) Environmental changes bridge evolutionary valleys
Science Advances 2:e1500921.

https://doi.org/10.1126/sciadv.1500921
- PubMed
- Google Scholar
(2012) Why transcription factor binding sites are ten nucleotides long
Genetics 192:973–985.

https://doi.org/10.1534/genetics.112.143370
- PubMed
- Google Scholar
(2022) Complex fitness landscape shapes variation in a hyperpolymorphic species
eLife 11:e76073.

https://doi.org/10.7554/eLife.76073
- PubMed
- Google Scholar
1. Stormo GD
(2000) DNA binding sites: representation and discovery
Bioinformatics 16:16–23.

https://doi.org/10.1093/bioinformatics/16.1.16
- PubMed
- Google Scholar
1. Stormo GD
2. Zhao Y
(2010) Determining the specificity of protein-DNA interactions
Nature Reviews. Genetics 11:751–760.

https://doi.org/10.1038/nrg2845
- PubMed
- Google Scholar
1. Struhl K
(1999) Fundamentally different logic of gene regulation in eukaryotes and prokaryotes
Cell 98:1–4.

https://doi.org/10.1016/S0092-8674(00)80599-1
- PubMed
- Google Scholar
1. Taylor LR
2. Provine WB
(1987) Sewall wright and evolutionary biology
The Journal of Animal Ecology 56:720.

https://doi.org/10.2307/5082
- Google Scholar
1. Taylor CF
2. Higgs PG
(2000) A population genetics model for multiple quantitative traits exhibiting pleiotropy and epistasis
Journal of Theoretical Biology 203:419–437.

https://doi.org/10.1006/jtbi.2000.1094
- PubMed
- Google Scholar
1. Trippe BL
2. Huang B
3. DeBenedictis EA
4. Coventry B
5. Bhattacharya N
6. Yang KK
7. Baker D
8. Crawford L
(2022) Randomized gates eliminate bias in sort‐seq assays
Protein Science 31:e4401.

https://doi.org/10.1002/pro.4401
- Google Scholar
1. Urtecho G
2. Tripp AD
3. Insigne KD
4. Kim H
5. Kosuri S
(2019) Systematic dissection of sequence elements controlling σ70 promoters using a genomically encoded multiplexed reporter assay in Escherichia coli
Biochemistry 58:1539–1551.

https://doi.org/10.1021/acs.biochem.7b01069
- PubMed
- Google Scholar
1. Vaishnav ED
2. de Boer CG
3. Molinet J
4. Yassour M
5. Fan L
6. Adiconis X
7. Thompson DA
8. Levin JZ
9. Cubillos FA
10. Regev A
(2022) The evolution, evolvability and engineering of gene regulatory DNA
Nature 603:455–463.

https://doi.org/10.1038/s41586-022-04506-6
- PubMed
- Google Scholar
1. Vermeij GJ
(2006) Historical contingency and the purported uniqueness of evolutionary innovations
PNAS 103:1804–1809.

https://doi.org/10.1073/pnas.0508724103
- PubMed
- Google Scholar
1. Visweswariah SS
2. Busby SJW
(2015) Evolution of bacterial transcription factors: how proteins take on new tasks, but do not always stop doing the old ones
Trends in Microbiology 23:463–467.

https://doi.org/10.1016/j.tim.2015.04.009
- PubMed
- Google Scholar
1. Wagner A
(2022) AI predicts the effectiveness and evolution of gene promoter sequences
Nature 603:399–400.

https://doi.org/10.1038/d41586-022-00384-0
- PubMed
- Google Scholar
1. Warren DJ
(2011) Preparation of highly efficient electrocompetent Escherichia coli using glycerol/mannitol density step centrifugation
Analytical Biochemistry 413:206–207.

https://doi.org/10.1016/j.ab.2011.02.036
- PubMed
- Google Scholar
1. Weeks R
2. Ostermeier M
(2023) Fitness and functional landscapes of the E. coli RNase III gene rnc
Molecular Biology and Evolution 40:msad047.

https://doi.org/10.1093/molbev/msad047
- PubMed
- Google Scholar
1. Weinberger E
(1990) Correlated and uncorrelated fitness landscapes and how to tell the difference
Biological Cybernetics 63:325–336.

https://doi.org/10.1007/BF00202749
- Google Scholar
(2014) The transcription factor titration effect dictates level of gene expression
Biophysical Journal 106:489A.

https://doi.org/10.1016/j.bpj.2013.11.4483
- Google Scholar
(2005) Perspective: sign epistasis and genetic constraint on evolutionary trajectories
Evolution; International Journal of Organic Evolution 59:1165–1174.

https://doi.org/10.1111/j.0014-3820.2005.tb01768.x
- PubMed
- Google Scholar
Preprint
(2024a) Entangled adaptive landscapes facilitate the evolution of gene regulation by exaptation
bioRxiv.

https://doi.org/10.1101/2024.11.10.620926
- Google Scholar
(2024b) The highly rugged yet navigable regulatory landscape of the bacterial transcription factor TetR
Nature Communications 15:10745.

https://doi.org/10.1038/s41467-024-54723-y
- PubMed
- Google Scholar
Book
1. Wickham H
(2016) Ggplot2: Elegant Graphics for Data Analysis
Springer-Verlag.

https://doi.org/10.1007/978-0-387-98141-3
- Google Scholar
(2004) Evolutionary changes in cis and trans gene regulation
Nature 430:85–88.

https://doi.org/10.1038/nature02698
- PubMed
- Google Scholar
(2006) The diversity challenge in directed protein evolution
Combinatorial Chemistry & High Throughput Screening 9:271–288.

https://doi.org/10.2174/138620706776843192
- PubMed
- Google Scholar
1. Wray GA
(2007) The evolutionary significance of cis-regulatory mutations
Nature Reviews. Genetics 8:206–216.

https://doi.org/10.1038/nrg2063
- PubMed
- Google Scholar
1. Wright S
(1931) Evolution in Mendelian Populations
Genetics 16:97–159.

https://doi.org/10.1093/genetics/16.2.97
- PubMed
- Google Scholar
Conference
1. Wright S
(1932)
The roles of mutation, inbreeding, crossbreeding and selection in evolution

Proc of the 6th International Congress of Genetics Preprint.
- Google Scholar
1. Wu NC
2. Dai L
3. Olson CA
4. Lloyd-Smith JO
5. Sun R
(2016) Adaptation in protein fitness landscapes is facilitated by indirect paths
eLife 5:e16965.

https://doi.org/10.7554/eLife.16965
- PubMed
- Google Scholar
1. Xie VC
2. Pu J
3. Metzger BPH
4. Thornton JW
5. Dickinson BC
(2021) Contingency and chance erase necessity in the experimental evolution of ancestral proteins
eLife 10:e67336.

https://doi.org/10.7554/eLife.67336
- PubMed
- Google Scholar
1. Yi X
2. Dean AM
(2019) Adaptive landscapes in the age of synthetic biology
Molecular Biology and Evolution 36:890–907.

https://doi.org/10.1093/molbev/msz004
- PubMed
- Google Scholar
1. Yona AH
2. Alm EJ
3. Gore J
(2018) Random sequences rapidly evolve into de novo promoters
Nature Communications 9:1530.

https://doi.org/10.1038/s41467-018-04026-w
- PubMed
- Google Scholar
1. Zhao Y
2. Ruan S
3. Pandey M
4. Stormo GD
(2012) Improved models for transcription factor binding site identification using nonindependent interactions
Genetics 191:781–790.

https://doi.org/10.1534/genetics.112.138685
- PubMed
- Google Scholar
1. Zhou T
2. Shen N
3. Yang L
4. Abe N
5. Horton J
6. Mann RS
7. Bussemaker HJ
8. Gordân R
9. Rohs R
(2015) Quantitative modeling of transcription factor binding specificities using DNA shape
PNAS 112:4654–4659.

https://doi.org/10.1073/pnas.1422023112
- PubMed
- Google Scholar

Article and author information

Author details

Cauã Antunes Westmann
1. Department of Evolutionary Biology and Environmental Studies, University of Zurich, Zürich, Switzerland
2. Swiss Institute of Bioinformatics, Quartier Sorge-Batiment Genopode, Lausanne, Switzerland
Contribution
Conceptualization, Data curation, Software, Formal analysis, Validation, Investigation, Visualization, Methodology, Writing – original draft, Writing – review and editing

For correspondence
caua.westmann@ieu.uzh.ch

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0003-1572-2279
Leander Goldbach
1. Department of Evolutionary Biology and Environmental Studies, University of Zurich, Zürich, Switzerland
2. Swiss Institute of Bioinformatics, Quartier Sorge-Batiment Genopode, Lausanne, Switzerland
Contribution
Software, Formal analysis

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0002-3816-1833
Andreas Wagner
1. Department of Evolutionary Biology and Environmental Studies, University of Zurich, Zürich, Switzerland
2. Swiss Institute of Bioinformatics, Quartier Sorge-Batiment Genopode, Lausanne, Switzerland
3. The Santa Fe Institute, Santa Fe, United States
Contribution
Conceptualization, Resources, Supervision, Funding acquisition, Investigation, Writing – original draft, Writing – review and editing

For correspondence
andreas.wagner@ieu.uzh.ch

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0003-4299-3840

Funding

Swiss National Science Foundation (310030_208174)

Cauã Antunes Westmann
Leander Goldbach

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

We extend our thanks to the UZH University Priority Research Program in Evolutionary Biology, the UZH flow cytometry facility, and the Functional Genomics Center Zurich for their technical support. We are especially thankful to Andrei Papkou for his guidance with computational analysis and for engaging in theoretical discussions.

Version history

Sent for peer review: November 11, 2024
Preprint posted: November 13, 2024
Reviewed Preprint version 1: February 24, 2025
Reviewed Preprint version 2: April 9, 2026
Version of Record published: July 21, 2026

Cite all versions

You can cite all versions using the DOI https://doi.org/10.7554/eLife.103774. This DOI represents all versions, and will always resolve to the latest one.

Copyright

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.