
==== Front
Mol Biol Evol
Mol Biol Evol
molbev
Molecular Biology and Evolution
0737-4038
1537-1719
Oxford University Press UK

39219319
10.1093/molbev/msae185
msae185
Discoveries
AcademicSubjects/SCI01130
AcademicSubjects/SCI01180
Promoters Constrain Evolution of Expression Levels of Essential Genes in Escherichia coli
https://orcid.org/0000-0002-7182-1637
Tsuru Saburo Universal Biology Institute, Graduate School of Science, The University of Tokyo, Bunkyo-ku, Tokyo 113-0033, Japan

https://orcid.org/0009-0002-9039-4658
Hatanaka Naoki Universal Biology Institute, Graduate School of Science, The University of Tokyo, Bunkyo-ku, Tokyo 113-0033, Japan

https://orcid.org/0000-0003-3554-4975
Furusawa Chikara Universal Biology Institute, Graduate School of Science, The University of Tokyo, Bunkyo-ku, Tokyo 113-0033, Japan
Department of Physics, Graduate School of Science, The University of Tokyo, Bunkyo-ku, Tokyo 113-0033, Japan
Center for Biosystems Dynamics Research (BDR), RIKEN, Suita, Osaka 565-0874, Japan

Xia Xuhua Associate Editor
Saburo Tsuru and Naoki Hatanaka are equally contributed.

Corresponding authors: Emails: saburotsuru@gmail.com; chikara.furusawa@riken.jp.
9 2024
02 9 2024
02 9 2024
41 9 msae18522 5 2024
31 7 2024
28 8 2024
17 9 2024
© The Author(s) 2024. Published by Oxford University Press on behalf of Society for Molecular Biology and Evolution.
2024
https://creativecommons.org/licenses/by-nc/4.0/ This is an Open Access article distributed under the terms of the Creative Commons Attribution-NonCommercial License (https://creativecommons.org/licenses/by-nc/4.0/), which permits non-commercial re-use, distribution, and reproduction in any medium, provided the original work is properly cited. For commercial re-use, please contact reprints@oup.com for reprints and translation rights for reprints. All other permissions can be obtained through our RightsLink service via the Permissions link on the article page on our site—for further information please contact journals.permissions@oup.com.

Abstract

Variability in expression levels in response to random genomic mutations varies among genes, influencing both the facilitation and constraint of phenotypic evolution in organisms. Despite its importance, both the underlying mechanisms and evolutionary origins of this variability remain largely unknown due to the mixed contributions of cis- and trans-acting elements. To address this issue, we focused on the mutational variability of cis-acting elements, that is, promoter regions, in Escherichia coli. Random mutations were introduced into the natural and synthetic promoters to generate mutant promoter libraries. By comparing the variance in promoter activity of these mutant libraries, we found no significant difference in mutational variability in promoter activity between promoter groups, suggesting the absence of a signature of natural selection for mutational robustness. In contrast, the promoters controlling essential genes exhibited a remarkable bias in mutational variability, with mutants displaying higher activities than the wild types being relatively rare compared to those with lower activities. Our evolutionary simulation on a rugged fitness landscape provided a rationale for this vulnerability. These findings suggest that past selection created nonuniform mutational variability in promoters biased toward lower activities of random mutants, which now constrains the future evolution of downstream essential genes toward higher expression levels.

evolutionary constraint
gene expression
mutant library
promoter
mutational robustness
Japan Society for the Promotion of Science 10.13039/501100001691 18H02427 22H05403 24K21985 17H06389 22K21344 Japan Science and Technology Agency 10.13039/501100002241
==== Body
pmcIntroduction

The rate and direction of phenotypic evolution largely depend on phenotypic variability, the tendency to vary phenotypically (Lynch and Walsh 1998; Walker 2007). It is widely recognized that organisms often exhibit nonuniform phenotypic variability, with certain variants arising more frequently than others in response to different genetic and environmental perturbations, even in the absence of natural selection (Darwin 1859; Waddington 1957; Smith et al. 1985; Noble et al. 2019). Therefore, it is crucial to elucidate the mechanisms underlying the bias in phenotypic variability to understand why phenotypic evolution proceeds as it does (Uller et al. 2018). Transcriptional variability of genes across different perturbations plays a significant role in phenotypic evolution because the evolution of many traits in organisms largely depends on the evolution of expression levels and their regulation (Wittkopp et al. 2004; Zheng et al. 2011). Interestingly, previous studies have identified differences in transcriptional variability from gene to gene in response to various perturbations (Denver et al. 2005; Rifkin et al. 2005; Landry et al. 2007; Tsuru and Furusawa 2024), where genes related to cellular growth and maintenance tend to exhibit lower transcriptional variability against random genetic perturbations to the genome in yeast (Landry et al. 2007) and bacteria (Tsuru and Furusawa 2024). These studies, supported by theoretical evidence (Kaneko 2009; Draghi and Whitlock 2012; Furusawa and Kaneko 2018), suggest that the transcriptional variability of genes is not a random occurrence devoid of biological significance but rather a product of natural selection.

What molecular mechanisms account for gene-to-gene differences in transcriptional variability? In principle, both cis- and trans-acting elements can influence transcriptional variability (Wittkopp et al. 2004; Tirosh et al. 2006; Landry et al. 2007; Payne and Wagner 2015; Tsuru and Furusawa 2024). Previous studies (Denver et al. 2005; Rifkin et al. 2005; Landry et al. 2007; McGuigan et al. 2014; Tsuru and Furusawa 2024) explored the transcriptional variability against different genetic perturbations using mutation accumulation (MA) experiments (Halligan and Keightley 2009), where independent mutants were generated by accumulating genomic mutations through genetic drift, and their transcriptome profiles were obtained under identical environmental conditions. Although MA experiments provide a powerful method to quantify transcriptional variability in the absence of strong natural selection, the indiscriminate accumulation of mutations on cis- and trans-acting elements (Tsuru et al. 2015) complicates the identification of the molecular mechanism underlying transcriptional variability. Accordingly, to overcome this difficulty, it is necessary to consider the contributions of cis- and trans-acting elements separately. Although such approaches have been extensively employed in yeast (Hornung et al. 2012; Metzger et al. 2015; Duveau et al. 2021), it remains largely unknown which elements practically contribute to gene-to-gene differences in transcriptional variability associated with functions of gene products.

Why have genes related to cellular growth and maintenance acquired low mutational variability? The evolutionary origin of mutational variability is debatable (de Visser et al. 2003). Mutational variability could result from direct natural selection for mutational robustness (Waddington 1942) or could be a byproduct of selection for robustness against environmental perturbations (Meiklejohn and Hartl 2002; de Visser et al. 2003). Several biological traits, such as metabolic flux (Ho and Zhang 2016) and cellular morphology (Ho and Zhang 2014), show signs of adaptive evolution through direct selection for mutational robustness, where traits highly related to cellular growth fitness are likely to acquire lower mutational variability through direct selection for mutational robustness. On the other hand, other studies have identified congruence in trait variability between environmental and genetic perturbations, such as in RNA structure (Szollosi and Derenyi 2009) and gene expression levels (Landry et al. 2007; Tsuru and Furusawa 2024), supporting the relevance of the byproduct of selection for environmental robustness. As an alternative scenario, mutational variability may spontaneously emerge with trait evolution, regardless of selection for robustness (Siegal and Bergman 2002; de Visser et al. 2003). To date, the most appropriate scenario for low transcriptional variability of genes related to cellular growth and maintenance remains largely unknown.

To address these questions, we focused on the mutational variability governed by cis-acting elements, that is, promoter regions, of genes in Escherichia coli. We created a mutant library for each promoter through random mutagenesis and analyzed the changes in promoter activity from the wild type to mutant using flow cytometry. To explore the relevant evolutionary scenarios for the mutational variability of promoters, we compared promoters controlling essential genes for cellular growth with those controlling nonessential genes. These groups were expected to experience different natural selection pressures on their transcriptional variability in the past. In addition, random sequences were used to construct synthetic promoters that experienced no selection for mutational robustness. Interestingly, when comparing the variance in promoter activity of the mutant libraries, these three promoter groups exhibited similar mutational variability in promoter activity, indicating the absence of signatures of natural selection for mutational robustness. In contrast, the promoters for essential genes showed a remarkable vulnerability to mutations compared to the other promoter groups, where variants with higher activity than the wild type arose less frequently than variants with lower activity. Our evolutionary simulation on a rugged fitness landscape explains the higher vulnerability of the promoters for essential genes. These results suggest that past selection created a bias in mutational variability in promoters, which now constrains the future evolution of downstream essential genes toward higher expression levels.

Results

Higher Expression Levels and Lower Transcriptional Variability of Essential Genes in E. coli

First, we analyzed the transcriptional variability of essential genes in E. coli in response to various genomic mutations (Fig. 1a). Previously (Tsuru and Furusawa 2024), we obtained transcriptome profiles of six independent mutants that had accumulated genomic mutations through genetic drift via the MA experiment (Tsuru et al. 2015), termed the Mut data set. Using this data set, we quantified the mean and standard deviation of the expression levels for each of thousands of genes (Fig. 1b). To compensate for the nonmonotonic mean dependency of the standard deviations, we calculated the vertical distance from the smoothed spline of the running median of the standard deviations, termed DMmut. The DMmut values reflect the mean transcriptional variability in response to different genetic perturbations (Tsuru and Furusawa 2024). We categorized thousands of genes into essential and nonessential genes for cell growth, following the definition of gene essentiality by Goodall et al. (2018). We confirmed that essential genes showed higher expression levels than nonessential genes across the different mutants (Fig. 1c). We also found that essential genes had lower DMmut than nonessential genes (Fig. 1d).

