On the incongruence of genotype-phenotype and fitness landscapes.

The mapping from genotype to phenotype to fitness typically involves multiple nonlinearities that can transform the effects of mutations. For example, mutations may contribute additively to a phenotype, but their effects on fitness may combine non-additively because selection favors a low or intermediate value of that phenotype. This can cause incongruence between the topographical properties of a fitness landscape and its underlying genotype-phenotype landscape. Yet, genotype-phenotype landscapes are often used as a proxy for fitness landscapes to study the dynamics and predictability of evolution. Here, we use theoretical models and empirical data on transcription factor-DNA interactions to systematically study the incongruence of genotype-phenotype and fitness landscapes when selection favors a low or intermediate phenotypic value. Using the theoretical models, we prove a number of fundamental results. For example, selection for low or intermediate phenotypic values does not change simple sign epistasis into reciprocal sign epistasis, implying that genotype-phenotype landscapes with only simple sign epistasis motifs will always give rise to single-peaked fitness landscapes under such selection. More broadly, we show that such selection tends to create fitness landscapes that are more rugged than the underlying genotype-phenotype landscape, but this increased ruggedness typically does not frustrate adaptive evolution because the local adaptive peaks in the fitness landscape tend to be nearly as tall as the global peak. Many of these results carry forward to the empirical genotype-phenotype landscapes, which may help to explain why low- and intermediate-affinity transcription factor-DNA interactions are so prevalent in eukaryotic gene regulation.

Introduction

Characterizing the relationship between genotype and phenotype is key to our understanding of evolution [1, 2]. For quantitative phenotypes, such as the expression level of a gene or the enzymatic activity of a protein, this relationship can be formalized as a genotype-phenotype landscape [3]. In such a landscape, genotypes represent coordinates in an abstract genotype space and their phenotype defines the elevation of each coordinate in this space [4]. The topographical properties of genotype-phenotype landscapes, such as their ruggedness, are influenced by epistasis [5]—non-additive interactions between mutations in their contribution to phenotype. These topographical properties have important evolutionary consequences, because they determine how mutation brings forth the phenotypic variation upon which selection acts [6, 7].

Technological advances are facilitating the construction and analysis of empirical genotype-phenotype landscapes at ever-increasing resolution, scale, and scope [8]. Example phenotypes include the enzymatic activity [9], binding affinity [10], allosteric profile [11], and fluorescence intensity of proteins [12], as well as exon inclusion levels [3], the expression levels of genes driven by regulatory elements [13], the expression patterns of gene regulatory circuits [14-16], and flux through metabolic pathways [17]. In analyzing the topographical properties of these landscapes and their evolutionary consequences, an assumption is often made that phenotype is a proxy for fitness [12, 18–24], thus rendering genotype-phenotype landscapes equivalent to fitness landscapes [25]. While this assumption may be justified under certain conditions, such as in directed protein evolution experiments [26], it is often the case that the relationship between phenotype and fitness is not so straightforward. For example, fitness may depend upon more phenotypes than those being assayed [27] or the relationship between phenotype and fitness may be inherently nonlinear, for example reflecting a tradeoff between the costs and benefits associated with a phenotype [28]. In the latter case, selection may favor a low or intermediate phenotypic value [29, 30]; e.g., an intermediate gene expression level, [31-33], enzyme efficiency [34] or protein production rate or activity [35]. Such non-linearities are a cause of epistasis [16, 36–38], and they can transform the effects of mutations as they map onto phenotype and fitness [37], thus rendering the topographical properties of a fitness landscape qualitatively different from those of its underlying genotype-phenotype landscape (Fig 1A). While we do not doubt that workers in the field are well aware that the topographical properties of a fitness landscape can differ from those of its underlying genotype-phenotype landscape, a systematic study of these differences is lacking.

Fig 1

Incongruence.

(A) Schematic illustration of how selection for an intermediate phenotypic value wopt can make a genotype-phenotype landscape incongruent with the resulting fitness landscape. (B) An additive three-locus, biallelic genotype-phenotype landscape with a single peak (gray filled circle; pgp = 1). One pairwise interaction is highlighted in blue. It exhibits no magnitude epistasis (ϵgp = 0) or sign epistasis. (C) The Gaussian phenotype-fitness map (Eq 1) is shown for three values of wopt (dashed line, wopt = 0; dotted line, wopt = 0.5; solid line, wopt = 1), with three values of σ shown for wopt = 0.5. (D) Applying the Gaussian phenotype-fitness map with wopt = 0 to the genotype-phenotype landscape results in a single-peaked fitness landscape (gray filled circle; pf = 1). The same pairwise interaction from (B) is highlighted in blue. It exhibits positive epistasis (ϵgp = 0.266), but no sign epistasis. (E) Applying the Gaussian phenotype-fitness map with wopt = 0.5 to the genotype-phenotype landscape results in a multi-peaked fitness landscape (gray filled circles; pf = 2). The same pairwise interaction from (B) is highlighted in red. It exhibits negative epistasis (ϵgp = −1.01), as well as reciprocal sign epistasis. (F) Applying the Gaussian phenotype-fitness map with wopt = 1 to the genotype-phenotype landscape results in a single-peaked fitness landscape (gray filled circle; pf = 1). The same pairwise interaction from (A) is highlighted in blue. It exhibits positive epistasis (ϵgp = 0.53), but no sign epistasis. In panels B, D-F, arrows point from genotypes with lower phenotypic or fitness values to genotypes with higher phenotypic or fitness values. The no sign epistasis motif is highlighted in blue and the reciprocal sign epistasis motif in red. Note the symmetry of landscapes in panels D and F.

Selection for low or intermediate phenotypic values is especially relevant to transcription factor-DNA interactions [39]. Transcription factors are sequence-specific DNA binding proteins that help regulate gene expression. They do so by binding DNA sequences (transcription factor binding sites) in regulatory regions such as promoters and enhancers to recruit or block the recruitment of RNA polymerase [40]. The regulatory effect of such a binding event depends in part on the affinity with which the DNA sequence is bound by the transcription factor [41, 42]. As such, binding affinity is an important molecular phenotype of transcription factor binding sites, upon which selection acts [43, 44]. While it is commonly assumed that selection increases binding affinity [21, 43–47], several lines of evidence suggest that selection for low or intermediate binding affinity also influences the evolution of transcription factor binding sites. For example, paralogous transcription factors often bind the same DNA sequences with high affinity, but different DNA sequences with low affinity [48, 49]. If an optimal gene expression pattern requires binding by just one of several transcription factor paralogs, such specificity can be achieved using low-affinity transcription factor binding sites, resulting in selection for low binding affinity [50]. Additional documented cases in which low-affinity binding sites play important regulatory roles include negative auto-regulation by high-copy number transcription factors in Escherichia coli [51], where high-affinity binding sites cause suboptimal noise suppression, and developmental patterning in Ciona intestinalis embryos [52, 53], where high-affinity binding sites cause deleterious ectopic gene expression patterns. Moreover, low-affinity binding sites are commonly observed in the regulatory portfolios of a diversity of organisms, including bacteria [51], yeast [54], fly [49, 55, 56], sea stars and sea urchins [57], as well as humans [58].

How incongruent are the topographies of genotype-phenotype and fitness landscapes when selection favors a low or intermediate phenotypic value? How does this depend on the ruggedness of the genotype-phenotype landscape? These are important questions, because the topography of a fitness landscape has implications for several evolutionary phenomena, including the evolution of genetic diversity [59], reproductive isolation [60], and sex [61], as well as the predictability of the evolutionary process itself [62]. How much we can learn about these phenomena from knowledge of a genotype-phenotype landscape depends on the genotype-phenotype landscape’s congruence with the fitness landscape. Despite decades of research on fitness landscapes [63, 64], these questions have not been addressed even in the context of classical theoretical models, such as Mt. Fuji [65], House-of-Cards [66], or NK landscapes [67]. They have also not been addressed in the context of biophysical models of genotype-phenotype landscapes or empirical genotype-phenotype landscapes. Here, we fill this knowledge gap by defining local and global measures of incongruence, which describe the topographical differences between a genotype-phenotype landscape and the corresponding fitness landscape when selection favors a low or intermediate phenotypic value. We use these measures to study incongruence in the context of the aforementioned theoretical models [65-67] and derive some fundamental results that are applicable to all empirical genotype-phenotype landscapes. We then consider the specific case of genotype-phenotype landscapes of transcription factor-DNA interactions, by first looking at an idealised biophysical model [68] and then taking a step further, by analysing 1,137 empirical genotype-phenotype landscapes, wherein genotypes are transcription factor binding sites and the phenotype is a measure of relative binding affinity [21]. We study transcription factor-DNA interactions because there is strong biological motivation for studying selection for low or intermediate binding affinity, as discussed above, and because a large number of empirical genotype-phenotype landscapes of transcription factor-DNA interactions are publicly available [21], thus facilitating the statistical analysis of their incongruence.

Results

We first present our measures of incongruence and the phenotype-to-fitness map used to incorporate the effect of selection for a low or intermediate phenotypic value. We use these measures and this map to study the incongruence of randomly generated genotype-phenotype landscapes, specifically two-locus biallelic landscapes and multi-locus biallelic landscapes. For the latter, we use NK landscapes [67], which include the corner cases of Mt. Fuji [65] (K = 0) and House-Of-Cards [66] (K = N − 1) landscapes (see Table 1 for a list of symbols). We then apply the principles learned from these model genotype-phenotype landscapes to genotype-phenotype landscapes of transcription factor-DNA interactions, first in the context of the mismatch model [68] (which enables us to study landscapes with more than two alleles per locus), and then in the context of empirical measurements of transcription factor-DNA interactions [21]. Finally, we study the evolutionary consequences of landscape incongruence.

Table 1

List of symbols.

w, wopt	Phenotypic value, Optimal phenotypic value
σ	Strength of selection
F(w)	Fitness value
ϵgp, ϵf	Epistasis in the genotype-phenotype landscape (gp) and fitness landscape (f) respectively
pgp, pf	Number of peaks in the genotype-phenotype landscape (gp) and fitness landscape (f) respectively
L	Length of the sequence
K	Ruggedness parameter of the NK model
a	Number of alleles at each locus
m	Number of mismatches in the mismatch model
〈l〉	Average length of adaptive walk
〈f〉	Mean fitness at equilibrium

Landscape incongruence

