RESEARCH ARTICLE 2016 © The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC). 10.1126/sciadv.1500921

MOLECULAR EVOLUTION

Environmental changes bridge evolutionary valleys Barrett Steinberg* and Marc Ostermeier†

INTRODUCTION The fitness landscape metaphor is often depicted as a three-dimensional hilly surface with distinct fitness peaks and valleys (Fig. 1). This representation flattens the high dimensionality of actual landscapes to aid our understanding of what a fitness landscape represents. Fitness landscapes depict the change in fitness along possible evolutionary pathways (1), and thus, the shape of the fitness landscape affects the dynamics and outcomes of evolution. Accordingly, much effort has gone into understanding the shape of the landscape and the degree to which ruggedness manifests in actual evolutionary scenarios (2–5). Ruggedness results from complex interactions between mutations, especially those exhibiting reciprocal sign epistasis (6). In rugged landscapes, local fitness maxima exist as “traps,” and therefore, the evolutionary endpoint of adaptation from a given starting point is not predetermined as it would be in a smooth landscape. The endpoint of each evolutionary pathway depends on selective history, landscape topology, and chance. However, one way in which evolution might escape traps and cross fitness valleys is through changes in the environment that remodel the landscape (1). Natural environments are intrinsically variable, and evolutionary pathways are contingent on past environmental interactions and selective history (7). Thus, the dynamics and complexity of selective pressures are likely to contribute to the routes of natural evolutionary pathways. There is a growing appreciation of the complex and important relationship between the environment and the fitness and epistatic effects of mutations (8, 9). Understanding the roles that environmental changes and landscape ruggedness play in adaptive pathways would allow insights into the predictability and potentiality of evolutionary pathways. How changes in environment or selection affect evolutionary outcomes is underexplored experimentally (8). One way to address this topic is to ask whether increasingly stringent positive selection is the best strategy for optimization. This strategy is ubiquitously used in directed Department of Chemical and Biomolecular Engineering, Johns Hopkins University, 3400 North Charles Street, Baltimore, MD 21218, USA. *Present address: Editas Medicine, 300 Third Street, First Floor, Cambridge, MA 02142, USA. †Corresponding author. E-mail: [email protected]

Steinberg and Ostermeier Sci. Adv. 2016; 2 : e1500921

22 January 2016

evolution experiments; however, on a rugged landscape, positive selection alone may quickly reach a local maximum, limiting evolutionary capacity. Studies of evolution under strong selection have indicated that evolution can follow few mutational pathways using positive selection, supporting the predictability of evolutionary endpoints (10). How evolution might cross fitness valleys remains an unresolved problem (1, 2). Models and simulations have suggested that random strategies, such as stochastic tunneling or recombination, may occasionally result in the crossing of fitness valleys (11, 12). Other studies have examined the effects of changing the starting point of evolution to arrive at different fitness maxima (13), but experiments rarely fundamentally change the nature of selection over the course of experimental evolution. Exceptions such as selection on alternating antibiotics (14) and a regimen of neutral drift followed by positive selection (15–21) have met with mixed success in evolving improved function.

RESULTS AND DISCUSSION Experimental design We wondered whether complex selection dynamics could facilitate the crossing of fitness valleys. Specifically, we examined whether environmental changes causing dynamic shifts in a protein fitness landscape would result in evolutionary pathways ending in distinct, potentially superior evolutionary endpoints. We designed an experiment in which environmental changes enabled selection toward sequences previously residing in a low-fitness valley and tested whether such an evolutionary path would arrive at distinct and superior outcomes (Figs. 1 and 2). We used Escherichia coli having a bandpass genetic circuit designed for tunable selection of b-lactam antibiotic resistance as a model system for systematically examining such a scenario (22). For cells with this gene circuit, b-lactam antibiotics such as ampicillin are required to induce the expression of tetracycline resistance and to do so even at sublethal levels of the b-lactam. When cells are plated in the presence of tetracycline and ampicillin, growth requires enough b-lactamase activity to degrade the ampicillin below lethal levels, but not too much so as to remove the signal that induces tetracycline resistance. This genetic circuit allowed us to 1 of 9

Downloaded from http://advances.sciencemag.org/ on March 30, 2016

In the basic fitness landscape metaphor for molecular evolution, evolutionary pathways are presumed to follow uphill steps of increasing fitness. How evolution can cross fitness valleys is an open question. One possibility is that environmental changes alter the fitness landscape such that low-fitness sequences reside on a hill in alternate environments. We experimentally test this hypothesis on the antibiotic resistance gene TEM-15 b-lactamase by comparing four evolutionary strategies shaped by environmental changes. The strategy that included initial steps of selecting for low antibiotic resistance (negative selection) produced superior alleles compared with the other three strategies. We comprehensively examined possible evolutionary pathways leading to one such highfitness allele and found that an initially deleterious mutation is key to the allele’s evolutionary history. This mutation is an initial gateway to an otherwise relatively inaccessible area of sequence space and participates in higher-order, positive epistasis with a number of neutral to slightly beneficial mutations. The ability of negative selection and environmental changes to provide access to novel fitness peaks has important implications for natural evolutionary mechanisms and applied directed evolution.

RESEARCH ARTICLE evolution experiments on TEM-1 that imposed positive selection for cefotaxime resistance found these same three mutations in the highestresistance alleles (13, 24–28). The GKTS allele, which adds the A42G mutation, has the highest resistance observed in these experiments (10, 27, 29), though other combinations of amino acids at these positions can confer slightly higher resistance (30). These directed evolution studies have used a typical methodology: mutation of the gene, positive selection for higher resistance, and iteration of these processes over several rounds.

B

