PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of ploscompComputational BiologyView this ArticleSubmit to PLoSGet E-mail AlertsContact UsPublic Library of Science (PLoS)
 
PLoS Comput Biol. 2010 July; 6(7): e1000848.
Published online 2010 July 8. doi:  10.1371/journal.pcbi.1000848
PMCID: PMC2900288

Redundancy and the Evolution of Cis-Regulatory Element Multiplicity

Philip E. Bourne, Editor

Abstract

The promoter regions of many genes contain multiple binding sites for the same transcription factor (TF). One possibility is that this multiplicity evolved through transitional forms showing redundant cis-regulation. To evaluate this hypothesis, we must disentangle the relative contributions of different evolutionary mechanisms to the evolution of binding site multiplicity. Here, we attempt to do this using a model of binding site evolution. Our model considers binding sequences and their interactions with TFs explicitly, and allows us to cast the evolution of gene networks into a neutral network framework. We then test some of the model's predictions using data from yeast. Analysis of the model suggested three candidate nonadaptive processes favoring the evolution of cis-regulatory element redundancy and multiplicity: neutral evolution in long promoters, recombination and TF promiscuity. We find that recombination rate is positively associated with binding site multiplicity in yeast. Our model also indicated that weak direct selection for multiplicity (partial redundancy) can play a major role in organisms with large populations. Our data suggest that selection for changes in gene expression level may have contributed to the evolution of multiple binding sites in yeast. We conclude that the evolution of cis-regulatory element redundancy and multiplicity is impacted by many aspects of the biology of an organism: both adaptive and nonadaptive processes, both changes in cis to binding sites and in trans to the TFs that interact with them, both the functional setting of the promoter and the population genetic context of the individuals carrying them.

Author Summary

TFs regulate gene expression by binding to specific sequences in the promoter regions of their target genes. Promoters often contain multiple copies of the same TF binding sites. How does this multiplicity evolve? One possibility is that individuals with multiple, redundant binding sites have higher fitness. However, nonadaptive processes are also likely to be important. Here, we develop a mathematical model of the evolution of TF binding sites to help us disentangle how different evolutionary mechanisms contribute to the evolution of binding site redundancy and multiplicity. We show that recombination is expected to promote the evolution of multiple binding sites. This prediction is corroborated by genome-wide data from yeast. Another important factor in the evolution of multiplicity predicted in our analysis is TF promiscuity, that is, the ability of a TF to bind to multiple sequences. In addition, our analysis indicated that direct selection can have large effects on the evolution of redundancy and multiplicity. Data from yeast identified selection for changes in expression level as a candidate mechanism for the evolution of multiple binding sites. We conclude that, although selection may play a major role in the evolution of multiplicity in regulatory regions, nonadaptive forces can also lead to high levels of multiplicity.

Introduction

Promoters frequently contain multiple functional regulatory elements [1]. For example, the regulatory region for stripe 2 of even-skipped (eve) of the fruit fly Drosophila melanogaster comprises 17 binding sites for four transcription factors (TFs), including five binding sites (B1–B5) for the activator bicoid (bcd) [2]. How does cis-regulatory element multiplicity evolve? There are three possibilities. First, perhaps “more is better” when it comes to TF binding sites. Multiple binding sites may cause changes in the level of gene expression or in its robustness against variation in TF concentrations [1], [3][5]. Second, multiplicity might be favored by selection, but independently of its functional consequences. For example, genotypes with many binding sites may be more likely to produce viable offspring after mutation or recombination with genotypes with fewer binding sites [6][9]. Third, cis-regulatory element multiplicity may arise by nonadaptive processes [9][11]. Stone and Wray [10] have shown that a population of An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e001.jpg diploid individuals could evolve two identical copies of a 6 base pair (bp) binding site in a 200-bp promoter every An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e002.jpg generations through random mutation and genetic drift alone. The intergenic regions of Saccharomyces cerevisiae are An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e003.jpg bp long on average, whereas those of multicellular eukaryotes can be orders of magnitude longer.

The common thread to all the evolutionary scenarios listed above is redundancy, the ability of structurally identical elements to contribute to the same function [12][16]. Redundancy is thought to be widespread in biological systems. In eukaryotes, a large proportion of genes are duplicates, and deletion of one copy often has little or no phenotypic effect because the other copy can compensate for the loss of function [17]. Functionality and redundancy are more difficult to establish for the case of multiple cis-regulatory elements [1]. The five bcd binding sites in eve the stripe 2 enhancer are not fully redundant because loss-of-function mutations to B1, B2 or B3 cause reduced eve stripe 2 expression and gain-of-function mutations to B4 and B5 lead to increased expression [2], [18]. However, redundancy was likely important in the evolution of these sites. When Ludwig and colleagues [3] compared the stripe 2 enhancers of different species of Drosophila, they found that some of them lacked the B3 site (Figure 1). This observation implies that the B3 site evolved recently in the lineage leading to the last common ancestor of D. melanogaster and D. simulans. Furthermore, the B3 site was probably redundant when it first appeared because the stripe 2 enhancers of three species lacking the B3 binding site were able to drive expression of a reporter gene in D. melanogaster embryos coincident with native eve stripe 2 (Figure 1). Thus, redundant transitional forms can, in principle, play an important role in the evolution of cis-regulatory element multiplicity [1], [19]. In this paper we develop a model of binding site evolution and use it to evaluate the plausibility of different scenarios for the evolution of cis-regulatory element redundancy and multiplicity. We then test predictions obtained from our model using data from yeast.

Figure 1
Evolution of the bcd binding sites in the eve stripe 2 enhancer in Drosophila.

Model

Here we introduce a model of binding site evolution. The model extends earlier phenomenological models of the evolution of cis-regulatory element redundancy [7], [9], [11], [13], [16] in that it considers binding sequences and their interactions with TFs explicitly, albeit in a simplified manner [4]. We also build upon recent attempts to apply the mutational network approach [20], [21] to the study of gene regulatory networks [22][24]. We then use our model to investigate the conditions favoring the evolution of multiple TF binding sites.

Gene regulation

A target gene has a promoter containing cis-regulatory sites for a number of TFs. An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e005.jpg denotes the jth binding site for An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e006.jpg. The TFs regulate expression of the target gene according to the following rules:

  1. An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e007.jpg binds preferentially to a canonical sequence An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e008.jpg of length n.
  2. The effect of a transcriptional activator An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e009.jpg on target gene expression through the jth binding site is given by An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e010.jpg, where An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e011.jpg is the expression level of An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e012.jpg, and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e013.jpg is a monotonically decreasing function of the number of mismatches An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e014.jpg between An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e015.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e016.jpg (i.e., the Hamming distance between the sequences), such that An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e017.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e018.jpg. See Figure S1 for examples of f functions.
  3. If An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e019.jpg is a repressor then An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e020.jpg is a monotonically increasing function of m, such that An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e021.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e022.jpg.
  4. The total effect of An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e023.jpg on gene expression is given by:
    equation image
    (1)
    An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e025.jpg will be positive for a transcriptional activator, and negative for a repressor.
  5. Target gene activity is a monotonically increasing function of An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e026.jpg.

