Genotype-phenotype landscapes exhibit composition bias in accessible mutational paths

We used protein-binding microarray data [46, 47] to construct empirical genotype-phenotype landscapes of transcription factor binding affinities for 746 transcription factors from 129 eukaryotic species, representing 48 distinct DNA binding domain structural classes (Methods; S1 Table). For each transcription factor, these data include an enrichment score (E-score)—a proxy for relative binding affinity—for all possible 32, 896 DNA sequences of length eight. We constructed one landscape per transcription factor, using only those DNA sequences that specifically bound the transcription factor, as indicated by an E-score exceeding 0.35 [22, 35, 36]. We represented each landscape as a genotype network, in which nodes are transcription factor binding sites and edges connect nodes if their corresponding sequences differ by a single point mutation (Fig 1b) [35]. For some transcription factors (∼ 37%), the genotype network fragmented into several disconnected components. When this occurred, we only considered the largest component, which always comprised more than 100 bound sequences. We refer to this as the dominant genotype network. We discarded non-dominant components because they usually comprised few sequences and rarely met our size requirement of 100 sequences. (S1 Fig). Each dominant genotype network formed the substrate of a genotype-phenotype landscape, whose surface was defined by the relative affinities (E-scores) of the network’s constituent binding sites [22]. We accounted for noise in the protein binding microarray data using a noise threshold δ, which allowed us to determine whether two DNA sequences differed in their binding affinity [22] (Methods). We refer to a mutation as accessible if it increases binding affinity more than δ, and we refer to a series of accessible mutations as an accessible mutational path (Fig 1c).

We developed a measure of composition bias on the unit interval, which we applied to entire landscapes and to accessible mutational paths (Methods). It is based on the null expectation that one transition occurs per every two transversions. When this measure equals 0.5, there is no composition bias. Values below 0.5 imply a composition bias toward transversions, whereas values above 0.5 imply a composition bias toward transitions. Fig 2a shows the distribution of composition bias across all 746 landscapes. When considering all mutations in a landscape, we observed variation in composition bias ranging from 0.19 to 0.75, but as a whole, this distribution did not differ from our null expectation of one transition per every two transversions (black distribution in Fig 2a, one-sample t test, t₇₄₅ = −0.4310, p = 0.66). In contrast, when we considered accessible mutational paths to the global peak, we found significant composition bias toward transversions (white distribution in Fig 2a, one-sample t test, t₇₄₅ = −14.3941, p < 10⁻⁴⁰). This bias was sometimes extreme. For example, for 68 transcription factors, fewer than 1 in 7 mutations on an accessible mutational path were transitions, more than a three-fold decrease relative to our null expectation. Such composition bias varied among transcription factors from different DNA binding domain structural classes (Fig 2b). For example, AT hook transcription factors had a strong bias toward transversions, because their binding sites are enriched for adenine and thymine bases, and thus have a very low GC content. At the opposite extreme were transcription factors with a basic helix-loop-helix domain (bHLH), which had a bias toward transitions, because the sequences they bind are more neutral in terms of GC content (average GC content = 0.61). Landscapes with either many high or many low GC content sequences were more likely to exhibit composition bias toward transversion mutations (S2 Fig). Overall, we observed more extreme bias toward transversions than toward transitions in accessible mutational paths (Fig 2a), a bias that became even more pronounced as we increased the noise parameter δ (S3 Fig). This is because transversions tend to cause larger changes in binding affinity [48], and thus have larger regulatory effects [49]—a phenomenon we observed mainly near the global peak (S4 Fig). Because these distributions deviate so strongly from our null expectation, we focus on the composition bias of accessible mutational paths in all subsequent analyses, and we refer to this simply as composition bias for brevity.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Accessible mutational paths exhibit composition bias toward transversions.

(a) The black bars show the distribution of composition bias across entire landscapes. The white bars show the distribution of composition bias across accessible mutational paths to the global peak, starting from the 10% of binding sites with the lowest affinities in each landscape. Gray indicates the overlap in the distributions. Data pertain to all 746 landscapes. (b) Composition bias of accessible mutational paths, with landscapes grouped by DNA binding domain structural class. Numbers in parentheses indicate the number of transcription factors per class in our dataset.

https://doi.org/10.1371/journal.pcbi.1008296.g002