Evolution experiments We subjected TEM-15 to eight rounds of experimental evolution using four different selection strategies: positive selection, neutral selection, negative selection, and oscillating selection (Fig. 2, B to D). As a control for evolution from low-resistance alleles using only positive selection, we also subjected TEM-1 to positive selection. The differences in the selection strategies occurred in the first three rounds of evolution and were followed by five rounds of positive selection to climb toward local fitness maxima. The positive selection experiments represent evolution in a fixed environment with a static landscape and are the typical directed evolution strategy. Because of the bandpass selection, our neutral selection avoids selection for alleles conferring resistance far above the starting allele, unlike previous directed evolution studies using neutral drift (15–21). During the negative selection steps in the oscillating and negative selection schemes, we set the level of resistance needed for growth above the background level of resistance exhibited by E. coli without a b-lactamase gene. We reasoned that mutations that completely deactivate a protein might be exceedingly common (such as nonsense mutations), whereas maintenance of above-background level of activity would better reflect the desired scenario that the gene is still required for cell viability and that it encodes a protein with some enzyme activity. We constructed five plasmids with distinct DNA barcodes identifying the five experiments to monitor and ensure against cross-contamination between the five different selection experiments. One round of evolution consisted of comprehensive codon mutagenesis on the plasmid-borne blactamase genes (isolated from the previous round), transformation of the newly mutated plasmid library into fresh E. coli cells, and selection on agar plates for the desired level of cefotaxime resistance (Fig. 2A). At least 106 library members were generated and subjected to selection in each round. Selections were tuned such that typically >100 and sometimes thousands of independent clones appeared at each step to limit the effects of evolutionary bottlenecking on the experimental outcome (data S1). Populations under negative and oscillating selection recovered from negative selection to resistance above that of TEM-15 after one round of positive selection (Fig. 2, D and E). This result supports the idea of the commonality of compensating mutations or suppressor mutations in restoring protein activity (31, 32) but might also reflect a reversal of deleterious mutations in combination with beneficial mutations.

C

Fig. 1. Simplified model representation of how evolution in alternating environments that modulate the fitness landscape can lead to the crossing of fitness valleys. The xy plane represents sequence space, whereas the z axis represents fitness. The solid red arrow indicates an evolutionary pathway following increasing positive selection pressure. (A) A sequence (circle) sits at an optimum separated from global optima by a fitness valley. (B) A change in the environment alters the landscape such that evolution brings the gene to a new sequence that previously resided in the valley. (C) Upon return to the original environment, the sequence resides in the valley but can now evolve to the global optimum, which lies uphill.

perform selections for b-lactam resistance within a narrow window of resistance (that is, one can select for resistance above a first threshold but below a second, higher threshold) (Fig. 2A). Thus, selections for low but nonzero resistance can be performed. We term such selections “negative” selections in the sense that high fitness (that is, the ability of the cells to grow) in such an environment is achieved only by alleles conferring low levels of b-lactam resistance. Similarly, neutral selection for maintenance of about the original level of resistance can be performed. These bandpass selections are tuned by the level of the b-lactam antibiotic in the environment and require the addition of tetracycline. Thus, these environments reshape the fitness landscape analogously to the depiction in Fig. 1, and our engineered bacteria are a suitable model system for examining evolution under landscape-altering environmental changes. The b-lactamase allele TEM-1 provides high antibiotic resistance to ampicillin but very low resistance to cefotaxime. Clinically isolated alleles of TEM-1 conferring increased cefotaxime resistance include TEM-15 (E104K/G238S) and TEM-52 (E104K/M182T/G238S), with the latter increasing resistance 4000-fold (23). Six independent directed Steinberg and Ostermeier Sci. Adv. 2016; 2 : e1500921

22 January 2016

Evolution outcomes At the end of the experiment, all strategies resulted in bacteria with a high cefotaxime resistance phenotype. The coding sequences of the genes from 47 of the most resistant colonies from each selection strategy were sequenced (data S2 to S6). All selection strategies showed sequence clustering, indicating that each strategy, for the most part, reached distinct areas of sequence space (Fig. 3A). This result emphasizes how much the environment and the nature of selective pressure can shape evolutionary 2 of 9

Downloaded from http://advances.sciencemag.org/ on March 30, 2016

A

RESEARCH ARTICLE B 10 2

10 0 TEM-15 TEM-1

10–2

0

2

4

6

Neutral Selected cefotaxime resistance (µg/ml)

Comprehensive saturation mutagenesis

C

Positive Selected cefotaxime resistance (µg/ml)

A

8

10 2

10 0 10–2

0

2

4

6

8

Round

Selection

E

Negative 10 2

10 0 10–2

0

2

4

Round

6

8

Oscillating 10 2

10 0 10–2

0

2

4

6

8

Round

Fig. 2. Schematic and course of experimental evolution on the b-lactamase gene. (A) General schematic. After a comprehensive saturation mutagenesis, selection proceeds in one of three manners depending on the environment. Positive selection (blue) enriches for sequences that provide resistance above a certain level. Neutral selection (yellow) enriches for sequences conferring about the same resistance as the starting gene. Negative selection (red) enriches for sequences conferring very low but above-background resistance. The latter two selection schemes are enabled by an engineered bandpass gene circuit present in the cell (22) and require the addition of the antibiotic tetracycline and different levels of cefotaxime, which constitute the difference in environment from the positive selection. (B to E) Course of the evolution experiments. Experiments differ in the first three rounds in which the bacteria underwent (B) positive selection, (C) neutral selection, (D) negative selection, or (E) oscillating selection. Each round represents both mutagenesis (line segments) and selection at the indicated level of cefotaxime (endpoints). Red circles indicate bandpass selections, and blue circles indicate positive selections. The dashed horizontal line indicates the resistance level conferred by TEM-15. Graphed resistance values indicate a lower limit on the resistance values of the entire population during selection (that is, the concentration at which the selection was performed).