Rules #2 and #3 are compatible with the two-state model for TF binding [4], [25][27]. Unless otherwise stated, our model deals with the evolution of the binding sites for a single transcriptional activator. For a discussion of how our model can be extended to repressors see ‘Generalizations and caveats’.

Functionality, multiplicity and redundancy

The target gene is considered functional if, given normal levels of expression of its transcriptional regulators An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e027.jpg, it is active above a threshold level, arbitrarily set to An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e028.jpg in this paper. Consider a promoter that contains An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e029.jpg binding sites for An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e030.jpg and is capable of sustaining gene function. A particular site An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e031.jpg is considered functional if binding to the site has an effect on gene expression, that is, if An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e032.jpg. Multiple binding sites (An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e033.jpg) are considered redundant if at least one of them can be deleted without affecting gene function. Full redundancy occurs when the viability of redundant and nonredundant genotypes is the same; partial redundancy occurs when the viability of redundant genotypes is higher than that of nonredundant ones [12], [15] (see also ‘Natural selection’ below). Note that, according to the above definitions, multiplicity does not imply redundancy (full or partial).

Mutation

In our model, the total effect of An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e034.jpg on the expression of a gene (An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e035.jpg, Equation 1) can change in three ways: a mutation in a binding site j that alters its An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e036.jpg (cis), a mutation in the coding sequence of An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e037.jpg that modifies the An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e038.jpg function directly (trans), and a change in the concentration of An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e039.jpg, An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e040.jpg. In the rest of the paper we consider only the first two types of evolution. We begin by considering the cis evolution of a single binding site.

One way to represent the evolution of a binding site is through its mutational network [21]. Two genotypes are connected in a mutational network if one genotype can be obtained from the other through a single mutation. For example, the sequences ACGCGC and ACGCAT are both connected to ACGCGT, but not to each other, in the mutational network of all possible DNA sequences of length An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e041.jpg base pairs (Figure 2A). If the mutation rate per base pair per generation is An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e042.jpg, then ACGCGT will mutate into ACGCGC with a probability An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e043.jpg. One difficulty with this approach is that even the relatively short sequences of TF binding sites (An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e044.jpg to 10 bp) define large mutational networks (e.g., the network of DNA sequences of length An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e045.jpg has An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e046.jpg sequences).

Figure 2
Condensed mutational network for a single binding site.

Given that binding site functionality in our model is determined by the number of mismatches m relative to the canonical sequence, we can simplify the mutational network of a binding site by collapsing all sequences with a given m. This grouping of genotypes is appropriate because every sequence contained within a given m class has exactly the same number of mutational neighbors, both within the m class (obtained by mutating already mismatched sites) and in the neighboring An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e048.jpg classes. A similar approach has been employed by others [4], [27]. We call the resulting network a condensed mutational network (Figure 2B). The m condensed genotypic class includes An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e049.jpg sequences. For example, if ACGCGT is the canonical sequence of the yeast TF Mbp1 [28], then both ACGCGC and ACGCAT belong to the An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e050.jpg condensed class. The condensed mutational network representation is extremely compact, implying a reduction from An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e051.jpg to An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e052.jpg states for a binding site of length n.

Now we introduce evolution over the condensed mutational network. Consider an infinite-sized population of asexual, haploid organisms. Each individual has a genotype at a binding site of length n, such that the population is distributed over a condensed mutational network. The state of the population is given by a vector of frequencies An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e053.jpg, where An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e054.jpg is the proportion of sequences in the population a Hamming distance m away from the canonical sequence. The population reproduces asexually, in discrete generations. Mutation causes the population to evolve according to the equation:

equation image
(2)

where An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e056.jpg is the state of the population at time t, and Q is the transition matrix such that An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e057.jpg is the probability that the offspring from an individual with i mismatches has j mismatches. Assuming that a sequence cannot acquire more than one mutation in a single generation (appropriate for realistic values of An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e058.jpg), the nonzero elements of row i of Q are given by: An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e059.jpg, An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e060.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e061.jpg. For example, the probability that the sequence ACGCGC from the An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e062.jpg condensed class will mutate into the canonical sequence ACGCGT (An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e063.jpg) is An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e064.jpg. The probability that it will mutate into a sequence from class An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e065.jpg is the probability that a mutation occurs at a site other than the already mismatched site: An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e066.jpg. And the probability that it will remain in the An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e067.jpg class is the probability that either no mutation occurs or that a mutation occurs at the mismatched site but does not result in the canonical sequence (CAn external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e068.jpgA or CAn external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e069.jpgG, but not CAn external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e070.jpgT): An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e071.jpg. The transition probabilities are the same for any other An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e072.jpg sequence, such as ACGCAT.

Natural selection

We introduce selection by assuming that target gene function is required for viability. The population can only occupy states within the viable portion of the condensed mutational network. Every generation, mutant genotypes may appear in the inviable part of the condensed mutational network, but they fail to reproduce. The evolutionary dynamics of the population can be described by restricting Equation 2 to the set of viable genotypes:

equation image
(3)

where ‘An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e074.jpg’ means entry-by-entry multiplication of the two vectors, An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e075.jpg is a vector of the viabilities of each genotypic class (1 and 0 correspond to viable and inviable, respectively) and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e076.jpg is the mean viability of the population.

Equations 2 and 3 can be generalized for genotypes with any number of binding sites K (Figures 3A, S2 and S3). This amounts to considering that each possible binding site defines an axis in a K-dimensional space, and that each point along that axis is the Hamming distance between the site and the canonical binding sequence of the corresponding TF. In the next four sections we consider the An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e077.jpg case in detail (see ‘Number of segregating binding sites’ for An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e078.jpg).

Figure 3
Condensed mutational networks for a promoter with An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e079.jpg binding sites (both with length An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e080.jpg).

Full redundancy

We begin by considering one of the simplest situations that can be represented in our model: an essential gene regulated by a single constitutively expressed activator An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e084.jpg. Gene function is required for viability. A binding site for this TF, An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e085.jpg, is functional only if it matches the canonical binding sequence exactly An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e086.jpg:

equation image
(4)

A single functional binding site is both necessary and sufficient to sustain gene function; additional functional binding sites are fully redundant.

The condensed mutational network for the case of An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e088.jpg sites is shown in Figure 3A. The axes represent the Hamming distances of each binding site relative to the canonical sequence. The viable portion of the condensed mutational network comprises the genotypic classes that have at least one functional site, and corresponds to the left and bottom edges of the mutational network in Figure 3A. There is only one redundant genotypic class: that possessing two functional binding sites (open circle in Figure 3A).