Our goal is to quantify the topographical differences between a fitness landscape and its underlying genotype-phenotype landscape, when selection favors a low or intermediate phenotypic value. To do so, we define measures of landscape incongruence, at both a local and a global scale. At a local scale, we quantify differences in pairwise epistasis amongst loci in the genotype-phenotype landscape relative to the same loci in the fitness landscape (Methods). We do so by classifying the type of magnitude epistasis (ϵ) as additive (i.e., no epistasis; ϵ = 0), positive (ϵ > 0), or negative (ϵ < 0). For a pair of loci, we then compare the type of epistasis in the genotype-phenotype landscape (ϵgp) to the type of epistasis in the fitness landscape (ϵf), and report whether this is the same in the two landscapes. We do the same for an additional classification of epistasis based on the absence or presence of sign epistasis. This results in three categories—no sign epistasis, simple sign epistasis, or reciprocal sign epistasis [69]. At a global scale, we quantify incongruence as the difference in the number of peaks in the fitness landscape (pf) relative to the genotype-phenotype landscape (pgp). Taken together, these local and global measures allow us to determine the extent to which selection for a low or intermediate phenotypic value increases or decreases the ruggedness of the fitness landscape, relative to the genotype-phenotype landscape.

Phenotype-to-fitness map

To study the effect of selection for low or intermediate phenotypic values, we use a Gaussian phenotype-to-fitness map F(w) centred around an optimal phenotypic value wopt,

The parameter σ determines the strength of selection by controlling how rapidly fitness decreases as the phenotype w deviates from the optimal phenotype wopt. Increasing σ decreases the strength of selection, because it broadens the fitness map around wopt and thus decreases the fitness differences between similar phenotypes. Similar maps are commonly used in evolutionary modeling frameworks, such as Fisher’s Geometric model [70, 71] and models of speciation [72], as well as in biophysical models of intermolecular interactions [29], including transcription factor-DNA interactions [46]. However, we emphasize that most of our measures of incongruence depend only on the rank ordering of phenotypic or fitness values, and are therefore independent of the exact shape of the phenotype-fitness map, so long as it is symmetric.

Fig 1B–1F illustrates the application of the phenotype-to-fitness map to a simple three-locus, biallelic genotype-phenotype landscape. This landscape is purely additive, and as a consequence, it exhibits no epistasis and has only one peak (Fig 1B). Applying the phenotype-to-fitness map (Fig 1C) to this landscape can change the amount and type of epistasis, as well as the location and number of peaks in the resulting fitness landscape, relative to the genotype-phenotype landscape (Fig 1D–1F). It can therefore cause incongruence between genotype-phenotype and fitness landscapes. Whereas this schematic and our analyses below pertain to a single phenotype, extending our model to multiple phenotypes is straightforward, as we later discuss.

Two-locus biallelic genotype-phenotype landscapes

We first study incongruence using the simplest form of genotype-phenotype landscape that is capable of exhibiting epistasis: a two-locus biallelic landscape. We represent genotypes as binary strings of length L = 2 and randomly assign a phenotype w to each genotype i, which we draw from a uniform distribution between 0 and 1. We then apply the Gaussian phenotype-fitness map (Eq 1) to generate the corresponding fitness landscape. We repeat this process 10,000 times for values of wopt ∈ [0, 1] (in increments of 0.01), and report the probability that the type of epistasis in the genotype-phenotype landscape is the same as in the fitness landscape. We first differentiate between no magnitude epistasis, positive epistasis, and negative epistasis, and then between no sign epistasis, simple sign epistasis, and reciprocal sign epistasis.

Incongruence in magnitude epistasis is highest when selection favors low phenotypic values

To determine whether the type of magnitude epistasis is the same in the genotype-phenotype landscape and fitness landscape, we calculate the product ϵgp ⋅ ϵf, which will be positive when the two landscapes have the same type of epistasis. We use this product, rather than a discrete categorization, to ensure analytical tractability. We obtain an analytical expression for ϵgp ⋅ ϵf by assuming σ to be large (S1 Appendix, Derivation 1): where w is the phenotypic value of genotype i with i ∈ {0, 1}2 and ϵ = w00w11 − w01w10, which is also known as multiplicative epistasis [73].

The probability that the type of epistasis is the same in the genotype-phenotype and fitness landscapes, P(ϵgp ⋅ ϵf > 0) can be computed using Monte-Carlo methods and yields results in good agreement with the randomly generated genotype-phenotype landscapes when σ is large. This is shown in Fig 2, where three trends are immediately apparent. First, for large σ, the probability that the type of epistasis is the same in the genotype-phenotype landscape and the fitness landscape increases as the optimal phenotype wopt increases, in agreement with the intuition that selection for large phenotypic values leaves the genotype-phenotype landscape mostly unchanged, except for the nonlinear rescaling introduced by the phenotype-fitness map. Second, even when wopt = 1, the probability that the type of epistasis is the same in the genotype-phenotype landscape and the fitness landscape is less than one. The reason is the nonlinear rescaling introduced by the phenotype-to-fitness map does not guarantee conservation of the type of magnitude epistasis, even though it does preserve the rank ordering of fitness values. Results obtained with the randomly generated genotype-phenotype landscapes show this effect becomes even more pronounced as σ decreases. This means that as the strength of selection for wopt increases, so does the likelihood of landscape incongruence. Finally, as σ decreases, the probability of retaining the type of epistasis from the genotype-phenotype landscape in the fitness landscape is not maximized at wopt = 1, but rather at a smaller wopt (e.g., at wopt ≈ 0.8 when σ = 0.01, Fig 2A). The reason is as σ decreases, the fitness function becomes extremely narrow and more phenotypic values are mapped to zero fitness (within computer precision i.e. any fitness <5.0 × 10−324 ≈ 0), resulting in cases where the genotype-phenotype landscape exhibits epistasis (i.e., ϵgp ≠ 0), but the fitness landscape does not (i.e., ϵf = 0), because all fitness values are zero. Fig 2B shows this is more likely to occur when σ is small and as wopt approaches its extreme values of 0 or 1. In these randomly generated landscapes, the probability of obtaining negative epistasis is the same as the probability of obtaining positive epistasis, and the probability of conserving the type of epistasis is independent of the type of epistasis in the genotype-phenotype landscape. In sum, these results show that selection for low or intermediate phenotypic values can modify the genotype-phenotype landscape, such that the resulting fitness landscape exhibits a different type of magnitude epistasis, and this effect is most pronounced when selection is strong and the optimal phenotypic value is low.

Fig 2

Local incongruence: Magnitude epistasis.

(A) The probability of retaining the type of epistasis, shown in relation to the optimal phenotype wopt. The black line shows the theoretical prediction and the dots show the results from randomly generated genotype-phenotype landscapes for different values of σ. The theoretical approximation agrees well with the results from randomly generated genotype-phenotype landscapes for large σ (i.e., σ ≥ 1). (B) The probability of observing zero magnitude epistasis in the fitness landscape, shown in relation to wopt, for different values of σ, which we selected to show the range of variation in the probability of observing zero epistasis.

Incongruence in sign epistasis is highest when selection favors intermediate phenotypic values

Next, we categorized landscapes as exhibiting no sign epistasis, simple sign epistasis, or reciprocal sign epistasis. These three motifs are shown in the center panel of Fig 3, where arrows point from genotypes with a lower phenotypic or fitness value to genotypes with a higher phenotypic or fitness value. Because the presence of sign epistasis only depends upon the partial ordering of the phenotypic values, we expect to retain the motif from the genotype-phenotype landscape in the fitness landscape as wopt → 1. We also expect to retain the motif as wopt → 0, because selecting for wopt = 0 simply flips all the arrows from the genotype-phenotype landscape in the fitness landscape, which does not change the categorization of the motif. Thus, we expect the probability of retaining the motif from the genotype-phenotype landscape in the fitness landscape to be “U” shaped and symmetric about wopt = 0.5. Fig 3 shows the probability of retaining or changing the motif from the genotype-phenotype landscape in the fitness landscape for 10,000 randomly generated two-locus biallelic landscapes, grouped according to the motif in the genotype-phenotype landscape. These results confirm the expected “U” shape of the probability of retaining the type of epistasis from the genotype-phenotype landscape in the fitness landscape as wopt is varied from 0 to 1.

Fig 3

Local incongruence: Sign epistasis.

The probability of retaining or changing the type of epistasis in the genotype-phenotype landscape, relative to the fitness landscape, shown in relation to the optimal phenotypic value wopt. Data are grouped based on whether the genotype-phenotype landscapes exhibits (A) no sign epistasis (blue), (B) simple sign epistasis (brown), or (C) reciprocal sign epistasis (red). The colours of the lines represent the type of epistasis in the resulting fitness landscape. These results are independent of σ, because they only depend on the rank ordering of fitness values. Notice the “U” shape of the probability of retaining the type of epistasis in each panel.

Simple sign epistasis cannot be modified into reciprocal sign epistasis

When the genotype-phenotype landscape has the no sign epistasis motif, selection for an intermediate phenotypic value can transform the landscape into the simple sign epistasis motif or the reciprocal sign epistasis motif (Fig 3A), in line with recent results on Fisher’s Geometric model [74]. When the genotype-phenotype landscape has the simple sign epistasis motif, selection for an intermediate phenotypic value can transform the landscape into a no sign epistasis motif, but not into a reciprocal sign epistasis motif (Fig 3B and S1 Appendix, Proof 1). Because reciprocal sign epistasis is a necessary condition for multiple peaks [5], this implies that genotype-phenotype landscapes with only simple sign epistasis motifs will always give rise to single peaked fitness landscapes, using the phenotype-fitness map considered here. Finally, when the genotype-phenotype landscape has the reciprocal sign epistasis motif, selection for an intermediate phenotypic value transforms the landscape into the no sign epistasis motif or the simple sign epistasis motif with equal probability (Fig 3C and S1 Appendix, Proof 2). Moreover, the probability of retaining the motif from the genotype-phenotype landscape in the fitness landscape is always higher than the probability of changing it when the genotype-phenotype landscape has the simple sign epistasis motif (Fig 3B) or the reciprocal sign epistasis motif (Fig 3C). This is not always true when the genotype-phenotype landscape has the no sign epistasis motif (Fig 3A), because at intermediate wopt the landscape is most likely to transform into the reciprocal sign epistasis motif. Taken together, these results show that selection for intermediate phenotypic values can modify genotype-phenotype landscapes with no sign epistasis into fitness landscapes with sign epistasis and vice versa.