outcomes. Sequences resulting from negative selection showed both the highest average genetic distance from TEM-15 (Fig. 3A and fig. S1) and the highest diversity (fig. S2). We discerned no apparent pattern in the location of mutations on the protein structure (fig. S3). We chose seven distinct sequences from each selection experiment (based on diversity within the tree) to subclone and analyze the effects of mutations on the coding sequence as opposed to those elsewhere on the plasmid. The coding sequences were ligated into a fresh plasmid backbone and transformed into fresh E. coli cells for resistance testing. All libraries produced alleles that conferred resistance higher than that conferred by GKTS, though an approximately eightfold drop in resistance upon subcloning of the genes indicated that mutations elsewhere on the plasmids contributed to the high resistance observed after round 8. We supposed that mutations elsewhere likely consisted of those that increased plasmid copy number or promoter strength. Strikingly, the negative selection strategy evolved alleles that conferred the highest resistance, with five of the seven alleles providing an eightfold resistance increase over that of GKTS (Fig. 3B). This result is interesting in light of theoretical predictions that a low starting fitness may dramatically enhance the chances for a population to reach the global maximum (33). Deep sequencing of populations during the course of the negative selection experiment We next sought insights into how negative selection produced this outcome. Deep sequencing of the population at each round of the negative selection regimen revealed a high expansion of genetic diversity during negative selection, consistent with our expectation that deleterious mutations are common (fig. S4). In contrast, during the subsequent positive Steinberg and Ostermeier Sci. Adv. 2016; 2 : e1500921

22 January 2016

selection steps, increasing purifying selection is apparent (fig. S4). In particular, the A184I mutation was highly enriched in the population (figs. S5 and S6). However, a large amount of diversity was maintained throughout the experiment. A total of 5836 unique sequences were found in deep sequencing data after the first positive selection, demonstrating that, despite rapidly changing modes of selection, we maintained a high diversity of alleles following the negative selection regimen. In addition, 85% (40 of 47) of the sequences from the final selection were unique, with an average difference of 6.1 AA mutations between unique mutants. Evolutionary pathway reconstruction using comprehensive fitness measurements To determine whether the environmental change resulting in negative selection was a key factor in the evolutionary history of an allele conferring a high resistance, we attempted to reconstruct all possible evolutionary intermediates between TEM-15 and one of the fittest alleles resulting from the negative selection regimen (BS-NEG-4) (fig. S7). This allele contained nine missense mutations relative to TEM-15 [A11S, A42G, L49M, Q90R, A184I, I208W, A224V, F230S, and R241F in Ambler notation (34)] and conferred an eightfold higher cefotaxime resistance than GKTS. Among the mutations in BS-NEG-4, A11S, A42G, L49M, and A224V have been reported in other TEM alleles, though none were present together in the same allele (35). We created a library designed to contain all 512 combinations of these nine mutations to determine the protein fitness (w) of all possible evolutionary intermediates in the evolution of BS-NEG-4 (assuming no mutation reversals or multiple mutations at the same amino acid position). Protein fitness quantifies the relative ability of 3 of 9

Downloaded from http://advances.sciencemag.org/ on March 30, 2016

β-Lactam resistance

D

Selected cefotaxime resistance (µg/ml)

Negative Neutral Positive

Selected cefotaxime resistance (µg/ml)

Selection coefficient

Round

RESEARCH ARTICLE A

Neutral Oscillating

Positive (TEM-15)

Negative 1. 0

Fold increase in cefotaxime resistance over TEM-1

B 14,000 12,000 10,000 8000 6000 4000 2000 0 Positive TEM-1

Positive

Neutral

Negative

Oscillating

TEM-15

Fig. 3. Diversity and conferred resistance of alleles after the eighth round of evolution. (A) Unrooted tree of Hamming distances between randomly selected alleles. Branches are color-coded by the evolution experiment as indicated. Scale bar indicates 1 unit of Hamming distance (that is, one missense mutation). (B) Increase in resistance conferred by seven alleles from each selection strategy as compared to TEM-1. Resistance was measured using the MIC assay. The seven alleles were selected from among those depicted in the tree to represent diverse sequences. The dashed line represents resistance conferred by the GKTS allele, whereas the TEM-15 allele conferred a 64-fold higher resistance relative to TEM-1.

the protein to confer a cefotaxime resistance phenotype: the higher the resistance, the higher the fitness. We determined the fitness of 98.8% (506 of 512) of all possible combinations of these nine mutations in a high-throughput manner involving bandpass selection and deep sequencing (data S7) (36). This method uses the bandpass circuit. Cells with a given level of cefotaxime resistance will grow best within a certain range of cefotaxime concentrations. A cell with higher resistance requires high concentrations, and a cell with lower resistance requires lower concentrations. Fitness is thus determined by plating the library on plates containing tetracycline and different cefotaxime concentrations, deep seSteinberg and Ostermeier Sci. Adv. 2016; 2 : e1500921

22 January 2016

Analysis of fitness landscapes and evolution pathways using selection-weighted attraction graphing We sought a method to visualize the high-dimensional space of the fitness landscape composed of the 512 potential evolutionary intermediates between TEM-15 and BS-NEG-4. For this purpose, we introduce selectionweighted attraction graphing (SWAG), a force-directed graphing method in which nodes represent sequences and edges connect nodes differing by only one missense mutation (Fig. 4). Forces between connected nodes are weighted proportionally to the selection strength, as determined by the nodes’ respective fitness values. These forces are applied on the xy plane, and the graph can be depicted in three dimensions by adding fitness as the z axis. In the graph, edges represent possible single steps in the evolutionary pathway, and the distance between connected nodes on the xy plane is inversely proportional to the likelihood of taking that step under positive selection. SWAG analysis of the intermediates between TEM-15 and BS-NEG-4 resulted in eight major clusters of sequences representing local fitness peaks (Fig. 4, A and B). Such clustering was robust to the weighting of selection strength in the forces. No clustering was observed when fitness values were randomly reassigned to the nodes. Clusters comprise groups of sequences with many interconnections (fig. S8) that tend to lead in succession to ever-increasing fitness values. The graph can be seen to have two “ranges” of four clusters each separated by a large gap devoid of nodes (Fig. 4). The two ranges differ by whether their sequences contain F230S (Fig. 4C). The range containing F230S lies far from the start point TEM15 and contains the endpoint BS-NEG-4. The proximal range contains TEM-15 and nearby nodes. Interrange edges primarily connect low-fitness alleles, typically with the lower fitness node lying in the F230S range. Although we cannot exclude the possibility that BS-NEG-4 could result from positive selection, we can say that the evolution of BS-NEG-4 is extremely unlikely under positive selection. Of the 9! = 362,880 possible pathways between TEM-15 and BS-NEG-4, we could analyze 347,616 (the missing pathways are due to six alleles not being frequent enough in the library to analyze fitness data on those alleles). Of the 347,616 pathways analyzed, only 16 involve increasing fitness as each mutation is added (Fig. 5). Furthermore, evolution under positive selection is more likely to terminate at intermediates that are local maxima than to reach BS-NEG-4 (Fig. 6). The SWAG analysis pictorially depicts this result. First, eight of the possible nine initial mutations (the exception being F230S) are located in the proximal range (fig. S9). Second, there are many more steps within a cluster than between clusters (fig. S8A). Finally, steps within a cluster tend to be fitter sequences than steps between clusters (fig. S8B), making intercluster steps even less likely. Thus, starting from TEM15, the most likely mutational pathways in a constant, positive selection environment progress first into one of the clusters of the proximal range and then preferentially climb within that cluster to local fitness maxima. 4 of 9