To test the sensitivity of these observations to our use of E-score as a proxy for relative binding affinity, we considered an alternative proxy—the median signal intensity Z-score—which is available for 713 of the 746 transcription factors in our dataset. We found a strong correlation between the composition bias calculated from the two scores (Pearson’s correlation coefficient r = 0.7212, p < 10⁻¹⁰) (S5 Fig), and although composition bias often varied quantitatively among the two scores, it did not often vary qualitatively. That is, it was uncommon for landscapes to switch from exhibiting a composition bias toward transversions to a composition bias toward transitions, or vice versa. There were only 83 transcription factors (12%) for which such switching occurred, and of these roughly half (∼ 48%) exhibited little to no composition bias when using E-scores (in the range (0.45, 0.55)). We therefore used E-scores for the rest of our analyses. In addition, we tested the sensitivity of our results to an alternative, less conservative definition of an accessible mutational path, which included any mutational steps that did not decrease binding affinity beyond some threshold (δ—our noise parameter). S6 Fig shows that this alternative definition of accessibility results in composition biases that are highly similar to those observed under our initial more conservative definition, which we therefore use for the rest of our analyses, unless otherwise noted.

Mutation bias interacts with composition bias to influence landscape navigability

We first explored how mutation bias and composition bias interact to influence adaptive evolution when the mutation supply is low and selection is strong. Under these population genetic conditions, only one mutation is present in the population at any time [50], which makes the process amenable to modeling as a random walk in genotype space (Methods). In this framework, each time step corresponds to the number of generations needed for a mutation to go to fixation, and the fixation probability of a mutant is proportional to its binding affinity and to the likelihood of that particular type of mutation occurring. The latter was determined by a mutation bias parameter α, which is defined similarly to our measure of composition bias (Methods). Specifically, when α = 0.5 there is no bias in mutation supply, values of α below 0.5 mean that the mutation supply is biased toward transversions, whereas values above 0.5 mean that the mutation supply is biased toward transitions.

As a measure of landscape navigability, we calculated the probability of reaching the global peak (Methods). Fig 3 shows how mutation bias and composition bias interact to influence landscape navigability. We found that mutation bias could either enhance or diminish landscape navigability, relative to an unbiased mutation supply, depending upon whether the bias in mutation supply aligned with the composition bias of the accessible mutational paths in the landscape (i.e., the two forms of bias were toward the same kind of mutations—transitions or transversions). This effect is clearly seen when comparing the value of mutation bias that maximizes the probability of reaching the global peak across the five panels of Fig 3 (dashed vertical lines), which group landscapes according to their composition bias. S7 Fig shows the correlation between the mutation bias that maximizes P_peak and composition bias, measured using two different definitions of accessible mutational paths. Landscapes with extreme composition bias exhibited increased sensitivity to mutation bias, relative to landscapes without composition bias. For example, for landscapes with a strong composition bias toward transversions (Fig 3a), the probability of reaching the global peak increased from 0.12 with a strong mutation bias toward transitions to 0.31 with a strong mutation bias toward transversions—a 2.6-fold increase. In contrast, for landscapes without composition bias (Fig 3c), only a 1.7-fold increase in the probability of reaching the global peak was obtained by varying the bias in mutation supply.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Mutation bias interacts with composition bias to influence landscape navigability.

The probability P_peak of reaching the global peak is shown for 19 different values of the mutation bias parameter α. The solid vertical lines indicate no bias in mutation supply (α = 0.5) and the dashed vertical lines indicate the value of α that maximizes P_peak. Landscapes are grouped based on their composition bias and the distribution of composition bias per panel is shown on top of each panel. The number of landscapes per panel is indicated is the bottom left corner.

https://doi.org/10.1371/journal.pcbi.1008296.g003

Overall, landscapes with strong composition bias were more navigable than those with intermediate or no bias, in that they had higher probabilities of reaching the global peak. We reasoned this is because landscapes with strong composition bias tend to contain fewer binding sites than those with intermediate or no bias (S8 Fig), and the probability of evolving to a landscape’s global peak decreases with the number of binding sites in the landscape (as shown in S9 Fig for an unbiased mutation supply). In addition, our finding that navigability will be enhanced when mutation bias and composition bias are aligned and diminished otherwise is insensitive to whether selection acts to increase or decrease binding affinity (S10 Fig). This is important because low-affinity binding sites often drive gene expression, for example in developmental enhancers [44, 45].

