Introduction

Gene regulation is critical to the development and homeostasis of living beings (Prud’homme et al., 2007). It is controlled by DNA sequences that encode motifs for transcription factors (TFs) to bind and either activate or repress transcription (Jacob and Monod, 1961; Seshasayee et al., 2011). Mutations to these motifs can modulate gene expression, leading to the evolution of new phenotypes (Fuqua et al., 2020; Prud’homme et al., 2007). Mutations can also create new regulatory sequences from non-regulatory DNA, a process called regulatory emergence (Rebeiz et al., 2011; Villar et al., 2015).

We define a motif as a DNA sequence computationally predicted to be compatible with TF binding. We specifically distinguish between a motif and a binding site, because binding sites refer to the location of where a TF binds in the genome, rather than a DNA sequence, and correspond to only a fraction of all motif locations in a genome (Avsec et al., 2021; Inukai et al., 2017; Krieger et al., 2022).

Although it is commonly thought that regulatory motifs influence both regulatory emergence and evolution, pertinent empirical evidence is sparse and conflicting. On the one hand, such evidence shows that motifs can facilitate promoter emergence and evolution. For example, recent experiments showed that a randomly synthesized DNA with a single motif for a particular TF favored the emergence of tissue-specific gene expression for that TF (Galupa et al., 2023). Likewise, transposable elements contain motifs compatible with the binding of particular TFs (Lynch et al., 2015; Villanueva-Cañas et al., 2019), and these motifs have biased gene regulatory networks towards utilizing these TFs in evolution (Lynch et al., 2015).

On the other hand, there is also evidence that regulatory motifs can impede promoter emergence and evolution. One reason for this is that different regulatory proteins can compete in their binding to shared motifs (Chauhan et al., 2022; Lagator et al., 2022). Another reason is that gene expression can be silenced when regulatory DNA contains multiple motifs for the same TF (Berman et al., 2002; Bykov et al., 2020; Loell et al., 2024; Perales et al., 2016). For example, in mice, retinal “silencers” contains clusters of motifs for the TF CRX, while retinal “enhancers” contain a mixture of CRX and other TF motifs (Loell et al., 2024). A third reason is that clusters of motifs for RNA polymerase binding can cause collisions between RNA polymerase molecules, leading to premature termination of transcription (Crampton et al., 2006; Hobson et al., 2012; Wang et al., 2023).

In prokaryotes the primary regulatory sequences are called promoters (Huerta and Collado-Vides, 2003; van Hijum et al., 2009). The canonical promoter recruits the RNA polymerase subunit, σ70 (Burgess et al., 1969; Huerta and Collado-Vides, 2003; Paget and Helmann, 2003), to transcribe target genes. These promoters consists of two promoter motifs called the -10 box and the -35 box (Huerta and Collado-Vides, 2003; Pribnow, 1975; Urtecho et al., 2019) whose sequences are very well characterized (consensus: TATAAT and TTGAAA, respectively)(Burr et al., 2000; Estrem et al., 1998; Ireland et al., 2020; Kinney et al., 2010; Urtecho et al., 2023, 2019; Warman et al., 2021). We also know how the DNA sequences of promoter motifs, as well as their distances (consensus 17±1 bp) influence gene expression (Huerta and Collado-Vides, 2003; Klein et al., 2021; Lagator et al., 2022; Mitchell et al., 2003; Tierrafría et al., 2022; Urtecho et al., 2019). Like in humans (Gotea et al., 2010) and invertebrates (Lifanov et al., 2003), most prokaryotic promoters contain multiple regulatory motifs (Huerta and Collado-Vides, 2003), and have multiple transcription start sites (Mendoza-Vargas et al., 2009). It is not clear whether these additional motifs influence promoter activity. If so, it may help to explain why it is difficult to predict promoter strength and activity (Lagator et al., 2022; Urtecho et al., 2023).

Recent work shows that promoters can emerge from promoter motifs or adjacent to them (Fuqua and Wagner, 2023). However, encoding -10 and -35 box motifs is insufficient to drive transcription. For instance, the E. coli genome contains clusters of -10 and -35 box motifs with little or no promoter activity. Such clusters are called promoter islands (Bykov et al., 2020; Panyukov and Ozoline, 2013; Purtov et al., 2014; Shavkunov et al., 2009). One study found that promoter activity can either emerge from a promoter island when a -35 box motif is mutated to match the consensus sequence, or counterintuitively, when -10 and -35 box motifs are destroyed (Bykov et al., 2020).

Thus, it is neither understood how additional promoter motifs influence promoter activity, nor is it clear when a promoter motif will facilitate or impede de-novo promoter emergence.

Here, we systematically study the ability of -10 and -35 box motifs to influence the evolution and emergence of promoters. To this end, we select 25 fragments from 25 promoter island sequences, with both regulatory and non-regulatory activity, and screen more than 240’000 mutant variants for promoter activity. We find that ∼67% of all promoter activity originates from -10 and -35 box motifs, but the likelihood of mutations creating new promoters does not correlate with the number of motifs. Additionally, mutations can create additional -10 and -35 box motifs to either increase or decrease promoter expression. These results highlight the complex interplay of promoter motifs throughout regulatory evolution and the emergence of new phenotypes.

Results

Parent sequences are enriched with motifs for -10 and -35 boxes