Not all genotypes within the viable portion of the condensed mutational network have the same reproductive value, defined as the proportion of viable offspring they produce. The reproductive value of condensed genotypic class i is given by: An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e089.jpg. The redundant genotype (0,0) has a reproductive value of An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e090.jpg because both of its mutational neighbors, (0,1) and (1,0), are viable. All other viable genotypes (with one functional binding site) have reproductive value An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e091.jpg (i.e., half of its mutational neighbors are inviable).

When the population reaches mutation-selection equilibrium it is not evenly distributed over all viable genotypes. Rather, the redundant genotype is approximately 2-fold overrepresented in the population, relative to other (nonredundant) genotypes (Figure 3B, squares). This finding is consistent with the prediction [20] that more highly connected genotypes in a neutral network should be overrepresented at equilibrium relative to a uniform distribution. But although the redundant genotype is overrepresented at equilibrium, redundancy cannot evolve easily in this model. That is because the (noncondensed) set of viable genotypes contains a single redundant genotype, but 8,190 nonredundant ones, 83% of which include at least 4 mismatches in the nonfunctional binding site. Thus, at equilibrium, the redundant genotype constitutes a miniscule proportion of the population (0.012%). This pattern is visible in the sum of the frequencies of all genotypes in a viable condensed genotypic class (Figure 3B, circles).

The model outlined above can be considered neutral with respect to redundancy because redundant and nonredundant genotypes have the same viability [9], [20]. The interaction between viability selection and the structure of the mutational network in this model does create indirect selection for multiplicity [20], but it is too weak to maintain a substantial proportion of redundant genotypes in the population. The results in this section are consistent with those obtained by Gerland and Hwa using a similar model [4]. In the next three sections, we build on this model by introducing different mechanisms independently, one at a time, and investigating how they affect the evolution of redundant genotypes.

Partial redundancy

Partial redundancy is thought to be more common than full redundancy [12], [15]. The presence of multiple binding sites might be advantageous if, for example, it changes the expression level of the target gene, or buffers expression against fluctuations in TF concentration [1], [3][5]. We model partial redundancy by setting the viabilities of redundant and nonredundant genotypes to An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e092.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e093.jpg, respectively. The equilibrium frequency of the redundant genotype increases with the strength of selection for redundancy (An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e094.jpg; Figure 4A). The effect of selection on redundancy undergoes a phase transition around the point where selection becomes strong relative to the rate of mutation from redundant to nonredundant genotypes (An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e095.jpg): the response to selection is small for weaker selection, but it increases sharply for stronger selection.

Figure 4
Effects of selection for multiplicity and recombination on the evolution of redundant cis-regulation.

Recombination

We incorporate recombination into our full redundancy model by taking into account the probability that each genotype has resulted from recombination between each available pair of genotypes (see Protocol S1 for details). Recombination is only allowed between binding sites, not within them. Recombination changes the evolutionary dynamics because it allows long steps across the mutational network. A modest amount of recombination between sites (An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e100.jpg) leads to the evolution of a high level of redundancy at mutation-recombination-selection equlibrium (Figure 4B). Lynch obtained similar results using a simpler model [9].

Our result can be understood by considering recombination between nonredundant genotypes containing different functional binding sites, that is, An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e101.jpg in Figure 3A. A recombination event between the sites produces two genotypes: one viable, with two functional binding sites (redundant), and another inviable, without any functional sites. In contrast, redundant genotypes always give rise to viable offspring, regardless of the kind of genotype they recombine with. This leads to strong selection against nonredundant genotypes. Stochastic simulations confirm that recombination can have a large effect on the evolution of redundancy even in finite populations (Figure S4).

TF promiscuity

Many TFs are promiscuous, that is, they can bind to several different sequences (Figure 1). Our model allows us to explore the implications of different levels and kinds of TF promiscuity for the evolution of redundancy (Figure S2B; also, see next section). We consider two different ways of increasing the promiscuity of the TF described in the basic full redundancy model (see Equation 4).

We begin by considering the “all or nothing” case where An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e102.jpg affects gene expression through a binding site An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e103.jpg according to the following relationship:

equation image
(5)

The binding site is functional only if its sequence differs from the canonical sequence of An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e105.jpg by no more than An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e106.jpg mismatches (Equation 4 is the special case for An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e107.jpg). Thus, the greater the value of An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e108.jpg, the more promiscuous the TF. Increasing An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e109.jpg expands both the size of the viable portion of the condensed mutational network and the number of redundant states. Figure 5 shows the viable portions of the condensed mutational networks for a stringent (a: An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e110.jpg) and a promiscuous TF (b: An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e111.jpg). The promiscuous TF evolves an equilibrium frequency of redundant genotypes two orders of magnitude greater than the stringent one.

Figure 5
TF promiscuity promotes the evolution of redundancy.

Generally, mismatches reduce the binding affinity and, therefore, the regulatory influence of a TF [1], [26], [29]. For example, in the D. melanogaster eve stripe 2 enhancer, the different bcd binding sites show different numbers of mismatches relative to the bcd consensus recognition sequence, inferred from in vitro binding assays [2], [30] (Figure 1). The deletion of binding sites with lower numbers of mismatches (B1 or B2, both with An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e116.jpg) result in much more severe reductions in stripe 2 expression, when compared to deletions in sites with higher numbers of mismatches (B4 or B5, both with An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e117.jpg; B3, with An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e118.jpg) [18], [31]. In addition, when the high-m sites B3–B5 were mutated into consensus sites (An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e119.jpg), they restored expression of a defective promoter lacking the B1 binding site [18]. To incorporate this type of “graded” TF promiscuity in our model, we defined An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e120.jpg as a decreasing function of An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e121.jpg (instead of a step-function as in Equation 5). Figure 5 shows two examples (c, d) that imply that graded promiscuity can promote the evolution of redundant cis-regulation more strongly than the all or nothing kind. The reason for this is that graded TF promiscuity can lead to the appearance of nonredundant genotypes containing multiple binding sites capable of sustaining gene function together but not in isolation (gray, Figure 5). These results show that nonredundant multiplicity can evolve from redundant transitional forms.

Number of segregating binding sites

Until now, our model has assumed that only alleles at An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e122.jpg binding sites segregate within a population at a given time. This assumption may not be met in reality. Long promoters provide the opportunity for more sites to arise by chance in a population [10] (Figure S2). Other factors that are expected to influence the number of segregating binding site alleles include the length of the site (n, Figure S2A), the match between the GC-content of the promoter and that of the canonical binding sequence [10], the promiscuity of the TF (m, Figure S2B; see previous section), the mutation rate and the population size.