Fig. 1. Transcriptional variability of essential genes in E. coli. a to d) Transcriptional profiles of six independent mutants of E. coli cultured under identical environmental conditions, termed the Mut data set. The mutants were created through the MA experiment to accumulate random genomic mutations by genetic drift. e to i) Transcriptional profiles of a single strain, MG1655, cultured under 76 different environmental conditions, termed the Env data set. a and e) Schematics of the data sets. b and f) Relationship between the means and standard deviations of mRNA expression levels for thousands of genes. The magenta and cyan lines represent running medians and their smoothed splines, respectively. Spearman's rank correlation coefficients (R) and P-values are shown. c and g) Mean expression levels of essential and nonessential genes. d and h) Transcriptional variability, represented as either DMmut or DMenv, of essential and nonessential genes. The DM values were calculated as the vertical distance of the standard deviations from the smoothed splines of the running median of standard deviations in b) and f). The lower and upper edges of the boxes in c), d), g), and h) represent the first (q1) and third (q3) quartiles, respectively. The horizontal lines in the boxes represent the medians (m). The whiskers from the boxes extend to the most extreme observed values inside inner fences, m ± 1.5 (q3 to q1). i) Relationship between mean expression levels and the shorter distance from oriC on the E. coli chromosome. Each dot represents a different gene.

Next, we explored the transcriptional variability of essential genes in response to different environmental perturbations (Fig. 1e and f) using a public data set (Sastry et al. 2019; Rychel et al. 2021) of the transcriptome profiles of a single strain of E. coli, MG1655, cultured under 76 different environmental conditions, termed the Env data set. Previously (Tsuru and Furusawa 2024), we quantified the mean-controlled standard deviations, termed DMenv, using a method similar to DMmut. The DMenv reflects transcriptional variability in response to different environmental perturbations. We confirmed that the essential genes exhibited higher expression levels (Fig. 1g) and lower variability (Fig. 1h) than the nonessential genes across different environmental conditions. Taken together, essential genes exhibited higher expression levels and lower transcriptional variability in response to different perturbations, which is similar to the characteristics of essential genes in yeast (Tirosh et al. 2006).

In E. coli, genes with many chromosomal copies or genes located at near the origin of replication (oriC) tend to have higher expression levels (Ying et al. 2014). However, essential genes shown in Fig. 1 are present as only single copies on the chromosome. In addition, these essential genes are widely spread on the chromosome and exhibited higher expression levels than neighboring nonessential genes over the chromosome (Fig. 1i). Therefore, gene dosage did not seem to explain the higher expression levels of essential genes. To explore whether the higher expression levels of essential genes are caused by chromosomal position or by higher-order chromosomal structures, we analyzed promoter-mediated gene expression from a plasmid based on the experimental data set obtained by Silander et al. (2012). Using flow cytometry, Silander et al. quantified promoter activity of plasmid copies of different promoter regions in E. coli, where green fluorescent protein (GFP) was transcriptionally fused to plasmid copies of promoters. Based on fluorescence intensity from GFP, we compared promoter activity between promoters controlling essential genes (termed essential promoters) and those controlling nonessential genes only (termed nonessential promoters) (supplementary fig. S1a, Supplementary Material online). We found that essential promoters tended to have higher activities than nonessential promoters in spite of plasmid copies. This result highlighted the importance of promoter regions in accounting for higher expression levels of essential genes independent of chromosomal positions or structures.

Like eukaryotes, bacteria often employ postreplicative DNA methylation for the epigenetic control of gene expression. In E. coli, DNA adenine methyltransferase plays a major role in the epigenetic control and the target GATC motifs are enriched in promoter regions (Oshima et al. 2002). Despite that adenine methylation at GATC sites in promoters can inhibit transcription for some genes (Casadesus and Low 2006), we found that differences in promoter activity between essential and nonessential promoters were almost independent of the presence or absence of GATC motifs in the promoter regions (supplementary fig. S1b, Supplementary Material online). These results implied that the higher expression levels of essential genes reflected higher promoter activities of essential promoters without known epigenetic controls.

Construction of Mutant Libraries of E. coli Promoters to Measure Mutational Effects

To explore the molecular mechanisms and signs of natural selection underlying the lower transcriptional variability of essential genes, we focused on the promoter regions of genes with different essentialities. In particular, we explored how promoter activity varies in response to new mutations in the promoter regions and how this variability relates to the essentiality of the downstream genes. To this end, we used 56 natural promoters from the E. coli promoter collection (Zaslaver et al. 2006). Importantly, to ensure an unbiased comparison of changes in promoter activities between different essentialities and to avoid measurement of promoter activity below the lower detection limit in flow cytometry, we avoided the use of natural promoters with very low activity that were slightly enriched in the promoters of nonessential genes. To achieve this, we selected natural promoters with moderate to high activity before mutagenesis. This selection was based on the known promoter activities of the collection as measured by Silander et al. (2012). Using the RegulonDB database (Tierrafria et al. 2022), we retrieved information regarding known transcriptional units regulated by these promoters. Based on the presence of essential genes in the transcription units, we categorized natural promoters into nonessential and essential promoters, termed Nes and Ess, respectively (Fig. 2a). Nonessential promoters control transcriptional units comprising only nonessential genes, whereas essential promoters regulate those that include at least a single essential gene. In addition to these natural promoters, we created synthetic promoters, termed Syn, as controls that experienced no or minimal selection for mutational robustness. To this end, we placed random sequences (∼200 bp) upstream of the gfp gene on a low-copy-number plasmid, pUA66-mCherry (Fig. 2b; supplementary fig. S2, Supplementary Material online). The library of mixed transformants containing this plasmid was analyzed using flow cytometry. The top 7% of cells showing higher expression levels of GFP (red shaded area in Fig. 2b) were sorted, and 16 unique clones were subsequently isolated. Importantly, these synthetic promoters experienced selection for their expression levels only once, and no longer experienced further mutagenesis or selection for expression level. Accordingly, they can be regarded as a control group that is almost free from selection for mutational robustness. Consequently, in the upstream of the GFP gene on the plasmid, we constructed three promoter groups, Nes, Ess, and Syn, comprising different sequences with low homology (supplementary fig. S3 and note S1, Supplementary Material online). The green fluorescence from GFP was used as a measure of promoter activity (Zaslaver et al. 2006; Silander et al. 2012). Despite the diversity in sequence, we confirmed that the mean activities of the three promoter groups were comparable to each other, as designed (Fig. 2c). We also confirmed that the obtained promoter activities were highly consistent with those measured by Silander et al. (2012) (Spearman's R = 0.94, P < 0.05; supplementary fig. S4a, Supplementary Material online). Compared to the Env data set, the activities of the natural promoters were positively correlated with known mRNA expression levels of the genes downstream of the chromosomal promoters in the wild-type strain (Spearman's R = 0.64; supplementary fig. S4b, Supplementary Material online), supporting that the wild-type promoter activities monitored by our method considerably reflect the activities of chromosomal promoters and their impact on the expression levels of downstream genes.

Fig. 2. Experimental design to measure expression levels of mutant promoter library. a) Natural promoters of E. coli were categorized into nonessential and essential promoters, termed Nes and Ess, respectively. Nonessential promoters control transcriptional units comprising only nonessential genes, while essential promoters regulate transcriptional units containing at least one essential gene. Gene essentiality was based on the requirement for cell growth, as defined by Goodall et al. (2018). The promoter regions were obtained from the E. coli promoter collection (Zaslaver et al. 2006). b) Synthetic promoters, termed Syn, were derived from a random oligo DNA pool. Transformants harboring a plasmid on which the random sequences were placed upstream of the gfp gene were analyzed using flow cytometry. The 16 unique transformants, satisfying a defined green fluorescence (shaded area), were isolated. The obtained synthetic promoters were named SP1–16. c) Natural and synthetic promoters were placed upstream of gfp on a plasmid. Transformants harboring each of these plasmids were analyzed using flow cytometry. Mean green fluorescence of each population is shown as small dots. The large dots and error bars represent the means and standard deviations within the groups, respectively. The number of promoters is indicated inside the panel. The P-value corresponds to ANOVA. d) Schematic showing construction of mutant library from a given promoter region. The promoter region was subjected to random mutagenesis by error-prone PCR. A mutant library was created through transformation to yield more than 100,000 genetically different mutants in a population. Unmutagenized plasmid was also transformed to create wild-type clones. e) Representative examples of log10-transformed green fluorescence intensity distributions of wild type and mutant library for the promoter region of the metG gene in E. coli (NPmetG). The means (vertical lines) and standard deviations (horizontal arrows) were calculated in each population.

For each of these promoters, we introduced mutations by error-prone PCR at a low frequency (6.3 mutations/kb on average; Fig. 2d) and integrated them into pUA66-mCherry to create a mutant library. Each mutant library had a genetic diversity of 105 or more genotypes. The promoter activity of each mutant library was measured by flow cytometry based on the green fluorescence derived from downstream gfp and was compared with the green fluorescence of the corresponding wild type (Fig. 2e). The introduction of mutations was limited to the promoter regions of the plasmid, while the chromosomal copies of the natural promoters were intact to avoid any confounding effects. The green fluorescence intensity was log10-transformed after subtracting the background (supplementary fig. S5, Supplementary Material online), and the mean and standard deviation of the GFP distribution were calculated (Fig. 2e; supplementary fig. S6, Supplementary Material online).

Mutational Variability in Promoter Activity was Similar between Different Promoter Groups