Promoter islands are ∼1kbps (kilo base pairs) long (Bykov et al., 2020; Panyukov and Ozoline, 2013; Purtov et al., 2014; Shavkunov et al., 2009), but for adequate mutational coverage, it is better to work with shorter sequences. We thus identified within 25 promoter islands, a region of 150 bp that showed the highest density of predicted -10 and -35 box motifs for further analysis (Fig 1A). To this end, we used position-weight matrices (PWMs) – computational tools to identify promoter motifs (Methods) (Hertz and Stormo, 1999). We refer to these 25 sequences as parent sequences (P1-25). On average, each parent contains ∼4.22-10 box motifs and ∼5.14 -35 box motifs on each of its two DNA strands (Fig S1).

Figure 1.

Inactive promoters vary widely in their potential to become active promoters

Next we created a mutant library of daughter sequences from the 25 parent sequences using an error-prone polymerase chain reaction (Fig 1B), and measured the extent to which these daughter sequences can drive gene expression using the well-established sort-seq procedure (Kinney et al., 2010; Peterman and Levine, 2016). Specifically, we cloned the library into a bidirectional fluorescence reporter plasmid (pMR1) to measure promoter activity on both DNA strands (top: GFP, bottom: RFP) (Westmann et al., 2018). We then transformed the plasmid library into E. coli (DH5α, Takara Japan), sorted the bacteria using fluorescence activated cell sorting (FACS) into eight fluorescence bins (none, weak, moderate, and strong fluorescence, for both GFP and RFP) (Fig 1C, see also Fig S2), and sequenced the inserts from the sorted fluorescence bins (see methods).

Our library comprised 245’639 unique daughter sequences, with an average of 4’913 daughter sequences per parent and orientation (top/bottom) and ∼2.5-point mutations per daughter (Fig S2). Each parent has 3×150 mutational neighbors that differ by a single mutation. We define coverage as the percentage of these 450 bases and positions that are present in at least one daughter sequence for each parent. (Daughters contain 1-10 point mutations.) The mean coverage was 93% (minimum: 80%; maximum: 100%) (Fig S3). The mutation library therefore encompasses the majority of neighboring mutations for the 25 parents.

We next quantified how strongly each daughter sequence could drive gene expression from the top and bottom strand by calculating a fluorescence score for both GFP and RFP expression, respectively. These scores are bounded between a lowest score of 1.0 arbitrary units (a.u.), indicating no expression, and a highest score of 4.0 a.u., indicating maximal expression (methods). We declared a daughter sequence to have an active promoter if its score exceeded 1.5 a.u., as this score lies at the boundary between no fluorescence and weak fluorescence based on the sort-seq bins (methods).

Some sequences in our library are unmutated parent sequences, which allowed us to determine that 8/25 of the parent sequences already encode active promoters before mutagenesis. Specifically, one parent drove expression on the top strand (P22), five did on the bottom strand (P6, P12, P13, P18, P21), and two drove expression on both strands (P19 and P20). See Fig S4 for fluorescence score distributions for each parent and its daughters, and Data S2 for all daughter sequence fluorescence scores.

We refer to the remaining 40 parent sequences and orientations as inactive. Although mutating each of the 40 inactive parent sequences could create a promoter, the likelihood P_new that a mutation creates an active promoter varies dramatically among parents. For each inactive parent and orientation, Fig 1D shows the percentage of active daughter sequences. The median P_new is 0.046 (std. ± 0.078), meaning that ∼4.6% of all mutants have promoter activity. The lowest P_new is 0.002 (P25) and the highest 0.41 (PI7), a 205-fold difference.

These P_new values do not significantly correlate with the number of daughter sequences, the GC-contents of the parents, and – most notably – with the number of -10 or -35 box motifs (least squares linear regression, p=0.529, 0.930, 0.852, 0.834 respectively) (Fig S5). Thus, the number of promoter motifs does not directly influence the likelihood of promoter emergence.

Promoters primarily emerge and evolve from specific subsets of -10 and -35 box motifs

Next we asked where both active and new promoters are located within the parent sequences. To this end, we first calculated for each parent and its respective daughters, the mutual information I_i(b, f) between the nucleotide identity at each position i and the corresponding fluorescence level f (Fig 2A). The essence of this calculation is to determine for every position i , the probability p_i(b) of that position being an A,T,C, or G and the probability of each sequence having a (rounded) fluorescence score of 1,2,3, or 4 a.u. p(f). Representing these probabilities as circles in a Venn-diagram, the overlap between the circles is the joint probability p_i(b, f) that position i has a specific base b AND a particular fluorescence score f. The larger this joint probability is compared to the individual probabilities, the more important the nucleotide identity at position i is for fluorescence, i.e., for promoter activity. This calculation can reveal positions where RNA polymerase and transcription factors bind (Belliveau et al., 2018; Ireland et al., 2020; Kinney et al., 2010), as well as DNA regions where promoters most readily emerge (Fuqua and Wagner, 2023). The calculation effectively converts all 245’639 data points into maps that reveal where active promoters either already exist or emerge in each of the parent sequences (see methods).

Figure 2.