The inferences about pairwise interactions can be carried forward to multi-locus biallelic landscapes because their genotype spaces, which are L-dimensional hypercubes, are composed of two-dimensional squares. Due to the adjacency of squares, in the three-locus case, the motifs of four out of the six squares are sufficient to determine the motifs of the rest of the squares, and for any L, the motifs of only 2 ⋅ (L − 1) squares are necessary to determine the motifs of all of the remaining squares. This is only a fraction 2/L of all the squares in the hypercube (because the total number of faces in an L-dimensional hypercube is ), which is clearly minuscule for large L. However, pairwise interactions are not sufficient to predict peak patterns, which may result from higher-order interactions [75]. We study these in the next section.

Multi-locus biallelic genotype-phenotype landscapes

We use the NK model [67] to study multi-locus biallelic genotype-phenotype landscapes (Methods). In this model, each locus in a genotype of length L epistatically interacts with K other loci (whereas N is typically used to denote the number of loci in this model, we use L for consistency with the rest of our text). As corner cases, this model includes Mt. Fuji landscapes [65] when K = 0 and House-of-Cards landscapes [66] when K = L − 1. For each combination of L and K, we use this model to randomly generate a genotype-phenotype landscape. We then apply the Gaussian phenotype-fitness map (Eq 1) to generate a fitness landscape. We repeat this process 10,000 times for wopt ∈ [0, 1] (in increments of 0.01), and report the average of the absolute change in the number of peaks, i.e., 〈|pf − pgp|〉, where pf is the number of peaks in the fitness landscape and pgp is the number of peaks in the genotype-phenotype landscape. This is our measure of global incongruence. We use the absolute value of the change in number of peaks so that we can average over many realisations of genotype-phenotype landscapes. However, since the sign of change is also important, we discuss that as well in the following sections.

Mt. Fuji landscapes

We begin with Mt. Fuji genotype-phenotype landscapes. Because these are single-peaked, selection for a low or intermediate phenotypic value can only maintain or increase the number of peaks from the genotype-phenotype landscape in the fitness landscape. Fig 4A shows this change in the number of peaks, and Fig 4B shows the probability that the number of peaks changes, in relation to wopt for landscapes with L = 2 to L = 8 loci. These trends are symmetric about wopt = 0.5, because Mt. Fuji landscapes are additive, so selecting for wopt = 0 is equivalent to selecting for wopt = 1 with regard to the change in the number of peaks. The reason is that selecting for wopt = 0 flips all of the arrows in the fitness landscape, relative to the genotype-phenotype landscape, which changes the location of the peak, but does not change the number of peaks. This is illustrated in Fig 1B and 1D. For a detailed explanation of the shape of the curves in Fig 4A, see S1 Appendix, Note 1. More obvious is the increase in the number of peaks, and the probability that the number of peaks increases, as L increases, the latter converging to one for all values of wopt except the extreme cases of wopt = 0 and wopt = 1. Note, however, that a high probability of increase in the number of peaks does not necessarily correspond to a high increase in the number of peaks, as can be seen from the different positions of the maxima in Fig 4A and 4B (S1 Appendix, Note 1). For large L (>10), the expected number of peaks in the fitness landscape increases exponentially with L [74]. In sum, these results show that selection for intermediate phenotypic values readily transforms Mt. Fuji genotype-phenotype landscapes, which are smooth and single-peaked, into rugged fitness landscapes, and that this effect is most pronounced for large L and intermediate wopt.

Fig 4

Global incongruence: Mt. Fuji and House-of-Cards genotype-phenotype landscapes.

The absolute change in the number of peaks and the probability that the number of peaks changes in the fitness landscape, relative to (A,B) Mt. Fuji and (C,D) House-of-Cards genotype-phenotype landscapes, shown in relation to wopt for L ∈ {2, 3..8}. These results are independent of σ, because they only depend on the rank ordering of fitness values.

House-of-Cards landscapes

We next study House-of-Cards genotype-phenotype landscapes. These landscapes are highly rugged, with an average of peaks [76], whereas the maximum possible number of peaks is 2. As such, selection for a low or intermediate phenotypic value can either increase or decrease the number of peaks in the fitness landscape, relative to the genotype-phenotype landscape. Fig 4C shows this change in the number of peaks, and Fig 4D shows the probability that the number of peaks changes, in relation to wopt for landscapes with L = 2 to L = 8 loci. The change in the number of peaks is symmetric about wopt = 0.5 for the two-locus case, where the number of peaks does not change as wopt → 0 or wopt → 1. However, this symmetry is lost for L > 2. The reason is that although the phenotype-fitness map flips all of the arrows in the fitness landscape when wopt = 0, relative to the genotype-phenotype landscape, this does not guarantee conservation of the number of peaks. The number of peaks is jointly determined by the adjacent faces of the hypercube and thus, only very specific changes in the directions of arrows in the genotype-phenotype landscape guarantees conservation of the number of peaks in the fitness landscape (S4 Fig). However, in contrast to Mt. Fuji genotype-phenotype landscapes, the magnitude of change in the number of peaks increases very little with L, despite an exponential increase in the maximum number of possible peaks. Moreover, the probability that the number of peaks changes is still less than one for large L. Finally, for large L, both the change in the number of peaks and the probability that the number of peaks changes are independent of wopt, so long as wopt is sufficiently less than one. This observation depends on the probability distribution used to generate these landscapes. Because the 2 phenotypes in the NK model are drawn from a uniform distribution, nearly the same number of these phenotypes will be close to the optimal phenotype wopt, so long as L is sufficiently large. Thus, averaging over all possible configurations of the genotype-phenotype landscape yields the same value of 〈|pf − pgp|〉 for every wopt. For sufficiently large L, this value is given by (S1 Appendix, Derivation 2):

So far we have focused on the absolute change 〈|pf − pgp|〉 in the number of peaks in the fitness landscape, relative to the genotype-phenotype landscape. For House-of-Cards genotype-phenotype landscapes, selection for a low or intermediate phenotypic value can either increase or decrease the number of peaks. We were therefore interested in finding out which outcome is more likely. While one might expect a decrease to be more likely, due to the extreme ruggedness of House-of-Cards genotype-phenotype landscapes, we find that the number of peaks is equally likely to increase or decrease (S1 Appendix, Proof 3). In sum, these results show that in House-of-Cards genotype-phenotype landscapes, selection for a low or intermediate phenotypic value increases or decreases the number of peaks in the fitness landscape with equal probability, and the severity as well as the probability of this change increases with L and is largely independent of wopt.

Global incongruence decreases as the ruggedness of the genotype-phenotype landscape increases

Finally, we study NK genotype-phenotype landscapes, which bridge the gap between Mt. Fuji and House-of-Cards landscapes in terms of ruggedness, as K increases from 0 to L − 1. Fig 5 shows the absolute change in the number of peaks in the fitness landscape relative to the genotype-phenotype landscape for genotypes of length L = 5 and L = 8, as K is increased from 0 to L − 1. Note the gradual transition from the trends observed for Mt. Fuji genotype-phenotype landscapes (Fig 4A) to those observed for House-of-Cards genotype-phenotype landscapes (Fig 4C) as K increases. From these trends, we conclude four principles of how the ruggedness of an NK genotype-phenotype landscape influences its incongruence with the fitness landscape. As an NK genotype-phenotype landscape becomes more rugged, incongruence (1) loses symmetry about wopt = 0.5, (2) becomes less sensitive to wopt, and (3) decreases in severity, at least in terms of the absolute change in the number of peaks. Finally, the probability of increasing the number of peaks is always greater than or equal to the probability of decreasing the number of peaks, with the equality holding for House-of-Cards genotype-phenotype landscapes. This last principle is both intuitive and informative—it tells us that on average, selection for low or intermediate values is more likely to increase the ruggedness of a fitness landscape, relative to the genotype-phenotype landscape. Thus, selection for low or intermediate values is more likely to break than to create phenotypic correlations between mutationally similar genotypes as they map onto fitness, rendering fitness landscapes more rugged than their underlying genotype-phenotype landscapes. In the subsequent sections we address whether and how these results apply to genotype-phenotype landscapes with more than two alleles per locus, specifically in the context of a biophysical model and experimental measurements of transcription factor-DNA interactions.

Fig 5

Global incongruence: NK genotype-phenotype landscapes.

The absolute change in the number of peaks in the fitness landscape, relative to the genotype-phenotype landscape, is shown in relation to wopt for genotypes of length (A) L = 5 and (C) L = 8, as K increases from zero to L − 1. The corresponding probability of change in the number of peaks is shown in relation to wopt for genotypes of length (B) L = 5 and (D) L = 8, as K increases from zero to L − 1. These results are independent of σ, because they only depend on the rank ordering of fitness values.

Genotype-phenotype landscapes of transcription factor-DNA interactions

Motivated by the common usage of low- and intermediate-affinity transcription factor binding sites in the regulatory portfolios of a diversity of organisms [49, 51, 54–57], we now study the incongruence of genotype-phenotype landscapes of transcription factor-DNA interactions and the corresponding fitness landscapes generated after selection for low or intermediate phenotypic values. In these landscapes, genotypes represent DNA sequences—transcription factor binding sites—and the phenotype of a DNA sequence is the affinity with which it binds a transcription factor [21]. Because the regulatory effects of transcription factor-DNA interactions are partly determined by binding affinity [41, 42] and mutations to transcription factor binding sites can alter binding affinity [48, 77], the topographies of genotype-phenotype landscapes of transcription factor-DNA interactions have important implications for the evolution of gene regulation [21]. We study these landscapes using both a biophysical model and experimental measurements of transcription factor-DNA interactions. We focus on transcription factor binding sites of length L = 8, because this is the length of the binding sites assayed by protein binding microarrays [77, 78]—the data used to construct the empirical genotype-phenotype landscapes of transcription factor-DNA interactions [21].

The mismatch model

We first study genotype-phenotype landscapes of transcription factor-DNA interactions generated using the so-called mismatch model [46, 47, 68]. The key assumption of this model is that the binding energy of a DNA sequence is a linear function of the number of mismatches between the sequence (genotype) and a transcription factor’s consensus sequence—the sequence it binds with the highest affinity. Further, each mismatch is assumed to have the same energetic cost and these costs combine additively to determine binding energy. This model results in a Mt. Fuji-like, permutation-invariant genotype-phenotype landscape, wherein the phenotype only depends on the number of differences between the genotype and the consensus sequence, but not on which loci in the genotype differ from the consensus sequence. Although this is a simplified model, it provides an opportunity to study the effects of having more than two alleles per locus and serves as a bridge to our analyses of empirical transcription factor-DNA interactions.