We confirmed that all mutant libraries showed larger variations than the corresponding wild-type clones, as designed (Fig. 3a), supporting that the difference in variations of GFP distributions between wild types and the mutant libraries reflected new genetic variations introduced by error-prone PCR rather than environmental changes or cell-to-cell heterogeneity. In addition, we also confirmed the well-known negative correlation between mean expression levels and standard deviations in GFP distributions regardless of the introduction of mutations or whether the promoters were natural or synthetic (Fig. 3b and c), which is consistent with previous observations (Zaslaver et al. 2006; Wolf et al. 2015). To compare different promoters, we compensated for the mean dependency of the standard deviations of both wild types and mutant libraries by subtracting the regression line obtained from wild types (Fig. 3b and c). We subsequently calculated MVprm, a measure of mutational variability in promoter activity, by subtracting the compensated standard deviation of wild types from that of mutant libraries. Unexpectedly, MVprm, showed similar values across the three promoter groups (Wilcoxon test, adjusted P > 0.05), even though each promoter might show unique values (Fig. 3d). These results indicated that mutational variability of the focal essential promoters showed little detectable signature of natural selection for mutational robustness despite their importance in cellular growth, even though such selection might be active. Moreover, MVprm showed no significant positive correlations with DMenv and DMmut of the downstream genes (Fig. 3e and f). The DM values were defined in Fig. 1. Using the RegulonDB and EcoCyc (Keseler et al. 2021) databases to retrieve the known direct regulatory relationships between genes, we counted the number of direct and unique regulations that each promoter receives from transcriptional regulators and explored the impact of the number of regulations on mutational variability. We found that a larger number of activatory transcriptional regulators facilitated mutational variability in promoter activity (Spearman's R = 0.31; Fig. 3g), probably due to a larger number of mutational target sites. The relationship between MVprm and the number of inhibitory transcriptional regulators was unclear because promoters with null regulation dominated the results (Fig. 3h). Despite that the number of activations was likely to make a difference in MVprm, the lack of significant correlations between MVprm and DM values suggests that the mutational variability of promoters provides no considerable explanation for the observed gene-to-gene differences in transcriptional variability of the downstream genes in response to different perturbations at the genome scale (Fig. 1). These results suggest that the gene-to-gene differences in transcriptional variability in response to genome-wide perturbations may reflect a certain trans-acting molecular mechanism, mainly other than the cis-acting one, consistent with previous observations (Tsuru and Furusawa 2024). This implication was also rational because promoter regions are shorter than trans-acting mutational sites for a single gene.

Fig. 3. Mutational variability in promoter activity. a) Comparison of the standard deviations in the GFP distributions between wild types and mutant libraries. The solid line represents y = x. b and c) Relationship between the mean and standard deviation of the GFP distributions for wild types (b) and mutant libraries (c). The solid lines represent a regression line for wild types and were used to compensate for the mean dependency of the standard deviations of both wild types and mutant libraries by subtraction. d) Mutational variability in promoter activity, termed MVprm, is shown for each promoter group. The MVprm was calculated by subtracting the compensated standard deviations of wild types from those of mutant libraries. The large dots and error bars are the means and standard deviations, respectively. Pairwise Wilcoxon test was examined, where ns represents no significance. e and f) The relationships between MVprm and the DM values (DMenv for e and DMmut for f). The DM values were defined in Fig. 1. g and h) Relationship between MVprm and the number of known activation g) and inhibition h) of natural promoters by transcriptional regulators. The relationship between promoters and known transcriptional regulators were based on RegulonDB (Tierrafria et al. 2022) and EcoCyc (Keseler et al. 2021). The gray, blue, and red dots in each panel represent synthetic (Syn), nonessential (Nes), and essential (Ess) promoters (defined in Fig. 2), respectively. Spearman's R and P-values are shown in a) to c) and e) to h).

Vulnerability in Essential Promotors against Mutations

Next, we explored the directionality of the mutational changes in promoter activity. Consistent with the low mutation rate in error-prone PCR, the mean expression levels of mutant libraries were similar to those of wild types (Spearman's R = 0.99; Fig. 4a). Nevertheless, we found that log10-fold changes in mean expression levels from wild types to mutant libraries followed a global tendency overexpression level, where the promoters with higher activities in wild types tended to decrease their activity in response to mutations (Spearman's R = −056; Fig. 4b). This tendency was retained even when the mean expression levels were calculated for the cells within the narrow gates with the fixed size overexpression levels (supplementary figs. S6 and S7, Supplementary Material online), supporting that the negative correlation was not a technical byproduct of our limitation of detecting a decrease in expression levels. This tendency may reflect bias in the mode, activation or inhibition, of transcriptional regulation (Kinney et al. 2010; Belliveau et al. 2018; Ireland et al. 2020). We hypothesized that promoters that receive more activation than inhibition might tend to reduce their activity in response to mutations. To test this possibility, we defined activation bias as the number of known activations subtracted by the number of inhibitions (supplementary table S1, Supplementary Material online). We found that the activation bias explained the observed mutational changes (Fig. 4c). Furthermore, the activation bias also explained the log10-fold changes after controlling the global trend using the regression line (Fig. 4d), indicating that a promoter becomes more vulnerable as it receives more activation among promoters with similar activities.

Fig. 4. Biased mutational changes in activity of essential promoters. a) Comparison of the means in the GFP distributions between wild types and mutant libraries. The solid line represents y = x. b) Relationship between log10-fold change of mean expression levels from wild types to mutant libraries and the expression levels of wild types. The solid line represents a regression line used to compensate for the mean dependency of the log10-fold changes by subtraction. c and d) Relationship between activation bias in natural promoters and the unnormalized c) and normalized d) log10-fold changes. The activation bias was calculated by subtracting the number of known inhibitions from those of activations by transcriptional regulators. Spearman's R and P-values are shown in a) to d). Gray, blue, and red dots in each panel represent synthetic (Syn), nonessential (Nes), and essential (Ess) promoters (defined in Fig. 2), respectively. e and f) Unnormalized (e) and normalized (f) log10-fold changes for each promoter group. The large dots and error bars are the means and standard deviations, respectively. Asterisks represent significance levels in pairwise Wilcoxon test (ns: P > 0.05; *: P ≤ 0.05; **: P ≤ 0.01). The P-values were adjusted using the BH method (Benjamini and Hochberg 1995).

Interestingly, synthetic and nonessential promoters showed similar mutational changes around zero, whereas essential promoters showed biased directionality in mutational changes toward negative values (Fig. 4e), even though the mean activities were almost the same among the three promoter groups (Fig. 2c). There were no significant differences in promoter length and GC content between the essential and nonessential promoters (Wilcoxon test, P > 0.05). In addition, there was no significant difference in the mutation rate in error-prone PCR between the essential and nonessential promoters (supplementary fig. S8a and table S2, Supplementary Material online). The mutational spectrum in error-prone PCR was quite similar between the two promoter groups (supplementary fig. S8b, Supplementary Material online). These results supported the notion that the observed bias was not derived from a technical bias. This bias in essential promoters was retained for the normalized log10-fold changes (Fig. 4f). Accordingly, essential promoters were more vulnerable than other groups in terms of promoter activity. In other words, this biased directionality indicated the low frequency of the variants with higher activities than the wild types in essential promoters. This biased variability might limit the transcriptional variability of essential genes in response to mutations in cis-regulatory regions, independently from transcriptional variability in response to genomic mutations (Figs. 1 and 3f).

Vulnerability in Essential Promotors against Mutations May Reflect Past Selection for Higher Activity

The difference in mutational vulnerability between essential and nonessential promoters (Fig. 4e and f) could reflect the difference in past selection pressures acting on higher expression levels between essential and nonessential genes. Based on the observation that essential genes tended to have higher expression levels than nonessential genes (Fig. 1), we hypothesized that the remarkable vulnerability of essential promoters reflects local optima in terms of promoter activity. That is, essential promoters might have evolved toward higher activity through past evolution and might have been stuck in local optima because of the high essentiality of downstream genes. In contrast, nonessential promoters might have experienced relatively relaxed selection for higher expression levels because of the low essentiality of downstream genes. To explore the relevance of our assumption, we compared the impact of knockdown on the growth fitness of E. coli between essential and nonessential genes using a public experimental data set reported by Hawkins et al. (2020). We found that knockdown of essential genes tended to be more deleterious for growth fitness than nonessential genes (supplementary fig. S9a, Supplementary Material online), supporting our assumption that essential genes are under stronger selection pressure toward higher expression levels than nonessential genes. To investigate the emergence of differences in vulnerability depending on selection pressure, we performed a numerical simulation of promoter evolution on a simple rugged fitness landscape (Fig. 5). For simplicity, we modeled promoter genotype as 10-bit binary numbers such as 1101011011. That is, our model promoters are composed of biallelic sites (0 or 1), which is simplified from real promoters comprising of quadallelic sites (A, T, G, and C). We assumed that the expression landscape is rugged and genotypes with higher expression levels are relatively rare, consistent with the previous experimental observations (Wolf et al. 2015; Yona et al. 2018). To satisfy this assumption, the mapping form genotype to expression level, equivalent to promoter activity, was defined by the rough Mt. Fuji (RMF) model, as illustrated in Fig. 5a and d. In short, the genotype with the highest expression level, g*, was first chosen at random. The expression levels of the other genotypes were decided as a decreasing linear function of the Hamming distance from g*. Subsequently, stochastic noise was added to generate ruggedness. Finally, expression levels were set to range from 0 (lower bound) to 1 (upper bound) as detailed in Materials and Methods section. Using this expression landscape, mutational changes were calculated for each genotype as the mean difference in the expression levels of all nearest neighbors with a single bit difference, or one Hamming distance, from the expression level of the focal genotype (Fig. 5b and e). This simple expression landscape captured the negative correlation between expression levels and mutational changes (Fig. 5e), which is consistent with the global trend observed in our experiment (Fig. 4b). This trend was likely due to the prevalence of local sinks (green dots) at lower expression levels and frequent local peaks (red dots) at higher expression levels (Fig. 5b and e). To clarify this point, we calculated the fraction of the nearest neighboring genotypes with higher expression levels for each genotype, termed uphill ratio. Local sinks tend to have higher uphill ratios, while local peaks tend to have lower uphill ratios. We confirmed that genotypes with higher uphill ratios dominated at lower expression levels, while genotypes with lower uphill ratios dominated at higher expression levels (Fig. 5f). Importantly, even at the same expression levels, local peaks tended to have negative mutational changes, whereas local sinks tended to show positive changes. The fitness landscape was then modeled by adding the basal fitness, Fb, to the expression levels multiplied by 1−Fb, where Fb ranged from 0 to 1 (Fig. 5c and g). That is, the expression level is under positive directional selection in our model, where the strength of selection pressure depends on Fb. A smaller Fb corresponds to a steeper fitness landscape, while a larger Fb results in a flatter fitness landscape. This simple model assumes a positive linear relationship between expression level and fitness (Fig. 5h), which was relevant to the experimental evidence in most cases (supplementary fig. S9b, Supplementary Material online). Based on the fact that knockdown of essential genes tends to be more deleterious for cell growth than nonessential genes (supplementary fig. S9a, Supplementary Material online), we assigned lower and higher values of Fb for essential and nonessential promoters, respectively. Thus, in our model, essential and nonessential promoters shared the same expression landscape but experienced different strengths of selection pressures. Using the constructed fitness landscape with different Fb values (0.3, 0.5, and 0.99), we performed an evolutionary simulation based on the standard Wright–Fisher model (Otto and Day 2007) with asexual haploid population without recombination and with reversible mutation. Initial genotypes were randomly generated for each simulation. To mimic ongoing evolution of real organisms, the evolutionary simulation was performed for relatively short generations. In practice, the most abundant mutants within the evolved populations after 30 generations were isolated and subjected to measurement of mutational changes (Fig. 5i and j). Despite the short generations, we confirmed that the evolved isolates in steeper fitness landscapes exhibited higher expression levels (Fig. 5k), which is consistent with the designed selection pressure. Comparing the mutational changes between different Fb values, we found that the mutants that evolved with smaller Fb showed negatively larger mutational changes, even at the same expression levels (Fig. 5j). These results suggest that higher activity and vulnerability in essential promoters are linked, and both can emerge simultaneously through stronger positive selection acting on higher expression levels of downstream essential genes.