As the number of segregating binding sites in the full redundancy model increases, the dimensionality of the model and the number of possible redundant genotypes also increase (Figures S2C and S3). The increase in the number of available redundant genotypes results in an increase in the total equilibrium frequency of these genotypes (Figure S2D). But although the number of redundant genotypes grows roughly exponentially with the number of segregating binding site alleles, the equilibrium frequency of redundant genotypes increases linearly, suggesting that the number of segregating binding site alleles has only a modest effect on the evolution of redundancy. The situation changes when the expected binding site copy number in a sequence is An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e123.jpg (e.g., points above the dashed line in Figures S2A and S2B). If that occurs, the maintenance of redundant cis-regulation does not require a selective explanation.

Generalizations and caveats

All the results derived above for an individual transcriptional activator can be generalized to two scenarios. First, to combinations of different transcriptional activators following similar rules. This would allow us to model the evolution of cis-regulatory element degeneracy (the equivalent of redundancy for elements that are structurally different [14]). Second, to transcriptional repressors (considered individually or in combination), where the function of the target gene is defined by its inactivity. Selection for decreased gene expression is expected to influence the evolution of the copy number of the binding sites of transcriptional repressors in the same way that selection for increased gene expression affects the evolution of the copy number of the binding sites of transcriptional activators. A major challenge for future work is to consider the simultaneous evolution of sites for activators and repressors in the same promoter.

Our model includes many simplifying assumptions, such as that the positions within a binding site influence TF binding uniformly and additively, and that TFs act additively through multiple binding sites. Additivity among the positions of a binding site appears to be a reasonable approximation [29], [32] and (together with uniformity) serves as the basis for the widely used two-state model [4], [25][27]. The other assumptions are not particularly realistic: synergistic effects among binding sites are commonplace [33], [34] and many TFs are not uniformly promiscuous (Table S1). The extent to which changing the assumptions of our model would modify our conclusions is not clear at present, and remains a fundamental question for future modeling.

Results

Cis-regulatory element multiplicity in yeast

To evaluate the level of regulatory multiplicity in the yeast genome, we have scanned all intergenic sequences depleted of nucleosomes [35] upstream of a single protein-coding gene (An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e124.jpg sequences, covering 8% of the genome) for 326 position weight matrix (PWM) models of 179 TFs from the literature [28], [36][38] (see Methods; Tables S1 and S2). In what follows, we analyse the 312 PWMs (corresponding to 176 TFs) predicted to have at least two binding sites in total. For simplicity, we refer to intergenic regions as promoters. On average, each promoter contained 0.08 binding sites of each PWM (standard deviation, An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e125.jpg).

We defined the amount of regulatory multiplicity (M) for a PWM as the proportion of promoters having at least one binding site that have two or more binding sites. On average, PWMs showed 7.1% multiplicity (An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e126.jpg; Table S3). The M measure of multiplicity is partly confounded with overall binding site copy number. To correct for this effect, we calculated the expected value of M for each PWM under the assumption that binding site copy number in a promoter region i of length An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e127.jpg (nucleosome depleted) is Poisson distributed with expectation An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e128.jpg, where An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e129.jpg is the observed number of binding sites in promoter j. Figure 6 shows that the observed cis-regulatory element multiplicity was approximately 40% higher than that expected under the null expectation (paired Wilcoxon test of the hypothesis that An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e130.jpg: An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e131.jpg).

Figure 6
The yeast genome shows excess cis-regulatory element multiplicity.

Evolutionary mechanisms

How did the excess multiplicity shown in Figure 6 evolve? Our model and the literature suggest three possibilities [1], [3][5], [9]: 1) recombination, 2) direct selection for increased robustness in gene expression, or 3) direct selection for increased gene expression. Each hypothesis makes a different prediction about promoters displaying binding site multiplicity: they should experience 1) higher recombination rates, or be upstream of genes showing 2) more robust expression patterns, or 3) higher expression levels (for activators; the opposite is expected for repressors). To test these hypotheses, we looked for genome-wide associations between cis-regulatory element multiplicity and a range of features of the promoters and the genes downstream of those promoters. Consider a genomic property An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e137.jpg (e.g., promoter length). For each PWM, we calculated an effect size An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e138.jpg, where An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e139.jpg (An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e140.jpg) is the mean of An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e141.jpg associated with promoters containing a single (multiple) binding site(s), and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e142.jpg is an unbiased estimate of the pooled standard deviation [39] (the effect size for binary traits, such as gene essentiality, was estimated as the arcsine transformed risk difference based on 2An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e143.jpg2 tables). We then combined the effect sizes and respective variances corresponding to different PWMs in a random-effects meta-analytic model. The results are summarized in Figure 7.

Figure 7
Cis-regulatory element multiplicity is associated with recombination and other genomic features.