To ensure that these results are comparable to the results on the empirical landscapes, we consider the negative of the binding energy as our phenotype, such that sequences that bind more strongly are assigned higher phenotypic values. We assume a phenotypic value of 1 for the genotype that is identical to the consensus sequence. For each mismatch between a genotype and the consensus sequence, we deduct a small positive value e, such that the phenotypic value of a genotype with m mismatches is A = 1 − m ⋅ e. Due to the permutation invariance, the genotype-phenotype landscape is highly degenerate, such that the number of genotypes N in a mismatch class m is distributed according to the asymmetric binomial distribution: where a is the number of alleles per locus (a = 4 for transcription factor binding sites, because they are DNA sequences). Fig 6A shows this distribution for transcription factor binding sites of length L = 8. Note that N is maximized at m = 6.

Fig 6

Global incongruence: Genotype-phenotype landscapes of transcription factor-DNA interactions.

Landscapes constructed using (A,B) the mismatch model and (C,D) experimental measurements from protein binding microarrays for 1,137 eukaryotic transcription factors. (A) The number of genotypes, shown in relation to mismatch class. (B) The absolute change in the number of peaks in the fitness landscape, relative to the genotype-phenotype landscape, shown in relation to the optimal mismatch class mopt. Labels indicate the number of genotypes per peak in the fitness landscape. Note the symmetry in the absolute change in the number of peaks around mismatch class mopt = 4, as well as the tripling of the number of genotypes per peak for each increment in mopt. The grey shaded circles are a schematic representation of the growing width of the peaks. (C) The number of genotypes per binding affinity class, where protein binding microarray E-scores are used as a proxy for relative binding affinity. Violin plots show the distribution, and box-and-whisker plots the 25–75% quartiles, across genotype-phenotype landscapes for the 1,137 transcription factors. Closed symbols and the dashed line denote the median of each distribution. (D) Violin plots of the distribution of the absolute change in the number of peaks in the fitness landscape, relative to the genotype-phenotype landscape, shown in relation to the optimal binding affinity wopt for σ = 0.15. The inset shows the mean absolute change in the number of peaks, in relation to wopt. The x-axes in (A,B) are arranged such that binding affinity increases when read from left to right, in qualitative agreement with the x-axes in (C,D). The results in panels A-C are independent of σ.

Fig 6B shows the incongruence between genotype-phenotype landscapes and fitness landscapes of transcription factor-DNA interactions constructed using the mismatch model, under selection for an optimal mismatch class mopt, reported in terms of the absolute change in the number of peaks in the fitness landscape, relative to the genotype-phenotype landscape. Based on our analysis of Mt. Fuji genotype-phenotype landscapes, we anticipated asymmetric incongruence about the mismatch class mopt = 6, because the distribution of the number of genotypes per mismatch class is asymmetric with a maximum at mismatch class m = 6 (Fig 6A) and we expected all genotypes in the optimal mismatch class mopt to be peaks. However, we observe symmetric incongruence about mopt = 4 (Fig 6B). This occurs because a > 2, which renders some genotypes in the same mismatch class mutational neighbours. Consequently, the peaks can be broad and include more than one genotype, thus resembling plateaus. Specifically, each genotype has (a − 1) − 1 mutational neighbours that are in the same mismatch class and are therefore part of the same peak. This leaves clusters of genotypes to be peaks in each mismatch class m, an expression that is maximized with m = 4 when L = 8, thus forcing symmetry in the absolute change in the number of peaks about mopt = 4. However, the width of the peaks increases as (a − 1), leading to a tripling of peak width for each increment in m (Fig 6B). Thus, even in this idealised genotype-phenotype landscape, features of empirical landscapes of TF-DNA interactions begin to emerge, such as broad peaks.

Empirical landscapes

We now study genotype-phenotype landscapes of transcription factor-DNA interactions generated using experimental data from protein binding microarrays [78] (Methods). For all possible DNA sequences of length L = 8, these data include an enrichment score (E-score) ranging from -0.5 to 0.5 that serves as a proxy for relative binding affinity, such that higher E-scores correspond to higher binding affinities [77, 78]. We have previously used these data to construct genotype-phenotype landscapes for 1,137 eukaryotic transcription factors, in which the surface of each landscape was defined by the E-score [21]. Due to limitations in the reproducibility of E-scores across microarray designs for genotypes that are bound non-specifically or with very low affinity [77, 78], each genotype-phenotype landscape only includes DNA sequences with an E-score exceeding a threshold of 0.35, which corresponds to a false discovery rate of 0.001 [77]. As shown in our previous work [21], these landscapes tend to exhibit little, if any, reciprocal sign epistasis and therefore comprise few peaks. As such, they bear resemblance to the genotype-phenotype landscapes constructed using the mismatch model. An important difference, however, is that genotypes in the lower half of the E-score range are omitted from each landscape due to the reproducibility issues mentioned above. Fig 6C shows the distributions of the number of genotypes across all 1,137 genotype-phenotype landscapes, grouped into six binding affinity classes. Whereas the lowest binding affinity class contains the most genotypes, we cannot determine if this is the true maximum, because we do not know what these distributions look like for lower binding affinity classes. However, assuming the energetic contribution of each binding site to be additive [68], we expect lower binding affinity classes to have fewer genotypes and the maximum to occur at an intermediate binding affinity class, as can be seen in the mismatch model and other models in literature [43].

Fig 6D shows the incongruence between these 1,137 empirical genotype-phenotype landscapes and their corresponding fitness landscapes, under selection for an optimal binding affinity wopt, reported in terms of the absolute change in the number of peaks in the fitness landscape, relative to the genotype-phenotype landscape for σ = 0.15. Interestingly, the effect of changing σ is small (S5B Fig) because increasing σ not only decreases the range of variation of fitness values, but also decreases the uncertainty in these values (Methods), leaving the number of peaks largely unchanged. The mean and variance in the absolute change in the number of peaks decreases as wopt increases to 0.5, in line with the intuition that selection for high wopt generates fitness landscapes that are topographically similar to the underlying genotype-phenotype landscape. More precisely, the percentage of landscapes that show a change in the number of peaks decreases from around 99% to 40%, as we go from wopt = 0.35 to wopt = 0.5. As anticipated from our results with additive genotype-phenotype landscapes (Mt. Fuji and the mismatch model), selection for low or intermediate phenotypic values is more likely to increase than to decrease the number of peaks, although the relative fraction of increase depends upon wopt—while around 99% of the landscapes show an increase in peaks for wopt = 0.35, this value decreases to around 40% for wopt = 0.5. Further, when we separately analyzed the single-peaked (≈ 66.75%) and multi-peaked (≈ 33.25%) genotype-phenotype landscapes, we found the single-peaked landscapes to be more incongruent (S5A Fig), in line with our results from the previous section. In sum, these results show that selection for low or intermediate phenotypic values tends to increase the ruggedness of a fitness landscape, relative to the underlying genotype-phenotype landscapes, rendering genotype-phenotype landscapes a poor proxy for fitness landscapes under such selection.

Mismatch model and empirical landscapes show different kinds of global incongruence

In contrast to genotype-phenotype landscapes constructed using the mismatch model, the incongruence of genotype-phenotype landscapes constructed using protein binding microarray data is highest for the lowest binding affinity class, rather than an intermediate class. There are three non-mutually exclusive explanations for this. First, these empirical landscapes are not purely additive [21], unlike the landscapes constructed with the mismatch model, so we do not expect perfect symmetry about an intermediate wopt. Second, as previously mentioned, protein binding microarray data do not capture the full range of binding affinity, so the lowest binding affinity class in our data (E-score = 0.35), which contains the most genotypes (Fig 6C), is unlikely to be the lowest binding affinity class, but rather an intermediate binding affinity class. Third, while the binding affinity of a sequence is highly correlated with the binding affinities of its mutational neighbors [79], this correlation is not perfect, so neighboring genotypes that are in the same mismatch class may not have sufficiently similar binding affinities to be considered part of the same peak, unlike in the mismatch model.

There are two additional differences between the incongruence of the empirical genotype-phenotype landscapes and those constructed using the mismatch model that are worth highlighting. First, in the empirical landscapes, the average height of the peaks is maximised when selecting for intermediate binding affinities and is lowest when selecting for the highest affinity (S6A Fig), whereas peak height is independent of mopt in the mismatch model. Second, in the empirical landscapes, peak width is maximized when selecting for the highest binding affinity (wopt = 0.5) (S6B Fig), whereas it is maximized when selecting for the lowest binding affinity (mopt = 8) in the mismatch model (Fig 6B). These two differences in incongruence are important for understanding evolutionary simulations on these landscapes, which are the focus of the next section.

Evolutionary consequences

A key finding of our analyses so far is that selection for low or intermediate phenotypic values is more likely to increase than to decrease the number of peaks in the fitness landscape, relative to the genotype-phenotype landscape. Since the ruggedness of a fitness landscape has important implications for a diversity of evolutionary phenomena [59-62], we now study the evolutionary consequences of this finding. We do so using two metrics: (1) the length 〈l〉 of a greedy adaptive walk, averaged over all possible genotypes as starting points, and (2) the mean population fitness at equilibrium 〈f〉 under deterministic mutation-selection dynamics (i.e., assuming an infinite population size). These metrics provide complementary information to measures of landscape ruggedness, such as the number of peaks. To see why, consider our results with the mismatch model, in which the fitness landscape had one or more global peaks of equal height (Fig 6A and 6B). We observed single-peaked fitness landscapes when the optimal mismatch class was mopt = 0 or mopt = 8, and the most rugged fitness landscapes when the optimal mismatch class was mopt = 4. We also observed a tripling of the number of genotypes per peak as mopt increased from 0 to 8, such that the landscape with mopt = 0 comprised a single peak with one genotype, the landscape with mopt = 4 comprised 70 peaks with 81 genotypes per peak, and the landscape with mopt = 8 comprised 1 peak with 6,561 genotypes. Which landscape topography is most conducive to adaptive evolutionary change?

NK landscapes