We identified ‘hotspots’ – regions of especially high mutual information – for each parent sequence (Fuqua and Wagner, 2023) (see methods). For example, two such hotspots exist in parent P19, which is an active promoter that drives expression from both DNA strands (Fig 2B). One hotspot overlaps with a -35 box motif and the other overlaps with a -10 box motif. These hotspots indicate that fluorescence changes when these DNA sequences mutate. Because these sequences additionally match a canonical promoter, they are the likely cause of the promoter activity in P19. Additionally, the number of hotspots that a sequence encodes also correlates modestly with the probability P_new (least squares linear regression, p=1.35×10^{- 4}, r² = 0.322) (Fig S5). In other words, the more hotspots a parent has, the more likely it is that mutations create a new promoter from it.

Across all 25 parent sequences and for both expression directions (top and bottom strand), we identified a total of 68 hotspots (Fig 2C). 29 hotspots occur in the parent sequences with promoter activity. They correspond to locations where mutations either increase or decrease promoter activity. ∼17% (5/29) of these hotspots overlap exclusively with -10 box motifs, ∼28% (8/29) with -35 box motifs, and ∼21% (6/29) overlap with both -10 and -35 box motifs. The remaining ∼34% of hotspots overlap with neither a -10 nor a -35 box motif (10/29).

The remaining 39 hotspots occur in inactive parent sequences. Mutations in these locations solely increase fluorescence. Of these 39 hotspots, ∼23% (9/39) overlap with -10 box motifs, ∼31% (12/39) with -35 box motifs, ∼13% with both a -10 and -35 box motifs (5/39), and ∼33% with neither motif (13/39) (Fig 2C-E). See Fig S6 for the mutual information plots for all parent sequences and Data S3 for a table with the hotspots and their coordinates.

Despite this extensive overlap between promoter motifs and hotspots, each parent contains additional -10 and -35 box motifs that do not overlap with hotspots. Of the 211 -10 box motifs, only 25 (∼12%) overlap with a hotspot. Similarly, of the 257 -35 box motifs, only 31 (∼12%) overlap with a hotspot. For example, the inactive parent P3 contains five -10 box motifs, but only one overlaps with a hotspot (Fig 2D). Similarly, the inactive P18 contains seven -10 and -35 box motifs each, but promoters only emerge in three hotspots (Fig 2E). In sum, the majority (∼67%) of new promoters evolve and emerge from motifs for -10 and -35 boxes, but the presence of a motif does not indicate that promoter activity can emerge or is encoded within the motif.

Promoters can emerge when mutations create motifs but not by destroying them

We hypothesized that promoters emerge when mutations create -10 and -35 box motifs. To test this hypothesis, we examined each hotspot and used PWMs to find out whether a subset of the daughter sequences gains a -10 or -35 box motif in the hotspot. We then asked whether gaining a -10 or -35 box motif leads to a significant increase in fluorescence that indicates the creation of a new promoter (Mann-Whitney U-test, Fig. S7 and methods)

We identified 5 and 4 hotspots in which gaining a -10 box motif or a -35 box motif creates a de-novo promoter by increasing fluorescence by more than +0.5 a.u. (Fig 3A). For ∼67% (6/9) of these hotspots, the locations of the newly created -10 and -35 box motifs overlap with preexisting -10 and -35 box motifs.

Figure 3.

As examples, we highlight two hotspots involving the creation of a -10 box. Parent P16 gains promoter activity when mutations create a -10 box downstream of a preexisting -35 box motif (Fig 3B), increasing median fluorescence by +1.99 a.u. (Fig 3C) (two-tailed MWU test, q=7.64×10^-32). We call the act of mutations creating a -10 box downstream of a preexisting -35 box motif Shiko emergence (Fuqua and Wagner, 2023). The term is a homage to the Sumo exercise Shiko, where the wrestler firmly plants their foot firmly on the floor, and lowers their other foot to a firm position. In contrast, parent P9 gains promoter activity when mutation creates a new -10 box that does not lie downstream of an existing -35 box motif (Fig 3D,E, median fluorescence gain +1.8 a.u., two-tailed MWU test, q=2.83×10^-173). See Fig S8 for the other hotspots in which -10 boxes are created.

We also highlight one hotspot where gaining a -35 box creates de-novo promoter activity. Specifically, parent P1 gains promoter activity when mutations create a -35 box upstream of a preexisting -10 box motif (Fig 3F), which increases fluorescence by +0.74 a.u. (Fig 3G, two-tailed MWU test, q=4.25x10^-42). This is also an example of Shiko Emergence (Shiko -35). See Fig S8 for the other hotspots in which -35 boxes are created. We did not find any hotspots where a -35 box is created without a downstream -10 box motif.

Prompted by an earlier report (Bykov et al., 2020), we also asked whether mutations that destroy -10 or -35 box motifs can increase promoter activity. No such mutations exist in our data (see Fig 3A).

To summarize, we found 9 hotspots where mutations created -10 and -35 boxes that lead to de-novo promoter emergence. When -10 boxes are gained, the majority (three of five) form downstream of a preexisting -35 box motif (Shiko -10 emergence). The remaining two -10 boxes do not form near a -35 box motif (gain -10, Fig 3H). Four of the four newly formed -35 boxes form upstream of a -10 box motif (Shiko -35). We did not find any hotspots where creating a -35 box without a downstream -10 box motif leads to promoter emergence. In the majority of the 9 hotspots (6 of 9), new -10 or -35 boxes overlap preexisting -10 and -35 box motifs.