Promoters with higher numbers of crossovers [40] showed significantly higher levels of binding site multiplicity (An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e147.jpg, An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e148.jpg; Figure 7), which is consistent with the recombination hypothesis. This association is explained in part by promoter length (Spearman's rank correlation coefficient: An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e149.jpg, An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e150.jpg). However, we believe that our data provide strong backing for the recombination hypothesis for three reasons. First, promoter length alone cannot explain the excess multiplicity illustrated in Figure 6 because it was considered in the calculation of An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e151.jpg. When the analysis was restricted to the subset of PWMs displaying excess multiplicity (An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e152.jpg), the effect size of crossover number was unchanged (An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e153.jpg, An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e154.jpg). Second, the effect size of the residual crossover number from a Poisson log-linear regression model with log-transformed promoter length as an explanatory variable decreased (An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e155.jpg), but remained statistically significant (An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e156.jpg, An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e157.jpg). Third, a measure of frequency of meiotic double-strand breaks (DSBs) per bp [41] was also elevated in promoters showing cis-regulatory element multiplicity (An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e158.jpg, An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e159.jpg; Figure 7).

Promoters showing binding site multiplicity tended to be upstream of genes showing low robustness in gene expression to various trans-perturbations [42] (all An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e160.jpg, An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e161.jpg, which contradicts the hypothesis that redundancy has evolved as a result of selection for robustness in gene expression. Although promoters with multiple binding sites were also more likely to contain a TATA box [43] (An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e162.jpg, An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e163.jpg), the results shown in Figure 7 did not change qualitatively when the analyses were repeated separately for genes with and without TATA boxes (not shown). Furthermore, multiplicity was not associated with protein expression noise [44] (An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e164.jpg, An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e165.jpg. Genes downstream of promoters with multiple binding sites tended to have higher expression levels (both protein and mRNA: An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e166.jpg, An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e167.jpg), which is consistent with the selection for expression hypothesis for activators, but not repressors.

Cis-regulatory element multiplicity was associated with several correlates of gene functionality (Figure 7). Promoters containing multiple sites tended to evolve more slowly [45] (divergence: An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e168.jpg, An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e169.jpg), and the genes downstream of these promoters tended to show higher levels of selective constraint [45] (An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e170.jpg: An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e171.jpg, An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e172.jpg) and to be involved in interactions with a greater number of other genes [46] (degree centrality: An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e173.jpg, An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e174.jpg). Genes with duplicates elsewhere in the genome were more likely to show binding site multiplicity [47] (An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e175.jpg, An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e176.jpg). Several gene ontology terms were significantly enriched in genes downstream of promoters containing multiple sites, including: plasma membrane, transporter activity, transcription regulator activity, DNA binding and transport (Table S5).

Discussion

Partial redundancy

Our mathematical model suggests that selection for multiple binding sites, that is, partial redundancy, can influence the evolution of cis-regulatory element multiplicity, provided that the redundant genotype has a selective advantage An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e177.jpg. Mutation rates per base pair in DNA-based organisms are of the order of An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e178.jpg [48]. Therefore, weak selection can play a major role in the evolution of cis-regulatory element multiplicity, provided that the effective population size is also large enough (An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e179.jpg) to render genetic drift negligible [4], [49].

We found a positive association in yeast between the presence of multiple binding sites for a TF and expression level of the downstream gene. This association is unlikely to have evolved neutrally or as a correlated response to the increases in multiplicity generated by recombination (see below) because gene expression patterns are under intense stabilizing selection [50], [51] and increases in gene expression are energetically costly [52]; also, the effect sizes of crossover number are not significantly correlated with those of either mRNA or protein abundance (both An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e180.jpg). Rather, the association between multiplicity and gene expression is consistent with an adaptive origin of cis-regulatory element multiplicity. Increasing the number of binding sites for transcriptional activators (inhibitors) in a promoter typically increases (decreases) gene expression [34], [53]. Since transcriptional activators are thought to be An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e181.jpg more common than repressors in yeast, and many TFs can perform either role [54], selection for different levels of expression of certain genes in certain environments could, over time, generate a positive association between cis-regulatory element multiplicity and expression level (provided that there is no strong overall bias towards selection for reduced expression). This adaptive scenario for the evolution of binding site multiplicity is consistent with the observation that functional TF binding sites have been frequently gained or lost in a lineage-specific manner among three closely related species of yeast [55], [56].

Our finding that cis-regulatory element multiplicity was associated with lower robustness in gene expression to various trans-perturbations contradicts the hypothesis that redundant genotypes benefit from being more robust [1], [5]. An earlier study reported a positive association between the number of binding sites for any TF—a possible correlate of both redundancy and degeneracy [14]—and variation in gene expression in yeast using different data from ours [57]. Our evidence is, of course, correlative: a more direct test would be to compare the robustness in expression of genes downstream of promoters containing multiple binding sites with that of the same genes with various combinations of sites mutated or deleted. Nevertheless, the observed relationships between multiplicity and robustness are also consistent with selection for changes in expression level. If the main consequence of gaining binding sites is to increase the effect of a TF on gene expression (An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e182.jpg in Equation 1), then changes in the levels of TFs, such as those caused by viable knockout mutations [42], are expected to lead to greater variance in these effects in promoters containing multiple sites, compared to promoters containing a single binding site.

Full redundancy

In addition, our model highlighted three candidate nonadaptive mechanisms for the evolution of cis-regulatory element multiplicity through fully redundant transitional forms. The first is the neutral evolution of multiple binding sites in long promoters. Such “trivial” redundancy is expected to occur in a long promoter if functional binding sites can occur over a large proportion of its length. Although the latter condition is difficult to evaluate in real organisms, intergenic regions longer than An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e183.jpg-bp are common in several mammals, including humans. Therefore, many mammalian promoters may be trivially redundant [9], [10]. This could explain the observation that approximately a third of human functional TF binding sites are not functional in rodents [58]. However, we do not expect that trivial redundancy played a dominant role in the evolution of multiplicity in organisms with relatively shorter promoters and larger populations, such as yeast.

The second mechanism is recombination. Based on our model we predict that recombination between binding sites on the order of An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e184.jpg will promote the evolution of cis-regulatory element redundancy. In yeast, a pair of sites 100-bp apart is expected to experience An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e185.jpg [48], [59]; if yeast only undergo sexual reproduction once every 1,000 asexual generations [60] we estimate An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e186.jpg, suggesting that this process has operated in yeast. We found a positive association between the presence of multiple binding sites for a TF and recombination rate in yeast. Estimates based on polymorphism data from 10 species of plants and animals [49] give An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e187.jpg bpAn external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e188.jpg, indicating that recombination is likely to be a powerful force in the evolution of cis-regulatory element multiplicity in other eukaryotes with relatively large populations. Our findings are in agreement with recent work showing that recombination selects for “mixable” genotypes [61], which leads to the evolution of higher mutational robustness [8], [62][64]. Our model predicts that redundant genotypes are robust to mutations in the binding sites, but this kind of mutational robustness does not imply robustness in the expression pattern of the downstream gene to trans-perturbations. In fact we found that cis-regulatory element multiplicity was associated with reduced robustness to perturbations in trans (see previous section).

The third nonadaptive mechanism indicated by our model is that increases in TF promiscuity promote the evolution of cis-regulatory element multiplicity. We could not test this prediction directly with our yeast data because we made an implicit assumption about the level of TF promiscuity when we scanned for binding sites. However, Bilu and Barkai [57], using a different data set from ours, reported that binding sites tended to be “fuzzier” (i.e., have lower PWM scores) when they appeared in promoter regions containing other binding sites for any TF. This observation is consistent with the prediction that graded TF promiscuity allows the existence of viable genotypes containing multiple binding sites, where each binding site is fuzzier than those found in viable genotypes containing fewer binding sites. Graded TF promiscuity is believed to be common [1], [26], [29], suggesting that multiplicity will often evolve through transitional forms showing redundant cis-regulation that then degenerate into nonredundant forms. If this evolutionary scenario is common, then lack of redundancy in extant genotypes containing multiple binding sites will be a poor indicator of whether or not its ancestral genotypes were redundant.

Conclusion

Our results suggest that redundant transitional forms can, indeed, play an important role in the evolution of cis-regulatory element multiplicity. Many aspects of the biology of an organism affect the evolution of redundancy and multiplicity: both adaptive and nonadaptive processes, both changes in cis to binding sites and in trans to the TFs that interact with them, both the functional setting of the promoter and the population genetic context of the individuals carrying them. Thus, understanding how gene networks evolve will require going beyond mere plausibility arguments into rigorous testing of specific mechanisms [4], [9]. We believe that the approach developed here provides a valuable framework to advance this research program.

Methods

Model analysis

The results reported in Figures 3 55 and were based on a deterministic version of the model (i.e., assuming infinite population size). The frequencies of different genotypic classes at mutation-selection or mutation-recombination-selection equilibrium were calculated by iterating populations for as long as necessary for genotypic class frequencies not to change by more than An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e189.jpg from one generation to the next.

TF binding site models

We used 326 PWMs summarizing the binding specificities of 179 putative yeast TFs reported in four studies [28], [36][38] (Table S1; see Protocol S1 for more details). Sequences scoring 95% or higher of the highest possible score for a given PWM were considered putative binding sites (on average, this allowed 1.23 mismatches, s.d. = 1.51). Each intergenic sequence was scanned with a PWM and its reverse complement and the number of matches were counted (simultaneous hits on exactly the same sequence and its reverse complement An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e190.jpg1 nucleotide were counted as a single hit; otherwise, binding sites overlapping over An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e191.jpg nucleotides were counted separately).

Other studies have attempted to distinguish between real and “impostor” binding sites by taking into account additional information, such as the degree of conservation of putative sites [36], [37]. We did not follow this approach because promoter sequence divergence is significantly correlated with many of the genomic features shown in Figure 7 (Table S4; see next section).

Genomic features

We calculated the following quantities for each intergenic region: 1) sequence length (including regions occupied by nucleosomes); 2) proportion of sequence occupied by nucleosomes [35]; 3) whether it contains a TATA box [43]. 4) GC content of the sequence; 5) a measure of the frequency of meiotic DSBs [41]; 6) proportion of nucleotides that differ between S. cerevisiae and S. paradoxus [45]. 7) number of crossover events [40]. We also calculated the following quantities for the gene downstream of these promoters: 1) three measures of robustness to trans-perturbations [42], derived from measurements of the variance in levels of gene expression (corrected for mean) across 167 viable knockout mutations (genetic), 30 wild isolates (genetic background), and 35 environments (environmental robustness); 2) essentiality, whether a homozygous knock-out of the gene was lethal [65], [66]; 3) whether the gene has a duplicate elsewhere in the genome [47]; 4) An external file that holds a picture, illustration, etc.
Object name is pcbi.1000848.e192.jpg, the ratio between the rates of nonsynonymous and synonymous site substitution based on the comparison between S. cerevisiae and S. paradoxus [45]; 5) degree centrality, the total number of interactions with other genes [46]; 6) protein expression noise [44]; 7) mRNA and 10) protein abundance [67], [68]. See Protocol S1 for more details.