Downloaded from http://advances.sciencemag.org/ on March 30, 2016

Positive (TEM-1)

quencing the b-lactamase gene on each plate, counting how many times each allele appears on each plate, and calculating the center of the concentration range at which the allele allows growth. This center concentration can be thought of as a proxy for the minimum inhibitory concentration (MIC) for cefotaxime. Fitness is expressed as this center concentration normalized to the center concentration of a reference (in this case, the resistance conferred by TEM-15). One mutation, F230S, was particularly deleterious (wF230S = 0.09 ± 0.01, with the fitness of TEM-15 set to a value of 1.0) and remained so even in combination with one to two other mutations, making it a highly likely starting point on the evolutionary pathway resulting in BS-NEG-4.

RESEARCH ARTICLE BS-NEG-4

A

B

BS-NEG-4

w TEM-15

D

+F230S

w

w –F230S

Fig. 4. SWAG of the possible evolutionary intermediates between TEM-15 and BS-NEG-4. Nodes are sequences, and edges connect nodes that differ by one missense mutation. Forces between connected nodes are weighted proportional to the selection strength, as determined by the node’s fitness values. (A) Two-dimensional SWAG. Edges are curved to aid visualization. (B to D) Three-dimensional SWAG has protein fitness as the z axis and is color-coded by (B) cluster, (C) whether the sequence contains F230S, and (D) protein fitness value.

However, when the environment changes such that low resistance is beneficial (that is, under negative selection), mutational steps to low-fitness values in the depths of the F230S range become advantageous (Fig. 4D). Then, upon a reversion back to the original environment (that is, back to positive selection), sequences tend to evolve up from the depths in the clusters of the F230S range, enabling access to BS-NEG-4. The environmental change can cause interrange bridges to form across gaps in the SWAG landscape that end in “gateway” sequences having F230S, allowing access to the sequence space of that range of peaks. Without this bridge, evolution of BS-NEG-4 is extremely unlikely. We conclude that evolution of BS-NEG-4 is extremely unlikely to occur under positive selection but can be accessed through a series of environmental changes that modulate the selective pressure on the gene. A deleterious gateway mutation allowing passage through a fitness valley via epistatic effects We next turned our attention to the role of F230S in the evolution of BSNEG-4. F230S is very deleterious to TEM-15 for cefotaxime resistance but is only mildly deleterious to ampicillin resistance in TEM-1 (36). The different effects of F230S on TEM-1 and TEM-15 might originate from a strong negative pairwise epistatic interaction with G238S, which Steinberg and Ostermeier Sci. Adv. 2016; 2 : e1500921

22 January 2016

lies at the opposite end of the same b-strand. Because we suspected F230S to be key in the evolution of this allele, we compared the epistatic and fitness effects of mutations of evolutionary intermediates containing F230S to intermediates lacking F230S (Fig. 7). Although alleles with F230S tend to have lower fitness, especially among intermediates with few mutations (Fig. 7A), several intermediates demonstrate higher fitness with the addition of F230S, including the global optimum among these sequences, which contains all but the Q90R and A184I mutations. We individually verified the MIC conferred by several of the fittest variants that contain F230S and observed up to a 32-fold improvement in resistance over that provided by GKTS (data S8). Seven of nine of the mutations of BS-NEG-4 are slightly beneficial as the sole mutation in TEM-15. However, the higher-order epistatic effect of combinations of these mutations in the absence of F230S is increasingly negative (Fig. 7B), consistent with previous observations that beneficial mutations usually interact with negative epistasis, creating diminishing returns (37–39). In contrast, higher-order epistatic interactions in intermediates containing F230S are almost always positive [252 of 256 (98%)] and become increasingly positive as mutations accumulate up to about four to five mutations, where the higher-order epistasis plateaus at e = +0.7. This epistasis value means that the average fitness of proteins having the 5 of 9

Downloaded from http://advances.sciencemag.org/ on March 30, 2016

C

TEM-15

RESEARCH ARTICLE A

B

Negative

Neutral

Positive

Fig. 5. Radial graph of all analyzed pathways from TEM-15 to BS-NEG-4. Each unit branch extending from the center indicates the gain of one of the nine amino acids from TEM-15 until the endpoint containing the full set of mutations is reached. (A) The line color indicates the log of the fitness change from the addition of the mutation. This graph pictorially shows that most pathways have steps with decreased fitness. (B) Same data as in (A), but with the color indicating the fitness relative to TEM-15 after the addition of the mutation. This graph pictorially shows that early steps in most pathways pass through mutants with fitness values less than that of TEM-15. Data is the same in A and B; only color coding is different.

pairwise epistasis between certain mutations (fig. S10). We suggest that F230S is more than just a key gateway mutation for traversing the fitness valley. F230S also serves to alleviate the negative interactions between beneficial mutations, creating the pathway to the final highactivity allele.

BS-NEG-4

w TEM-15

Fig. 6. Three-dimensional SWAG of all evolutionary pathways in which fitness increases monotonically with each additional mutation. Fitness is represented by height on the z axis. Color indicates clusters. Each node is sized by the likelihood of arrival, as measured by PageRank algorithm weighted by selection strength between nodes.