Gaining boxes can increase or decrease promoter activity (by more than one mechanism)

We next asked how gaining and losing -10 and -35 box motifs can change expression driven by active parent sequences. To this end, we repeated the previous analysis but for the active parents.

We first identified 5 and 3 hotspots in which losing a -10 box motif or -35 box motif decreases expression (Fig S9A). These particular hotspots revealed where promoters are located in each parent. For example, for parent P12 (Fig S9), losing a -35 box decreases fluorescence by -1.57 a.u. (Fig S9) (two-tailed MWU test, q=6.28×10^-67), and losing a -10 box decreases fluorescence by -1.41 a.u. (two-tailed MWU test, q=6.28×10^-91). Based on these data and the close spacing (18 bp) of these two motifs, we can infer that they constitute the binding sites for the promoter in P12 (Fig 4A). See Fig. S9 for additional examples. We did not find any hotspots in which destroying a -10 box motif or a -35 box motif increases expression (see Fig S9A).

Figure 4.

We also found 4 and 1 hotspot(s) in which gaining a -10 or -35 box further increases expression. More surprisingly, in 4 and 6 hotspots gaining a -10 or -35 box decreases expression (n= 4+1+4+6, 15 hotspots with increased or decreased expression) (see Fig S9A). To understand how this is possible, we analyzed each parent sequence where this occurs in detail. This analysis allowed us to group such events into three categories. That is, mutation can 1) create a tandem motif that lowers expression, 2) shift -10 and -35 boxes to modulate expression, or 3) create additional -10 and -35 box motifs that destroy the original promoters and lower expression (Fig 4). Next we briefly describe examples of each category.

For 2 of the 15 hotspots (∼13%), mutations create a tandem motif, where two -10 or -35 box motifs occur in head-to-tail orientation. For example, when mutations of the aforementioned parent sequence P12 create a -35 box motif upstream of the existing -35 box, median fluorescence decreases by 1.1 a.u. (two-tailed MWU test, q=3.35×10^-5) (Fig 4B). Similarly, when mutations in P12 create a -10 box motif upstream of the mapped -10 box, median fluorescence also decreases by 0.61 a.u. (two-tailed MWU test, q=4.36×10^-4) (Fig 4C).

For 6 of the 15 hotspots (40%), mutations shift the location of -10 and -35 boxes. For example, mutations can shift the -10 box of parent P12 +1 bp away from the -35 box, which increases median expression by +0.78 a.u. (Fig 4D) (two-tailed MWU test, q=5.85×10^-13). See Fig S10 for additional hotspots where such shifting occurs (Fig 4E).

For 3 of the 15 hotspots (20%), we find that mutations can create -10 or -35 box motifs overlapping -10 and -35 boxes from the putative promoters we mapped in Figure S9. For example, for parent P22, where our mutational data had allowed us to map the active -35 and -10 boxes (Fig 4F), mutations can change the -35 box into a -10 box motif, which decreases expression by 0.90 a.u. (Fig 4G) (two-tailed MWU test, q=7.29×10^-7). Thus, while the loss in expression appears to be caused by a gain of a motif, this gain also entails of a loss of an existing, active box (Fig 4H). See Fig S10 for additional examples.

The remaining 4 of 15 hotspots (∼27%) are located in P20 (Fig. S11), which encodes a bidirectional promoter. However, because our mutational data did not allow us to confidently map this promoter on either strand of P20, we cannot resolve why gaining -10 or -35 box motifs may increase expression (in 2 of the 4 P20 hotspots) or decrease expression (in the remaining 2 hotspots, Fig. S11) of this parent sequence.

Discussion

We found that sequences enriched with -10 and -35 box motifs without promoter activity vary dramatically (more than 200-fold) in the probability P_new that mutations create active promoters in them (P_new = 0.002 to 0.41, see Fig 1D). When trying to understand the reasons for this variation, we found that P_new is not significantly correlated with the number of existing -10 or -35 box motifs (least squares linear regression, p=0.852, 0.834 respectively). In addition, hotspots of promoter emergence coincide with only ∼12% of the -10 and -35 box motifs. However, P_new correlates with the number of mutual information hotspots (least squares linear regression, p=1.35×10^-4, r² = 0.322, see Fig S5). For instance, three parental sequences with high P_new (P8, P9, and P16, see Fig 1D) encode multiple hotspots (3,3, and 5 hotspots respectively, see Fig S6). In each of these hotspots, promoters emerge when mutations create new -10 and -35 boxes (see Fig 3 and Fig S8). Sequences with the lowest emergence probability, such as the bottom strand of P25, contain no hotspots (Fig S6). Therefore, promoter emergence is best explained by how readily mutations create additional promoter motifs, not by the number of preexisting -10 and -35 box motifs. Additionally, ∼78% (7/9 hotspots) of emergence hotspots, gain new -10 and -35 boxes adjacent to corresponding -10 and -35 box motifs (Shiko Emergence, see Fig 3H), suggesting that the location of new and preexisting motifs also influences promoter emergence.