Software

The model was analysed using Mathematica 6 (http://www.wolfram.com/mathematica/). Sequence and statistical analyses were done using R 2.9.0 (http://www.r-project.org/) and Bioconductor 2.4 (http://www.bioconductor.org/).

Supporting Information

Protocol S1

Supplementary Methods. Sections: modeling recombination; yeast data; software.

(0.09 MB PDF)

Figure S1

Examples of f functions consistent with the viable portions of the condensed mutational networks in Figure 5. The dashed line indicates the threshold for driving gene expression.

(0.06 MB PDF)

Figure S2

Redundancy is more likely to evolve if there are more segregating binding site alleles. (A) Expected number of exact matches to canonical binding sequences of different lengths (n) in promoters of different lengths (L). (B) Expected number of matches to an 8-bp canonical binding sequence allowing for different numbers of mismatches (m) in promoters of different L. The value of m models different levels of TF promiscuity. In (A) and (B) values are means and 95% confidence intervals of 10 independent sets of 104 random sequences with the same average GC content as yeast intergenic regions (except for n  = 8, m  = 0 and L≤200, where 60 sets of sequences were used). Dashed lines mark an expected number of 2 binding sites. (C) Number of redundant genotypes and (D) total equilibrium frequency of redundant genotypes for different numbers of segregating binding sites (K). For K  = 2, the model is that shown in Figure 2. See Figure S3 for K  = 3.

(0.09 MB PDF)

Figure S3

Condensed mutational networks for a promoter with K  = 3 binding sites (all with length n  = 6). (A) Diagram of gene with three binding sites. (B) Condensed mutational network. Axes represent the numbers of mismatches of each binding site relative to the canonical sequence. Each node represents a genotypic class. As in Figure 2, the promoter regulates an essential gene such that at least one canonical binding site is required for activity. The nodes shown in black define the viable portion of the condensed mutational network. The nodes in gray represent inviable genotypes. (C) Shows only the viable portion of the condensed mutational network. The genotypes highlighted in gray are redundant.

(0.17 MB PDF)

Figure S4

Stochastic simulations of the effect of recombination. Populations of different sizes (N) are initialized at mutation-selection equilibrium. (A) r/μ  = 0.1, (B) r/μ  = 1, and (C) r/μ  = 10. In all cases, we used μ  = 0.1, an unrealistically high value. Values are medians of the frequencies of redundant genotypes for 500 replicate populations. In populations of both sizes redundancy evolves quickly, but is then lost by drift. Dotted lines show the deterministic expectation (see Figure 3B).

(0.08 MB PDF)

Table S1

Binding site PWM models used in our study. Values are the length (8.7 ± 2.6 bp, mean ± standard deviation), GC content (0.53 ± 0.19) and mean information content (I) per position (1.33 ± 0.29) of each PWM (after processing as described in Protocol S1). The value of I can vary between 0 and 2, and is a measure of the energy contribution of a position to TF binding. Each PWM is summarized by its canonical sequence: “.” indicates a position with I  = 0; “[ / ]” indicates bases with the same weight at a given position. Letters in parentheses after TF names indicate the study from which we took the PWM data: B, Badis et al. (2008); H, Harbison et al. (2004); M, MacIsaac et al. (2006); Zhu et al. (2009).

(0.06 MB PDF)

Table S2

Binding site PWM models not considered in our study. These PWMs were excluded because they were almost identical to the PWMs listed in the ‘Equivalent’ column, shown in Table S1. Letters in parentheses after TF names indicate the study from which we took the PWM data: B, Badis et al. (2008); H, Harbison et al. (2004); M, MacIsaac et al. (2006); Zhu et al. (2009).

(0.04 MB PDF)

Table S3

Data used to construct Figure 6. The data are sorted by decreasing M Obs/M Exp. The second and third columns show the Total number of binding sites revealed in a scan across the number of promoter regions in the ‘Prom’ column (the numbers vary because the minimum promoter length considered in each scan is twice the length of the PWM). Letters in parentheses after TF names indicate the study from which we took the PWM data: B, Badis et al. (2008); H, Harbison et al. (2004); M, MacIsaac et al. (2006); Zhu et al. (2009).

(0.05 MB PDF)

Table S4

Matrix of correlations between genomic features. Values are Spearman's rank correlation coefficients (ρ). Data used for the cluster analysis in Figure 7.

(0.04 MB PDF)

Table S5

Association between multiplicity and GO slim terms from each domain. Significance levels after correction for multiple comparisons using the Holm method: *, P<0.01; **, P<0.001; ***, P<0.0001.

(0.06 MB PDF)

Acknowledgments

We thank G. Balázsi, T. Cooper, S. Proulx, C. Ray and M. Siegal for helpful discussions. S. Proulx shared the yeast gene expression robustness data.

Footnotes

The authors have declared that no competing interests exist.

This work was supported in part by grants from the National Science Foundation (EF-0742803) and the James S. McDonnell Foundation to RBRA. TP was supported in part by a fellowship from the Fundação para a Ciência e a Tecnologia (Portugal). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

References

1. Wray GA, Hahn MW, Amores A, Balhoff JP, Pizer M, et al. The evolution of transcriptional regulation in eukaryotes. Mol Biol Evol. 2003;20:1377–1419. [PubMed]
2. Stanojevic D, Small S, Levine M. Regulation of a segmentation stripe by overlapping activators and repressors in the Drosophila embryo. Science. 1991;254:1385–1387. [PubMed]
3. Ludwig MZ, Patel NH, Kreitman M. Functional analysis of eve stripe 2 enhancer evolution in Drosophila: rules governing conservation and change. Development. 1998;125:949–958. [PubMed]
4. Gerland U, Hwa T. On the selection and evolution of regulatory DNA motifs. J Mol Evol. 2002;55:386–400. [PubMed]
5. O'Meara MM, Bigelow H, Flibotte S, Etchberger JF, Moerman DG, et al. Cis-regulatory mutations in the Caenorhabditis elegans homeobox gene locus cog-1 affect neuronal development. Genetics. 2009;181:1679–1686. [PubMed]
6. Wagner A. Redundant gene functions and natural selection. J Evol Biol. 1999;12:1–16.
7. Proulx SR, Phillips PC. The opportunity for canalization and the evolution of genetic networks. Am Nat. 2005;165:147–162. [PubMed]
8. Azevedo RBR, Lohaus R, Srinivasan S, Dang KK, Burch CL. Sexual reproduction selects for robustness and negative epistasis in artificial gene networks. Nature. 2006;440:87–90. [PubMed]
9. Lynch M. The evolution of genetic networks by non-adaptive processes. Nat Rev Genet. 2007;8:803–813. [PubMed]
10. Stone JR, Wray GA. Rapid evolution of cis-regulatory sequences via local point mutations. Mol Biol Evol. 2001;18:1764–1770. [PubMed]
11. Wagner GP, Otto W, Lynch V, Stadler PF. A stochastic model for the evolution of transcription factor binding site abundance. J Theor Biol. 2007;247:544–53. [PubMed]
12. Thomas JH. Thinking about genetic redundancy. Tr Genet. 1993;9:395–399. [PubMed]
13. Nowak MA, Boerlijst MC, Cooke J, Maynard Smith J. Evolution of genetic redundancy. Nature. 1997;388:167–171. [PubMed]
14. Edelman GM, Gally JA. Degeneracy and complexity in biological systems. Proc Natl Acad Sci U S A. 2001;98:13763–8. [PubMed]
15. Wagner A. Distributed robustness versus redundancy as causes of mutational robustness. BioEssays. 2005;27:176–188. [PubMed]
16. Frank SA. Evolutionary dynamics of redundant regulatory control. J Theor Biol. 2008;255:64–68. [PMC free article] [PubMed]
17. Gu ZL, Steinmetz LM, Gu X, Scharfe C, Davis RW, et al. Role of duplicate genes in genetic robustness against null mutations. Nature. 2003;421:63–66. [PubMed]
18. Arnosti DN, Barolo S, Levine M, Small S. The eve stripe 2 enhancer employs multiple modes of transcriptional synergy. Development. 1996;122:205–14. [PubMed]
19. Dermitzakis ET, Bergman CM, Clark AG. Tracing the evolutionary history of Drosophila regulatory regions with models that identify transcription factor binding sites. Mol Biol Evol. 2003;20:703–714. [PubMed]
20. van Nimwegen E, Crutchfield JP, Huynen M. Neutral evolution of mutational robustness. Proc Natl Acad Sci U S A. 1999;96:9716–9720. [PubMed]
21. Cowperthwaite MC, Meyers LA. How mutational networks shape evolution: Lessons from RNA models. Annu Rev Ecol Evol Syst. 2007;38:203–230.
22. Wagner A. Circuit topology and the evolution of robustness in two-gene circadian oscillators. Proc Natl Acad Sci U S A. 2005;102:11775–11780. [PubMed]
23. Ciliberti S, Martin O, Wagner A. Robustness can evolve gradually in complex regulatory gene networks with varying topology. PLoS Comput Biol. 2007;3:e15. [PubMed]
24. Azevedo RBR, Lohaus R, Paixão T. Networking networks. Evol Dev. 2008;10:514–515. [PubMed]
25. Berg O, von Hippel P. Selection of DNA binding sites by regulatory proteins. statistical-mechanical theory and application to operators and promoters. J Mol Biol. 1987;193:723–750. [PubMed]
26. Gerland U, Moroz J, Hwa T. Physical constraints and functional characteristics of transcription factor-DNA interaction. Proc Natl Acad Sci U S A. 2002;99:12015–12020. [PubMed]
27. Berg J, Willmann S, Lässig M. Adaptive evolution of transcription factor binding sites. BMC Evol Biol. 2004;4:42. [PMC free article] [PubMed]
28. Badis G, Chan ET, van Bakel H, Pena-Castillo L, Tillo D, et al. A library of yeast transcription factor motifs reveals a widespread function for Rsc3 in targeting nucleosome exclusion at promoters. Mol Cell. 2008;32:878–887. [PMC free article] [PubMed]
29. Stormo G, Fields D. Specificity, free energy and information content in protein–DNA interactions. Tr Biochem Sci. 1998;23:109–113. [PubMed]
30. Driever W, Thoma G, Nüsslein-Volhard C. Determination of spatial domains of zygotic gene expression in the Drosophila embryo by the affinity of binding sites for the bicoid morphogen. Nature. 1989;340:363–367. [PubMed]
31. Small S, Blair A, Levine M. Regulation of even-skipped stripe 2 in the Drosophila embryo. EMBO J. 1992;11:4047–57. [PubMed]
32. Benos PV, Bulyk ML, Stormo GD. Additivity in protein-DNA interactions: how good an approximation is it? Nucl Acids Res. 2002;30:4442–4451. [PMC free article] [PubMed]
33. Giniger E, Ptashne M. Cooperative DNA binding of the yeast transcriptional activator GAL4. Proc Natl Acad Sci U S A. 1988;85:382–386. [PubMed]
34. Anderson GM, Freytag SO. Synergistic activation of a human promoter in vivo by transcription factor Sp1. Mol Cell Biol. 1991;11:1935–1943. [PMC free article] [PubMed]
35. Lee W, Tillo D, Bray N, Morse RH, Davis RW, et al. A high-resolution atlas of nucleosome occupancy in yeast. Nat Genet. 2007;39:1235–1244. [PubMed]
36. Harbison CT, Gordon DB, Lee TI, Rinaldi NJ, Macisaac KD, et al. Transcriptional regulatory code of a eukaryotic genome. Nature. 2004;431:99–104. [PMC free article] [PubMed]
37. MacIsaac KD, Wang T, Gordon DB, Gifford DK, Stormo GD, et al. An improved map of conserved regulatory sites for Saccharomyces cerevisiae. BMC Bioinformatics. 2006;7:113. [PMC free article] [PubMed]
38. Zhu C, Byers KJRP, McCord RP, Shi Z, Berger MF, et al. High-resolution DNA-binding specificity analysis of yeast transcription factors. Genome Res. 2009;19:556–566. [PubMed]
39. Hedges L, Olkin I. Statistical Methods for Meta-Analysis. Orlando, FL: Academic Press; 1985.
40. Mancera E, Bourgon R, Brozzi A, Huber W, Steinmetz LM. High-resolution mapping of meiotic crossovers and non-crossovers in yeast. Nature. 2008;454:479–485. [PMC free article] [PubMed]
41. Buhler C, Borde V, Lichten M. Mapping meiotic single-strand DNA reveals a new landscape of DNA double-strand breaks in Saccharomyces cerevisiae. PLoS Biol. 2007;5:e324. [PubMed]
42. Proulx SR, Nuzhdin S, Promislow DEL. Direct selection on genetic robustness revealed in the yeast transcriptome. PLoS ONE. 2007;2:e911. [PMC free article] [PubMed]
43. Basehoar AD, Zanton SJ, Pugh BF. Identification and distinct regulation of yeast TATA box-containing genes. Cell. 2004;116:699–709. [PubMed]
44. Newman JRS, Ghaemmaghami S, Ihmels J, Breslow DK, Noble M, et al. Single-cell proteomic analysis of S. cerevisiae reveals the architecture of biological noise. Nature. 2006;441:840–846. [PubMed]
45. Kellis M, Patterson N, Endrizzi M, Birren B, Lander ES. Sequencing and comparison of yeast species to identify genes and regulatory elements. Nature. 2003;423:241–254. [PubMed]
46. Lee I, Li Z, Marcotte E. An improved, bias-reduced probabilistic functional gene network of baker's yeast, Saccharomyces cerevisiae. PLoS ONE. 2007;2:e988. [PMC free article] [PubMed]
47. Kellis M, Birren BW, Lander ES. Proof and evolutionary analysis of ancient genome duplication in the yeast Saccharomyces cerevisiae. Nature. 2004;428:617–24. [PubMed]
48. Drake JW, Charlesworth B, Charlesworth D, Crow JF. Rates of spontaneous mutation. Genetics. 1998;148:1667–86. [PubMed]
49. Lynch M. The Origins of Genome Architecture. Sunderland, Mass.: Sinauer Associates; 2007.
50. Denver DR, Morris K, Streelman JT, Kim SK, Lynch M, et al. The transcriptional consequences of mutation and natural selection in Caenorhabditis elegans. Nat Genet. 2005;37:544–8. [PubMed]
51. 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:220–223. [PubMed]
52. Wagner A. Energy constraints on the evolution of gene expression. Mol Biol Evol. 2005;22:1365–1374. [PubMed]
53. Kusch T, Storck T, Walldorf U, Reuter R. Brachyury proteins regulate target genes through modular binding sites in a cooperative fashion. Genes Dev. 2002;16:518–529. [PubMed]
54. Guelzim N, Bottani S, Bourgine P, Képès F. Topological and causal structure of the yeast transcriptional regulatory network. Nat Genet. 2002;31:60–63. [PubMed]
55. Borneman AR, Gianoulis TA, Zhang ZD, Yu H, Rozowsky J, et al. Divergence of transcription factor binding sites across related yeast species. Science. 2007;317:815–819. [PubMed]
56. Doniger SW, Fay JC. Frequent gain and loss of functional transcription factor binding sites. PLoS Comput Biol. 2007;3:e99. [PMC free article] [PubMed]
57. Bilu Y, Barkai N. The design of transcription-factor binding sites is affected by combinatorial regulation. Genome Biology. 2005;6:R103. [PMC free article] [PubMed]
58. Dermitzakis ET, Clark AG. Evolution of transcription factor binding sites in mammalian gene regulatory regions: Conservation and turnover. Mol Biol Evol. 2002;19:1114–1121. [PubMed]
59. Kaback DB, Steensma HY, de Jonge P. Enhanced meiotic recombination on the smallest chromosome of saccharomyces cerevisiae. Proc Natl Acad Sci U S A. 1989;86:3694–3698. [PubMed]
60. Tsai IJ, Bensasson D, Burt A, Koufopanou V. Population genomics of the wild yeast Saccharomyces paradoxus: Quantifying the life cycle. Proc Natl Acad Sci U S A. 2008;105:4957–4962. [PubMed]
61. Livnat A, Papadimitriou C, Dushoff J, Feldman MW. A mixability theory for the role of sex in evolution. Proc Natl Acad Sci U S A. 2008;105:19803–19808. [PubMed]
62. Gardner A, Kalinka AT. Recombination and the evolution of mutational robustness. J Theor Biol. 2006;241:707–715. [PubMed]
63. Kim K, Fernandes V. Effects of ploidy and recombination on evolution of robustness in a model of the segment polarity network. PLoS Comput Biol. 2009;5:e1000296. [PMC free article] [PubMed]
64. Misevic D, Ofria C, Lenski RE. Sexual reproduction reshapes the genetic architecture of digital organisms. Proc R Soc B. 2006;273:457–464. [PMC free article] [PubMed]
65. Giaever G, Chu A, Ni L, Connelly C, Riles L, et al. Functional profiling of the Saccharomyces cerevisiae genome. Nature. 2002;418:387–391. [PubMed]
66. Winzeler E, Shoemaker D, Astromoff A, Liang H, Anderson K, et al. Functional characterization of the S. cerevisiae genome by gene deletion and parallel analysis. Science. 1999;285:901. [PubMed]
67. Holstege FCP, Jennings EG, Wyrick JJ, Lee TI, Hengartner CJ, et al. Dissecting the regulatory circuitry of a eukaryotic genome. Cell. 1998;95:717–728. [PubMed]
68. Ghaemmaghami S, Huh W, Bower K, Howson RW, Belle A, et al. Global analysis of protein expression in yeast. Nature. 2003;425:737–741. [PubMed]
69. Drosophila 12 Genomes Consortium. Evolution of genes and genomes on the Drosophila phylogeny. Nature. 2007;450:203–218. [PubMed]

Articles from PLoS Computational Biology are provided here courtesy of Public Library of Science