On NK landscapes, the length of the greedy adaptive walk 〈l〉 is minimised when wopt = 0.5 for K < L − 1, whereas for K = L − 1 (House-of-Cards landscapes), 〈l〉 is independent of wopt. Moreover, the change in mean fitness at equilibrium is always positive, increases with L and K, and tends to be maximized at wopt = 0.5 (S1 Appendix, Note 2). In sum, these results show that on biallelic landscapes, selection for an intermediate phenotypic value decreases the length of a greedy adaptive walk and increases mean fitness at equilibrium, despite increasing the overall ruggedness of the fitness landscape, relative to the genotype-phenotype landscape. Whereas previous work has shown that peak accessibility increases with alphabet cardinality (a) due to the existence of indirect paths [10, 80], we show below that despite this increased accessibility, results qualitatively similar to the biallelic case (a = 2) also hold for the multi-allelic case of the mismatch model and of the empirical landscapes of TF-DNA interactions (a = 4).

Mismatch model

In the fitness landscapes generated using the mismatch model, the length of the greedy adaptive walk averaged over all starting genotypes is given by

Fig 7A shows this expression for L = 8 and a = 4. Due to the additive nature and degeneracy of this genotype-phenotype landscape, the fitness landscapes have monotonically decreasing fitness values as the mutational distance from any peak increases. Therefore, the length of any individual greedy adaptive walk is simply the absolute difference of the mismatch class of the starting genotype and that of the optimal mismatch class (i.e., |mopt − m|). As such, the average length of the greedy adaptive walk is minimized when selecting for mopt = 6 (Fig 7A), because this maximizes the number of genotypes in adaptive peaks (Fig 6B; 28 global peaks × 729 genotypes per peak). If instead, we chose to start the greedy walk from only non-peak genotypes, the walk would still be minimised for mopt = 6, because this class has the maximum number of genotypes that are Hamming distance one away from the peak genotypes. Recall that in each fitness landscape generated upon selection for mopt, all the peaks are of the same height and therefore, the greedy walk always terminates on a global peak. This is in contrast with the NK landscapes, and as we will see below, the empirical landscapes.

Fig 7

Rugged fitness landscapes need not impede adaptation.

The average length of an adaptive walk 〈l〉 and the change in mean population fitness at equilibrium 〈f〉 is shown for landscapes of transcription factor-DNA interactions generated using (A,C) the mismatch model for e = 0.05 and (B,D) protein binding microarray data. In (B,D), violin plots show the distribution, and box-and-whisker plots the 25–75% quartiles, across the 1,137 empirical landscapes for each optimal binding affinity wopt. The large variability of 〈l〉 at intermediate and high wopt is a consequence of the random diffusion of the population on non-peak plateaus, which results in longer walks. In (C), μ = 0.1, σ = 1 while in (D) μ = 0.1, σ = 0.15.

For each mopt, we next calculated the change in mean fitness at equilibrium under deterministic mutation-selection dynamics, relative to the fitness landscape generated for mopt = 0. To do so, we exploit the permutation-invariance of these landscapes to group genotypes into a lower-dimensional state space defined by mismatch class (Methods). Specifically, we construct a transition matrix that defines the probability that a genotype from one mismatch class mutates into another, based on the frequency and fitness of genotypes in each mismatch class. We iterate this matrix until we reach steady state, which is guaranteed by the Frobenius-Perron theorem to be independent of initial conditions, because the matrix is irreducible [81]. Fig 7C shows the change in mean fitness at equilibrium in relation to the optimal mismatch class mopt = 0. Our first observation is that the change in mean fitness is always positive for mopt > 0 when μ = 0.1 and σ = 1, even though selection for such mismatch classes always causes an increase in the number of peaks in the fitness landscape, relative to the underlying genotype-phenotype landscape (Fig 6B). Our second observation is that the change in mean fitness relative to mopt is always unimodal. For example, selection for mopt = 6 maximizes the change in mean fitness when μ = 0.1 and σ = 1 (Fig 7C). More generally, the mopt that maximizes the change in mean fitness depends on the interplay between the mutation rate μ and the strength of selection σ. See S10 Fig for the phase diagram. When selection is strong or mutation is weak, mismatch class mopt = 8 maximizes the change in mean fitness. As the strength of selection decreases or the mutation rate increases, the mismatch class that maximizes the change in mean fitness decreases from mopt = 8 to mopt = 7 and then to mopt = 6, where it remains for weak selection or high mutation rates. This is because when σ is small (i.e., selection is strong), it is costly to step down from a peak and therefore, selecting for the class with the broadest peak (mopt = 8) leads to the highest equilibrium mean population fitness. As σ increases and selection becomes weaker, it is no longer as costly to step down from a peak and therefore, selecting for the class with the maximum number of genotypes in peaks (mopt = 6) maximizes mean fitness at equilibrium. While this effect is reminiscent of the “survival of the flattest” [82] phenomenon, the crucial difference is that in the mismatch model, the peaks always have the same height and thus, there is no trade-off between peak height and width. Another way of altering the strength of selection is by changing e, the energetic cost of a mismatch. Larger e corresponds to stronger selection and the phase diagram changes accordingly (S11 Fig).

The empirical genotype-phenotype landscapes are topologically more complex [79, 83] than the genotype-phenotype landscapes generated with the mismatch model, which are regular graphs (i.e., every genotype has (a − 1) ⋅ L mutational neighbors). We therefore used simulations to calculate the average length of a greedy adaptive walk 〈l〉 in these empirical landscapes [84], initiating the walks from all non-peak genotypes in the fitness landscape. Moreover, each binding affinity measurement (E-score) in the empirical landscapes is associated with a noise threshold that is used to determine whether two genotypes truly differ in phenotype (Methods). This noise threshold can cause landscapes to have large non-peak plateaus, in which many genotypes have indistinguishable fitness. We therefore modified the greedy adaptive walk such that when a non-peak plateau was encountered, we chose a random mutational neighbor of indistinguishable fitness for the next step in the walk, disallowing reversion mutations. We repeated this process until the plateau was traversed and a sequence with higher fitness was reached. Finally, we terminated the walk when a peak sequence was reached. Therefore, the walk was primarily a deterministic greedy walk, with some stochasticity due to the non-peak plateaus.

Fig 7B shows 〈l〉, averaged over 100 simulations of the adaptive walk from each initial condition, in relation to wopt. In contrast to the mismatch model, 〈l〉 is shortest for wopt = 0.35 and wopt = 0.5 and slightly higher for intermediate wopt. This is because, as in the mismatch model, 〈l〉 is correlated with the total number of genotypes in peaks, which depends upon both the number of peaks and their widths. Whereas selecting for wopt = 0.35 leads to the largest number of peaks (Fig 6D), selecting for wopt = 0.5 leads to the broadest peaks (S6B Fig), thus explaining the minimisation of 〈l〉 when selecting for these extreme phenotypes. We note that in the absence of plateaus (i.e., when the noise threshold is zero), 〈l〉 increases monotonically with wopt, because in this case, 〈l〉 is inversely correlated with the number of peaks in the fitness landscape, which decreases monotonically with wopt. Regardless of the noise threshold, when selecting for low or intermediate phenotypic values, most walks terminate at a local, rather than the global fitness peak. However, these local peaks tend to be nearly as high as the global peak, especially when selecting for intermediate phenotypic values (S9 Fig). These results hold for all noise thresholds, and are therefore not a consequence of the existence of plateaus in the genotype-phenotype landscapes (see S1 Appendix, Note 2 for explanation).

Next, we simulated deterministic mutation selection dynamics (Methods). For each wopt, Fig 7D shows the change in mean fitness at equilibrium for μ = 0.1, σ = 0.15, relative to the fitness landscape generated after selecting for wopt = 0.5. As in the mismatch model, the change in mean fitness tends to be positive, despite the increase in the number of peaks caused by selection for low or intermediate phenotypic values (Fig 6D). Moreover, the mean change in fitness at equilibrium is maximized when selection favors an intermediate phenotypic value (wopt = 0.425 for μ = 0.1, σ = 0.15). This can be explained by the average peak heights and widths that occur when selecting for intermediate phenotypic values (S6A and S6B Fig for σ = 0.15). Exactly which wopt maximizes the change in mean fitness depends on the interplay between the mutation rate μ and the strength of selection σ (S12 Fig), converging on 0.425 for large μ and σ, similar to the convergence seen at mopt = 6 for the mismatch model.

Discussion

Non-linear relationships between phenotype and fitness have been observed in a diversity of biological systems [28, 31, 33, 35, 85, 86], often reflecting a trade-off between the costs and benefits of a particular phenotype, such as antibiotic resistance [87, 88]. It is well established that these non-linearities, either in the genotype-phenotype or phenotype-fitness map are a cause of epistasis [36–38, 74]. Moreover, when mutations have epistatic interactions in their contribution to phenotype, a non-linear phenotype-fitness map can change the form of these interactions from negative to positive, or vice versa [16], as well as introduce or remove sign epistasis [85]. Our work complements these empirical observations and expands upon previous theoretical work [74, 89], by systematically quantifying how and how often selection for a low or intermediate phenotypic value introduces or removes epistasis in the fitness landscape. Specifically, we show that the probability of changing the type of magnitude epistasis (e.g., positive to negative) is highest when selecting for low phenotypic values and the probability of introducing or removing sign epistasis is highest when selecting for intermediate phenotypic values. Further, we show that the simple sign epistasis motif cannot be converted into reciprocal sign epistasis, implying that genotype-phenotype landscapes with only simple sign epistasis motifs will remain single peaked and globally congruent to their corresponding fitness landscapes.

Another key finding of our analysis is that selection for low or intermediate phenotypic values is more likely to increase than to decrease the number of peaks, with the probability of the two types of change being equal only in House-of-Cards genotype-phenotype landscapes. This means that additive genotype-phenotype landscapes will tend to be incongruent with their fitness landscapes, whereas rugged genotype-phenotype landscapes will not. While increased landscape ruggedness is typically thought to frustrate the evolutionary process, because it limits the amount of adaptive phenotypic variation mutation can bring forth [4], our evolutionary simulations show this need not be the case. Specifically, we find that the rugged fitness landscapes caused by selection for low or intermediate phenotypic values comprise local adaptive peaks that are nearly as tall as the global adaptive peak. Moreover, these local peaks tend to be accessible from throughout the landscape via a small number of sequential mutations that monotonically increase fitness. As a result, the mean population fitness at equilibrium is almost always higher when selecting for low or intermediate phenotypic values than when selecting for a high phenotypic value.