-10 and -35 box motifs, however, influence where promoters emerge. We found that the majority (∼67% ) of all mutual information hotspots overlap with -10 or -35 box motifs (see Fig 2C). This can partially be explained by our finding that mutations create additional -10 and -35 boxes upstream or downstream of -35 and -10 box motifs respectively, to build promoters (see Fig 3). Additionally, mutations within the promoter motifs can shift their locations, change a -35 box into a -10 box motif, or vice versa (see Fig 3 and Fig S8). This finding may help to explain why a previous study concluded that destroying promoter motifs creates promoter activity (Bykov et al., 2020). Specifically, destroying a -10 or -35 box motif can create a de-novo promoter, but only if this destruction simultaneously creates a new -10 or -35 box (see Fig 3A).

When asking how mutations change the activity of already active parent sequences, we found that their ability to shift the locations of -10 and -35 boxes is especially important. For example, in 6 hotspots, mutations that shift -10 and -35 boxes either increase or decrease expression (see Fig 4E). This shifting also explains why ∼67% of hotspots in active sequences overlap with existing -10 and -35 box motifs. This finding is consistent with previous studies, which demonstrate that the distance between -10 and -35 boxes is an integral component of promoter strength, and that deviations from the 17±1 bp spacing can result in lower expression (Klein et al., 2021; Lagator et al., 2022; Urtecho et al., 2019).

Promoter activity can be further modulated by the creation of tandem motifs in parent P12. Such tandem motifs do not abolish gene expression but decrease it by about one half (Fig 4B and Fig 4C). One possible explanation for this halving is that RNA polymerase can bind to both parts of the tandem motif (i.e. both -10 boxes or both -35 boxes), but only binding to one of these two yields promoter activity. To understand the biophysical mechanism behind this phenomenon is a task for future work. Regardless of this mechanism, the observation underlines that DNA sequences flanking a promoter can be important for protein binding (Mitchell et al., 2003), as has also been observed in eukaryotic transcriptional regulation (Afek et al., 2014; Chiu et al., 2022; Farley et al., 2015; Gordân et al., 2013). Collectively, the shifting of -10 and -35 boxes, and the creation of tandem motifs demonstrate how mutations can readily create promoters with a wide dynamic range of activity. This ability facilitates natural selection’s fine-tuning of a promoter to precise expression levels, within a small sequence search space that is easily accessible by few DNA mutations.

Overall, our study demonstrates that mutations to -10 and -35 box motifs are the primary paths to create new promoters and to modulate the activity of existing promoters. On the one hand, this is exciting for studying de-novo genes and the evolution of novel gene expression, because it suggests that promoter emergence and evolution may be predictable (Lässig et al., 2017; Stern and Orgogozo, 2008). On the other hand, our observation that promoter motifs influence where new regulatory sequences can emerge, leads to many exciting insights about regulatory emergence throughout the genome, plasmids, and transposable elements.

For instance, new enhancer sequences in eukaryotes primarily emerged from specific families of transposable elements (Villar et al., 2015; Yue et al., 2014). Furthermore, transposons encoding particular motifs influence which TFs are utilized in a gene regulatory network (Lynch et al., 2015). Collectively, observations like these suggest that the identity of a DNA sequence heavily influences how gene regulation emerges and evolves, and thereby which novel phenotypes are most likely to arise by chance.

Data availability

Data S1 is a list of the DNA sequences for the primers and promoter island parent sequences in this study. Data stored as an Excel spreadsheet.

Data S2 is a dataframe with all of the daughter sequences and their respective fluorescence scores after cleaning the dataset. Data stored as a csv file.

Data S3 is a dataframe with each identified hotspot and whether it overlaps with a -10 or -35 box. Data stored as a csv file.

Data S4 is a dataframe with the output from the gain-loss analysis for both the top and bottom strand. Data stored as a csv file.

Data S5 is an Excel spreadsheet to rapidly reproduce the main figures (to the best of Excel’s capabilities).

A Jupyter notebook with all of the python scripts to recreate the analysis and figures in this study, as well as Data S1-S5 are available at the Github repository: https://github.com/tfuqua95/promoter_islands

The sequence reads are available at the Sequence Read Archive (SRA) with accession number: PRJNA1071572: https://www.ncbi.nlm.nih.gov/bioproject/1071572

Methods

Position weight matrices (PWMs)

To create PWMs for -10 and -35 box motifs, we obtained a list of -10 and -35 box sequences from Regulon DB (Tierrafría et al., 2022), and converted them into PWMs using Biopython.motifs from the Biopython package (Cock et al., 2009) (v1.79). We deliberately assumed a uniform nucleotide composition (25% A, T, C, or G) in this process, because doing so allowed us to test PWM prediction performance in a variety of different sequence backgrounds. For every input sequence, a PWM returns an output score in bits. The higher this score is, the more likely it is that the input sequence can be bound by the protein the PWM represents.

To classify whether an output score is strong enough to be classified as a -10 or -35 box motif, we used the standard practice of setting a threshold equal to the information content of the PWM (Hertz and Stormo, 1999). The -10 box in our study has an information content of 3.98 bits and the -35 box of 3.39 bits. We classified input sequences with scores at least this high as -10 or -35 box motifs, respectively. In Fig S9 we also plot “low-affinity” -10 and -35 box motifs. For these motifs, we halved the thresholds (1.99 and 1.70 bits, respectively).

Amplifying parent sequences and plasmid MR1 with PCR