Mutation bias interacts with composition bias to influence the predictability of evolution

When the mutation supply is low, a bias in mutation supply can influence an evolving population’s adaptive trajectory through a genotype-phenotype landscape, making some mutational paths more likely than others [25, 26]. Mutation bias may therefore influence the predictability of evolution. We next explored how mutation bias interacts with composition bias to influence the predictability of evolution, quantifying predictability using a measure of path entropy [51] (Methods). This measure takes on low values when an evolving population tends to take few mutational paths to the global peak, each with high probability. It takes on high values when an evolving population tends to take many mutational paths to the global peak, each with low probability. It is therefore inversely related to the predictability of evolution. S11 Fig shows the relationship between path entropy and the mutation bias parameter α for 50 randomly chosen landscapes.

As expected, we found that mutation bias is considerably more likely to increase the predictability of evolution than to decrease it, in line with the notion of mutation bias as an orienting factor in evolution [25]. Across all landscapes, there was at least one mutation bias value that decreased path entropy relative to when there was no mutation bias (Fig 4a), and on average 73% of mutation bias values decreased path entropy relative to when there was no mutation bias (S12 Fig). Moreover, a mutation bias toward transitions minimized path entropy more often than one toward transversions (Fig 4a). Mutation bias therefore readily increases the predictability of evolution. Fig 4b shows how the misalignment between mutation bias and composition bias usually leads to the minimization of path entropy, thus increasing predictability. For landscapes with little to no composition bias, a strong mutation bias toward transitions was more likely to minimize path entropy than a strong mutation bias toward transversions. This is because in landscapes with no composition bias there are twice as many transversions than transitions, so a mutation bias toward transitions represents a greater evolutionary constraint in such cases. To determine the extent to which entropy changes in response to mutation bias, we calculated the ratio of the entropy observed in the absence of mutation bias to the minimum entropy caused by mutation bias (entropy_{no bias}/entropy_min), for each landscape. Path entropy decreased by an average of approximately 2-fold and 3-fold in the most extreme mutation biases toward transversions and transitions, respectively (Fig 4c).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Mutation bias interacts with composition bias to impact the predictability of evolution.

(a) Distribution of the mutation bias parameter α that minimizes path entropy for each landscape. (b) Mutation bias parameter that minimizes path entropy, shown in relation to composition bias. (c) Relative entropy change (entropy_{no bias}/entropy_min), shown in relation to the mutation bias parameter α that minimizes path entropy. (d) Distribution of the mutation bias parameter α that maximizes path entropy for each landscape. (e) Mutation bias parameter that maximizes path entropy, shown in relation to composition bias. (f) Relative entropy change (entropy_max/entropy_{no bias}), shown in relation to the mutation bias parameter α that maximizes path entropy. Data in all panels pertain to all 746 landscapes.

https://doi.org/10.1371/journal.pcbi.1008296.g004

Perhaps counter-intuitively, we found that mutation bias can also decrease the predictability of evolution. For 693 landscapes (93%), there was at least one mutation bias value that increased path entropy relative to when there was no mutation bias (Fig 4d), and on average across all 746 landscapes, 27% of mutation bias values increased path entropy relative to when there was no mutation bias (S12 Fig). To explore this further, we determined the mutation bias parameter α that maximized path entropy for each landscape. Fig 4d shows the distribution of this parameter, which varied across the full range of α values considered. We hypothesized that the value of α that maximizes path entropy will be positively correlated with composition bias, such that path entropy will be maximized when these biases are aligned. Fig 4e reveals this positive correlation, which is statistically significant, but weak in explanatory power (Spearman’s rank correlation coefficient ρ = 0.13, p < 10⁻⁴). The reason is that our measure of composition bias does not capture how mutations are distributed across accessible mutational paths, which strongly influences the value of α that maximizes path entropy (S13 Fig). In addition, Fig 4f shows the relative change in entropy (entropy_max/entropy_{no bias}) in relation to the mutation bias parameter α that maximized path entropy. It reveals that path entropy increased by up to 7-fold, and that the largest increases were associated with the most extreme biases in mutation supply, either toward transitions or toward transversions. In these cases, evolution became less predictable because an evolving population could traverse a greater diversity of accessible mutational paths. In sum, these analyses reveal that mutation bias and composition bias interact to influence the predictability of evolution, in most cases increasing predictability, but in many others decreasing it. Thus, whether mutation bias acts as an orienting or dispersive factor in evolution depends upon the prevalence and type of composition bias in the landscape.

