Edited by Richard E. Lenski, Michigan State University, East Lansing, MI, and approved October 6, 2015 (received for review May 26, 2015)
Suboptimal fitness peaks are generally recognized as causing evolutionary stasis. Here, we show that these constraints can be overcome in an adaptive manner by reconstructing mutational trajectories for a transcription factor and its DNA binding site in variable environments. Cross-environmental tradeoffs, typically associated with evolutionary limitations, are an essential enabling component of this evolutionary mechanism. Our results underscore the importance of characterizing environmental dependencies when studying genetic interactions and provide the clearest indication so far that environmental variability can accelerate evolution hampered by stasis in constant conditions. Given that environmental variations and tradeoffs are ubiquitous, this evolutionary mechanism may be relevant to a wide range of genetically constrained phenotypes and major evolutionary transitions.
Epistatic interactions can frustrate and shape evolutionary change. Indeed, phenotypes may fail to evolve when essential mutations are only accessible through positive selection if they are fixed simultaneously. How environmental variability affects such constraints is poorly understood. Here, we studied genetic constraints in fixed and fluctuating environments using the Escherichia coli lac operon as a model system for genotype–environment interactions. We found that, in different fixed environments, all trajectories that were reconstructed by applying point mutations within the transcription factor–operator interface became trapped at suboptima, where no additional improvements were possible. Paradoxically, repeated switching between these same environments allows unconstrained adaptation by continuous improvements. This evolutionary mode is explained by pervasive cross-environmental tradeoffs that reposition the peaks in such a way that trapped genotypes can repeatedly climb ascending slopes and hence, escape adaptive stasis. Using a Markov approach, we developed a mathematical framework to quantify the landscape-crossing rates and show that this ratchet-like adaptive mechanism is robust in a wide spectrum of fluctuating environments. Overall, this study shows that genetic constraints can be overcome by environmental change and that cross-environmental tradeoffs do not necessarily impede but also, can facilitate adaptive evolution. Because tradeoffs and environmental variability are ubiquitous in nature, we speculate this evolutionary mode to be of general relevance.
It is widely believed that epistatic interactions can direct evolutionary change (1⇓⇓⇓⇓⇓–7). Epistasis has been implicated in shaping RNA (8) and protein (4, 6, 7, 9) sequences, sensing (5) and translation (10) functions, and developmental programs (11) and speciation (12⇓–14). Phenotypes may be difficult to evolve not because they are impossible biochemically or physically, but because essential mutations are mutually dependent and must be fixed together to be selected positively (5, 15⇓–17). How such genetic constraints can be overcome has been considered previously: population expansion or subdivision can limit negative selection and maintain less fit phenotypes (18, 19), and large populations and long waiting times can enable the joint fixation of multiple mutations (20), whereas recombination can join mutant alleles (21⇓–23). Other mechanisms include drift (24⇓–26), partial penetrance (27), and nonheritable lifetime plasticity (28, 29). However, how the constraining effects of such genetic interactions are affected by environmental variability remains poorly understood. It has been shown that mutational effects (30⇓⇓–33) and epistasis itself (34, 35) can depend on the environment, that bacterial resistance evolution can be contingent on the rate of antibiotic increase (36), and that adaptation in silico can be accelerated by environmental change (37⇓⇓–40). These observations suggest that the effects of environmental variability may go beyond merely producing variable selective pressures that favor certain phenotypes but also, could be involved in controlling phenotype accessibility and stasis.
To investigate how environmental variability affects genetic constraints, we focused on a model system for genotype–environment interactions, the lac regulatory system of Escherichia coli. Its physiology has been studied extensively: in the presence of lactose, expression of the lac genes allows E. coli cells to import and metabolize lactose, whereas in the absence of lactose, repression of these genes limits physiological costs (41, 42). The ability to regulate lac expression relies on the binding of the lac repressor to the lac operator DNA sequence upstream of the coding region (Fig. 1A). We surmised that the coevolution of such protein–DNA interfaces could be severely constrained by epistatic interactions, such that some genotypes are inaccessible by positive selection in single-mutation steps. In lock–key recognition, mutating either lock or key is expected to lead to recognition loss (2, 15). At the same time, mutating both lock and key may produce a different, better-matching pair. Indeed, the lac transcription factor phylogeny suggests extensive historic adaptation of the repressor–operator interface and reveals multiple homologous repressors that bind specifically to their cognate operator (43, 44). Furthermore, the lock–key recognition of the lac regulatory system is a highly specific function that is confined to a limited number of residues, which restricts the range of adaptive solutions. Indeed, extensive mutational analysis of the lac repressor–operator interface has shown that just two repressor residues and four operator bases control binding specificity (Fig. 1 A and B) (45, 46). Moreover, in contrast to phenotypes that are affected by many unknown mutations, one can identify the genetic interactions between mutations in controlling residues that are key to genetic constraint.
Repression and expression ability of lac repressor–operator mutants. (A) Schematic representation of the E. coli lac system. β-Galactosidase (LacZ) and the lac permease (LacY) are coregulated by the repressor LacI. Expression is induced by IPTG. Red lines correspond to mutated positions. (B) The multimeric lac repressor in green bound to its operator DNA. Red indicates mutated positions responsible for specific repressor–operator binding. Y and Q are the mutated amino acid residues on positions 17 and 18 in the DNA binding helix of the lac repressor, and 4 and 4′ g–c and 5 and 5′ t–a are the mutated base pairs in the operator DNA. We note that the genotype represented here is YQ:tggt. (C) Characterization of 64 lac repressor–operator variants. The starting and final sequences are indicated by white and black circles, respectively. R is the inverse of the measured expression level in the absence of IPTG. E is the measured expression level in the presence of IPTG (Methods).
We mutated sites in the lac operator DNA and the multimeric lac repressor that control binding specificity (Fig. 1 A and B) and measured the expression level of the downstream lac genes in two contrasting conditions (Methods). We quantified the ability to repress the lac genes, R, as the inverse of the measured lac expression level in the absence of inducing ligand (Methods and Fig. 1C). The ability to express the lac genes, E, was quantified by the measured lac expression level in the presence of ligand (Fig. 1C). Note that the repression ability (R) is, thus, not the inverse of the expression ability (E). We constructed four lac repressor–operator variants that have been predicted (45⇓⇓–48) to display binding: PK:agga, PS:acca, MK:acca, and YQ:tggt, where the first two letters indicate the controlling repressor residues and the last four letters indicate the controlling operator bases (Fig. 1 A and B). We focused on the latter two variants, because they displayed substantial fold changes between the induced and noninduced expression levels (R•E was 6 and 55, respectively), with E being approximately equal but R being about 20-fold lower for MK:acca. The MK:acca genotype is, thus, able to regulate lac expression but can improve repression ability by mutating the repressor (MK to YQ) and operator (acca to tggt).
We investigated the interaction between noncognate pairs by swapping around the two operators. The ability to repress (R) was found to be low (i.e., expression in the absence of ligand was high) for MK:tggt and YQ:acca (100- to 200-fold lower than for the cognate pair YQ:tggt). These data were consistent with the reciprocal sign epistasis hypothesized for lock–key interactions: changing either of the binding partners alone leads to binding loss, but changing the other partner as well restores it. This notion was supported by the overall expression levels for MK:tggt and YQ:acca, which were high and unresponsive to ligand (R•E = 1). However, although the presence of reciprocal sign epistasis is required, it is not sufficient to constrain phenotypes on suboptima (15). Indeed, the repressor and operator modifications both involve multiple mutations, and their one by one fixation in particular order (1) could confer continuous improvements in repression ability.
To test the accessibility of trajectories considering all possible orders of all essential mutations, we constructed the remaining intermediate genotypes between MK:acca and YQ:tggt. In total, 6! = 720 direct trajectories can be taken along the 26 = 64 genotypes. Whether a mutation is positively selected, and thus accessible, depends on the sign rather than the magnitude of the associated fitness change. Fitness and phenotypic changes have the same sign when the phenotype–fitness relation is monotonic. In the environment without ligand, for instance, a trajectory is then not accessible when it contains mutations that decrease R. Below, we first consider this case of monotonic phenotype fitness relations, and later, we relax this assumption to also consider a range of nonmonotonic relations.
Analysis showed that all trajectories contain depressions in both R and E (Fig. 2 and Fig. S1). The depressions are at least two mutations wide and peak at a width of five mutations, whereas the involved decrease is at least 3-fold and reaches up to ∼100-fold. Thus, none of the trajectories to YQ:tggt are accessible by fixing mutations one by one in either of two environments. Although this analysis concerns only direct trajectories (i.e., without mutational reversions), allowing for reversions did not open up accessible trajectories in either of the environments (SI Results, sections S6). Overall, these data indicate that higher-order genetic interactions (i.e., epistasis involving multiple mutations) limit optimization of the lac regulatory phenotype in each of two environments.
Genetic constraints in constant environments. (A) The valley width is the number of mutations required to increase repression or expression ability above the previous suboptimum in a given mutational trajectory. For each trajectory, the widest valley is tabulated. The horizontal bar indicates one mutation. (B) Valley width in repression ability for all direct 720 mutational trajectories from MK:acca to YQ:tggt. (C) Valley width in expression ability. (D) The valley depth is the fold decrease in repression or expression ability within the valleys of a given trajectory. For each trajectory, the deepest valley is tabulated. Note that the widest and the deepest valleys for a given trajectory do not necessarily coincide. (E) Valley depth in repression ability for all direct 720 mutational trajectories from MK:acca to YQ:tggt. (F) Valley depth in expression ability. Welch’s t tests were performed to determine the statistical significance of differences in measured repression and expression ability values (Methods).
Tradeoffs in repression and expression over longer sequence space. (A) Valley depth in both environments for all direct 720 trajectories. Valley depth is the largest decrease along a mutational trajectory. (B) Valley width in both environments for all direct 720 trajectories. Valley width equals the number of mutations required to achieve a level that is greater than or equal to the previous suboptimum. Points of corresponding trajectories are depicted with a small random offset for visualization purposes. Selectively accessible trajectories are indicated in black. One does not observe selectively accessible trajectories that have a small (A) valley depth or (B) valley width in both repression and expression abilities.
How does environmental variability affect these constraints? We first explore this question with individual trajectories starting with MK:acca. For instance, R can be increased through an operator mutation (MK:acca to MK:tcca) in the environment without ligand but then remains trapped, because the other mutations yield no additional improvements (Fig. 3 A and B and Fig. S2). However, switching to the other environment opens up various trajectories that increase E, such as the repressor mutation MK:tcca to MQ:tcca. After additional increase in E (MQ:tcca to MQ:tgca), the system becomes trapped again on a suboptimum. Concomitantly, note the low R for this genotype (Fig. S2A), which indicates an overall weakened binding between these LacI and operator variants. Switching back to the first environment now allows escape by a compensatory mutation that provides access to YQ:tggt by reestablishing binding and increasing R. We found that a significant fraction of the direct trajectories (21%) becomes accessible in this manner (Figs. S1 and S3). Accessibility was afforded by diverse patterns of environmental change, although not by all (Table S1). Interestingly, the number of accessible trajectories starting in the environment E [with isopropyl-β-d-1-thiogalactopyranoside (IPTG)] outnumbered the accessible trajectories starting in the environment R (without IPTG; 92 vs. 56), although the main improvement is made in R. Overall, these findings indicate that mutational pathways that fail to confer gradual optimization in either constant condition can do so when alternating between these same conditions.
Escape from genetic constraint in fluctuating environments. (A) Mutational trajectory accessible by continuous improvements in a changing environment corresponding to B and C. Red indicates mutated position. Forward arrows indicate mutations conferring increases in repression or expression ability, and backward arrows indicate decreasing or neutral steps. The two-color bars represent changing environments that confer continuous improvements. Without environmental changes, the system would be trapped at MK:tcca and MQ:tgca, where no additional improvements are possible in the current environment. (B) Expression and repression ability along the trajectory indicated in A. Data are represented as means (n = 3 or 4) ± SEMs. The trajectory starts at MK:acca (white circles) and ends at YQ:tggt (large black circles). This trajectory contains a valley in (Upper Right) repression ability (blue line) and (Lower) expression ability (orange line). (Upper Left) Blue lines indicate mutations that confer improvements in repression ability, and orange lines indicate mutations that confer improvements in expression ability. The gray area indicates the envelope of all trajectories. (C) Schematic representation of repression and expression ability along the trajectory indicated in A in the variable environment. Semitransparent lines indicate inaccessible mutations. Gray dotted lines indicate environmental changes allowing escape from suboptima.
Lock–key interactions trap genetic trajectories. Lock–key binding repressor–operator pairs constrain the trajectory depicted in Fig. 3 A and B (black trajectory). The black and cyan traces together form a valley in both (A) repression ability and (B) expression ability. The tcca operator is a low-expression operator in both environments (Fig. 1C and Fig. S5), which initially retains the repression ability at a high level. After this operator is mutated to tcga or tgca, the repression ability drops and is only restored to a higher repression ability when a new operator–repressor binding pair is formed at YQ:tgga. Cross-environmental tradeoffs can readily be observed by comparing A with B; initially, high levels of repression ability correspond to low levels of expression ability. The mutation to YK:tcca (in the cyan dashed trajectory) is neutral in repression ability according to the statistical test, and it is, therefore, considered not selectively accessible. Means are represented ± SEM. The same situation occurs for the mutation to YQ:tgga (in the cyan trajectory) in expression ability.
Genetic constraint in fluctuating environment circumvented by an escape route. (A) In the opaque trajectory, continuous positive selection on all mutations from MK:acca (large white circles) to YQ:tggt (large black circles) is prohibited, because the second mutation (purple) is detrimental in both environments. A mutational escape route (transparent line), bypassing the genetic constraint, is provided by a mutation that is beneficial in the absence of inducer (black arrows). The accessible (transparent) and nonaccessible (opaque) trajectories coalesce at Hamming distance 5. Dashed lines indicate detrimental effects. Means are represented ± SEM. (B) Schematic representation of the changes in repression and expression abilities vs. the number of mutations from MK:acca, with increasing levels of both repression and expression abilities on the vertical axis.
Environmental patterns allowing accessibility
The above findings cannot be understood from environmental interactions that alter the magnitude of mutational effects, because they would affect only the depth of constraining valleys. Rather, they indicate the importance of cross-environmental tradeoffs, in which increases in R occur at the expense of decreases in E, and vice versa, increases in E lead to decreases in R (Fig. 3 B and C). Such tradeoffs are also referred to as GxE interactions (49) and can give rise to higher-order GxGxE (34) and sign environmental epistasis (36). In the lac repressor–operator system, cross-environmental tradeoffs between E and R were pervasive (55% of all mutational steps in the direct trajectories) (Fig. S4) and can be understood mechanistically. For example, a low but significant level of repression can be maintained in the presence of inducer through residual binding (50). We found for several genotypes (22 genotypes) that the induced expression level is significantly lower than the highest measured level for the involved operator (Fig. 1C), consistent with residual binding of induced repressors reducing expression. Mutations that increase (decrease) the overall repressor–operator affinity in both environments will increase (decrease) both the repression ability without inducer as well as the residual repression with inducer, leading to opposite effects on R and E and hence, to cross-environmental tradeoffs (Figs. S4 and S5). These tradeoffs have consequences for the relations between constraints in different conditions. We found multiple local optima for each of two environments (3 in R and 13 in E), but none coincided at the same genotype. This feature allows trajectories to repeatedly surf ascending slopes and hence, traverse valleys in a ratchet-like manner: when trapped on a local optimum, the system can wait for an environmental change that enables repositioning on a new ascending slope.
Tradeoffs between expression and repression abilities. Histogram of the changes in mean repression and mean expression ability over all mutational steps in all direct 720 trajectories, with the effect on R on the horizontal axis and the effect on E on the vertical axis. The color gradient indicates fractions of mutations occurring at increasing Hamming distance from the starting genotype, with lighter green corresponding to larger Hamming distances. Mutations in the upper left and lower right quadrants indicate tradeoffs, whereas mutations falling in the upper right quadrant are accessible in both environments. Mutations in the lower left quadrant signal inaccessibility in both environments (34). Most mutations fall in either the lower right or the upper left quadrant, which indicates that a large amount of mutational steps exhibit cross-environmental tradeoff.
Expression measurements and absolute values. Mean expression levels of all operator–repressor combinations measured as the developed fluorescence in the fluorescein di-β-d-galactopyranoside expression-level assay (Methods) in either the absence or presence of the inducer IPTG. The expression level is expressed in arbitrary units. Computation of the population transition probabilities.
To assess the robustness of this evolutionary mode for different environmental conditions, we extended an evolution model based on a fixed environment Markov process (51) to include environmental fluctuations (Methods, Fig. 4 A–C, SI Results, sections S1–S6, and Figs. S6 and S7). We considered a discrete time Moran process in the strong selection weak mutation regime (52), in which mutational reversions are allowed and trajectories can be of arbitrary length (Methods) (53). Here, we first assume fitness and selection to be proportional to phenotype (E and R) and then, test different nonlinear and nonmonotonic phenotype–fitness relations (SI Results, section S7 and Fig. S8). Consistent with the observed constraint in fixed conditions (Fig. 2), we found that the rate to evolve to YQ:tggt from MK:acca (crossing rate kc) is null for either constant environment (Fig. 4D, fraction of time in environment E of 0 or 1). However, kc is consistently above zero when the environmental fluctuation rate kf is lower than the mutation rate km (Fig. 4D, blue and green lines) and maximized when kf = km, consistent with previous related work (39). That kc is maximal when kf = km can be understood as follows: for kf >> km, there is an effective averaging over the two environments, resulting in a constrained condition, whereas for kf << km, the waiting time for an environment-triggered escape is long enough to allow mutational escape in one of the two environments. Environmentally triggered escape is, thus, found to be robust to changes in the ratio between the times spent in the two environments (Fig. 4D).
Landscape crossing in stochastically alternating environments. (A) Environmental fluctuations and occurrence of mutations (crosses) (SI Methods and SI Results, section S1). Environments R and E refer to the environment selecting for repression or expression ability, respectively. (B) Schematic representation of genotype space. Large white and gray circles are start and end genotypes, respectively, of mutational trajectories. Arrows indicate increasing repression ability (R; blue) or expression ability (E; orange); arrow thickness reflects magnitude and hence, transition probability (SI Methods and SI Results, section S2). Shadowed arrows indicate one possible path of continuous improvement from the initial to the final genotypes. The structure of the space is schematic and does not reflect the actual system. (C) Schematic depiction of the Markov chain method for computing crossing rates. The probability vector lists all N genotypes, with the gray scale indicating the probability of populating a genotype at a given indicated time. Initially, only the beginning genotype is populated. The N × N environment-dependent transition probability matrices (colored squares) reflect the arrows in B: a matrix entry at position i,j indicates the transition probability from genotype i to genotype j. Each matrix multiplication yields a novel probability of genotype occupancy after a mutation occurred in a given environment. This illustration is schematic: we use an (infinite time limit) analytical solution for this process, considering a range of possible scenarios of environmental fluctuations. (D) Crossing rate kc as a function of the fraction of time spent in each environment for different environmental fluctuation rates kf. The unit of time is the time between two mutations. The red line indicates environment dwell time <<1, meaning that the environment fluctuates much faster than the time between mutations. The top green line to bottom blue line indicates environment dwell times of c = 1, 2, 5, 10, 20, and 50 (i.e., decreasing kf or decreasing frequency of environmental fluctuation) (SI Methods and SI Results, sections S3–S6). The crossing rate is the inverse of the mean number of mutations that are necessary to cross the landscape. The absolute maximum crossing rate is 0.17 (6−1; corresponding to six mutations being fixed).
(A) Diagram for the tree of events relative to a mutation + fixation/extinction step in a given environment and their associated probabilities: gi is the isogenic population with genotype gi,
Timeline of events for the Markov chain computation. Example of a temporal sequence of alternating environments for a scenario with an approximate 50%–50% time fraction between environment R and E (
Robustness of landscape crossing for nonlinear cost–benefit tradeoffs. (Left) b and zc are cost–benefit parameters used to convert relative expression levels into relative fitness (see Eq. S10 and Eq. S11). (Center) Relative fitness in each environment as a function of the measured average relative expression level for each genotype and for the b and zc values on the left. Dots indicate experimentally measured average relative expression levels, and lines indicate Eqs. S10 and S11 plotted as a guide for the eyes. (Right) Crossing rate kc in a variable environment scenario computed as previously described (environment dwell time
The above result remained valid for nonlinear phenotype–fitness relations, reflecting, for instance, expression costs that make the dependence weaker for higher E (SI Results, section S7 and Fig. S8, rows 1 and 2) (41, 54). The escape mechanism also remained robust for nonmonotonic phenotype–fitness relations, in which fitness first increases with increasing E and then decreases (Fig. S8, row 3). We found that the mechanism broke down when E increases do not confer fitness increases at all, reflecting a scenario in which lac expression costs outweighs the benefits in the inducing environment (Fig. S8, row 4). However, the latter is not observed experimentally (41, 42).
To further probe the requirements for crossing multipeaked landscapes, we considered computationally generated random landscapes. In constant conditions, mutational trajectories become trapped on these landscapes as expected, and hence, the probability to find selectively accessible paths is low (Fig. S9, absence of cross-environmental tradeoff). We found that, when we progressively increase the level of cross-environmental tradeoff (as quantified by the fraction of mutations having opposite effects in the two environments), the probability of finding accessible paths increases until practically all of the landscapes that were generated contained accessible paths in fluctuating conditions at a tradeoff level of 0.5 or more (Fig. S9). We note that this value was 0.47 for the experimentally determined landscape. Overall, these results indicated that cross-environmental tradeoffs are a central ingredient for overcoming constraining genetic interactions in fluctuating environments by positive Darwinian evolution.
Effect of cross-environmental tradeoff level between two alternating environments on the crossing ability of randomly generated landscapes. We consider a landscape over six loci of two alleles each (i.e., on a hypercube of dimension 6 and 64 genotypes). For each tradeoff level (i.e., the fraction of mutations having opposite effects in the two environments), random fitness values for each genotype and each environment are drawn to satisfy the tradeoff level between environments. The crossing ability is determined by computing for each tradeoff level the probability p over n = 1,000 (pairs of) landscapes of containing at least one selectively accessible path from the starting genotype to the final genotype in fluctuating conditions. Error bars represent the binomial error
To investigate how environmental variability affects genetic constraints within evolutionary trajectories, we have systematically mapped the genotype–phenotype landscape spanning two matching pairs of transcription factors and their DNA binding sites. Consistent with theoretical considerations on the evolution of molecular recognition, we find that the resulting landscapes in different contrasting environments are highly rugged: none of the mutational trajectories between the two pairs are selectively accessible by fixing one mutation at a time. Landscapes ruggedness has been studied previously for phenotypes that involve multiple possible genes and spontaneously evolved mutations (5, 16, 21, 55, 56). Here, we studied a phenotype that has a well-understood genetic basis, is restricted to one physical site, and hence, is controlled by a reduced set of mutations. The data show that epistatic interactions between these controlling mutations give rise to distinct peaks within genotype space.
Rugged landscapes constitute a specific type of genetic constraint that is distinct from, for instance, the presence of sign epistasis, which imposes specific fixation order but does allow accessible paths (1). Theoretical investigations have considered how escape from entrapment on suboptimal peaks is enabled; they addressed limiting negative selection (18, 19), simultaneously fixing mutations (20), recombination (21⇓–23), and drift (24⇓–26). These mechanisms do not permit adaptive evolution by positive selection of single mutations. Here, we find that such gradual adaptive evolution by positive selection is possible when considering temporal alternations between different environments, with cross-environmental tradeoffs as a second essential ingredient. These tradeoffs are pervasive in the quantified landscapes. They allow the environmental changes to displace local fitness peaks in such a way that the system can be ratcheted in an adaptive manner through fitness valleys by repeatedly climbing locally ascending slopes. We note that, at high mutation rates, such as for instance those observed in viral evolution, multiple variants may arise before one sweeps through the population, resulting in clonal interference and competition between variants. Although in this case, the sweep time rather than the time between mutations becomes the relevant timescale, the tradeoff ratchet mechanism remains valid in this regime.
Crossing fitness valleys by positive selection may seem contradictory. However, one must distinguish between the global fitness of the population and the relative fitness differences between individuals. When trapped, the relative fitness decreases conferred by mutations lead to negative selection within the population. In contrast, the global fitness decreases conferred by environmental changes affect the entire population and do not lead to negative selection. These global fitness decreases open up possibilities for mutations that confer relative fitness gains again. These repeated gains do not come completely for free, because the mechanism does not only drive forward but also backward: after a global fitness decrease induced by environmental change, the positively selected mutations that opened up point away from the ancestor but also, back to it. The occurrence of backward mutations does not cancel the positive selection of this evolutionary mode but merely reduces its efficacy. Thus, environmental changes can serve to overcome fitness decreases that constrain genotypes in constant environments.
The role of tradeoffs in facilitating adaptive change is notable, because they are typically associated with adaptive constraint. Tradeoffs can be rationalized mechanistically for regulatory phenotypes that rely on lock–key molecular recognition, which has been shown to be important genetic material for evolutionary novelties (57, 58). Evolutionary optimization in the presence of tradeoffs can be seen as a sequence of compensatory mutations; recognition loss caused by mutating one binding partner can be restored by complementary mutations in the other binding partner. Note that tradeoffs are present more generally in nature (59, 60) along with environmental fluctuations (61), and hence, the reported adaptive mechanism is not necessarily limited to regulatory phenotypes.
Organisms are known to occupy remarkable wide ranges of contrasting environments. Our findings suggest that the regulatory systems that enable these wide environmental niches are constrained by pervasive (reciprocal) sign epistatic interactions, resulting in rugged landscapes. An inability to optimize regulatory responses to environmental change can result in stasis and the emergence of environmental specialists. However, we find that the proposed mode of adaptation can overcome these constraints and facilitate the evolution of generalists that can occupy diverse niches. Interestingly, such adaptability comes at the cost of weakened environmental robustness of desirable genotypes. Indeed, the stabilizing selection of such a genotype at the top of a fitness peak can be disrupted by tradeoffs and environmental change, which together, can displace peaks. Some environmental robustness may alternatively be afforded by regions in genotype space that are flat in both environments. However, except for a “neutral chain” (YK:agga-YQ:agga-MQ:agga-MQ:tgga-MQ:tcga), no other neutral region was observed in our landscape.
Genetic constraints are commonly regarded as a key factor in adaptive stasis and major evolutionary transitions, such as for instance the evolution of sex, multicellularity or symbiosis, and eusociality (62, 63). Overcoming genetic constraints by tradeoffs in a fluctuating environment could be relevant to these transitions given its few requirements and their ubiquitous presence. Additionally, this evolutionary mode may have implications for clinically or biotechnologically relevant problems, such as the treatment of infections with multidrug protocols (64⇓–66) and the evolutionary engineering of antibodies (67). It will be of interest to test whether this adaptive mode affects the dynamics of experimental evolution in variable environments.
We used E. coli strain BW23473 (68) F−, Δ(argF-lac)169, ΔuidA3::pir+, recA1, rpoS396(Am), endA9(del-ins)::FRT, rph-1, hsdR514, rob-1, and creC510 (Yale Coli Genetic Stock Collection 7837) for all experiments.
We used a single-copy plasmid that can be induced to a multicopy plasmid for cloning purposes (pETcoco-2; Novagen) to mimic natural expression levels. The plac, lacZ, and lacY WT sequences were derived from E. coli strain MG1655 (Yale Coli Genetic Stock Collection 6300). To facilitate the insertion of mutants (see below), the plac region was amplified until the ninth codon in the coding sequence of lacZ, and a BmtI restriction site was introduced, whereas the WT amino acid sequence was preserved. The remaining lacZ and lacY sequence was amplified and inserted in the polylinker sequence of pETcoco-2 to make pMdV53. Additional information is in SI Methods.
plac operator variants differed in the 4, 5, 5′, and 4′ nucleotide residues as depicted in Fig. 1 B and C. This O1 operator was made palindromic by the deletion of the central c–g base pair (69). In LacI, amino acid residues 17 and 18 were altered using optimal codons for the amino acid variants. Additional information is in SI Methods.
All growth and expression measurements were performed in Mops EZ Rich Defined Medium (Teknova) with 0.2% glucose as the carbon source and supplemented with 1 mM thiamine-HCl, 1 mM uracil, and appropriate antibiotics. For expression ability measurements, the medium was supplemented with 1 mM IPTG.
Cultures were grown at 37 °C in PerkinElmer Victor3 and VictorX3 plate readers with 200 µL medium per well in a black clear-bottom 96-well plate (NUNC 165305). Cells were fixed after the cultures had reached an OD in the plate reader of at least 0.015 and at most 0.07 by adding 20 µL fluorescein di-β-d-galactopyranoside fixation solution [109 µM fluorescein di-β-d-galactopyranoside (MarkerGene Technologies Inc.), 0.15% formaldehyde, 0.04% DMSO in deionized water]. Fluorescence development was measured every 8 min (excitation = 480 nm and emission = 535 nm), and the OD was measured at 600 nm. When cells were not induced with IPTG during growth, 1 mM IPTG was added to each well immediately before or after fixation, because the assay is sensitive to the amount of IPTG present in the medium (7). Analysis of the fluorescence to quantify the LacZ expression level was as described in ref. 7. Additional information is in SI Methods.
For each genotype, expression measurements were performed in either triplicates or quadruplicates. Expression levels measured in one 96-well plate were normalized to a control genotype on that plate. Throughout all of the analysis, whenever the values between two genotypes were compared (either expression or repression ability values or relative fitness values obtained through the cost–benefit treatment), a two-tailed Welch’s t test was performed to test the null hypothesis that the two genotypes do not have significantly different values (hypothesis of neutrality). The null hypothesis was rejected if the p value obtained from the test was lower than a significant level of 0.05. If the null hypothesis was not rejected, the mutations distinguishing the two genotypes were considered neutral (e.g., we would take the condition
We modeled the dynamics of evolution in fluctuating environment using a discrete time Moran process in the regime of strong selection weak mutation and through a heterogeneous Markov process. Additional information is in SI Methods.
We used a single-copy plasmid that can be induced to a multicopy plasmid for cloning purposes (pETcoco-2; Novagen) to mimic natural expression levels; lacI was deleted from the pETcoco-2 and replaced with the T0 and Trev terminators from the SacI-SpeI fragment of pInv-110 (70). The introduced placI was derived from the overexpressing placIq sequence from pTrc99A (71). By site-directed mutagenesis, the t at −35 was replaced with the c nucleotide (Quick Change Site-Directed Mutagenesis Kit; Qiagen), leading to native LacI expression levels. The auxiliary operator O3, located at the end of the coding sequence of the WT lacI, has a low affinity for LacI [about 300-fold lower than O1 (72)]. It is involved in weak autorepression of the lac repressor (73) and was, hence, deleted while preserving the LacI amino acid sequence (7); placI combined with lacI from pRD007 (which lacks the auxiliary O3 operator sequence) (7) and a spacer sequence derived from Drosophila Kinesin (pCascade5) by overlap PCR (Pfu Ultra; Invitrogen) were inserted in pETcoco-2. The plac, lacZ, and lacY WT sequences were derived from Escherichia coli strain MG1655 (Yale Coli Genetic Stock Collection 6300); lacA was not incorporated in the assayed operon, because it decreases the active concentration of IPTG in the cell (74, 75). To facilitate the insertion of mutants (see below), the plac region was amplified until the ninth codon in the coding sequence of lacZ (Phusion Polymerase; Finnzymes; New England Biolabs), and a BmtI restriction site was introduced, whereas the WT amino acid sequence was preserved. This fragment was inserted in the pETcoco-2 vector backbone. The remaining lacZ and lacY sequence was amplified (Phusion Polymerase; Finnzymes; New England Biolabs) and inserted in the polylinker sequence of pETcoco-2 to make pMdV53. All restriction enzymes and T4 ligase used were obtained from New England Biolabs, and all primers were obtained from MWG Eurofins.
lacI and plac operator mutants were constructed separately; plac operator variants differed in the 4, 5, 5′, and 4′ nucleotide residues as depicted in Fig. 1 B and C. This O1 operator was made palindromic by the deletion of the central c–g base pair (69). In LacI, amino acid residues 17 and 18 were altered using optimal codons for the amino acid variants. Half-sides of each lacI or plac operator variant were constructed using Pfu Ultra Polymerase (Invitrogen) with one flanking primer and one primer harboring the desired nucleotide variation. After gel purification (Gel Extraction Kit; Qiagen) and ExoI Treatment (New England Biolabs), the half-sides with the desired nucleotide variation in the region of overlap were amplified by overlap PCR to obtain full-length lacI or plac operator variants. Combinations of lacI-spacer-plac mutants were constructed by overlap PCR (Pfu Ultra; Invitrogen), inserted in the pMdV53 backbone, and verified by Sanger sequencing.
Cultures were grown at 37 °C in PerkinElmer Victor3 and VictorX3 Plate Readers with 200 µL medium per well in a black clear-bottom 96-well plate (NUNC 165305). OD at 600 nm was recorded every 4 min, and every 29 min, 9 µL sterile water was added to each well to counteract evaporation. Between measurements, the well plate was shaken at double orbit with a diameter of 2 mm. Cells were fixed after the cultures had reached an OD in the plate reader of at least 0.015 and at most 0.07 by adding 20 µL fluorescein di-β-d-galactopyranoside fixation solution [109 µM fluorescein di-β-d-galactopyranoside (MarkerGene Technologies Inc.), 0.15% formaldehyde, 0.04% DMSO in deionized water]. Fluorescence development was measured every 8 min (excitation = 480 nm and emission = 535 nm) along with the OD at 600 nm. Shaking and dispensing conditions were as described above. When cells were not induced with IPTG during growth, 1 mM IPTG was added to each well immediately before or after fixation, because the assay is sensitive to the amount of IPTG present in the medium (7). Analysis of the fluorescence to quantify the LacZ expression level was as described in ref. 7.
We modeled the dynamics of evolution in fluctuating environment using a discrete time Moran process in the regime of strong selection weak mutation (SSWM) and through a heterogeneous Markov process.
Briefly, for each environment [without IPTG (environment R) or with IPTG (environment E)], a
We compute the landscape crossing rate for different scenarios of temporally alternating environments, where the frequency of the environmental changes and the time fraction between environments are taken as parameters. To this aim, we determine the mean first passage time to genotype YQ:tggt for a population evolving from genotype MK:acca and experiencing alternately two types of selective environments: one being selective for the repression ability (environment R: absence of IPTG) and one being selective for the expression ability (environment E: presence of IPTG).
We consider an infinite-sized population following a discrete time Moran process in the regime of SSWM (52, 77). In this regime, the time for fixation (or extinction) of a newly appeared mutant allele in the population is far less than the time between mutations. The population is, therefore, essentially assumed to be isogenic.
Here, we formulate a heterogeneous Markov process that models the evolution of the population by a probability density diffusion process through a dynamic adaptive landscape. We build on similar approaches for fixed environments (51, 53). Below, we explain the construction of the model.
The population evolves in a landscape of 26 = 64 genotypes subject to the following dynamics, where the unit of time is taken as the time for a new mutant allele to appear in the population.
At time 0, the population is assumed to carry the starting genotype g1 = MK:acca.
At each time step, a single-point mutant genotype gj appears in the population of genotype gi, with j ≠ i, and this mutation goes to either fixation or extinction, with probabilities that depend on fj/fi, the relative fitness value of the newly appeared genotype with respect to the fitness of other individuals in the population (Fig. S6).
For each genotype, the fitness depends on the environment. Environment R selects for repression ability, whereas environment E selects for expression ability. One can consider various ways in which fitness depends on the measured phenotype. Here, we first consider a direct equivalence between fitness and ability (
For each environment, a specific Markov matrix is computed that contains the transition probabilities between isogenic states of the population. In the Moran model, the general expressions for the fixation probabilities of a genotype gj of average fitness
The Markov process and the associated Markov matrices describe the evolution dynamics on the genotype space by the population submitted to the corresponding transition probabilities. Our model allows for reversal of previously fixed mutation (e.g., after reversal of the environment and/or in a modified genetic background). As such, the considered evolutionary process is not unidirectional but can be a back and forth evolution on the landscape, despite neglecting deleterious and neutral substitutions. Therefore, computing the evolution rate from genotype MK:acca = g1 to genotype YQ:tggt = g64 is equivalent to computing the mean first passage time of the population to the genotype YQ:tggt.
For a Markov process described by a single Markov matrix M, the mean first passage time of the system to a specific state sr can be computed by (i) setting sr as an absorbing state, i.e., fixing
The components Nij of N give the expected number of time steps in which the system is found in state sj before it is absorbed starting from state si. The mean first passage time to state sr assuming si to be the starting state is then simply given by
To compute the mean first passage time to genotype YQ:tggt in alternating environments, we need to compute the mean first passage time in the case of a Markov process described by the alternating application of the two transition matrices MR and ME and their associated transient state submatrices
To assess different frequencies of environmental changes and different values of fraction of time spent in each environment, we define the probability w to be in environment R and a dwell time c being the number of successive mutation + fixation/extinction cycles occurring before any new possible change in environment is allowed (i.e., the transition matrix applied is allowed to change only each time a number c of mutation + fixation/extinction cycles has occurred) (Fig. 4 A–C and Figs. S6 and S7).
The fundamental matrix
When the environment fluctuates much faster than the time between mutations, the isogenic population
Analyzing the accessibility of all of the paths individually (i.e., not only the 6! = 720 direct paths but also, indirect paths that allow for mutation reversal) is computationally intractable. For instance, for two genotypes distant from just five mutations, there is a total of 18,651,552,840 simple paths (79). Nevertheless, computing the landscape crossing rate kc (SI Results, sections S1–S4) indicates if the landscape is crossable, taking indirect paths into account in the analysis.
Fig. 4D shows that, for fractions of time w = 0 or 1 (i.e., for fixed environment of either R or E), the landscape crossing rates kc are identically null in all cases. This result shows that the landscape is not crossable in a fixed environment of either R or E and that, therefore, all indirect paths also contain a valley.
Fitness (or growth rate) might not be a linear or even monotonic function of expression levels. First, the benefit associated with the expression of a beneficial gene can level off beyond a certain expression level. Second, expression of a beneficial gene can also come at a cost that increases with the level of expression of the gene. To test the robustness of our results, we analyzed the landscape crossing rates in the context of a cost–benefit scenario.
To explore a range of cost–benefit scenarios, we used a parameterized version of a cost–benefit function applied previously to the lac operon (41). In the environment that selects for the expression ability (environment E), we define the following relative fitness
b [a benefit coefficient; b = benefit/(benefit + cost)],
zc (a critical relative expression level characterizing the nonlinearity of the benefit), and
z (the relative expression level of the corresponding genotype).
The term
In the environment selecting for repression ability (environment R), there is no benefit associated with expression, and the benefit term of Eq. S10 drops. The relative fitness then simply becomes a linear decreasing function of the relative expression level:
For given values of b and zc, the fitness values
Fig. S8 shows the result for four different sets of the b and zc parameters, exploring increasingly extreme cost–benefit scenarios.
This analysis illustrates that, except for extreme cases (uniformly null landscape crossing rate) (Fig. S8, row 4), landscape crossing stays feasible, and the mechanism can sustain saturation in benefit (Fig. S8, rows 1 and 2) or highly nonlinear cost–benefit relationships with moderate fitness decrease (Fig. S8, row 3).
We thank D. M. Weinreich, J. A. G. M. de Visser, T. Paixão, J. Polechová, T. Friedlander, and A. E. Mayo for reading and commenting on earlier versions of the manuscript and B. Houchmandzadeh, O. Rivoire, and M. Hemery for discussions and suggestions on the Markov computation. Furthermore, we thank F. J. Poelwijk for sharing plasmid pCascade5 and pRD007 and Y. Yokobayashi for sharing plasmid pINV-110. We also thank the anonymous reviewers for remarks on the initial version of the manuscript. This work is part of the research program of the Foundation for Fundamental Research on Matter, which is part of the Netherlands Organization for Scientific Research (NWO). M.G.J.d.V. was (partially) funded by NWO Earth and Life Sciences (ALW), project 863.14.015.
Author contributions: M.G.J.d.V., A.D., and S.J.T. designed research; M.G.J.d.V., A.D., and V.S. performed research; M.G.J.d.V., A.D., and S.J.T. analyzed data; and M.G.J.d.V., A.D., and S.J.T. wrote the paper.
The authors declare no conflict of interest.
This article is a PNAS Direct Submission.
This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1510282112/-/DCSupplemental.