Fig. 5. Numerical simulation of promoter evolution on a ragged fitness landscape. a to c) Schematics of landscapes of expression levels a and b) and fitness c). a) Expression landscape was constructed by the RMF model. b) Mutational change of each genotype was calculated as the mean of differential expression levels of the nearest neighboring genotypes (arrow heads, mutant types) from the focal genotype (colored circles, wild type). The nearest neighbors differ from the focal genotype by one Hamming distance. The color bar represents the fraction of the genotypes with higher expression levels than the focal genotype among the nearest neighbors, termed uphill ratio. c) Fitness landscape was modeled by adding basal fitness (Fb, horizontal line) to the expression landscape and normalizing it. d and e) Representatives of expression landscape d) and mutational change e). Genotype was modeled as 10-bit binary numbers. Hamming distance represents the distance from the genotype with the global optima in expression level and fitness. f) Relationship between expression level and uphill ratio. The uphill ratio was equally divided into three classes (low, mid, and high). Expression level was also equally divided into ten classes. Frequency of each class of uphill ratio was calculated as the number of corresponding genotypes divided by the total number of genotypes within the class. g and h) Representatives of fitness landscape g) and relationship between fitness and expression level h) (Fb = 0.3). i) Mutational changes of the evolved isolates obtained through evolutionary simulation using the standard Wright–Fisher model. The basal fitness was varied by mimicking strong positive (Fb = 0.3, left), moderate positive (Fb = 0.5, middle), and nearly neutral selection (Fb = 0.99, right) acting upon higher expression levels. Ten independent expression/fitness landscapes were generated for each condition of basal fitness. In each fitness landscape, 600 independent evolutionary simulations were examined where initial genotypes were randomly selected. Dots represent the most abundant mutants within the population after 30 generations in each run. Solid lines represent the running medians. Dashed lines are guides for the eye. j) Enlarged figure of the running medians in i). k) Expression level of the evolved isolates. Asterisks represent significance levels of the adjusted P-value (the BH method) in the Wilcoxon test (****: P ≤ 0.0001).

Despite that our simple linear model basically agreed with the known relationship between expression level and growth fitness in most cases (supplementary fig. S9, Supplementary Material online), an alternative model might be suitable for some genes such as yidC (supplementary fig. S9b, Supplementary Material online). These genes exhibited almost no changes in the growth fitness in response to small knockdown level, while growth fitness gradually decreases with increases in knockdown level beyond certain knockdown levels. To deal with such biphasic cases, we constructed the linear-plateau model (supplementary fig. S10, Supplementary Material online). In short, the linear-plateau model has two phases: a positive linear phase and a flat plateau phase, where the joint point of the two phases is characterized by Ep ranging from 0 to 1 (supplementary fig. S10d, Supplementary Material online). A larger Ep corresponds to a higher optimal expression level. Using the linear-plateau model with different Ep values, we performed an evolutionary simulation by following the same procedure as the linear model (supplementary fig. S10e and f, Supplementary Material online). We found that the evolved isolates under the larger Ep exhibited higher expression levels (supplementary fig. S10g, Supplementary Material online) and negatively larger mutational changes, even at the same expression levels (supplementary fig. S10f, Supplementary Material online). These results seemed to be relevant by considering neutral evolution prevalent in the fitness landscapes with the smaller Ep. These results highlighted an impact of optimal expression level on activity and vulnerability in promoters, suggesting that higher optimal expression levels could also contribute to higher activity and mutational vulnerability in essential promoters. Overall, our findings imply that both higher optimal expression and stronger directional selection could contribute to mutational vulnerability in expression levels for essential genes. Although the expression of essential promoters tended to decrease due to mutations (Fig. 4e and f), a few essential promoters with lower mean expression levels exhibited neutral to positive log10-fold changes in expression response to mutations. A prominent example was the essential promoter controlling tilS, which exhibited the lowest expression level among wild-type essential groups (supplementary table S2, Supplementary Material online). This promoter showed a positive value in log10-fold change in response to mutations, indicating the frequent occurrence of mutants with higher expression levels than wild types. These essential promoters might suggest weaker directional selection or lower optimal expression levels for the downstream essential genes.

Consistent results were also obtained using the NK model (Kauffman and Weinberger 1989) for rugged expression landscapes (supplementary figs. S11 and S12, Supplementary Material online, for the linear and linear-plateau models, respectively). The NK model, similar to the RMF model, frequently provided sinks (genotypes with higher uphill ratios) at lower expression levels and peaks (genotypes with lower uphill ratios) at higher expression levels (supplementary fig. S11c, Supplementary Material online). However, compared to the RMF model (Fig. 5f), sinks were relatively rare at middle to high expression levels (beyond 0.3 in expression level). Interestingly, peaks were also rare at middle expression levels (0.3 to 0.6 in expression level) in the NK model. That is, genotypes with moderate uphill ratios were relatively frequent at middle to high expression levels. The low frequency of sinks and peaks indicates fewer genotypes with higher or lower mutational changes. Consequently, the difference in mutational changes between different genotypes at the same expression levels was relatively small in the NK model. Therefore, the difference in mutational changes between different Fb (or Ep) values in the NK model was small relative to the RMF model (supplementary figs. S11g and S12f, Supplementary Material online). Despite these minor differences, both expression landscapes showed consistent outcomes for the linear (Fig. 5j and k for the RMF model and supplementary fig. S11g and h, Supplementary Material online for the NK model, respectively) and linear-plateau models (supplementary fig. S10f and g, Supplementary Material online for the RMF model and supplementary fig. S12f and g, Supplementary Material online for the NK model, respectively). These results support the broad applicability of our key assumptions underlying rugged expression landscapes, suggesting the robustness of the outcomes derived from our simple models.

Discussion

Transcriptional variability differs between genes. Using public transcriptome data sets for E. coli, we found that essential genes exhibited lower transcriptional variability across different mutants and environmental conditions than nonessential genes. To explore the underlying mechanism, we focused on the promoter regions. We quantified the changes in promoter activity changes in response to mutations by constructing mutant libraries of natural and synthetic promoters. As a result, essential promoters did not exhibit a significant difference in mutational variability compared with nonessential and synthetic promoters. Unexpectedly, we found that essential promoters showed remarkable vulnerability to mutations, where variants with higher activities than wild type were rare and those with lower activity were abundant. These results imply that essential promoters have evolved to reach local optima in terms of promoter activity through past selections for higher expression levels of downstream essential genes. Our numerical evolutionary simulation confirmed the relevance of this scenario. Overall, our results suggest the evolutionary impact of essentiality on the directionality of the mutational variability of promoters and imply that the directional bias in mutational variability of promoter regions might constrain further evolution of expression levels of essential genes toward higher expression levels.

There are three well-known evolutionary scenarios for the emergence of mutational robustness in many biological systems (de Visser et al. 2003), namely adaptive (Waddington 1942), congruent (Meiklejohn and Hartl 2002), and intrinsic (Siegal and Bergman 2002) scenarios. The adaptive scenarios predict that essential traits have evolved to reduce mutational variability and minimize the emergence of harmful mutants through direct selection for mutational robustness. The congruent scenarios predict that mutational robustness is a byproduct of evolution toward low variability against environmental perturbations. Contrary to these two scenarios assuming adaptive evolution toward low variability, the intrinsic scenarios predict that the mutational robustness of the traits emerges spontaneously and inherently along with the emergence of traits without selection against variability. In the case of promoter activity, an intrinsic scenario corresponds to the mutational robustness of promoters emerging spontaneously with the emergence of promoter activity without selection for mutational robustness, which corresponds to the mutational variability of synthetic promoters. Our results demonstrated that both adaptive and congruent scenarios failed to explain the observed mutational variability in natural promoters because essential promoters showed similar mutational variability to the other promoters, and the natural promoters exhibited no significant positive correlation with the transcriptional variability of downstream genes against different environmental perturbations. We note that the adaptive scenario might explain the higher expression levels, not transcriptional variability, of essential genes (Fig. 1c and g), which might minimize the risk of mutational knockdown for essential genes to a lethal level. In contrast, the intrinsic scenario explained the observed mutational variability because the natural promoters showed similar mutational variability to that of synthetic promoters. Thus, the mutational variability of essential promoters is unlikely to be a consequence of adaptive evolution toward robustness against genetic and environmental perturbations.

Why did essential promoters show no less mutational variability than other promoters, despite the fact that the reduction in expression levels of essential genes is more likely to be harmful to cell growth than nonessential genes (Cui et al. 2018)? One possible reason is that promoter sequences with significant mutational robustness are rare, although several possible mechanisms underlying mutational robustness in promoter regions are hypothesized (Payne and Wagner 2015). In other words, the expression landscape of the promoters might be rugged everywhere as we modeled in Fig. 5a. In fact, synthetic promoters derived from different random sequences achieved higher promoter activities equivalent to natural promoters, which was consistent with previous observations (Wolf et al. 2015; Yona et al. 2018). A second possible reason is that selection for higher expression levels on a rugged fitness landscape (Fig. 5c) prevents the essential promoters at local optima from traversing the expression landscape to reach genotypes with higher mutational robustness. Consistently, essential promoters showed biased mutational effects, where most mutations in essential promoters tended to reduce promoter activity (Fig. 4e and f). Interestingly, this scenario agreed with the predictions from a machine-learning model based on an empirical expression landscape of yeast promoters (Vaishnav et al. 2022). Thus, our results suggest several plausible hypotheses that may be shared by different organisms. Ideally, future studies should explore evolutionary scenarios for diverse promoters that were not tested in this study and validate the hypothetical landscape of promoters more thoroughly. The increase of the sample set of promoters would also compensate for our small sample sets and would enhance power to detect a potential subtle difference in mutational variability between different promoter groups.