Mutation bias influences the distribution of polymorphic populations in genotype-phenotype landscapes

We next explored how mutation bias influences adaptive evolution when the mutation supply is high and selection is strong. Under these population genetic conditions, multiple mutations coexist in the population and compete for fixation. Because this process is challenging to model analytically, we used computer simulations of a Wright-Fisher model of evolutionary dynamics (Methods). As in our previous analyses, we used the probability of reaching the global peak as a measure of landscape navigability. We found that mutation bias has no effect on navigability in this “pick-the-winner” regime (S14 Fig), as expected on theoretical grounds [25]. However, we reasoned that mutation bias may still influence other properties of an evolving population, specifically those that depend upon the population’s distribution in the landscape. To explore this possibility, we used a measure called an overlap coefficient, which quantifies the similarity of two populations as the proportion of individuals that are common to both (Methods). This coefficient takes on its minimum value of 0 when there are no individuals in common between two populations; it takes on its maximum value of 1 when both populations are identical, having the same individuals in the same proportions. We applied this measure to pairs of populations after they had evolved for 1000 generations, reaching steady state (S15 Fig). As a baseline for comparison, we first calculated the overlap coefficient for pairs of replicate populations. That is, pairs of populations with identical initial conditions, but with different random number generator seeds (Fig 5a). This allowed us to assess how different we expect two evolved populations to be at steady state, due solely to the stochastic nature of the evolutionary simulations. For replicate populations, the overlap coefficient ranged from 0.912 to 1, with a median of 0.976 and a 75th percentile of 0.817 (Fig 5b). This indicates that while the stochastic nature of the evolutionary simulation can cause large changes in a polymorphic population’s distribution in a landscape, it usually does not. Replicate populations tend to converge on highly similar distributions. In contrast, when we calculated the overlap coefficient for pairs of evolved populations with identical initial conditions (including identical random number generator seeds), but different values of the mutation bias parameter, we observed far less overlap (two-sample Kolmogorov-Smirnov test, D = 0.2695, p < 10⁻³¹). Specifically, the overlap coefficient ranged from 0.692 to 1, with a median of 0.908 and a 75th percentile of 0.274 (Fig 5b). This indicates that mutation bias often has a strong influence on an evolved population’s distribution in a genotype-phenotype landscape. The strength of this influence depends upon how different the mutation bias values are in the populations being compared (Fig 5c), with larger differences corresponding to larger changes in population distribution, a trend that also holds for infinitely large populations (S16 Fig).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Evolving polymorphic populations are more sensitive to changes in mutation bias than to the stochastic nature of the evolutionary simulations.

(a) Schematic figure of our experimental design. For each landscape and combination of population size and mutation rate (Nμ), we considered 10 replicates for each of 10 different initial conditions and 19 values of the mutation bias parameter α. Importantly, we used the replicate number to seed the random number generator of each evolutionary simulation, facilitating the comparison of variation across replicates versus across the mutation bias parameter α. For example, the matrix elements indicated in gray contain the information necessary to compare the effects of the mutation bias parameter with the stochasticity of the evolutionary simulations, for one initial condition. (b) Overlap coefficient for pairs of evolved populations that differ in random number generator seed (“Replicates”) or in mutation bias parameter (“Mutation bias”). Notches indicate medians, whiskers indicate the 25th and 75th percentiles, and cross symbols indicate outliers. (c) Overlap coefficient for pairs of evolved populations, shown in relation to the difference in their mutation bias parameters.

https://doi.org/10.1371/journal.pcbi.1008296.g005

Mutation bias and composition bias interact to influence the evolution of genetic diversity and mutational robustness in polymorphic populations