F230S mutation is five times higher than expected based on the effects of the individual mutations in TEM-15. We emphasize that this strong positive epistasis with F230S only occurs in the higher-order epistasis terms. The benefits of F230S are not so apparent when one examines the effects of F230S on pairwise epistasis among the other eight mutations, though F230S tends to have a slightly positive effect on the Steinberg and Ostermeier Sci. Adv. 2016; 2 : e1500921

22 January 2016

Conclusions Our results are important for applied directed evolution experiments, which almost exclusively use positive selection as a strategy. In our study, a selection regimen including negative selection generated the highest activity variants analyzed and outperformed strict positive selection. This result may appear counterintuitive to protein engineers but could prove useful in escaping plateaus in protein functionality. A few directed evolution studies have demonstrated the advantage of neutral drift in the evolution of activity toward alternate substrates (15–21), but our results are the first to show that a regimen involving negative selection can rival or even, perhaps, outperform other strategies, including strict positive selection and neutral drift. Although the neutral selection regimen performed here also accesses alternate fitness maxima, negative selection reached more distant sequences with superior properties on average. Thus, negative selection may be a more effective means of generating highly diverse and effective alternative starting points for directed evolution experiments. Our study lends further credence to models of rough fitness landscapes shaped by epistasis. SWAG enabled us to visualize and gain insights into potential evolutionary pathways using network analysis. Our technique creates statistical inference from evolutionary data that can be 6 of 9

Downloaded from http://advances.sciencemag.org/ on March 30, 2016

Fitness

RESEARCH ARTICLE A

F230S absent F230S present

8

Fitness of each mutant

of the isopropyl-b-D-thiogalactopyranoside (IPTG)–inducible tac promoter and confers spectinomycin resistance. Induction of b-lactamase expression was accomplished using 300 mM IPTG and the maintenance antibiotic spectinomycin (50 mg/ml).

9

7 6 5 4 3 2 1 0

0

2

4

6

8

10

Epistasis of each group of mutants

B

1.2

F230S absent F230S present

1 0.8 0.6 0.4 0.2 0 –0.2 –0.4 –0.6 –0.8

0

2

4

6

8

10

Number of mutations

Fig. 7. Protein fitness and epistasis along possible evolutionary pathways between TEM-15 and BS-NEG-4. (A) Protein fitness values of all possible evolutionary intermediates. (B) N-order epistasis between all mutations as a function of the number of mutations. Proteins containing F230S are noted in red, whereas all other proteins are shown in blue. Points are jittered for clarity. Larger points indicate the mean for each data set, with bars extending to the SD.

understood intuitively. We demonstrate SWAG visualization of a landscape with an evolutionary valley and multiple fitness peaks. SWAG landscape visualization should be broadly useful for depicting fitness landscapes and evolutionary pathways. Our experiments address the predictions of Wright (1) regarding the importance of environmental change in evolutionary processes. We present evidence that environmental changes that alter the fitness landscape can extend access to areas of sequence space that are difficult to explore in a constant environment. In our particular example, this access resulted in superior evolutionary outcomes. Whether our result is the exception or the rule will require testing of other genes and systems.

MATERIALS AND METHODS Plasmids All experiments used the expression plasmid pSkunk3 with a 9–base pair (bp) barcode indicating selection type downstream of the b-lactamase gene. The barcodes did not alter the level of antibiotic resistance conferred by the plasmid. The plasmid pSkunk3 contains the gene under the control Steinberg and Ostermeier Sci. Adv. 2016; 2 : e1500921

22 January 2016

Selection For bandpass experiments, plasmid libraries were isolated from the NEB 5-alpha cells and transformed into E. coli SNO301 in the presence of plasmid pTS40 (22). Bandpass selections were performed, as described previously (22, 36), on solid LB agar plates with spectinomycin (25 mg/ml), chloramphenicol (25 mg/ml), and 300 mM IPTG at 37°C for 24 hours [except for tetracycline (10 mg/ml) and various concentrations of cefotaxime]. Selections under non-bandpass conditions were performed in NEB 5-alpha cells on LB broth agar containing 300 mM IPTG, spectinomycin (25 mg/m), and various concentrations of cefotaxime at 37°C for 24 hours. Roughly 106 colony-forming units were plated for each selection experiment. Cell stocks were not previously exposed to b-lactam antibiotics before selection. At the end of the experiment, the 47 alleles were picked for sequencing based on colony size on plates, with antibiotics at the final selection level (that is, the largest 47 colonies were picked, first from the highest-concentration plate and then from the next highest-concentration plate, as necessary to reach 47). Screening and prevention of cross-contamination Between rounds, the 9-bp barcodes on each plasmid were used to monitor for contamination via a polymerase chain reaction (PCR)– based screen. Primers were designed to specifically amplify DNA containing each barcode to diagnose the presence of cross-contamination between the selection experiments. In the event of contamination, the desired fraction of the contaminated library was amplified with the appropriate barcodes and ligated into fresh vector to eliminate the contaminants. Deep sequencing of populations during evolution experiments Bacteria were recovered from plates with LB broth containing 10% glycerol and pelleted by centrifugation. Plasmid DNA was isolated from the pelleted cells. Amplicons were prepared for deep sequencing via PCR from each plasmid prep with the addition of a 3-bp barcode on either side of the amplicons indicating the round from which each amplicon was isolated. PacBio deep sequencing was performed on these amplicons using a three-pass circular consensus. Custom MATLAB scripts described previously (36) were used to align, analyze, and quantify reads and mutation composition. Reads were filtered for quality score, length, and quantity of insertions and deletions. Each read was then aligned to the reference gene to identify barcodes and mutations. 7 of 9

Downloaded from http://advances.sciencemag.org/ on March 30, 2016

Number of mutations

Mutagenesis PFunkel mutagenesis was performed as described (36, 40). Briefly, plasmids were transformed into E. coli CJ236 to prepare single-stranded template DNA containing uracil, using filamentous phage R408. Mutagenic oligonucleotides with degenerate codon NNN were targeted to every codon in the b-lactamase gene. Following mutagenesis, all reaction libraries were transformed into electrocompetent NEB 5-alpha cells (New England Biolabs) to prepare a naïve library for each round. All naïve libraries exceeded 106 transformants following mutagenesis.