In this study, we focused only on mutations in promoter regions to explore the mechanism underlying the lower transcriptional variability in essential genes. As a result, we found that mutants with higher promoter activity were relatively rare for essential promoters. However, this result does not exclude the contribution of any other mechanism underlying the observed transcriptional variability in response to genome-wide perturbations (Fig. 1). Notably, trans-acting mechanisms, such as transcriptional regulators, may also play a role. As it is, the previous studies reported that genes exhibiting greater transcriptional variability across diverse environmental and genetic perturbations tend to be regulated by a larger number of unique transcriptional regulators in E. coli (Tsuru and Furusawa 2024) and yeast (Landry et al. 2007). A higher number of regulators provide more trans-mutational target sites for the regulated genes, thereby enhancing transcriptional variability in response to genomic perturbations. Based on this evidence, it is plausible that essential genes have less transcriptional regulation. To support this hypothesis, we used data from RegulonDB and EcoCyc to compare the number of regulations between essential and nonessential genes (supplementary fig. S13, Supplementary Material online). First, we found a smaller fraction of essential genes (90 of 248 genes, 36%) lacking known transcriptional regulators than nonessential genes (1976 of 3995, 49%) (supplementary fig. S13a, Supplementary Material online). This finding aligns with the positive activation bias values observed (Fig. 4c and d), probably attributable to the rich activation by sigma factors aimed at enhancing expression levels. On the other hand, among the genes regulated by at least one transcriptional regulator, we observed that essential genes received fewer transcriptional regulations than nonessential genes (supplementary fig. S13b, Supplementary Material online), supporting this hypothesis. This observation is consistent with previous findings in yeast (Bilu and Barkai 2005). Thus, our study underscores the dual contributions of both cis- and trans-acting mechanisms including their epistatic interactions in constraining the transcriptional variability of essential genes.

Materials and Methods

Bacterial Strains and Medium

We used E. coli DH5α electro-cells (Takara Bio, Japan, 9027) to construct the pUA66-mCherry plasmid and its derivatives. LB medium (LB Broth, Miller, BD Difco, USA, 244620) supplemented with 25-µg/mL kanamycin (kanamycin sulfate, Wako, Japan, 115-00342), referred to as LB/Km medium, served as the selective medium. The SOC medium, derived from SOB medium (SOB Medium, BD Difco, USA, 244310) supplemented with 22 mM Glucose (Wako, Japan, 049-31165), was used for recovery culture posttransformation. The E. coli K-12 MG1655 strain served as the host for the constructed plasmids and was employed for measuring promoter activity via flow cytometry. Bacterial cells for flow cytometry were cultured in a synthetic medium, a derivative of M9 minimal medium (pH 7.2), supplemented with 47.7 mM Na2HPO4 (disodium hydrogenphosphate 12-Water, Wako, Japan, 196-02835), 22.0 mM KH2PO4 (Wako, Japan, 164-22635), 8.55 mM NaCl (Wako, Japan, 191-01665), 9.35 mM NH4Cl (Wako, Japan, 014-03005), 11 mM Glucose, 1 mM MgSO4 (magnesium sulfate heptahydrate, Wako, Japan, 138-00415), 50 µM CaCl2 (Wako, Japan, 038-24985), and 3.8 µM thiamin (thiamin hydrochloride, Wako, Japan, 201-00852), with supplementation of 25-µg/mL kanamycin when required.

Selection of Natural Promoters

The E. coli promoter collection (Zaslaver et al. 2006) (Horizon Discovery Ltd., UK, PEC3877) was used to obtain the promoter regions of the natural promoters in E. coli. Initially, we filtered out the promoters with low activities that were near our detection limit based on the promoter activities of this library measured by Silander et al. (2012). This filtration excluded 12% and 17% of essential and nonessential promoters, respectively. Subsequently, 115 promoters were randomly selected from this subset library carrying the pUA66 plasmid. The DNA sequences of these promoter regions were verified using Sanger sequencing (Azenta Life Sciences, USA) and compared with the reference genome sequence of E. coli K-12 MG1655 (Blattner et al. 1997) (GenBank ID: U00096). Subsequently, 63 of the 115 strains without mutations in their promoter sequences were identified, and plasmids were extracted from these colonies. Among these, 56 strains were randomly selected for the subsequent experiments (supplementary table S2, Supplementary Material online). Each promoter was denoted by the “NP” (natural promoter), followed by the name of the ORF present in the immediate downstream sequence in the genome. For example, NPmetG represents the promoter regions upstream of the metG gene on the chromosome.

Plasmid Construction