We next asked how mutation bias interacts with composition bias to influence the evolution of genetic diversity. We reasoned that when mutation bias is aligned with composition bias, evolving populations will be less constrained in their exploration of the landscape and will therefore accumulate greater genetic diversity. To explore this possibility, we measured the genetic diversity of populations at steady state using Shannon’s diversity index (Methods). This measure takes on its maximum value of 1 when the population comprises all possible individuals in equal proportions. For DNA sequences, this means that all four bases are equally likely to appear at all positions in the sequence. The measure takes on its minimum value of 0 when the population comprises N copies of just a single individual. Fig 6a–6e shows that mutation bias can either increase or decrease genetic diversity, relative to an unbiased mutation supply, and depending on whether mutation bias aligns with composition bias. This effect was most pronounced in landscapes with strong composition bias, either toward transversions (Fig 6a) or transitions (Fig 6e), and when the mutation supply was high. For example, in landscapes with a strong composition bias toward transversions (Fig 6a), a bias in mutation supply could change genetic diversity 1.9-fold when Nμ = 50, but had almost no effect when Nμ = 5. In these cases, genetic diversity could reach levels higher than those observed on landscapes with little to no bias (Fig 6c). Calculating genetic diversity using nucleotide diversity π results in similar patterns (S17 Fig). Specifically, diversity increases when mutation bias aligns with composition bias and decreases otherwise. To further illustrate how mutation bias and composition bias interact to influence population diversity, we show in S18 Fig the allele frequency spectra of populations evolved on two different landscapes, which we chose as illustrative examples because of their strong composition bias toward transversions (S18a Fig) and toward transitions (S18b Fig). As expected, the allele frequency spectra reflect the sequence bias of the mutation process, with more transition polymorphisms present when mutation is biased toward transitions and more transversion polymorphisms present when mutation is biased toward transversions. These examples also illustrate two different ways in which the interaction between mutation bias and composition bias influence the frequency of mutations in the population. In S18a Fig, polymorphisms are present, but they are rare, and the sequence that evolves to the highest frequency is the same as in the absence of mutation bias. In contrast, in S18b Fig some transition polymorphisms come to dominate the population, such that the sequence that evolves to the highest frequency is not the same as in the absence of mutation bias.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. Mutation bias interacts with composition bias to influence the evolution of genetic diversity and mutational robustness.

(a-e) The average genetic diversity and (f-j) mutational robustness of evolved populations at steady state is shown for 19 different values of the mutation bias parameter α, and for each of three different values of mutation supply Nμ (see legend). The solid vertical lines indicate no bias in mutation supply (α = 0.5). Landscapes are grouped based on their composition bias and the distribution of composition bias per panel is shown on top of each panel as in Fig 3.

https://doi.org/10.1371/journal.pcbi.1008296.g006

We then explored the potential consequences of these changes in genetic diversity. First, we characterized how mutation bias and composition bias interact to influence the mutational robustness of binding sites in evolved populations at steady state. We quantified the mutational robustness of an individual binding site as the fraction of all possible mutations to that site that created another site that was also part of the landscape [35, 37]. The mutational robustness of a population of binding sites was then simply calculated as the average mutational robustness of its constituent sites. Fig 6f–6j shows that mutation bias can increase or decrease the mutational robustness of an evolved population, relative to an unbiased mutation supply, especially in landscapes with strong composition bias, either toward transversions (Fig 6f) or transitions (Fig 6j). In contrast, for landscapes with little to no composition bias (Fig 6h), only the most extreme values of mutation bias influenced mutational robustness. For landscapes with strong composition bias, the changes in mutational robustness caused by mutation bias mirrored the changes observed for genetic diversity (compare Fig 6a and 6e to Fig 6f anf 6j). In these landscapes, therefore, a decrease in genetic diversity was associated with a decrease in mutational robustness. Moreover, we observed that populations evolving on landscapes with composition bias tended to be less mutationally robust at steady state. The reason is that landscapes with composition bias tended to comprise genotypes with fewer mutational neighbors, relative to landscapes without composition bias (S19 Fig). When mutation bias aligned with composition bias, those genotypes with more mutational neighbors were more likely to evolve, because the mutation bias oriented the population toward those genotypes.

Mutation bias influences the evolvability of polymorphic populations