RESEARCH ARTICLE

Phylogenetic analysis We used the Hamming distance between any two given sequences at the amino acid level to construct a distance matrix using custom scripts in MATLAB (fig. S1). Next, we applied a neighbor-joining algorithm in MATLAB to generate the tree from all combined sequences and cluster branches. The distance tree was visualized using FigTree (41). Selection-weighted attraction graphing SWAG network analysis was accomplished using Gephi (42) and calculated fitness values. Forces between nodes were weighted by the calculated relative selection strength between nodes (sr, the ratio of fitnesses wN and wN+1 between node N and node N + 1, which differ by one amino acid mutation) sr ¼

wN þ1 wN

ð1Þ

Nodes that differed by only one amino acid mutation were connected. The ForceAtlas2 algorithm (43) was used to generate graphs. Clustering was accomplished using the modularity algorithm weighted by selection strength. PageRank (44) was implemented and also weighted by selection strength. Steinberg and Ostermeier Sci. Adv. 2016; 2 : e1500921

22 January 2016

Epistasis calculations Pairwise epistasis between mutation A with fitness wA and mutation B with fitness wB was calculated as


eA•B ¼ log10

wAB wo wA wB



ð2Þ

in which wo is the wild-type fitness (45). Higher-order epistasis of order N with mutations i,j,k,…,n was calculated using Eq. 3 ei;j;k;…;n ¼ log10

wi;j;k;…;n wNo −1 Pni wi

! ð3Þ

Among three mutations A, B, and C, the amount by which a third mutation changes the epistasis is the same regardless of which mutation is considered the third mutation De3 ¼ eA•BjC − eA•B ¼ eA•CjB − eA•C ¼ eB•CjA − eB•C

ð4Þ

in which ei•j|k refers to the pairwise epistasis between i and j in the context of an allele containing mutation k. This calculation was used in fig. S10. MIC testing MIC testing was performed in NEB 5-alpha cells by growing 1 ml of each culture to be tested overnight in LB broth in the presence of spectinomycin (50 mg/ml). Saturated growths were diluted 1:100 and grown until above OD600 (optical density at 600 nm) = 0.3. All samples were then diluted back to OD600 = 0.003, and 1 ml was spotted onto LB broth agar plates with 300 mM IPTG and various concentrations of cefotaxime in twofold increments. Growth was evaluated after 24 hours of incubation at 37°C. The MIC was the lowest concentration of cefotaxime at which no growth was visible.

SUPPLEMENTARY MATERIALS Supplementary material for this article is available at http://advances.sciencemag.org/cgi/ content/full/2/1/e1500921/DC1 Fig. S1. Heat map of the matrix of Hamming similarity between unique proteins found from 47 of the most resistant alleles from each selection regimen following round 8. Fig. S2. Comparison of the differences in sequences obtained using the different selection regimen. Fig. S3. Structural mapping of mutations found in each selection scheme following round 8. Fig. S4. Heat map of the frequency of mutations observed during the negative evolution regimen, as analyzed by deep sequencing. Fig. S5. Stacked frequency distribution of mutations found by deep sequencing at each step in the negative evolution regimen as a function of codon position over eight rounds. Fig. S6. Stacked distribution of mutations found in the top 47 alleles from all selection strategies following round 8 as a function of codon position in b-lactamase. Fig. S7. Structural mapping of mutations found in the BS-NEG-4 allele, which confers a high resistance to cefotaxime. Fig. S8. Pathways between clusters found in SWAG analysis. Fig. S9. Two-dimensional SWAG landscape of pathways between TEM-15 and BS-NEG-4 with the nine single mutants indicated. Fig. S10. Change in pairwise epistasis (De3) between two mutations on the pathway from TEM15 to BS-NEG-4 in response to the addition of F230S. Data S1. Evolutionary history of selection schemes. Data S2. Sequences resulting from the positive selection regimen (TEM-1). Data S3. Sequences resulting from the positive selection regimen (TEM-15). Data S4. Sequences resulting from the neutral selection regimen. Data S5. Sequences resulting from the negative selection regimen. Data S6. Sequences resulting from the oscillating selection regimen. Data S7. Fitness and epistasis values for BS-NEG-4 pathways. Data S8. Additional MIC data.

8 of 9

Downloaded from http://advances.sciencemag.org/ on March 30, 2016

Deep sequencing to determine protein fitness Fitness values were determined essentially as previously described (36) and summarized here in a way that includes the minor differences. Bandpass selection selects for cells with activity at or near a particular level, as determined by the cefotaxime concentration. Thus, plating at a variety of different cefotaxime concentrations partitions the library into bins of relative fitness levels. Thus, by deep sequencing the plasmid DNA to determine the frequency of each allele on each plate, the fitness for each mutant can be determined. To determine fitness values, we performed selections over 11 plates containing doubling levels of cefotaxime (from 0.01 to 10.24 mg/ml) under bandpass conditions [300 mM IPTG, spectinomycin and chloramphenicol (25 mg/ml), and tetracycline (10 mg/ml) on LB broth agar at 37°C for 24 hours]. Selected libraries were swept with LB broth with 10% glycerol, cells were pelleted, and plasmid DNA was isolated. Amplicons were generated directly from each plasmid DNA prep. During PCR, 3-bp barcodes indicating the cefotaxime concentration on which each sequence was found were added at both ends. After amplifying the library from each plate, PacBio deep sequencing data were generated for each selection condition. We used custom Matlab scripts to align, analyze, and quantify reads and amino acid mutation composition. Reads were filtered for quality score (reported average quality score > 30, or average reported probability of error = 10−3), length (length of read less than 1100 bp but more than 930 bp), and quality of alignment to the reference gene (entirety of reference gene aligned to read) and the barcode (perfect match accepted only). We tabulated allele counts for each plate. In our previous work (36), the counts were used directly in the fitness calculation. Here, we adjusted the counts to reflect the desired proportion of sequencing reads from each plate as determined by the frequency of colonies appearing on that plate. Next, for each protein sequence, the plate with the highest counts was found, and the counts in that plate and in the surrounding two plates (five in total) were used to calculate the final fitness measurement as before (36). Only fitnesses for protein sequences with five or more counts were considered. Fitness values (w) were defined relative to the TEM-15 fitness.