Finally, while there have been several attempts at investigating genotype-phenotype-fitness landscapes in the past [74, 90, 91]—some models have been very specific to the system of interest and others are agnostic to any mechanistic details [92]. We tried to bridge this gap, by applying a Fisher’s Geometric model-like phenotype-fitness function to biophysically motivated and empirically determined genotype-phenotype landscapes. Further, our results on the mismatch model and the 1,137 landscapes of TF-DNA interactions may help to explain the prevalence of low- and intermediate-affinity binding sites in the control of gene expression. Prior work has suggested an entropic argument [93]: As with certain RNA secondary structures [94] or regulatory circuit motifs [95], low- and intermediate-affinity binding sites appear more frequently simply because they are more findable. That is, because a transcription factor binds more distinct DNA sequences with low or intermediate affinity than with high affinity, low- and intermediate-affinity binding sites are more likely to evolve to control gene expression. Our work complements this arrival of the frequent argument [96] by showing that low- and intermediate-affinity binding sites are not only more likely to arise de novo due to their increased frequency, selection for such sites also generates fitness landscapes that are more conducive to adaptation—in terms of increased mean fitness at equilibrium and decreased average length of adaptive walks, than fitness landscapes that were generated by selection for high affinity binding sites.

We made several modeling assumptions to simplify our analyses of transcription-factor DNA interactions, the relaxation of which may open new avenues for future research. First, we assumed a single nonlinearity in the relationship between genotype, phenotype, and fitness. As noted by Domingo et al. [37], from transcription to RNA processing, translation, and protein folding and all the way up to protein activity and cellular fitness, there are many layers of biological organization where the effects of a mutation can be transformed. To date, our knowledge of how such nonlinearities combine to modify genotype-phenotype landscapes is based on a small number of experimental studies (e.g., ref [16]). A systematic theoretical analysis could build off the work presented here, for example by incorporating the sigmoidal relationship between binding site occupancy and binding affinity in modeling transcription factor-DNA interactions [97], such that fitness depends nonlinearly on occupancy, rather than affinity. Other nonlinearities, such as those caused by transcription factor cooperativity [97], could also be included.

Second, we assumed that selection acts directly on a single phenotype—binding affinity. While this assumption is common in models of the evolution of transcription factor binding sites [21, 43, 46] and is supported by empirical data [43, 44], the relationship between binding affinity and fitness is not so direct, because it is modulated by gene expression. Gene expression depends on a variety of factors, including the presence, arrangement, and affinities of binding sites for other competing or cooperating transcription factors in promoters and enhancers [98], as well as local sequence context [99], chromatin context [100], DNA methylation [101], and local transcription factor concentrations [102]. Existing modeling frameworks that relate the architecture of entire regulatory regions to gene expression patterns may provide a path forward [103], facilitating the study of landscape incongruence when fitness depends upon the multitude of molecular phenotypes characteristic of eukaryotic gene regulation. Alternatively, our modeling framework could be extended to include multiple phenotypes by defining fitness in terms of the differences between a vector of phenotypes and a vector of optimal phenotypes, rather than the scalars considered here. Incongruence could then be quantified between the fitness landscape generated by selecting for the highest phenotypic value of each phenotype and that generated by selecting for a combination of low and intermediate values of the phenotypes. Our results correspond to a special case of this extended model, wherein all phenotypes except one are exactly attuned to their optimal values.

Third, we assumed a Gaussian phenotype-fitness map, which is commonly employed in a diversity of modeling frameworks [29, 70, 72, 104], including those for transcription factor-DNA interactions [46]. Alternative symmetric phenotype-fitness maps (e.g., ref. [105]) will only affect our results quantitatively, because many of our findings, such as the changes in sign epistasis motifs and number of peaks, only depend on the rank ordering of fitness values. However, we expect asymmetric phenotype-fitness maps, such as those uncovered in experimental studies of biological systems such as viruses [30] and yeast [106], to affect our results qualitatively. For example, in our analyses of NK landscapes, we often observed symmetries in incongruence around an intermediate phenotypic value. These symmetries will almost certainly be lost. Understanding how asymmetric phenotype-fitness maps affect the incongruence of genotype-phenotype landscape is therefore an outstanding challenge.

Finally, our study may open new lines of research on dynamic genotype-phenotype and fitness landscapes [107-113]. For example, whereas we studied selection for a fixed phenotypic optimum, this optimum may in fact change in space or in time. Our results imply that even gradual changes in the phenotypic optimum may lead to abrupt changes in fitness landscape topography, which may have implications for an evolving population’s ability to track this optimum and thus for population persistence and extinction. Moreover, because our measures of incongruence can be applied to any pair of landscapes so long as they are defined over the same set of genotypes, they are also applicable whenever a phenotype is mapped non-linearly to another phenotype. Ideally, we would be able to make inferences about phenotypic architecture based on the topographical properties of higher-level phenotypic or fitness landscapes—as was previously done for an antibiotic resistance phenotype [114]. However, because the phenotype-fitness map we study is not invertible, we can only make such inferences in limited cases. For instance, when the fitness landscape is single-peaked, we can make the probabilistic inference that the underlying genotype-phenotype landscape is also likely to be single-peaked, because the phenotype-fitness map is more likely to increase than to decrease the number of peaks. In contrast, when the fitness landscape has multiple peaks, we can only infer that the underlying genotype-phenotype landscape does not solely comprise simple sign epistasis motifs. Beyond that, we cannot infer the topographical properties of the underlying genotype-phenotype landscape. It could be smooth or rugged. Additionally, our measures may shed light on the kinds of topographical alterations induced by fluctuating environmental factors, such as DNA methylation [101] or the presence of protein partners [115] on transcription factor-DNA interactions. As our ability to experimentally interrogate such complexities in the relationship between genotype, phenotype, and fitness continues to improve, we anticipate a sharpened focus on landscape dynamics and their implications for the evolutionary process.

Methods

Epistasis

In a genotype-phenotype landscape, we classify the type of magnitude epistasis between a pair of loci using the following linear combination of phenotypic values: where w represents the phenotype of genotype i ∈ {0, 1}2. When ϵgp = 0, there is no magnitude epistasis, because the phenotypic effects of the two alleles combine additively; when ϵgp > 0, there is positive epistasis, because the phenotypic effects of the two alleles are greater than expected based on their individual phenotypic effects; when ϵgp < 0, there is negative epistasis, because the phenotypic effects of the two alleles are less than expected based on their individual phenotypic effects.

Analogously, in the fitness landscape, we classify the type of magnitude epistasis between a pair of loci using the following linear combination of fitness values: where F represents the fitness of the corresponding genotype i. As in the genotype-phenotype landscape, ϵf = 0, ϵf > 0, and ϵf < 0 indicate additive, positive, and negative epistatic interactions among loci, respectively.

We use ϵgp and ϵf to calculate the fraction of pairs of loci that have the same type of epistasis in the genotype-phenotype landscape and the fitness landscape, which we determine as the product of ϵgp and ϵf. This is because the type of epistasis is the same in the two landscapes when ϵgp ⋅ ϵf > 0. While it is theoretically possible for both ϵgp and ϵf to be 0, in which case the type of epistasis would be the same in the two landscapes yet the condition ϵgp ⋅ ϵf > 0 would not be satisfied, this never happens in practice.

NK Landscapes

We constructed the NK landscapes using the adjacency neighbourhood scheme, wherein each locus i of a genotype of length L, interacts with K adjacent loci to the right of locus i and 0 ≤ K ≤ L − 1. We used periodic boundary conditions, such that the L-th locus interacts with the first K loci, and so on.

The phenotype w(τ) of genotype τ is computed as the sum of the individual contributions of all loci, each of which depends on K other interacting loci: where represents the contribution of the i-th locus, which depends on K other loci . We drew the contributions of each of the 2 possible configurations from a uniform distribution between 0 and 1. Finally, we re-scaled the phenotypic values by subtracting the minimum and dividing by the maximum, such that they fell between 0 and 1.

Empirical genotype-phenotype landscapes of transcription factor-DNA interactions

We studied the incongruence of 1,137 genotype-phenotype landscapes of transcription factor-DNA interactions. The procedure for constructing these landscapes has been described elsewhere [21]. In brief, each landscape corresponds to a single transcription factor, the genotypes it contains represent DNA sequences of length L = 8 that specifically bind the transcription factor, and the phenotype of each sequence is a quantitative proxy for relative binding affinity, which defines the surface of the landscape. These phenotypes are reported as enrichment scores (E-scores) derived from protein binding microarrays. In each landscape, two genotypes are considered mutational neighbors if they differ by a single small mutation, specifically a point mutation or a small indel that shifts an entire contiguous binding site by a single base [79]. We performed all analyses on the dominant genotype network (i.e., the largest connected component of the network), which comprises the vast majority of genotypes in each landscape [21].

We used the Genonets Python package (version 0.31) to characterize the topographical properties of the empirical genotype-phenotype landscapes [116]. Specifically, we used this package to compute the number of peaks per landscape and the number of genotypes in the peaks of each landscape. These calculations rely on a noise threshold δ, which is used to determine whether two genotypes actually differ in phenotype. For each transcription factor, we used the value of δ reported in ref. [21], which was derived from a comparison of binding affinity measurements across two protein binding microarray designs.

While characterizing the topographies of the resulting fitness landscapes, we also had to transform δ following the rules of error propagation. Accordingly, the noise in fitness values, dF, depends on the noise in the phenotypic values dw = δ as follows:

To compute the number of peaks and their widths in the fitness landscapes, we adapted Genonets to specify different noise threshold values (dF) for each genotype.

Mutation-selection dynamics for the mismatch model

We grouped genotypes according to their mismatch class m and iterated a series of selection and mutation steps until the population reached an equilibrium distribution. In each selection step, the frequency of each mismatch class X was scaled by its fitness F and then normalized: where is the frequency of the mismatch class m after the selection step.

The selection step was followed by a mutation step, in which genotypes could either mutate within their mismatch class or mutate to an adjacent mismatch class. The total mutation probability is μ, so with probability 1 − μ, the genotype does not mutate. The frequency of each mismatch class was thus updated as where is the frequency of the mismatch class m after the mutation step and μ is the mutation probability. The mutation step implicitly accounts for the different number of genotypes per mismatch class. Finally, X(t + 1) was set to and t was incremented. This process was repeated until equilibrium.

Mutation-selection dynamics for the empirical landscapes