Because mutational robustness is a cause of evolvability [21, 52, 53], we explored whether mutation bias and composition bias interact to influence evolvability—defined in this context as the ability of mutation to bring forth new binding phenotypes (Methods) [35]. To test this, we calculated the average number of transcription factors that bind the one-mutant neighbors of any individual in the population at steady state, and computed the difference between these averages for populations with a strong mutation bias toward transversions (α = 0.05), relative to an unbiased mutation supply (α = 0.5), as well as toward transitions (α = 0.95), relative to an unbiased mutation supply. S20 Fig shows how a bias in mutation supply affected the evolvability of polymorphic populations for 128 transcription factors from Arabidopsis thaliana and 128 transcription factors from Mus musculus, the two species with the most transcription factors in our dataset. In approximately 70% of the landscapes, mutation bias had no influence on evolvability. However, in the remaining 30% of landscapes, mutation bias could increase or decrease the number of transcription factors accessible via point mutation. In A. thaliana, this ranged from plus or minus 6 transcription factors (∼ 5%), whereas in M. musculus this ranged from plus 12 to minus 7 transcription factors (between ∼ 5% and ∼ 9%). These observations suggest that mutation bias is capable of orienting evolving populations both toward and away from more evolvable regions of genotype-phenotype landscapes. However, this effect is apparently independent of composition bias. This may seem counterintuitive at first glance, because the interaction between mutation bias and composition bias influences both genetic diversity and mutational robustness, two properties that facilitate evolvability [53]. However, the global peaks of the landscapes considered here are sufficiently narrow that any increase of diversity or robustness within them is not sufficient to cause a change in the number of phenotypes accessible via mutation. Said differently, the fraction of sequence space covered by these peaks is too small to facilitate mutational access to the landscapes of other transcription factors. We speculate that in landscapes with broader, mesa-like peaks, mutation bias and composition bias could interact to influence evolvability.

Data

We constructed genotype-phenotype landscapes using data from protein binding microarrays, which we downloaded from the UniPROBE [46] and CIS-BP [47] databases. These data include proxies for the relative binding affinity of a transcription factor to all possible (4⁸ − 4⁴)/2 + 4⁴ = 32,896 eight-nucleotide, double-stranded DNA sequences (transcription factor binding sites). These proxies include an enrichment score (E-score), which for each sequence is a function of the fluorescence intensities of a subset of probes that contain the sequence and a subset of probes that do not contain that sequence [66], and a Z-score, which for each sequence is the difference between the logarithm of the median fluorescence intensities of a subset of probes that contain the sequence and the logarithm of the median fluorescence intensities of all probes, reported in units of standard deviation. We considered a sequence as specifically bound by a transcription factor if its E-score exceeded 0.35, following previous work [22, 35, 37, 67]. We included a landscape in our dataset if its dominant genotype network comprised at least 100 bound sequences. According to these criteria, our dataset included 746 transcription factors, representing 129 eukaryotic species, and 48 DNA-binding domain structural classes. Details are provided in S1 Table.

Constructing and analyzing genotype-phenotype landscapes

For each transcription factor, we used the Genonets Server [68] to construct a genotype-phenotype landscape from the set of sequences that specifically bound the factor (i.e., with an E-score > 0.35). We represented each such sequence as a node in a genotype network, and connected nodes with edges if their corresponding sequences differed in a single point mutation. Note that we did not consider indels in our definition of a mutation, like we did in our previous work [22, 35–37]. The reason is that we were interested in understanding the influence of a form of composition bias defined by point mutations (transitions vs. transversions), and we were concerned that the inclusion of indel-based edges would confound our analyses. For genotype networks that were fragmented into multiple components, we only considered the largest component, which we call the dominant genotype network.

Each dominant genotype network served as the substrate of a genotype-phenotype landscape, the surface of which was defined by relative binding affinity. For our main analyses, this was captured by the E-score, whereas for some of our sensitivity analyses, this was captured by the Z-score. We studied accessible mutational paths to the global peaks of these landscapes. An accessible mutational path comprises edges that each confer an increase in binding affinity greater than the noise threshold parameter δ [22]. For each transcription factor, we calculated δ as the residual standard error of a linear regression between the affinity values of all bound sequences from the two replicate protein binding microarrays [22]. Thus, each transcription factor had its own δ, which reflects the noise in the replicated measurements for that particular transcription factor.

Mutation bias and composition bias