RESEARCH ARTICLE REFERENCES AND NOTES

Steinberg and Ostermeier Sci. Adv. 2016; 2 : e1500921

22 January 2016

Acknowledgments: We thank the Johns Hopkins Deep Sequencing and Microarray Core Facility for performing the PacBio deep sequencing. Funding: This work was supported by the National Science Foundation (grants DEB-0950939 and CBET-1402101 to M.O.). Author contributions: B.S. and M.O. conceived the approach, analyzed the data, and wrote the paper. B.S. performed all experiments. Competing interests: M.O. is co-inventor on patent applications for the mutagenesis technique. This intellectual property has been licensed to Revolve Biotechnologies, for which M.O. serves as scientific advisor and has stock options. B.S. cofounded Revolve Biotechnologies and holds equity in the company. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors. Submitted 10 July 2015 Accepted 12 November 2015 Published 22 January 2016 10.1126/sciadv.1500921 Citation: B. Steinberg, M. Ostermeier, Environmental changes bridge evolutionary valleys. Sci. Adv. 2, e1500921 (2016).

9 of 9

Downloaded from http://advances.sciencemag.org/ on March 30, 2016

1. S. Wright, The roles of mutation, inbreeding, crossbreeding and selection in evolution. Proc. 6th Int. Congr. Genet. 1, 356–366 (1932). 2. S. Gavrilets, Fitness Landscapes and the Origin of Species (Princeton Univ. Press, Princeton, NJ, 2004). 3. J. A. G. M. de Visser, J. Krug, Empirical fitness landscapes and the predictability of evolution. Nat. Rev. Genet. 15, 480–490 (2014). 4. M. S. Breen, C. Kemena, P. K. Vlasov, C. Notredame, F. A. Kondrashov, Epistasis as the primary factor in molecular evolution. Nature 490, 535–538 (2012). 5. A. I. Podgornaia, M. T. Laub, Pervasive degeneracy and epistasis in a protein-protein interface. Science 347, 673–677 (2015). 6. F. J. Poelwijk, S. Tănase-Nicola, D. J. Kiviet, S. J. Tans, Reciprocal sign epistasis is a necessary condition for multi-peaked fitness landscapes. J. Theor. Biol. 272, 141–144 (2011). 7. S. J. Gould, Wonderful Life: The Burgess Shale and the Nature of History (W. W. Norton & Company, New York, 1990). 8. K. M. Taute, S. Gude, P. Nghe, S. J. Tans, Evolutionary constraints in variable environments, from proteins to networks. Trends Genet. 30, 192–198 (2014). 9. V. Mustonen, M. Lässig, From fitness landscapes to seascapes: Non-equilibrium dynamics of selection and adaptation. Trends Genet. 25, 111–119 (2009). 10. D. M. Weinreich, N. F. Delaney, M. A. DePristo, D. L. Hartl, Darwinian evolution can follow only very few mutational paths to fitter proteins. Science 312, 111–114 (2006). 11. D. M. Weinreich, L. Chao, Rapid evolutionary escape by large populations from local fitness peaks is likely in nature. Evolution 59, 1175–1182 (2005). 12. C. S. Gokhale, Y. Iwasa, M. A. Nowak, A. Traulsen, The pace of evolution across fitness valleys. J. Theor. Biol. 259, 613–620 (2009). 13. M. L. M. Salverda E. Dellus, F. A. Gorter, A. J. M. Debets, J. van der Oost, R. F. Hoekstra, D. S. Tawfik, J. A. G. M. de Visser, Initial mutations direct alternative pathways of protein evolution. PLOS Genet. 7, e1001321 (2011). 14. M. F. Schenk, S. Witte, M. L. M. Salverda, B. Koopmanschap, J. Krug, J. A. G. M. de Visser, Role of pleiotropy during adaptation of TEM‐1 b‐lactamase to two novel antibiotics. Evol. Appl. 8, 248–260 (2014). 15. J. D. Bloom, P. A. Romero, Z. Lu, F. H. Arnold, Neutral genetic drift can alter promiscuous protein functions, potentially aiding functional evolution. Biol. Direct 2, 17 (2007). 16. J. D. Bloom, F. H. Arnold, In the light of directed evolution: Pathways of adaptive protein evolution. Proc. Natl. Acad. Sci. U.S.A. 106, 9995–10000 (2009). 17. W. S. Smith, J. R. Hale, C. Neylon, Applying neutral drift to the directed molecular evolution of a b-glucuronidase into a b-galactosidase: Two different evolutionary pathways lead to the same variant. BMC Res. Notes 4, 138 (2011). 18. K. L. Petrie, G. F. Joyce, Limits of neutral drift: Lessons from the in vitro evolution of two ribozymes. J. Mol. Evol. 79, 75–90 (2014). 19. R. D. Gupta, D. S. Tawfik, Directed enzyme evolution via small and effective neutral drift libraries. Nat. Methods 5, 939–942 (2008). 20. S. Bershtein, K. Goldin, D. S. Tawfik, Intense neutral drifts yield robust and evolvable consensus proteins. J. Mol. Biol. 379, 1029–1044 (2008). 21. G. Amitai, R. D. Gupta, D. S. Tawfik, Latent evolutionary potentials under the neutral mutational drift of an enzyme. HFSP J. 1, 67–78 (2007). 22. T. Sohka, R. A. Heins, R. M. Phelan, J. M. Greisler, C. A. Townsend, M. Ostermeier, An externally tunable bacterial band-pass filter. Proc. Natl. Acad. Sci. U.S.A. 106, 10135–10140 (2009). 23. C. Poyart, P. Mugnier, G. Quesne, P. Berche, P. Trieu-Cuot, A novel extended-spectrum TEM-type b-lactamase (TEM-52) associated with decreased susceptibility to moxalactam in Klebsiella pneumoniae. Antimicrob. Agents Chemother. 42, 108–113 (1998). 24. M. Barlow, B. G. Hall, Predicting evolutionary potential: In vitro evolution accurately reproduces natural evolution of the TEM b-lactamase. Genetics 160, 823–832 (2002). 25. G. Kopsidas, R. K. Carman, E. L. Stutt, A. Raicevic, A. S. Roberts, M.-A. V. Siomos, N. Dobric, L. Pontes-Braz, G. Coia, RNA mutagenesis yields highly diverse mRNA libraries for in vitro protein evolution. BMC Biotechnol. 7, 18 (2007). 26. M. C. Orencia, J. S. Yoon, J. E. Ness, W. P. C. Stemmer, R. C. Stevens, Predicting the emergence of antibiotic resistance by directed evolution and structural analysis. Nat. Struct. Biol. 8, 238–242 (2001). 27. W. P. C. Stemmer, Rapid evolution of a protein in vitro by DNA shuffling. Nature 370, 389–391 (1994).