Plasmid pUA66 was modified by replacing the promoter regions of gfp, flanked by the BamHI and XhoI restriction sites, with a red fluorescence protein (mCherry) expression cassette. This cassette comprised dual terminators RNAI and TSAL (Reynolds et al. 1992), promoter PLtetO-1 (Lutz 1997), Shine-Dalgarno sequence SD8 (Ringquist et al. 1992), and codon-optimized mCherry (Dunlop et al. 2008) coding sequences. The terminators were obtained from pNS2-σVL (Dunlop et al. 2008) (Addgene, Plasmid #26756), and the assembled sequences of the other elements were commercially synthesized (Fasmac Co., Ltd., Japan). The assembled cassette was inserted between the BamHI and XhoI recognition sequences of the pUA66 plasmid in the direction opposite to that of gfp (supplementary fig. S2, Supplementary Material online), yielding plasmid pUA66-mCherry. PCR amplification was performed using KOD One PCR Master Mix Blue (Toyobo, Japan, KMM-201) and pUA66-mCherry as a template with primers A (GCTTTCCAGGGATCCTCTAGATTTAAGAAGG) and B (GAGCAGTACTCTCGAGGTGAAGACGAAAGG) to include all sequences except for the mCherry expression cassette. The PCR product was column purified using the FastGene Gel/PCR Extraction Kit (NIPPON Genetics, Japan, FG-91302). The purified PCR product was double digested with BamHI (Takara Bio, Japan, 1010A) and XhoI (Takara Bio, Japan, 1094A) overnight at 37 °C, followed by column purification. The digested product is referred to as the vector.

Construction of Promoter Library Comprising of Random Sequences

The construction of the synthetic promoters was based on a previous study (Wolf et al. 2015). Single-stranded oligonucleotides (CCTTTCGTCTTCACCTCGAG- (N200)- GGATCCTCTAGATTTAAGAAGG) were synthesized by Eurofins Genomics (Ebersberg, Germany). The 5′ and 3′ ends of the oligos contained XhoI and BamHI restriction sites, respectively, and were homologous to primers C (CCTTTCGTCTTCACCTCGA) and D (CCTTCTTAAATCTAGAGGATCC), respectively (supplementary fig. S2, Supplementary Material online). A 200-base region of random nucleotides was present between them. PCR amplification using primers C and D was performed using a mixture of single-stranded oligos as a template to create double-stranded DNA. The PCR products were column purified and double digested using BamHI and XhoI for 3 h at 37 °C, followed by column purification. Thus, an insert for the following ligation reaction was obtained. The ligation reaction was performed using T4 DNA Ligase (Takara Bio, Japan, 2011A) overnight at 16 °C involved mixing the vector and inserting at a molar ratio of 1:3. The resulting ligation product solution was column purified and DNA was extracted using 20 µL of sterile water. All purified ligation product solutions were added to 50 µL of DH5α electro-cells and electroporated at 1,800 V (Eporator, Eppendroph, Germany, 4309000027). Transformed cells were incubated at 37 °C for 1 h in 1 mL of SOC medium, followed by overnight incubation at 37 °C in 1 mL of LB medium supplemented with 2 µL of 25-mg/mL kanamycin solution.

Selection of Synthetic Promoters Using Flow Cytometry

Green and red fluorescence, as well as forward scatter of the cells in the DH5α-transformed random promoter library, was measured using a flow cytometer (FACSAria III, BD, USA). Gates were set for the green and red channels to exclude cells with low GFP and high mCherry fluorescence, respectively. Subsequently, 10,000 cells without red fluorescence but with green fluorescence within the selected range (Fig. 2b) were collected by cell sorting, according to the manufacturer's instructions. The sorted cells were spread on LB/Km agar plates and incubated overnight at 37 °C. Bacterial colonies were randomly picked, and their inserted synthetic promoter regions were sequenced by Sanger sequencing. Consequently, 16 synthetic promoters without improper mutations at the restriction sites were identified and named SP1–16, (supplementary table S1, Supplementary Material online).

Random Mutagenesis Using Error-Prone PCR

Error-prone PCR was conducted using 0.1 ng of each promoter region as a template, with 35 cycles of amplification performed using the GeneMorph II Random Mutagenesis Kit (Agilent Technologies, USA, 200550) and primers C and D. The resulting PCR products underwent gel purification and were double digested with BamHI and XhoI for 3 h at 37 °C, followed by column purification to create mutant inserts. Ligation was performed overnight at 16 °C using T4 DNA Ligase, with the vector and mutant inserts mixed at a molar ratio of 1:3. The ligation reaction solution was column purified, and DNA was extracted using 20 µL of sterile water.

Construction of Mutant Promoter Library

The purified ligation product solution was added to 50 µL of DH5α electro-cells and electroporated at 1,800 V. Transformed cells were cultured at 37 °C for 1 h in 1 mL of SOC medium, and a portion of the culture was plated on LB/Km agar medium for overnight growth at 37 °C. Bacterial colonies were counted to determine the transformation efficiency. If the transformation efficiency exceeded 105 cells/reaction, the next step was performed; otherwise, the process was repeated. The transformed culture was diluted with 5 mL of LB/Km medium and grown overnight. Plasmid extraction from an overnight culture was performed using the FastGene Plasmid Mini Kit (NIPPON Genetics, Japan, FG-90502). The obtained plasmid solution was added to 50 µL of MG1655 electrocompetent cells and electroporated at 1,800 V. Transformed cells were cultured at 37 °C for 1 h in 1 mL of SOC medium, and a portion of the solution was plated on LB/Km agar medium for overnight growth at 37 °C to count the number of colonies and determine transformation efficiency. If the transformation efficiency exceeded 105 cells/reaction, the next step was performed; otherwise, the process was repeated. A sterile 2-mL 96-well deep well plate (Greiner Bio-One, Austria, 780271) was used to add 1 mL of the synthetic medium to each well, followed by the addition of 10 µL of the overnight culture containing the MG1655 transformants. The bacterial cells were grown at 37 °C with shaking at 800 rpm until confluent growth was achieved. Subsequently, another plate was prepared with 1 mL of fresh synthetic medium added to each well, and 10 µL of the confluently grown cultures were diluted into each well. The cultures were grown for 2.5 h at 37 °C with shaking at 800 rpm and were used as mutant libraries for flow cytometry.

Verification of PCR Random Mutagenesis

The number of mutations introduced into the promoters by error-prone PCR was determined using Sanger sequencing. Three to four colonies were randomly selected from six randomly selected mutant libraries, including essential and nonessential promoters (supplementary table S2, Supplementary Material online). Colony PCR was performed using KOD One PCR Master Mix Blue with primers E (TTACTTTGCAGGGCTTCCCAA) and F (CCAGTCTTTCGACTGAGCCT) to amplify the regions including the promoter regions (supplementary fig. S2, Supplementary Material online). PCR products were column purified and verified by Sanger sequencing. The average mutation rate was calculated to be 6.2 mutations/kb. Mutation rates and mutational spectrum for each promoter group were shown in supplementary fig. S8, Supplementary Material online.

Construction of Wild-Type Promoter Control Group

To ensure accurate measurement of GFP fluorescence from cells with wild-type promoters, wild-type plasmids were transformed into MG1655 cells by electroporation. Transformed cultures were plated on LB/Km agar medium and incubated overnight, and the resulting colonies were cultured in LB/Km medium. Cell cultures of these control groups in the synthetic medium were prepared and used as wild types for flow cytometry.

Flow Cytometry

Bacterial cells were dispensed in M9 buffer prior to analysis using a flow cytometer (BD, USA, FACSAria III). Fluorescence intensity from GFP (GFP FI) and mCherry (mCherry FI) was measured using a 488-nm laser with 515- to 545-nm emission filter and a 561-nm laser with 563- to 588-nm emission filter, respectively. The following PMT voltage settings were used: forward scatter (FSC), 300; side scatter (SSC), 340; green fluorescence, 700; and red fluorescence, 600. The events exceeding the defined thresholds of FSC (200) and SSC (200) were recorded. A total of 1,000,000 events were measured per library. The MG1655 strain was used as a control for autofluorescence of GFP FI.

Data Analysis of Flow Cytometry

FCS files were converted to CSV format using the flowCore (Ellis et al. 2024) package in R (R Core Team 2023). Subsequent analyses were conducted using custom R scripts. To compare green fluorescence among cells in similar physiological states, a narrow gate for the forward scatter channel was set to capture 40% of the events, including the most frequent values for MG1655 (supplementary fig. S5a, Supplementary Material online). The top 0.5% value of the green fluorescence intensity distribution of MG1655 within the gate was determined as the autofluorescence intensity (supplementary fig. S5b, Supplementary Material online). For the events of wild types and mutant libraries, the autofluorescence intensity was subtracted from the green fluorescence intensity of the events within the gate. Events with a negative green fluorescence intensity after subtracting the autofluorescence intensity were excluded from the analysis. In addition, events with a red fluorescence intensity greater than 103 [a.u.] were regarded as the cells carrying pUA66-mCherry and were excluded from the analysis. Clonal cell populations exhibit huge cell-to-cell heterogeneity in gene expression even under identical environmental conditions, resulting in long-tailed or log-normal shaped distributions in protein numbers (Tsuru et al. 2009; Taniguchi et al. 2010). This stochastic nature often results in rare cells exhibiting very highly expression levels within cell populations, which greatly influences the variance or standard deviation of protein number distributions. To ensure robust measurement of variation in fluorescent distributions, the green fluorescence intensity after subtracting the autofluorescence intensity was log10-transformed and used to calculate the means and standard deviations. We note that this log-transformation might limit the detection of possible small differences in mutational variability in promoter activities between different promoter groups.

Classification of Promoter Based on Essentiality of Transcribed Unit

Gene essentiality was based on a previous study (Goodall et al. 2018), which identified 248 genes designated as essential across different data sets. In bacteria, clusters of genes are controlled by the same promoter sequences and are transcribed as single mRNAs, termed transcription units. Promoters controlling at least one essential gene in the transcription units were defined as essential promoters (Fig. 2a). The known transcription units were retrieved from the RegulonDB (Tierrafria et al. 2022) database. Of the 56 promoters, 18 were classified as essential and 38 as nonessential as listed in supplementary table S1, Supplementary Material online.

Estimation of Transcriptional Variability

Transcriptional variability of each gene in E. coli in response to different genetic perturbations was estimated using known transcriptome profiles obtained from six independent mutants derived from MG1655 cultured under a single environmental condition, termed the Mut data set (Tsuru and Furusawa 2024). These mutants were subjected to MA experiments to accumulate genomic mutations via genetic drift (Halligan and Keightley 2009; Tsuru et al. 2015). In addition, to quantify the transcriptional variability in response to different environmental conditions, known transcriptome profiles of MG1655 cultured under 76 different environmental conditions (Sastry et al. 2019; Rychel et al. 2021; Lamoureux et al. 2021), termed the Env data set (Tsuru and Furusawa 2024), were utilized. The expression levels were determined by log2-transformation and quantile normalization. The means and standard deviations in expression levels for each gene across different genetic perturbations and environmental conditions were quantified. To compensate for the mean dependency of the standard deviations (Fig. 1b and f), the DM values for each gene were calculated as the vertical distance of each standard deviation from a smooth running median of the standard deviations for each data set, as detailed previously (Tsuru and Furusawa 2024). The resulting DM values were termed DMmut and DMenv for the Mut and Env data sets, respectively. The DM values were used as measures of transcriptional variability in response to genetic and environmental perturbations.

Phylogenetic Analysis of Promoter Sequence of Length 73 Bases

Synthetic promoters underwent a blastn (Altschul et al. 1990) search to confirm that there was no considerable similarity to any natural sequences (supplementary note S1, Supplementary Material online). Posterior sequences of 73 bases were aligned using the msa (Bodenhofer et al. 2015) package in R (Clustal Omega), where 73 bases represented the shortest length of a synthetic promoter. A phylogenetic tree was constructed using the neighbor-joining method.

Phylogenetic Analysis of Core Elements of sigma70 Promoters

We identified sigma70 promoters from our natural promoters using RegulonDB. The core elements comprising the −10/−35 elements of these promoters were concatenated and subjected to phylogenetic analysis. The Levenshtein distance (Levenshtein 1966) between the focal elements and the consensus elements, which consist of TTGACA and TATAAT, was calculated. A phylogenetic tree was constructed as previously described.

Identification of the Number of Transcriptional Regulators Regulating the Natural Promoters

Reported regulatory interactions between genes were compiled using all interactions with experimental evidence from RegulonDB and EcoCyc (Keseler et al. 2021) for both transcription factors and sigma factors. The mode of regulation was classified as either activation or inhibition, and the mode of regulation by sigma factors was classified as activation.

Numerical Simulation of Promoter Evolution

The promoter genotype, g, was modeled by a bit string comprising either zeros or ones of length N = 10 (e.g. 1110001100). The expression level, E, was defined by the RMF model (Aita et al. 2000; Neidhart et al. 2014) as:

E(g)=−cD(g,g*)+η(g),

where c is a scaling coefficient equal to one in this study, g* is the genotype with the highest expression level before adding noise, D(g,g*) is the Hamming distance between g and g*, and η is an independent and identically distributed random variable. A standard normal distribution with a mean equal to 0 and a standard deviation equal to 1 for η was used to generate a rugged expression landscape. To construct an expression landscape, g* was first chosen at random from the genotype space. E(g) was subsequently calculated using the above equation. The expression landscape was then normalized using the minimal and maximal values of E(g) to range from 0 to 1, as follows:

E^(g)=E(g)−min(E(g))max(E(g))−min(E(g)).

The fitness of a given genotype g was defined as:

F(g)=(1−Fb)E^(g)+Fb,

where 1 ≥ Fb ≥ 0. Fb represents the basal fitness when the expression level is 0 and determines the strength of the selection pressure acting on the expression level. Different values of Fb (0.3, 0.5, and 0.99) were used to demonstrate the evolution under strong positive selection, moderate positive selection, and nearly neutral evolution, respectively. It was assumed that essential genes underwent strong positive selection on their expression levels, while nonessential genes underwent relatively moderate selection for their expression levels, based on the empirical fact that knockdown of essential genes was more harmful for cell growth than nonessential genes (Li et al. 2016). Ten independent fitness landscapes were constructed for each Fb. Evolutionary simulation on the fitness landscapes was performed using the standard Wright–Fisher model (Otto and Day 2007) assuming a haploid asexual population with a population size of 105 and a mutation rate of 10−3 mutations per genotype per generation. The initial genotypes were randomly selected for each simulation, and 600 independent simulations were examined for each fitness landscape. For each generation, mutants with genotypes differing by one Hamming distance from the current population were generated at the mutation rate. The frequency of each genotype within the population was updated based on the fitness landscape. Subsequently, a new population for the next generation was formed by sampling mutants, with replacement, from the current population using a binomial distribution with the given population size. The most frequent mutants in a population after 30 generations were isolated in each simulation, resulting in 6,000 mutants for each Fb. Mutational changes in expression levels were calculated for each isolate by averaging the expression levels of all nearest neighboring genotypes that differed by one Hamming distance from the focal isolate. Custom R scripts were used for the numerical simulations, modified from scripts constructed in a previous study (Obolski et al. 2018; Song and Zhang 2021).

Data Visualization

All figures were generated in R and the following packages were used. Phylogenetic trees of the promoters were visualized using the ape (Paradis and Schliep 2019) and ggtree (Yu et al. 2017) packages. Other illustrative plots were generated using the ggplot2 (Wickham 2016) and ggpubr (Kassambara 2023) packages.

Supplementary Material

msae185_Supplementary_Data

Acknowledgments

The authors thank Hiroyo Koike for the technical assistance of Sanger sequencing. The authors also appreciate Dr. Olin Silander and Dr. Luise Wolf for fruitful technical advice on the construction of a synthetic promoter library. The pNS2-σVL plasmid was a gift from Dr. Michael Elowitz (Addgene, plasmid # 26756).

Supplementary Material

Supplementary material is available at Molecular Biology and Evolution online.

Funding

This study was supported by the Japan Society for the Promotion of Science (JSPS) KAKENHI (18H02427, 22H05403, and 24K21985 to S.T.; 17H06389 to C.F. and S.T.; and 22K21344 and 24H01798 to C.F.) and the Japan Science and Technology Agency (JST) ERATO (JPMJER1902 to S.T. and C.F.).

Data Availability

Custom R scripts used in the manuscript are available at GitHub (https://github.com/tsuruubi/Variability_PromoterAct_E_coli). The raw fcs data are available at Zenodo (wild types for https://doi.org/10.5281/zenodo.11108555; mutant libraries for https://doi.org/10.5281/zenodo.11108786).
==== Refs
References

Aita  T, Uchiyama  H, Inaoka  T, Nakajima  M, Kokubo  T, Husimi  Y. Analysis of a local fitness landscape with a model of the rough Mt. Fuji-type landscape: application to prolyl endopeptidase and thermolysin. Biopolymers. 2000:54 (1 ):64–79. 10.1002/(SICI)1097-0282(200007)54:1<64::AID-BIP70>3.0.CO;2-R.10799982
Altschul  SF, Gish  W, Miller  W, Myers  EW, Lipman  DJ. Basic local alignment search tool. J Mol Biol.  1990:215 (3 ):403–410. 10.1016/S0022-2836(05)80360-2.2231712
Belliveau  NM, Barnes  SL, Ireland  WT, Jones  DL, Sweredoski  MJ, Moradian  A, Hess  S, Kinney  JB, Phillips  R. Systematic approach for dissecting the molecular mechanisms of transcriptional regulation in bacteria. Proc Natl Acad Sci U S A. 2018:115 (21 ):E4796–E4805. 10.1073/pnas.1722055115.29728462
Benjamini  Y, Hochberg  Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc Series B (Method). 1995:57 (1 ):289–300. 10.1111/j.2517-6161.1995.tb02031.x.
Bilu  Y, Barkai  N. The design of transcription-factor binding sites is affected by combinatorial regulation. Genome Biol. 2005:6 (12 ):R103. 10.1186/gb-2005-6-12-r103.16356266
Blattner  FR, Plunkett  G, Bloch  CA, Perna  NT, Burland  V, Riley  M, Collado-Vides  J, Glasner  JD, Rode  CK, Mayhew  GF, et al  The complete genome sequence of Escherichia coli K-12. Science. 1997:277 (5331 ):1453–1462. 10.1126/science.277.5331.1453.9278503
Bodenhofer  U, Bonatesta  E, Horejs-Kainrath  C, Hochreiter  S. MSA: an R package for multiple sequence alignment. Bioinformatics. 2015:31 (24 ):3997–3999. 10.1093/bioinformatics/btv494.26315911
Casadesus  J, Low  D. Epigenetic gene regulation in the bacterial world. Microbiol Mol Biol Rev. 2006:70 (3 ):830–856. 10.1128/MMBR.00016-06.16959970
Cui  L, Vigouroux  A, Rousset  F, Varet  H, Khanna  V, Bikard  D. A CRISPRi screen in E. coli reveals sequence-specific toxicity of dCas9. Nat Commun.  2018:9 (1 ):1912. 10.1038/s41467-018-04209-5.29765036
Darwin  C . On the origin of species by means of natural selection, or preservation of favoured races in the struggle for life. London: John Murray; 1859.
Denver  DR, Morris  K, Streelman  JT, Kim  SK, Lynch  M, Thomas  WK. The transcriptional consequences of mutation and natural selection in Caenorhabditis elegans. Nat Genet. 2005:37 (5 ):544–548. 10.1038/ng1554.15852004
de Visser  JA, Hermisson  J, Wagner  GP, Ancel Meyers  L, Bagheri-Chaichian  H, Blanchard  JL, Chao  L, Cheverud  JM, Elena  SF, Fontana  W, et al  Perspective: evolution and detection of genetic robustness. Evolution. 2003:57 (9 ):1959–1972. 10.1111/j.0014-3820.2003.tb00377.x.14575319
Draghi  JA, Whitlock  MC. Phenotypic plasticity facilitates mutational variance, genetic variance, and evolvability along the major axis of environmental variation. Evolution. 2012:66 (9 ):2891–2902. 10.1111/j.1558-5646.2012.01649.x.22946810
Dunlop  MJ, Cox  RS, Levine  JH, Murray  RM, Elowitz  MB. Regulatory activity revealed by dynamic correlations in gene expression noise. Nat Genet.  2008:40 (12 ):1493–1498. 10.1038/ng.281.19029898
Duveau  F, Vande Zande  P, Metzger  BP, Diaz  CJ, Walker  EA, Tryban  S, Siddiq  MA, Yang  B, Wittkopp  PJ. Mutational sources of trans-regulatory variation affecting gene expression in Saccharomyces cerevisiae. eLife. 2021:10 :e67806. 10.7554/eLife.67806.34463616
Ellis  B, Haaland  P, Hahne  F, Meur  L, NolwennGopalakrishnan  N, Spidlen  J, Jiang  M, Finak  G, Granjeaud  S. flowCore: flowCore: basic structures for flow cytometry data. R package version 2.14.2. [2024; accessed 2024 Apr]. https://bioconductor.org/packages/flowCore/.
Furusawa  C, Kaneko  K. Formation of dominant mode by evolution in biological systems. Phys Rev E. 2018:97 (4 ):042410. 10.1103/PhysRevE.97.042410.29758752
Goodall  ECA, Robinson  A, Johnston  IG, Jabbari  S, Turner  KA, Cunningham  AF, Lund  PA, Cole  JA, Henderson  IR. The essential genome of Escherichia coli K-12. mBio. 2018:9 (1 ):1–18. 10.1128/mBio.02096-17.
Halligan  DL, Keightley  PD. Spontaneous mutation accumulation studies in evolutionary genetics. Annu Rev Ecol Evol Syst.  2009:40 (1 ):151–172. 10.1146/annurev.ecolsys.39.110707.173437.
Hawkins  JS, Silvis  MR, Koo  BM, Peters  JM, Osadnik  H, Jost  M, Hearne  CC, Weissman  JS, Todor  H, Gross  CA. Mismatch-CRISPRi reveals the co-varying expression-fitness relationships of essential genes in Escherichia coli and Bacillus subtilis. Cell Syst. 2020:11 (5 ):523–535.e9. 10.1016/j.cels.2020.09.009.33080209
Ho  WC, Zhang  J. The genotype–phenotype map of yeast complex traits: basic parameters and the role of natural selection. Mol Biol Evol. 2014:31 (6 ):1568–1580. 10.1093/molbev/msu131.24723420
Ho  WC, Zhang  J. Adaptive genetic robustness of Escherichia coli metabolic fluxes. Mol Biol Evol.  2016:33 (5 ):1164–1176. 10.1093/molbev/msw002.26733489
Hornung  G, Oren  M, Barkai  N. Nucleosome organization affects the sensitivity of gene expression to promoter mutations. Mol Cell. 2012:46 (3 ):362–368. 10.1016/j.molcel.2012.02.019.22464732
Ireland  WT, Beeler  SM, Flores-Bautista  E, McCarty  NS, Roschinger  T, Belliveau  NM, Sweredoski  MJ, Moradian  A, Kinney  JB, Phillips  R. Deciphering the regulatory genome of Escherichia coli, one hundred promoters at a time. eLife. 2020:9 :e55308. 10.7554/eLife.55308.32955440
Kaneko  K . Relationship among phenotypic plasticity, phenotypic fluctuations, robustness, and evolvability; Waddington's legacy revisited under the spirit of Einstein. J Biosci. 2009:34 (4 ):529–542. 10.1007/s12038-009-0072-9.19920339
Kassambara  A.  ggpubr: ‘ggplot2’ based publication ready plots. R package version 0.6.0. [2023; accessed 2024 Apr]. https://rpkgs.datanovia.com/ggpubr/.
Kauffman  SA, Weinberger  ED. The NK model of rugged fitness landscapes and its application to maturation of the immune response. J Theor Biol. 1989:141 (2 ):211–245. 10.1016/S0022-5193(89)80019-0.2632988
Keseler  IM, Gama-Castro  S, Mackie  A, Billington  R, Bonavides-Martinez  C, Caspi  R, Kothari  A, Krummenacker  M, Midford  PE, Muniz-Rascado  L, et al  The EcoCyc database in 2021. Front Microbiol. 2021:12 :711077. 10.3389/fmicb.2021.711077.34394059
Kinney  JB, Murugan  A, Callan  CG  Jr, Cox  EC. Using deep sequencing to characterize the biophysical mechanism of a transcriptional regulatory sequence. Proc Natl Acad Sci U S A. 2010:107 (20 ):9158–9163. 10.1073/pnas.1004290107.20439748
Lamoureux  CR, Decker  KT, Sastry  AV, Rychel  K, Gao  Y, McConn  JL, Zielinski  DC, Palsson  BO.  2021. A multi-scale transcriptional regulatory network knowledge base for Escherichia coli. bioRxiv 439047. 10.1101/2021.04.08.439047, 11 April 2021, preprint: not peer reviewed.
Landry  CR, Lemos  B, Rifkin  SA, Dickinson  WJ, Hartl  DL. Genetic properties influencing the evolvability of gene expression. Science. 2007:317 (5834 ):118–121. 10.1126/science.1140247.17525304
Levenshtein  VI . Binary codes capable of correcting deletions, insertions, and reversals. Soviet physics doklady. 1966:10 (8 ):707–710.
Li  XT, Jun  Y, Erickstad  MJ, Brown  SD, Parks  A, Court  DL, Jun  S. tCRISPRi: tunable and reversible, one-step control of gene expression. Sci Rep. 2016:6 (1 ):39076. 10.1038/srep39076.27996021
Lutz  R . 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 Res.  1997:25 (6 ):1203–1210. 10.1093/nar/25.6.1203.9092630
Lynch  M, Walsh  B. Genetics and analysis of quantitative traits. Sunderland (MA): Sinauer; 1998.
McGuigan  K, Collet  JM, McGraw  EA, Ye  YH, Allen  SL, Chenoweth  SF, Blows  MW. The nature and extent of mutational pleiotropy in gene expression of male Drosophila serrata. Genetics. 2014:196 (3 ):911–921. 10.1534/genetics.114.161232.24402375
Meiklejohn  CD, Hartl  DL. A single mode of canalization. Trends Ecol Evol.  2002:17 (10 ):468–473. 10.1016/S0169-5347(02)02596-X.
Metzger  BP, Yuan  DC, Gruber  JD, Duveau  F, Wittkopp  PJ. Selection on noise constrains variation in a eukaryotic promoter. Nature. 2015:521 (7552 ):344–347. 10.1038/nature14244.25778704
Neidhart  J, Szendro  IG, Krug  J. Adaptation in tunably rugged fitness landscapes: the rough mount fuji model. Genetics. 2014:198 (2 ):699–721. 10.1534/genetics.114.167668.25123507
Noble  DWA, Radersma  R, Uller  T. Plastic responses to novel environments are biased towards phenotype dimensions with high additive genetic variation. Proc Natl Acad Sci U S A. 2019:116 (27 ):13452–13461. 10.1073/pnas.1821066116.31217289
Obolski  U, Ram  Y, Hadany  L. Key issues review: evolution on rugged adaptive landscapes. Rep Prog Phys. 2018:81 (1 ):012602. 10.1088/1361-6633/aa94d4.29051394
Oshima  T, Wada  C, Kawagoe  Y, Ara  T, Maeda  M, Masuda  Y, Hiraga  S, Mori  H. Genome-wide analysis of deoxyadenosine methyltransferase-mediated control of gene expression in Escherichia coli. Mol Microbiol. 2002:45 (3 ):673–695. 10.1046/j.1365-2958.2002.03037.x.12139615
Otto  SP, Day  T. A biologist's guide to mathematical modeling in ecology and evolution. Princeton (NJ): Princeton University Press; 2007.
Paradis  E, Schliep  K. ape 5.0: an environment for modern phylogenetics and evolutionary analyses in R. Bioinformatics. 2019:35 (3 ):526–528. 10.1093/bioinformatics/bty633.30016406
Payne  JL, Wagner  A. Mechanisms of mutational robustness in transcriptional regulation. Front Genet. 2015:6 :322. 10.3389/fgene.2015.00322.26579194
R Core Team . R: a language and environment for statistical computing. Vienna (Austria). [2023; accessed 2024 Apr]. https://www.R-project.org/.
Reynolds  R, Bermúdez-Cruz  RM, Chamberlin  MJ. Parameters affecting transcription termination by Escherichia coli RNA polymerase. J Mol Biol.  1992:224 (1 ):31–51. 10.1016/0022-2836(92)90574-4.1372365
Rifkin  SA, Houle  D, Kim  J, White  KP. A mutation accumulation assay reveals a broad capacity for rapid evolution of gene expression. Nature. 2005:438 (7065 ):220–223. 10.1038/nature04114.16281035
Ringquist  S, Shinedling  S, Barrick  D, Green  L, Binkley  J, Stormo  GD, Gold  L. Translation initiation in Escherichia coli: sequences within the ribosome-binding site. Mol Microbiol.  1992:6 (9 ):1219–1229. 10.1111/j.1365-2958.1992.tb01561.x.1375310
Rychel  K, Decker  K, Sastry  AV, Phaneuf  PV, Poudel  S, Palsson  BO. iModulonDB: a knowledgebase of microbial transcriptional regulation derived from machine learning. Nucleic Acids Res. 2021:49 (D1 ):D112–D120. 10.1093/nar/gkaa810.33045728
Sastry  AV, Gao  Y, Szubin  R, Hefner  Y, Xu  S, Kim  D, Choudhary  KS, Yang  L, King  ZA, Palsson  BO. The Escherichia coli transcriptome mostly consists of independently regulated modules. Nat Commun. 2019:10 (1 ):5536. 10.1038/s41467-019-13483-w.31797920
Siegal  ML, Bergman  A. Waddington's canalization revisited: developmental stability and evolution. Proc Natl Acad Sci U S A. 2002:99 (16 ):10528–10532. 10.1073/pnas.102303999.12082173
Silander  OK, Nikolic  N, Zaslaver  A, Bren  A, Kikoin  I, Alon  U, Ackermann  M. A genome-wide analysis of promoter-mediated phenotypic noise in Escherichia coli. PLoS Genet. 2012:8 (1 ):e1002443. 10.1371/journal.pgen.1002443.22275871
Smith  JM, Burian  R, Kauffman  S, Alberch  P, Campbell  J, Goodwin  B, Lande  R, Raup  D, Wolpert  L. Developmental constraints and evolution: a perspective from the mountain lake conference on development and evolution. Q Rev Biol.  1985:60 (3 ):265–287. 10.1086/414425.
Song  S, Zhang  J. Unbiased inference of the fitness landscape ruggedness from imprecise fitness estimates. Evolution. 2021:75 (11 ):2658–2671. 10.1111/evo.14363.34554581
Szollosi  GJ, Derenyi  I. Congruent evolution of genetic and environmental robustness in micro-RNA. Mol Biol Evol. 2009:26 (4 ):867–874. 10.1093/molbev/msp008.19168567
Taniguchi  Y, Choi  PJ, Li  GW, Chen  H, Babu  M, Hearn  J, Emili  A, Xie  XS. Quantifying E. coli proteome and transcriptome with single-molecule sensitivity in single cells. Science. 2010:329 (5991 ):533–538. 10.1126/science.1188308.20671182
Tierrafria  VH, Rioualen  C, Salgado  H, Lara  P, Gama-Castro  S, Lally  P, Gomez-Romero  L, Pena-Loredo  P, Lopez-Almazo  AG, Alarcon-Carranza  G, et al  RegulonDB 11.0: comprehensive high-throughput datasets on transcriptional regulation in Escherichia coli K-12. Microb Genom. 2022:8 (5 ):mgen000833. 10.1099/mgen.0.000833.35584008
Tirosh  I, Weinberger  A, Carmi  M, Barkai  N. A genetic signature of interspecies variations in gene expression. Nat Genet. 2006:38 (7 ):830–834. 10.1038/ng1819.16783381
Tsuru  S, Furusawa  C. Genetic properties underlying transcriptional variability in Escherichia coli. bioRxiv 589659. 10.1101/2024.04.15.589659, 16 April 2024, preprint: not peer reviewed.
Tsuru  S, Ichinose  J, Kashiwagi  A, Ying  BW, Kaneko  K, Yomo  T. Noisy cell growth rate leads to fluctuating protein concentration in bacteria. Phys Biol. 2009:6 (3 ):036015. 10.1088/1478-3975/6/3/036015.19567940
Tsuru  S, Ishizawa  Y, Shibai  A, Takahashi  Y, Motooka  D, Nakamura  S, Yomo  T. Genomic confirmation of nutrient-dependent mutability of mutators in Escherichia coli. Genes Cells. 2015:20 (12 ):972–981. 10.1111/gtc.12300.26414389
Uller  T, Moczek  AP, Watson  RA, Brakefield  PM, Laland  KN. Developmental bias and evolution: a regulatory network perspective. Genetics. 2018:209 (4 ):949–966. 10.1534/genetics.118.300995.30049818
Vaishnav  ED, de Boer  CG, Molinet  J, Yassour  M, Fan  L, Adiconis  X, Thompson  DA, Levin  JZ, Cubillos  FA, Regev  A. The evolution, evolvability and engineering of gene regulatory DNA. Nature. 2022:603 (7901 ):455–463. 10.1038/s41586-022-04506-6.35264797
Waddington  CH . Canalization of development and the inheritance of acquired characters. Nature. 1942:150 (3811 ):563–565. 10.1038/150563a0.
Waddington  CH . The strategy of the genes. London: George Allen & Unwin; 1957.
Walker  JA . A general model of functional constraints on phenotypic evolution. Am Nat. 2007:170 (5 ):681–689. 10.1086/521957.17926290
Wickham  H. ggplot2: elegant graphics for data analysis. [2016; accessed 2024 Apr]. https://ggplot2.tidyverse.org.
Wittkopp  PJ, Haerum  BK, Clark  AG. Evolutionary changes in cis and trans gene regulation. Nature. 2004:430 (6995 ):85–88. 10.1038/nature02698.15229602
Wolf  L, Silander  OK, van Nimwegen  E. Expression noise facilitates the evolution of gene regulation. eLife. 2015:4 :1–48. 10.7554/eLife.05856.
Ying  BW, Tsuru  S, Seno  S, Matsuda  H, Yomo  T. Gene expression scaled by distance to the genome replication site. Mol Biosyst. 2014:10 (3 ):375–379. 10.1039/C3MB70254E.24336896
Yona  AH, Alm  EJ, Gore  J. Random sequences rapidly evolve into de novo promoters. Nat Commun. 2018:9 (1 ):1530. 10.1038/s41467-018-04026-w.29670097
Yu  G, Smith  DK, Zhu  H, Guan  Y, Lam  TT-Y. Ggtree: an R package for visualization and annotation of phylogenetic trees with their covariates and other associated data. Methods Ecol Evol.  2017:8 (1 ):28–36. 10.1111/2041-210X.12628.
Zaslaver  A, Bren  A, Ronen  M, Itzkovitz  S, Kikoin  I, Shavit  S, Liebermeister  W, Surette  MG, Alon  U. A comprehensive library of fluorescent transcriptional reporters for Escherichia coli. Nat Methods.  2006:3 (8 ):623–628. 10.1038/nmeth895.16862137
Zheng  W, Gianoulis  TA, Karczewski  KJ, Zhao  H, Snyder  M. Regulatory variation within and between species. Annu Rev Genomics Hum Genet. 2011:12 (1 ):327–346. 10.1146/annurev-genom-082908-150139.21721942