We report both mutation bias and composition bias relative to the null expectation that one transition occurs for every two transversions. Letting Ti and Tv represent mutation rates of transitions and transversions, respectively, we define mutation bias as (1) A mutation bias of α = 0.5 corresponds to the null expectation of one transition per two transversions. Values below 0.5 mean there are more transversions than expected under the null, whereas values above 0.5 mean there are more transitions than expected under the null.

Composition bias was measured in the same way, except with Ti and Tv representing the number of transitions and transversions in a landscape, or in an accessible mutational path to the global peak of a landscape.

Origin-fixation model of evolutionary dynamics

We used an origin-fixation model to study evolutionary dynamics when the mutation supply is low (Nμ << 1). This was implemented using Markov chains, a memoryless process that gives the jumping probability from one genotype i to another genotype j using the matrix (2) where f_i,j is the relative difference in binding affinity b (3) and (4)

Then in general, the probability of going from any state to another state in a Markov chain given by the matrix P (Eq (2)) after t steps is (5)

Wright-Fisher model of evolutionary dynamics

We carried out simulations of a Wright-Fisher model to study evolutionary dynamics when the mutation supply is high (Nμ > 1). Each simulation was initialized with a monomorphic population comprising N copies of the same sequence, chosen from the bottom 10% of binding affinity values in the landscape. In each generation t, N sequences were chosen from the population at generation t − 1 with replacement and with a probability that was linearly proportional to binding affinity. Mutations were introduced to these sequences at a rate μ per sequence per generation with a mutation bias α. For each of the 746 landscapes, we performed 15 replicate simulations for each initial condition, using 19 linearly spaced mutation bias values between 0.05 and 0.95, and 3 mutation supply values (Nμ ∈ {5, 20, 50}). Each simulation ran for 1000 generations, which was sufficient to ensure that the population had reached steady state.

Landscape navigability

As a measure of landscape navigability, we calculated the probability of reaching the global peak, starting from the 10% of sequences in the dominant genotype network with the lowest binding affinity. For low mutation supply, this was calculated as the average probability of going from the initial sequences to the global peak using Eq (5) after t = 1000 steps. For high mutation supply, this was calculated as the fraction of simulations per landscape in which at least 50% of the population reached the global peak.

Path entropy

PathMAN (Path Matrix Algorithm for Networks), is a publicly available Python script that efficiently calculates path statistics of a given Markov process [51]. We employed PathMAN to calculate the Shannon’s entropy of the path distribution, which accounts for the predictability of the process. Low entropy means that few paths with large probability dominate the process, while large entropy means that several low-probability paths contribute. We calculated the path entropy for 19 different values of the mutation bias parameter within the range [0.05,0.95], in order to find the mutation bias parameters α_min and α_max that minimize and maximize the path entropy for each landscape.

Overlap coefficient

The overlap coefficient between two different polymorphic populations A and B was calculated as (6) where A and B are multisets—sets that permit multiple instances of an element. The cardinality of such multisets is defined as (7) where the number of occurrences of the element x in the multiset is indicated by the multiplicity function m(x).

Then C is the multiset defined as C = A ∩ B, with multiplicity function (8) For example, if A = {1, 1, 2, 2, 2, 3} and B = {1, 2, 2, 4}, then C = {1, 2, 2} and the overlap coefficient is O_A,B = 0.75.

Quasispecies dynamics of infinite populations

Since the matrix in Eq (2) is non-negative and connected, the Perron-Frobenius theorem for non-negative matrices applies [69]. Hence, the steady-state distribution of an infinite size population on a genotype network is determined by the eigenvector associated to the largest eigenvalue of the matrix in Eq (2). Per each landscape, the eigenvectors were computed numerically for 19 different values of mutation bias parameter α within the range [0.05,0.95].

Genetic diversity

We measured the diversity of a population as the average Shannon’s diversity over all genotypes, normalized by the maximum diversity per landscape log₂(n): (9) where n is the number of sequences in the landscape, p_i is the fraction of the steady-state population that is at sequence i.

Evolvability

We quantified the evolvability of a transcription factor’s binding sites as follows. First, we determined the set of binding sites that have evolved at steady state in our Wright-Fisher simulations. Then we enumerated the set of DNA sequences that differed by one mutation from any of these binding sites, but were not themselves part of the focal transcription factor’s landscape. Finally, we determined the fraction of transcription factors in our dataset that these one-mutant neighbors bound. This fraction was our measure of evolvability.