28. M. Zaccolo, E. Gherardi, The effect of high-frequency random mutagenesis on in vitro protein evolution: A study on TEM-1 b-lactamase. J. Mol. Biol. 285, 775–783 (1999). 29. M. K. Singh, B. N. Dominy, The evolution of cefotaximase activity in the TEM b-lactamase. J. Mol. Biol. 415, 205–220 (2012). 30. E. Firnberg, M. Ostermeier, The genetic code constrains yet facilitates Darwinian evolution. Nucleic Acids Res. 41, 7420–7428 (2013). 31. K. M. Brown, M. S. Costanzo, W. Xu, S. Roy, E. R. Lozovsky, D. L. Hartl, Compensatory mutations restore fitness during the evolution of dihydrofolate reductase. Mol. Biol. Evol. 27, 2682–2690 (2010). 32. N. Tokuriki, D. S. Tawfik, Stability effects of mutations and protein evolvability. Curr. Opin. Struct. Biol. 19, 596–604 (2009). 33. P. Hegarty, A. Martinsson, On the existence of accessible paths in various models of fitness landscapes. Ann. Appl. Probab. 24, 1375–1395 (2014). 34. R. P. Ambler, A. F. W. Coulson, J. M. Frère, J. M. Ghuysen, B. Joris, M. Forsman, R. C. Levesque, G. Tiraby, S. G. Waley, A standard numbering scheme for the class-A b-lactamases. Biochem. J. 276, 269–270 (1991). 35. Q. K. Thai, F. Bös, J. Pleiss, The lactamase engineering database: A critical survey of TEM sequences in public databases. BMC Genomics 10, 390 (2009). 36. E. Firnberg, J. W. Labonte, J. J. Gray, M. Ostermeier, A comprehensive, high-resolution map of a gene's fitness landscape. Mol. Biol. Evol. 31, 1581–1592 (2014). 37. M. F. Schenk, I. G. Szendro, M. L. M. Salverda, J. Krug, J. A. G. M. de Visser, Patterns of epistasis between beneficial mutations in an antibiotic resistance gene. Mol. Biol. Evol. 30, 1779–1787 (2013). 38. H.-H. Chou, H.-C. Chiu, N. F. Delaney, D. Segrè, C. J. Marx, Diminishing returns epistasis among beneficial mutations decelerates adaptation. Science 332, 1190–1192 (2011). 39. R. C. MacLean, G. G. Perron, A. Gardner, Diminishing returns from beneficial mutations and pervasive epistasis shape the fitness landscape for rifampicin resistance in Pseudomonas aeruginosa. Genetics 186, 1345–1354 (2010). 40. E. Firnberg, M. Ostermeier, PFunkel: Efficient, expansive, user-defined mutagenesis. PLOS One 7, e52031 (2012). 41. A. Rambaut, A. Drummond, FigTree (tree figure drawing tool 2012); http://tree.bio.ed.ac. uk/software/figtree/ 42. M. Bastian, S. Heymann, M. Jacomy, Gephi: An open source software for exploring and manipulating networks, International AAAI Conference on Weblogs and Social Media, San Jose, CA, 2009. 43. M. Jacomy, T. Venturini, S. Heymann, M. Bastian, ForceAtlas2, a continuous graph layout algorithm for handy network visualization designed for the Gephi software. PLOS One 9, e98679 (2014). 44. L. Page, S. Brin, R. Motwani, T. Winograd, The PageRank Citation Ranking: Bringing Order to the Web (Stanford InfoLab, Stanford, 1999). 45. B. Østman, A. Hintze, C. Adami, Impact of epistasis and pleiotropy on evolutionary adaptation. Proc. Biol. Sci. 279, 247–256 (2012).

Environmental changes bridge evolutionary valleys Barrett Steinberg and Marc Ostermeier (January 22, 2016) Sci Adv 2016, 2:. doi: 10.1126/sciadv.1500921

This article is publisher under a Creative Commons license. The specific license under which this article is published is noted on the first page. For articles published under CC BY licenses, you may freely distribute, adapt, or reuse the article, including for commercial purposes, provided you give proper attribution.

The following resources related to this article are available online at http://advances.sciencemag.org. (This information is current as of March 30, 2016): Updated information and services, including high-resolution figures, can be found in the online version of this article at: http://advances.sciencemag.org/content/2/1/e1500921.full Supporting Online Material can be found at: http://advances.sciencemag.org/content/suppl/2016/01/19/2.1.e1500921.DC1 This article cites 40 articles,14 of which you can be accessed free: http://advances.sciencemag.org/content/2/1/e1500921#BIBL

Science Advances (ISSN 2375-2548) publishes new articles weekly. The journal is published by the American Association for the Advancement of Science (AAAS), 1200 New York Avenue NW, Washington, DC 20005. Copyright is held by the Authors unless stated otherwise. AAAS is the exclusive licensee. The title Science Advances is a registered trademark of AAAS

Downloaded from http://advances.sciencemag.org/ on March 30, 2016

For articles published under CC BY-NC licenses, you may distribute, adapt, or reuse the article for non-commerical purposes. Commercial use requires prior permission from the American Association for the Advancement of Science (AAAS). You may request permission by clicking here.