The empirical landscapes comprise far fewer genotypes than the mismatch landscape. We therefore defined the state space of the empirical landscapes in terms of the individual genotypes in each landscape, merging each genotype with its reverse complement [79]. The mutational neighbors of each genotype were those that differed by a single small mutation, namely a point mutation or an indel that shifted an entire contiguous binding site by a single base [79]. In this analysis, we did not account for uncertainties in the fitness values, because these do not influence our results in this evolutionary regime (see below).

The recursion relation for mutation-selection dynamics in discrete time is where is the vector of genotype frequencies, x is the frequency of the ith genotype, μ is the probability of a single point mutation or a small indel mutation, L is the length of the genotypes, d(τ, τ) is the minimum of the mutational distance between genotypes τ and τ and between genotypes τ and the reverse compliment of τ, f is the fitness of the ith genotype, and is the mean population fitness at time t.

We linearized the dynamics by substituting , yielding from which we retrieved the normalized genotype frequencies with . In matrix form, the dynamics are where the mutation matrix M has elements and the selection matrix S is a diagonal matrix with S = f, where f is the fitness of sequence τ. We calculated the eigenvector corresponding to the largest eigenvalue of M ⋅ S to determine the equilibrium state of the dynamics, which is guaranteed by the Frobenius-Perron theorem to be unique and stable.

Because the fitness values have some uncertainty around them, we additionally performed a sensitivity analysis by adding Gaussian noise to the fitness values, with standard deviation equal to 1/3 ⋅ dF, such that 99.7% of the sampled fitness values fell within F±δ. Our results are robust to these perturbations, with the exception of two parameter combinations—σ = 0.1, μ = 0.001 and σ = 0.1, μ = 0.0025. For small σ, dF is large and therefore such sensitivity is expected. However, the quantitative changes in our results were small, leading us to conclude that fitness uncertainties do not significantly influence our results in this evolutionary regime.

Derivations, proofs and notes.

(PDF)

Click here for additional data file.

Simple sign epistasis motif cannot be changed into a reciprocal sign epistasis motif.

A simple sign epistasis motif (shown in brown) cannot be changed into a reciprocal sign epistasis motif, for any wopt. The no sign epistasis motif is shown in blue. The remaining colours represent the neighbourhood of the phenotypic values (w) corresponding to each genotype (g), where i ∈ {1, 2, 3, 4} and the genotypes are labelled in descending order, such that the genotype with the highest phenotypic value is labelled 1, and so on. The four bottom arrows point to the transformed motifs when wopt belongs to the neighbourhood from which the arrow emanates.

(PDF)

Click here for additional data file.

Reciprocal sign epistasis motif can be changed into the no sign epistasis motif and the simple sign epistasis motif with equal probability.

Selection for wopt transforms the reciprocal sign epistasis motif (shown in red) into the no sign epistasis motif (shown in blue) and the simple sign epistasis motif (shown in brown) with equal probability. The neighbourhood of each phenotypic value is shown in a different colour. For the no sign epistasis motif to emerge (top), the fitness values need to be “separated”, while for the simple sign epistasis motif to emerge (bottom), the fitness values need to be “interspersed”.

(PDF)

Click here for additional data file.

Evolutionary consequences of landscape ruggedness in NK landscapes.

(A,C) The average length of a greedy adaptive walk and (B,D) the change in mean fitness at equilibrium, relative to fitness at equilibrium when selecting for wopt = 1, are shown for (A,B) L = 5 and (C,D) L = 8. In (B,D), μ = 0.1 and σ = 0.5.

(PDF)

Click here for additional data file.

Selection for wopt = 0 can change the number and location of peaks.

Selection for wopt = 0 does not change the type of epistasis motif for any “square” in the fitness landscape, relative to the genotype-phenotype landscape, yet it can change the number and location of peaks. To understand how, consider that any two adjacent faces of the hypercube (e.g., gray faces above) are sufficient to determine whether the genotypes on their common edge are peaks or not. After selecting for wopt = 0, the type of epistasis motif does not change in the adjacent faces, yet the number and location of the peaks on their common edge does change. Peak genotypes are shown in red. Arrows point from lower to higher phenotypic or fitness values.

(PDF)

Click here for additional data file.

Sensitivity analysis for empirical landscapes.

Mean absolute change in the number of peaks in the fitness landscape relative to the genotype-phenotype landscape, shown in relation to wopt, for (A) single-peaked vs. multi-peaked genotype-phenotype landscapes and (B) three values of σ.

(PDF)

Click here for additional data file.

Average peak height and width for empirical landscapes.

Average peak (A) height and (B) width for 1,137 empirical landscapes, shown in relation to the optimal binding affinity wopt. Violin plots show the distribution across the landscapes for each wopt. Data include both local and global peaks.

(PDF)

Click here for additional data file.

Greedy adaptive walks on NK landscapes.

The average height of peaks reached by greedy adaptive walks in NK landscapes with σ = 0.5, shown in relation to wopt, for (A) L = 5 and (B) L = 8, with K = 0…L − 1.

(PDF)

Click here for additional data file.

Change in mean fitness at equilibrium for NK landscapes.

Change in mean fitness at equilibrium for NK landscapes with (A,B) L = 5 and (C,D) L = 8, shown in relation to wopt, for different values of σ and μ.

(PDF)

Click here for additional data file.

Greedy adaptive walks on empirical landscapes.

(A) Local peaks on which the adaptive walks terminate tend to be nearly as tall as the global peaks in the 1,137 empirical landscapes. Violin plots show the distribution of the fractional height of local peaks reached by greedy adaptive walks, relative to the height of the global peak, for each optimal binding affinity wopt. (B) Violin plots show the distribution of the fraction of walks terminating on the global peak, for each optimal binding affinity wopt.

(PDF)

Click here for additional data file.

Change in mean fitness at equilibrium for mismatch model landscapes.

The mismatch class mopt that maximizes mean fitness at equilibrium is shown in relation to the strength of selection σ and the mutation rate μ. The three smaller panels show mean fitness at equilibrium in relation to mopt for three combinations of σ and μ. The value of mopt that maximizes mean fitness is indicated with a gray rectangle.

(PDF)

Click here for additional data file.

Sensitivity analysis for the change in mean fitness at equilibrium for mismatch model landscapes.

The mismatch class mopt that maximizes mean fitness at equilibrium is shown in relation to the strength of selection σ and the mutation rate μ for three different values of the mismatch penalty e: (A) e = 0.025, (B) e = 0.05 and (C) e = 0.1.

(PDF)

Click here for additional data file.

Change in mean fitness at equilibrium for empirical landscapes.

The binding affinity wopt that maximizes mean fitness at equilibrium is shown in relation to the strength of selection σ and the logarithm of the mutation rate (log μ) for the 1,137 empirical landscapes. The three smaller panels show the distributions of mean fitness at equilibrium as violin plots, in relation to wopt, for three combinations of σ and μ. Box-and-whisker plots show the 25–75% quartiles. The value of wopt that maximizes mean fitness is indicated. These results are robust to perturbations of the fitness values (Methods), with the exception of two parameter combinations—σ = 0.1, μ = 0.001 and σ = 0.1, μ = 0.0025. For these parameter combinations, the binding affinity wopt that maximizes mean fitness at equilibrium changes from 0.35 to 0.385 (on average) and 0.375 respectively.

(PDF)

Click here for additional data file.

20 Jun 2022

Dear Prof. Payne,

Thank you very much for submitting your manuscript "On the incongruence of genotype-phenotype and fitness landscapes" for consideration at PLOS Computational Biology.

As with all papers reviewed by the journal, your manuscript was reviewed by members of the editorial board and by several independent reviewers. In light of the reviews (below this email), we would like to invite the resubmission of a significantly-revised version that takes into account the reviewers' comments.

While the reviewers found your work interesting, they made several suggestions for clarifications and raised additional questions. These points appear valid and should be addressed carefully, which will require extensive revisions. If you decide to submit a revised version, please respond to all reviewer comments, especially those of reviewer #2, and take them into account when preparing the revised manuscript.

We cannot make any decision about publication until we have seen the revised manuscript and your response to the reviewers' comments. Your revised manuscript is also likely to be sent to reviewers for further evaluation.

When you are ready to resubmit, please upload the following:

[1] A letter containing a detailed list of your responses to the review comments and a description of the changes you have made in the manuscript. Please note while forming your response, if your article is accepted, you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out.

[2] Two versions of the revised manuscript: one with either highlights or tracked changes denoting where the text has been changed; the other a clean version (uploaded as the manuscript file).

Important additional instructions are given below your reviewer comments.

Please prepare and submit your revised manuscript within 60 days. If you anticipate any delay, please let us know the expected resubmission date by replying to this email. Please note that revised manuscripts received after the 60-day due date may require evaluation and peer review similar to newly submitted manuscripts.

Thank you again for your submission. We hope that our editorial process has been constructive so far, and we welcome your feedback at any time. Please don't hesitate to contact us if you have any questions or comments.

Sincerely,

Tobias Bollenbach

Associate Editor

PLOS Computational Biology

Ville Mustonen

Deputy Editor

PLOS Computational Biology

***********************

While the reviewers found your work interesting, they made several suggestions for clarifications and raised additional questions. These points appear valid and should be addressed carefully, which will require extensive revisions. If you decide to submit a revised version, please respond to all reviewer comments, especially those of reviewer #2, and take them into account when preparing the revised manuscript.

Reviewer's Responses to Questions

Comments to the Authors:

Please note here if the review is uploaded as an attachment.

Reviewer #1: The manuscript “On the incongruence of genotype-phenotype and fitness landscapes” investigates how the ruggedness of fitness landscape is related to the underlying genotype-phenotype landscape when different selective pressures are applied. By using multisite biallelic genotype-phenotype landscapes and Gaussian phenotype-to-fitness map centered around phenotypic optimum, authors derive analytical relations for local and global properties of resulting fitness landscapes. The results are verified computationally for genotype-phenotype landscapes in the NK model and in the experimentally obtained dataset containing binding affinity of TF to DNA. The important finding is that selection for low and intermediate phenotypic values significantly increases ruggedness of the landscape, but interestingly the evolutionary accessibility of global peak is not affected.