We acquired the genomic coordinates for the promoter island sequences from (Panyukov and Ozoline, 2013). We used a polymerase chain (PCR) reaction to amplify 78 parent sequences from the genome, using Q5 polymerase (NEB, USA product #M0491). The reaction amplifies each parent sequence from the genome, and additionally concatenates short sequences to the 5’ and 3’ ends of the inserts (5’-GGCTGAATTC…insert…GGATCCTTGC-3’). These overhangs allowed us to carry out the error-prone PCR with a single set of primers. See Data S1 for the list of primers and Fig S12 for an overview of the PCR.

To amplify the parent sequences from the genome, we performed each PCR in a 50 uL reaction volume with the following reagents: 1 uL of each primer (2× primers, 2uL total, concentration 100 uMol), 5 uL Q5 reaction buffer, 1 uL 10 mM dNTPs (Thermo Scientific, USA, product #R0191), 1 uL of the genomic DNA extracted using a Wizard Genomic DNA Purification Kit according to the manufacturer’s instructions (A1120, Promega USA), 1 uL Q5 high-fidelity DNA polymerase, and 40 uL of water. We carried out the PCR in a thermal cycler (C1000 Touch Thermal Cycler, Bio-Rad, USA) with 30 cycles at 55°C for annealing and 72°C for primer extension, both for 30 seconds. To purify the PCR products, we separated them on a 1% agarose gel, excised them, and used a Qiagen QIAquick Gel Purification Kit (Qiagen, Germany, product #28706) following the manufacturer’s instructions. To verify the purification’s success and estimate the final concentrations of each product, we used a Nanodrop One spectrophotometer (Thermo Scientific, USA).

To create the mutagenesis library of parent sequences, we pooled the 78 PCR products such that each product had the same concentration in the pool. We then performed a single error-prone polymerase chain reaction on the pooled products to create a mutagenesis library. To this end, we used GoTaq polymerase (Promega, USA, product #M3001). Specifically, we combined 1 uL of two PCR primers complementary to both the pMR1 plasmid and to the sequences concatenated to the ends of the initial PCR (1 uL forward and 1 uL reverse primer, each at a 100 uMol l^-1) (see Data S1 and Fig S12), 10 uL of GoTaq reaction buffer, 1 uL 10 mM dNTPs (Thermo Scientific, USA, product #R0191), 1 uL of pooled DNA, and 1 uL of GoTaq polymerase. We added water to a total of 98 uL, then added 2 uL of 15 mMol MnCl₂ to the reaction. We then excised and purified the product as described for the parent sequences.

We additionally used Q5 polymerase to create and amplify linear copies of pMR1 (see Data S1 for primer sequences and Fig S12). We carried out the necessary PCR as described for the parent sequences but for an extension time of 2 minutes and 30 seconds.

Molecular cloning

We cloned the error-prone PCR library into pMR1 at the EcoRI and BamHI sites using NEBuilder (New England Biolabs, USA, product #E2621) as follows: 100 ng of linear pMR1, 25 ng insert, 5 uL NEBuilder mastermix, and water to 10 uL total volume. To increase the reaction efficiency, we added 1 uL of T4 DNA ligase buffer (New England Biolabs, USA, product #B0202S). Using a thermal cycler (C1000 Touch Thermal Cycler, Bio-Rad, USA) we incubated each assembly reaction for 1 hour at 50°C.

Subsequently, we electroporated our product directly into E. coli DH5α electrocompetent cells (Takara, Japan, product #9027) by adding 2 uL of (unpurified) product to 2× 50 uL aliquots of E. coli (100 uL total) in a 2 mm electroporation cuvette (Cell Projects, England, product #EP-202), and electroporated the cells using a Bio-Rad MicroPulser (Bio-Rad, USA). For recovery, we then added 1 mL of Super Optimal Broth with Catabolite Repression Medium (SOC, provided with the electrocompetent cells (Takara, Japan, product #9027), shaking the cells at 230 RPM at 37°C for 1.5 hours (Infors HT, Switzerland, Multitron).

We then plated 5 uL of the bacterial culture onto a petri dish containing LB-agar supplemented with 100 ug/ml of chloramphenicol, and incubated the plate at 37°C overnight. In parallel, we added 9 mL of LB-chloramphenicol (100 ug/mL) to the remaining 995 uL of bacteria (in SOC), and incubated the resulting culture overnight. The following morning, we manually counted colonies on one quarter of the plate, and estimated the total number of clones (cells) in the original electroporated culture by multiplying this value by 4×200.

Sort-seq controls

We synthesized three control promoter sequences (Data S1) as fiducial markers for defining the fluorescence sorting gates for sort-seq. The insert DNA sequences were chemically synthesized (Integrated DNA Technologies, USA). We then amplified these promoters, cloned them into plasmid MR1, and transformed E.coli with them as described in “Molecular Cloning.” The first is the promoter iGem bba_j23110 on the top DNA strand, which drives moderate GFP expression from plasmid pMR1. The second is the reverse complement of this promoter, which drives moderate RFP expression from the bottom DNA strand of pMR1. The third and final (negative) control is the wild-type pMR1 plasmid without any insert in its multiple cloning sequence.

Sort-Seq

To perform sort-seq, we first added 100 uL of the glycerol stock of our mutagenesis library to 1 mL of LB-chloramphenicol (100 ug/ml). We then allowed the resulting culture to shake at 230 RPM and 37°C overnight (∼16 hours) on an Infors HT (Multitron, Switzerland) shaking incubator. Subsequently, we centrifuged the culture and washed the cell pellet twice with Dulbeco’s Phosphate Buffered Saline (PBS) (Sigma, USA, D8537). We then resuspended the pellet in PBS.

To detect fluorescence caused by GFP and RFP expression, we measured green fluorescence using a 488 nm laser (FITC-H) at 750 volts, and red fluorescence with a 633 nm laser (PE-H) at 510 volts on an Aria III cell sorter (BD Biosciences, USA). We used the following fluorescence gates for sort-seq (Figure S2).

RFP bin #1 (none, i.e., no fluorescence): We defined a minimum boundary to not sort cell debris, salts, empty droplets, and other impurities, as the larger of two PE-H values. These values are the minimum PE-H of the negative control (empty pMR1 plasmid) and the minimum PE-H of the opposite fluorophore positive control (GFP in this case). The upper boundary of this bin is the highest PE-H value detected for these same controls.

RFP bin #4 (high): We defined the lower boundary of this bin using the mean fluorescence level of the positive RFP control. This bin does not have an upper bound, because it encompasses the highest levels of fluorescence.

RFP bins #2 and #3 (low and moderate): the lower bound of bin #2 begins at the upper bound of bin #1. The upper bound of bin #3 ends at the lower bound of bin #4. The upper bound of bin #2, which is identical to the lower bound of bin #3, equals the average of the lower boundary of bin #4 and the upper bound of bin #1.

We defined the bins for GFP analogously, but with GFP and RFP controls swapped, with the following exception: Because the GFP positive control produces a bimodal FITC-H distribution, we defined the lower bound of bin #4 for green fluorescence as the peak of the higher mode of this distribution. See Fig S2 for an illustration of the bins.

We sorted the library using two consecutive rounds of sorting on consecutive days. On the first day, we sorted 2’160’000, 71’545, 28’399, and 22’720 cells into GFP bins none, low, moderate, and high, respectively; and 2’160’000, 540’000, 101’636, and 76’898 cells into RFP bins none, low, moderate, and high, respectively. After sorting, we added 1 mL of SOC medium (Sigma, USA, product #CMR002K) to the sorted cultures, incubating the cells at 37°C and shaking for two hours at 230 RPM (Infors HT, Switzerland, Multitron). Following this recovery period, we added 9 mL of LB-Chloramphenicol (100 ug/ml chloramphenicol) to each culture and incubated them overnight (∼16 hours) under the same conditions.

On day two, we once again washed the cultures and re-sorted cells in each culture into the same fluorescence bin. That is, we sorted cells from GFP-none into GFP-none, GFP-low into GFP-low, and so forth. This resorting step ensures that each genotype fluoresces within the same boundaries as on the previous day, and reduces the likelihood of propagating sorting errors into subsequent analyses. The resulting number of cells for each bin were as follows: GFP-none: 6’480’000, GFP-low: 216’000, GFP-moderate: 90’000, GFP-high: 69’000, RFP-none: 6’480’000, RFP-low: 540’000, RFP-moderate: 306’000, RFP-high 231’000 cells. We allowed the cells to recover as described at the end of day one.

After ∼16 hours of overnight culture, we created glycerol stocks for each culture (8 cultures total) as described in “Molecular Cloning.” From the remaining culture volume, we isolated the plasmids using a Qiagen QIAprep Spin Miniprep kit (Qiagen, Germany, product #27104) as described by the manufacturer. We amplified the plasmid inserts using PCR as described in “Molecular Cloning”, using primers (Data S1) that add a unique multiplexing barcode onto the 5’-end of each PCR product from each bin. In addition, we isolated the inserts from the unsorted library, and endowed them too with a unique barcode. This procedure yielded 9 (=8+1) sequence samples, each with a unique barcode. We then pooled equimolar amounts of these samples, and sequenced them using Illumina paired-end sequencing (Eurofins GmbH, Germany).

Processing sequencing reads

Using Flash2 (Magoč and Salzberg, 2011), we merged the paired-end reads that resulted from Illumina sequencing. Each paired-end read contains palindromic cut-sites for 5’-EcoRI (GAATTC) and 3’-BamHI (GGATCC). To orient all of the paired-end reads in the same direction, we searched each paired-end read for both cut-sites, discarding any sequence that did not encode both, and took the reverse complement of the read if the BamHI site was upstream of the EcoRI site. We then identified the tagmented barcode sequence upstream of the EcoRI site, classifying from which fluorescence bin the read originated. We then cropped the sequences to remove the barcodes, BamHI, and EcoRI sites, and counted the total number of reads within every bin. We refer to each of these unique paired-end read sequences as daughter sequences.

To focus on point mutations rather than insertions and deletions, we removed daughter sequences with a length different from the wild-type 150 bp. To minimize the influence of sequencing errors in our dataset, we removed daughter sequences that did not occur at least once in the unsorted library, and daughter sequences with less than thirty reads in the entire dataset.

To map each daughter sequence to its respective template parent sequence, we calculated the Hamming distance (i.e. the number of nucleotide differences between two sequences) between each daughter sequence and the 78 parent sequences, and identified the parent as the sequence with the smallest Hamming distance. To minimize our likelihood of mapping the daughter sequences to the wrong parent sequence, we removed daughter sequences from the dataset for which the closest Hamming distance was greater than 10. Because the error-prone PCR created different numbers of daughter sequences per parent sequence, we also removed all parent and their respective daughter sequences if there were less than 2’000 unique daughter sequences. After these calculations and filtration steps, our dataset contained 245’639 daughter sequences derived from 25 parent sequences.

To calculate fluorescence scores (F ) we used in the next analysis step, we calculated (F ) using Equation (1):

Here, f is an integer index corresponding to each of the four fluorescence bins, i.e., f = 1,2,3,4 corresponds to no fluorescence, low, moderate, and high fluorescence, respectively. Reads_f is the number of times the sequence was identified in each fluorescence bin f. We calculated F for both green and red fluorescence. Each score is in arbitrary units (a.u.) of fluorescence.

See Data S2 for a csv file with the daughter sequences and their respective fluorescence values, which we analyzed further with the code on the GitHub repository: https://github.com/tfuqua95/promoter_islands.

Mutual information

To calculate mutual information between nucleotides and fluorescence scores, we first randomly split the dataset into three equal-sized pseudo-replicates (r1, r2, and r3) to minimize the effects of sampling biases as well as sorting errors when calculating mutual information. For each pseudo-replicate, we calculated the mutual information I_i(b, f) between the nucleotide identity b and fluorescence score f score at each position i (1<i<150), of every daughter sequence using equation (2), as defined in ref. (Kinney et al., 2010):

In the left term on the right-hand side of equation (2), b indexes each nucleotide base (b = A,T,C,G) and f denotes the fluorescence score (see Processing sequencing reads) rounded to the nearest whole number ( f=1,2,3,4 ). The expression p_i(b) denotes the relative frequency of base b (A, T, C, or G) at position i . p(f)) is the relative frequency of a daughter sequence with a (rounded) fluorescence score of 1,2,3, or 4. p_i(b, f) is the observed joint probability, i.e., the relative frequency of a daughter sequence having a score of f and base b at position i.

The right term on the right-hand side of equation (2) is a correction term for differences in degrees of freedom (9 degrees of freedom), where n_b = 4 is the number of possible bases n_f = 4 is the number of possible (rounded) fluorescence scor es, and N is the size of the library. The fewer degrees of freedom there are and the larger the library is, the smaller the correction term.

For the mutual information plots shown in the text, we computed a Gaussian filter with the ndimage (alpha = 1) function in scipy (v1.8.1). In these plots, the solid line is always the average mutual information between r1, r2, and r3 at each position, and the lightly shaded region indicates ± 1 standard deviation between r1, r2, and r3 at each position.

Identifying mutual information hotspots

We calculated the mutual information for each parent and its daughters (see Mutual Information), and used the following procedure to automatically identify mutual information hotspots. First, we subtracted the standard deviation from the mean mutual information at each nucleotide position i, and truncated all values smaller than 0.0005 bits to 0 bits. We refer to the resulting quantity as adjusted mutual information. We then established the 90^th percentile of this adjusted mutual information at each nucleotide position i as a threshold to identify mutual information hotspots. Specifically, we identified mutual information hotspots as nucleotide positions where the adjusted mutual information 1) lies above this threshold, and 2) is greater than the mutual information at the positions i-1 and i+1. See Data S3 for hotspot coordinates.

Analysis of motif gain and loss

To identify hotspots in which gaining and losing a -10 or -35 box associates with a significant change in gene expression (fluorescence), we performed the following procedure. First, we drew a 6 bp sliding window starting at the beginning of each parent sequence. Within this window, we calculated the respective -10 and -35 PWM score for every daughter sequence (see Fig S7A). For each daughter sequence, if the binding score was greater than a designated threshold (-35 box = 3.39 bits, -10 box = 3.98 bits, see “Position Weight Matrices”), we added the fluorescence score to a “Motif (+)” list. If the PWM score was below this threshold, we added the score to a “Motif (-)” list (see Fig S7B).

If each list contained more than 10 values, we then tested the null hypothesis that the fluorescence scores have the same central tendency, using a two-sided Mann-Whitney U test (mannwhitneyu function from scipy.stats, v1.8.1). If the p-value was less than 0.05 and was located within ±3 bp from a hotspot, we recorded the value in Data S4. We then advanced the 6 bp window in steps of one nucleotide position until it had reached the final end of the parent sequence, and repeated this calculation at each location of the sliding window. See Fig S7.

We performed this analysis for both GFP and RFP fluorescence scores on their respective DNA strands for both -10 and -35 boxes and for all 25 parent sequences (Data S4). To account for multiple hypothesis testing, we calculated corrected q-values using a Benjamini-Hochberg correction (false-discovery rate 0.05).

For subsequent analyses, we only focused on significant associations found within ± 3 bp from the peak of each hotspot (see “Identifying promoter emergence hotspots” and “Mutual information”), and associations where the median fluorescence difference between the lists box present and box absent is larger than the absolute value of 0.5 arbitrary units (a.u.).

Acknowledgements

The European Molecular Biology Organization (EMBO) supports TF with a postdoctoral fellowship (ALTF 963-2021). AW and his group is supported by the Swiss National Science Foundation (grants 310030_208174).

We would like to thank the Wagner lab for all discussions and conversations, even if they were not always pertinent to science. We additionally thank Baxter for inspiring Tim with the analysis pipeline during a painfully slow and cold morning walk.