The question of how different selective pressures modify both level of epistasis and number of adaptive peaks when genotype-phenotype landscape is mapped to the fitness landscape is a very interesting one with potential implications in both theoretical and experimental studies of biological evolution. The manuscript provides a very valuable contribution in terms of identifying different factors affecting local and global measures of ruggedness when mapping genotype-phenotype landscape into fitness landscape. I highly appreciate supporting computational results with analytical calculations that could guide our intuition about system behavior. The paper is well written, technically sound and results are well presented. My only major concern, is how the multi-allelic landscape is mapped to bi-allelic landscape, and whether this mapping is not biasing the evolutionary accessibility of resulting fitness landscape. Still, I expect this concern could be addressed. In summary, I find that this work deserves a broad dissemination and will be of interest for the PLoS Computational Biology audience.

Major issue

1. The multi-allelic genotype-phenotype landscape that mimics energy landscape of TF-DNA interactions is transformed into bi-allelic landscape through the threshold model. The threshold model results in genotype-phenotype landscape with large plateaus corresponding to low and intermediate phenotypic values. If now one phenotypic value corresponding to one of these plateaus is selected as most optimal it is not surprising that although the landscape is very rugged the optimal peaks can be reached with very short adaptive walks.

However, in more realistic approach the energy landscape of TF-DNA would likely not have these plateaus, but rather a spread of phenotypic values (see Mustonen et al. 2008, cited as [43]), with empirical distribution of these energies resembling rather a continuous bulk of non-specific bindings and additional peak for specific bindings. My objection with the current mismatch model is that applying a fitness map to the idealized plateau-like landscape, could be qualitatively different to applying the same fitness map to a landscape in which plateaus are not present. In practice one could introduce noise to phenotypic values in plateau-like landscape, or consider more fine-grained distinction of binding energies to introduce larger number of discrete phenotypic levels. The question would be whether fitness map will still result in fitness landscape with short adaptive walks toward the global fitness optimum?

I would appreciate at least a discussion of such case, as my impression is that current conclusion about no change in accessibility of adaptive peaks is simply a consequence of the plateau-like mapping introduced by threshold model.

Minor issues

2. lines 167-168, please give explicitly the cut-off for zero fitness, instead of “within computer precision”.

3. Eq 5 and lines 436-439, the properties of accessible pathways in multi-allelic landscapes were studied in Zagorski et al., PLoS Comp Bio 2016, indicating that increasing number of alleles would increase the average length of mutational pathways towards global optimum. Does it hold for eq. 5?

4. figure 7 caption, there are no “Closed symbols and dashed line” in figures 7B and 7D. Please remove or correct figures.

Reviewer #2: The topic of the paper are the changes in landscape properties that occur when a genotype-phenotype landscape is combined with a nonlinear phenotype-fitness map. Specifically, the authors consider a one-dimensional real-valued phenotype that is mapped to fitness via a nonlinear function with a unique optimal phenotype value. The primary motivation for the work is the perception that "genotype-phenotype landscapes are often used as a proxy for fitness landscapes", and that this may be potentially misleading. While I agree that this problem exists to some extent (mainly because organismal fitness is oftentimes a complex, multidimensional construct), I feel that it is also something of a straw man, since most workers in the field are very well aware of the distinction between the two concepts.

For the same reason, I'm not really happy with the title of the paper (though I know that it has already been modified

compared to the bioRxiv version of the manuscript). The authors address an important problem, viz., the transformation

properties of multidimensional phenotypic epistasis under nonlinear phenotype-fitness maps, and I think this could be better

reflected in the title. For the specific case where phenotypic epistasis is absent (and therefore the nonlinear phenotype-fitness

map is the sole source of epistasis for fitness) this problem has been addressed extensively (in the context of Fisher's geometric

model and elsewhere), and the authors cite a number of studies in this direction, but the generalization to epistatic

genotype-phenotype landscapes is novel and timely.

I am generally in favor of publication of the manuscript in PLoS Computationial Biology, but ask the authors to address the

following specific points prior to acceptance.

1. Related to the previous remark, I think the authors should emphasize that their analysis is not restricted to nonlinear

phenotype-fitness maps, but applies whenever one phenotype is mapped nonlinearly to another. In such a situation, the analysis of the resulting landscape can provide information about the underlying phenotypic architecture. As an example, Zwart et al. [Heredity 2018] concluded from an analysis of a genotype-resistance landscape for synonymous mutations in an antibiotic resistance enzyme that the underlying phenotype cannot be one-dimensional.

2. One of the most remarkable results of the paper is that simple sign epistasis motifs (SSE) cannot be transformed into

reciprocal sign epistasis (RSE) under a unimodal phenotype-fitness map. Unless I am mistaken, this has an even more remarkable (and quite counterintuitive) corollary. A well-known result by Poelwijk et al. (JTB 2011) states that a multipeaked fitness landscape has to contain at least one instance of RSE. Together with the new result, this would then imply that a genotype-phenotype landscape where all 2-faces have SSE remains single-peaked under arbitrary unimodal phenotype-fitness transformations [it is straightforward to see that such genotype-phenotype landscapes can be constructed for any L]. If the authors agree with this conclusion, I think it is worth mentioning.

3. Conceptually, it is not clear to me why the authors use the absolute change in the number of peaks as a measure of "incongruence", rather than just the different p_f - p_gp. Isn't is important to know whether the number of peaks increases or decreases?

4. Figure 4A and lines 225-227: I cannot see any maxima in the figure for L \\geq 4, and in fact I think the argument in the SM

(Note 1) is flawed. When the optimal phenotype is between two integer multiples of 1/L, there will be peaks in both layers of the cube and there is no reason why the total number of peaks should decline.

5. As it stands, Eq.(3) makes no sense, because the right hand side diverges with L. This should be rewritten to make it into a

proper limit statement.

6. Figure 7B: Where does the enormous variability of the adaptive walk length come from?

7. Supplementary figure 3: There are analytic results for the mean length of greedy adaptive walks in the HoC and NK-models

(e.g., Orr, JTB 2003; Nowak and Krug, J. Stat. Mech. 2015). In particular, the greedy walk length for the HoC model is

e-1 \\approx 1.7 for large L, which is smaller than the results shown in Fig. S3A and C. Do the authors possibly use a nonstandard definition of walk length where a 'zeroth' step is included?

**********

Have the authors made all data and (if applicable) computational code underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data and code underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data and code should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data or code —e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: Yes

Reviewer #2: Yes

Figure Files:

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, . PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email us at .

Data Requirements:

Please note that, as a condition of publication, PLOS' data policy requires that you make available all data used to draw the conclusions outlined in your manuscript. Data must be deposited in an appropriate repository, included within the body of the manuscript, or uploaded as supporting information. This includes all numerical values that were used to generate graphs, histograms etc.. For an example in PLOS Biology see here: http://www.plosbiology.org/article/info%3Adoi%2F10.1371%2Fjournal.pbio.1001908#s5.

Reproducibility:

To enhance the reproducibility of your results, we recommend that you deposit your laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. Additionally, PLOS ONE offers an option to publish peer-reviewed clinical study protocols. Read more information on sharing protocols at https://plos.org/protocols?utm_medium=editorial-email&utm_source=authorletters&utm_campaign=protocols

11 Aug 2022

Submitted filename: ResponseLetter_Final.pdf

Click here for additional data file.

30 Aug 2022

Dear Prof. Payne,

We are pleased to inform you that your manuscript 'On the incongruence of genotype-phenotype and fitness landscapes' has been provisionally accepted for publication in PLOS Computational Biology.

Before your manuscript can be formally accepted you will need to complete some formatting changes, which you will receive in a follow up email. A member of our team will be in touch with a set of requests.

Please note that your manuscript will not be scheduled for publication until you have made the required changes, so a swift response is appreciated.

IMPORTANT: The editorial review process is now complete. PLOS will only permit corrections to spelling, formatting or significant scientific errors from this point onwards. Requests for major changes, or any which affect the scientific understanding of your work, will cause delays to the publication date of your manuscript.

Should you, your institution's press office or the journal office choose to press release your paper, you will automatically be opted out of early publication. We ask that you notify us now if you or your institution is planning to press release the article. All press must be co-ordinated with PLOS.

Thank you again for supporting Open Access publishing; we are looking forward to publishing your work in PLOS Computational Biology.

Best regards,

Tobias Bollenbach

Academic Editor

PLOS Computational Biology

Ville Mustonen

Section Editor

PLOS Computational Biology

***********************************************************

Reviewer's Responses to Questions

Comments to the Authors:

Please note here if the review is uploaded as an attachment.

Reviewer #1: In the revised manuscript authors have addressed all my major and minor concerns. In particular, authors performed simulations of landscapes with fine-grained distinction of phenotypic values, to validate whether accessibility of such landscapes is affected in a different way than in the landscapes with plateau regions. As a result, authors find that there are quantitative differences in the difference of peaks, the length of adaptive walks, and fraction of walks that reach global peak, but the overall trends are qualitatively the same for fine-grained and plateau-like landscapes. This was illustrated in the reviewer figures and addressed in the text of the revised manuscript.

I find the revised manuscript to provide a very valuable contribution to systems and computational biology communities and strongly support its publication in PLoS Computational Biology.

Reviewer #2: The authors have fully addressed the comments in my previous report. I recommend publication of the manuscript in its present form.

**********

Have the authors made all data and (if applicable) computational code underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data and code underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data and code should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data or code —e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: Yes

Reviewer #2: Yes

**********

PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: Yes: Joachim Krug

13 Sep 2022

PCOMPBIOL-D-22-00693R1

On the incongruence of genotype-phenotype and fitness landscapes

Dear Dr Payne,

I am pleased to inform you that your manuscript has been formally accepted for publication in PLOS Computational Biology. Your manuscript is now with our production department and you will be notified of the publication date in due course.

The corresponding author will soon be receiving a typeset proof for review, to ensure errors have not been introduced during production. Please review the PDF proof of your manuscript carefully, as this is the last chance to correct any errors. Please note that major changes, or those which affect the scientific understanding of the work, will likely cause delays to the publication date of your manuscript.

Soon after your final files are uploaded, unless you have opted out, the early version of your manuscript will be published online. The date of the early version will be your article's publication date. The final article will be published to the same URL, and all versions of the paper will be accessible to readers.

Thank you again for supporting PLOS Computational Biology and open-access publishing. We are looking forward to publishing your work!

With kind regards,

Zsofi Zombor

PLOS Computational Biology | Carlyle House, Carlyle Road, Cambridge CB4 3DN | United Kingdom ploscompbiol@plos.org | Phone +44 (0) 1223-442824 | ploscompbiol.org | @PLOSCompBiol