Systematic mappings of the effects of protein mutations are becoming increasingly popular. Unexpectedly, these experiments often find that proteins are tolerant to most amino acid substitutions, including substitutions in positions that are highly conserved in nature. To obtain a more realistic distribution of the effects of protein mutations, we applied a laboratory drift comprising 17 rounds of random mutagenesis and selection of M. HaeIII, a DNA methyltransferase. During this drift, multiple mutations gradually accumulated. Deep sequencing of the drifted gene ensembles allowed determination of the relative effects of all possible single nucleotide mutations. Despite being averaged across many different genetic backgrounds, about 67% of all nonsynonymous, missense mutations were evidently deleterious, and an additional 16% were likely to be deleterious. In the early generations, the frequency of most deleterious mutations remained high. However, by the 17th generation, their frequency was consistently reduced, and those remaining were accepted alongside compensatory mutations. The tolerance to mutations measured in this laboratory drift correlated with sequence exchanges seen in M. HaeIII’s natural orthologs. The biophysical constraints dictating purging in nature and in this laboratory drift also seemed to overlap. Our experiment therefore provides an improved method for measuring the effects of protein mutations that more closely replicates the natural evolutionary forces, and thereby a more realistic view of the mutational space of proteins.
The ability to reliably measure and predict the effects of amino acid mutations in proteins is of fundamental importance to protein engineering and design and for understanding protein evolution and human genetic variation. Data regarding the effects of individual mutations originate from two major sources: sequence analysis of natural proteins, and laboratory experiments. Phylogenetic analyses enable insights regarding how protein sequences diverge [1,2] and what dictates the purging of mutations [3,4]. Protein phylogenies also allow us to predict whether a given mutation might be deleterious or neutral, assuming that the fitness effects of mutations correlate with their occurrence in orthologous sequences (reviewed in [5–7]). The algorithms are relatively accurate in predicting disease-causing mutations [8,9]. However, many nonsynonymous SNPs predicted to have a deleterious effect are not clearly associated with a disease phenotype [10], either because they are rare [11] or because a deleterious effect in a single gene often results in no phenotype at the organismal level [12]. Indeed, the effects of mutations at the isolated protein and organismal levels do not necessarily overlap. Predictors may also fail in assigning deleterious effects to mutations in highly conserved sites that when mutated experimentally appear to be neutral [13]. Exhaustive datasets listing the effects of all mutations within a given gene/protein, independently of organismal effects, would therefore greatly improve prediction [14]. Systematic experimental mappings of the effects of mutations within one given gene/protein are therefore crucial for understanding protein evolution, as well as an attractive resource for improving predictions [15,16] and for refining protein design algorithms [17]. Experiments that systematically map the effects of mutations in a given protein are generally conducted through either saturation mutagenesis, using NNS codons to diversify individual sites [18–21], or random mutagenesis along the entire gene using error-prone replication [22]. In both cases, the diversified gene repertoires are subjected to selection that purges deleterious mutations, and then sequenced to identify which mutations are tolerated. Recently, advanced gene synthesis technologies and deep sequencing have yielded exhaustive mappings (for examples see [23–38]; reviewed in [16]). However, although deep mutational scanning provides a powerful means of studying protein structure-function, there remain challenges that are yet to be tackled [16]. Foremost, the relevance of the results of laboratory experimental mappings our understanding of natural protein evolution may be limited. Specifically, there is a disagreement between the trends indicated by experimental mappings versus natural protein diversity. In silico analyses of natural protein diversities suggest that the vast majority of mutations are deleterious [39–43]. Given enough drift, mutations at other sites enable the acceptance of certain deleterious substitutions. But at a background of a given sequence, most substitutions would result in the loss of configuration stability and/or function [41,43–45]. Experimental mappings, however, mostly portray a different picture–the majority of mutations are tolerated (for examples, see [25–27,33,36,46,47]). Accordingly, a poor correlation between the acceptance of mutations in the laboratory and the occurrence of the same exchanges in natural orthologs of the studied protein has been noted [25–27,33,36]. For example, positions that are 75–90% conserved in Hsp90 tolerated a range of amino acids some of which are not seen in any ortholog [25]. However, at lower expression levels, these mutations did reveal deleterious fitness effects, thus indicating that the sensitivity of the experimental system is a key parameter [48]. The comparison of results from different experimental mappings is also problematic. The experiments not only address different proteins, but also apply different mutagenesis strategies and methods of determining the effects of mutations. In some cases the measured effects of mutations relate to growth of the host organism (e. g. , antibiotics resistance) and in others to the biochemical function of the targeted protein in isolation (e. g. levels of fluorescence, or of DNA methylation, as applied here). Nonetheless, the disagreement between tolerance in the laboratory and occurrence amongst natural sequence raises several questions. Does the absence of a given exchange within natural orthologs indicate its deleterious fitness effect, or does the sparse and sporadic sampling of natural sequences prevents reliable prediction? Do laboratory experimental setups adequately reproduce the constraints that shape protein sequences in nature, or do tolerance or acceptance of a mutation in the laboratory have limited relevance to the evolutionary history/future of a protein. Finally, obtaining a realistic distribution of the fitness effects (DFE) of protein mutations remains a worthy goal [36,46,49,50]. To address the above questions, and obtain a more realistic distribution of fitness effects of protein mutations, we have set up a laboratory system that better mimics the manner by which protein sequences diverge in nature. To this end, we performed 17 iterative rounds of random mutagenesis and purifying selection. This laboratory experiment does not address crucial elements of natural drifts (mutation rates, population sizes, and organismal fitness demands). It does, nonetheless, mimic the process of prolonged accumulation of mutations under purifying selection to maintain the protein’s structural and functional integrity (hence the term ‘neutral drift’). As a model, we used a bacterial DNA methyltransferase, M. HaeIII, which can be readily placed under purifying selection in the laboratory. At different rounds along this prolonged drift, the ensembles of gene variants that survived the purifying selection were subjected to deep sequencing. The naïve, unselected mutational repertoire was similarly sequenced. This enabled us to determine the frequency of occurrence, and hence the relative fitness effects of all single nucleotide mutations in M. HaeIII. As described in the following pages, our results differed from those of other experimental mappings in several key respects.
M. HaeIII is a DNA methyltransferase isolated from Haemophilus aegyptius. Being part of the bacterial restriction-modification system, this enzyme selectively methylates GGCC DNA sequences, and thereby protects DNA from digestion by the cognate endonuclease, HaeIII. Sequence specific methylation-restriction offers a facile way of performing laboratory evolution. As described in earlier works [51,52] (S1 Fig), M. HaeIII' s open reading frame was randomly mutated using PCR with an error-prone polymerase, cloned into an expression plasmid and transformed to E. coli. In each bacterium, the encoding plasmid is methylated, or not, depending on whether the M. HaeIII gene variant it encodes is properly folded and functional. Following ‘expression’ of the plasmid encoded M. HaeIII variants in individual bacteria, the plasmids were pooled and digested with HaeIII, thus selecting for M. HaeIII’s native specificity [51,52]). The starting point for these experiments was a variant of M. HaeIII optimized for soluble and functional expression in E. coli. This variant carried four consensus substitutions replacing the amino acid in M. HaeIII with the one that dominates in all M. HaeIII’s orthologs [50]. For the sake of simplicity we refer to this variant as wild-type M. HaeIII. These consensus substitutions are likely to have a stabilizing, compensatory effect, and spontaneously accumulate in accelerated, laboratory drifts. They may thus allow a larger variety of deleterious mutations to be accepted, especially during the first rounds of mutagenesis [53]. In Haemophilus aegyptius, M. HaeIII is under a strong and constitutive selection pressure imposed by the presence of the cognate restriction enzyme HaeIII–a DNase that would cause chromosomal breaks unless the genome is methylated at all HaeIII sites. The HaeIII restriction-modification system is naturally encoded by single copy chromosomal genes [54,55]. In our experimental system, M. HaeIII was encoded by a multi-copy plasmid (~400 copies per cell). To avoid unrealistic enzyme doses, expression was driven from a tightly controlled promoter with no induction. Although M. HaeIII’s levels in Haemophilus aegyptius are unknown, its expression level in the E. coli cells of our experimental setup is extremely low (a similar plasmid showed no detectible GFP signal when inducer levels were ≤20 μg/ml [56], and we used no inducer). This basal expression level was nonetheless sufficient to enable wild-type M. HaeIII to methylate all GGCC sites, not only in the encoding plasmid, but also within the E. coli host' s chromosome, as is the case with natural methyltransferases [51]. M. HaeIII underwent 17 rounds of random mutagenesis, at an average mutational rate of 2. 2±1. 6 nucleotide mutations per gene per generation followed by purifying selection (i. e. , digestion of the encoding plasmids with HaeIII nuclease). To avoid false positives due to mutations in GGCC sites, the applied M. HaeIII’s coding sequence (ORF) contained no GGCC sites whist the encoding plasmid contained 14 such sites including three sites within the antibiotic selection marker. Each round, selection was repeated three times (i. e. , repeated isolation of plasmid DNA pool from the grown bacteria, digestion with HaeIII, and transformation into E. coli). Subsequently, the drifted M. HaeIII’s ORF was mutagenized and recloned into a fresh plasmid for the next round of selection. We ensured the same level of selection pressure and the absence of bottlenecks throughout: ≥105 independent transformants were passed to the next round (effective population size, Ne > 105). The drifting M. HaeIII thus met the conditions that essentially eliminate the possibility of mutations fixing by chance (1/Ne <10−5). Mutations that were enriched are therefore likely to have provided a selective advantage, most typically, as shown below, a compensatory effect. By the 17th round, the drifted genes carried on average of 18±1. 6 mutations per gene in total, and 9. 6±0. 7 nonsynonymous mutations per gene (determined in parallel by deep sequencing and conventional Sanger sequencing of the full length ORFs of randomly chosen genes; see ‘Dynamics of the laboratory drift‘ below). The mutational spectrum of the unselected, naïve gene pool (dubbed G0), and of the pools after three (G3), seven (G7) and seventeen (G17) rounds of selection, were analyzed by Illumina high-throughput sequencing. Mutations were identified using a script that aligned all codon triplets to the reference gene, wild-type M. HaeIII (S1 File). The background frequency that stems from the Illumina sequencing PCRs and sequencing errors was determined using the sequencing data for the region upstream of the randomly mutated ORF of M. HaeIII (the N-terminal fused His tag that was part of the encoding vector and was hence not subjected to mutagenesis). In this manner, the frequencies of all single, double and triple nucleotide mutations were determined by the fraction of sequence contigs that carry a mutation out of all contigs that covered the respective position. The amino acid mutational frequencies were subsequently determined by summing up the frequencies of all codon triplets that yield a given amino acid (See S1, S2 and S3 Files). In theory, M. HaeIII’s sequence space includes 6,580 possible amino acid mutations; i. e. , 329 positions, each mutated to all other 19 different amino acids or to a stop codon (Table 1). However, the immediate mutational space originates from single nucleotide mutations. Subsequent nucleotide mutations within the same codon were found (dubbed double and triple mutations; Table 1), but at very low frequencies (in G17—an average of 0. 052% and 0. 003% of nonsynonymous double and triple mutations). However, these double and triple mutations only appeared at later stages, and after many other positions had changed due to single nucleotide mutations. Single nucleotide missense mutations also dominate polymorphism, and thus, our analysis focused on their effects. We thus examined all 1,957 possible missense single point mutations, namely all amino acid exchanges accessible by single nucleotide mutations; Table 1). The effects of stop codons were also examined as described in ‘Tolerance of nonsense mutations’. All possible single nucleotide mutations were detected in the unselected G0 library at the raw data level– 329X9 = 2,961 possible single nucleotide mutated codons, that in turn yield 1,957 possible single nucleotide amino acid mutations (raw data provided as S2 File). However, 77 mutations at G0 were observed at lower than background frequencies. That a mutation is observed under what we defined as the background frequency is not necessarily an indication that it did not occur. The background frequency was derived from averaging the frequencies for relatively few positions compared to the measured ones (20 positions versus 329) and thus it is conceivable that a small fraction of the latter (<1%) will deviate from this average. Indeed, out of the 77 mutations, 35 were detected in the later, selected rounds, and some were even enriched. The remaining 42 mutations were also observed with under background frequencies in later rounds. We therefore assume that they were strongly purged and assigned them as eliminated mutations. Overall, our analysis related to the complete set of 1,957 single nucleotide missense mutations with >98% (1,915/1,957) of these being covered with complete confidence. The spectrum of mutations covered by our experiment was dictated by the genetic code, M. HaeIII’s DNA sequence, and by the nucleotide substitution matrix that underlined our mutagenesis protocol. Although we used an engineered, error-prone DNA polymerase, the obtained spectrum of mutations was similar to that naturally observed in E. coli. Specifically, a transition/transversion ratio of ~1. 3 was observed in our naïve repertoire (G0) similar to what has been observed in the comparison of closely related E. coli genomes (0. 91 or 1. 3, [57,58], S1 Table). The variability in mutation frequencies along M. HaeIII positions in the unselected G0 library was relatively high (1. 07 ± 0. 24% mutations/position). Thus, mutation frequencies varied not only by the type of base substitution (e. g. transitions, transversion; S2A Fig), but also according to the position of the mutated base along M. HaeIII’s gene. To verify that this variability is not the outcome of limited sampling in G0 (the naïve repertoire that underwent only one round of mutagenesis) we compared the frequencies of synonymous mutations in the unselected library, G0, and in the selected one, G3. As expected, synonymous mutations were under relatively weak selection (detailed below) and thus their frequencies, certainly within the early rounds, largely reflect the rate of mutagenesis. Indeed, the frequencies of synonymous mutations in G0 and in G3 were highly correlated (R = 0. 9, S2B and S2E Fig). By G17, the correlation was still significant although weaker indicating some degree of selection on synonymous mutations (R = 0. 6; S2C Fig). The observed frequencies in the unselected library, f (G0), therefore appear to provide a reliable measure for the positional rates of occurrence of mutations in all 17 mutagenesis steps of the drift. However, for the 77 mutations (out of 1,957) with lower than background frequencies in G0 (Table 1), the rate of the occurrence, f (G0), was based on the base substitution table derived from G0 (S2A Fig). Mutations were retained, purged or enriched in each round of our experiment. The change in frequency along the drift therefore reflects the effects of selection per each mutation or, as defined here, their relative fitness effects (Wrel). The frequency of a given mutation in a given round (f (Gn) ) is dictated by its relative fitness (Wrel), and relates to the frequency of this mutation in the previous round, f (Gn-1) plus the frequency of re-occurrence at round n. For example, the frequencies of neutral mutations (Wrel = 1) are essentially equal to their cumulative rate of occurrence (f (Gn) ~ n f (G0) ). Conversely, the frequencies of deleterious mutations (Wrel < 1) decrease from round to another, in an exponential manner, and their observed frequency is lower than expected from their rate of occurrence (f (Gn) < n f (G0) ). The opposite applies for beneficial mutations (Wrel > 1). However, since the genes in our drifting ensembles contained multiple mutations, and the applied sequencing approach does not reveal the specific mutational composition of individual genes, the Wrel values measured here relate to the effect of a given mutation at the background of many different genetic compositions. For better and for worse, the measured Wrel values therefore represent an average that ignores epistatic interactions between mutations. This averaging has obvious drawbacks, and may cause biases due to hitchhiking and clonal interference (e. g. a highly deleterious mutation would result in every other mutation on the same gene having a low Wrel) However, under our experimental setup, new mutations are reintroduced in each round of mutagenesis, allowing multiple resampling of the effects of each given mutation at the background of many different mutations, and thus reducing the probability of hitchhiking. Indeed, as indicated below, there are clear indications that hitchhiking and clonal interference did not bias the observed Wrel values. We also note that, in general, allele frequencies, and thereby fitness effects of mutations, are measured in populations comprising individuals with different genetic backgrounds with certain caveats [59–61]. Foremost, the number of sequenced alleles needs to be in the order of thousands [59]–a demand that is amply met in our experiment. Thus, if a mutation is on average purged (Wrel <<1), we can conclude that it has deleterious effects on M. HaeIII’s structure and/or function independently of the specific genetic background. We used the following model to derive the relative fitness effect, Wrel, from the mutational frequencies observed in the selected versus unselected libraries (see Methods for details). Following the first round: f (G1) =f (G0) ⋅Wrel, For the subsequent rounds, the observed mutational frequency, f (Gn), is derived from the mutations inherited from the previous round f (Gn-1) plus the mutations newly incorporated in this round. The latter corresponds to the frequency of this mutation in the naïve, unselected ensemble, f (G0), as discussed above: f (Gn) =[f (Gn−1) +f (G0) ]⋅Wrel (1) Eq (1) corresponds to a geometrical series that has no closed solution. Thus, to derive the Wrel values of each mutation, we calculated the expected frequency ratio (f (Gn) f (G0) ) for a series of discrete Wrel values from absolutely deleterious (Wrel = 0) to highly beneficial (Wrel = 3. 5; see Methods and S3 Fig). In this manner, each of the 1,957 amino acid mutations measured by the deep sequencing (each derived from the respective single nucleotide mutation; Table 1), were assigned a Wrel value. We used the variability in the relative fitness effects of synonymous mutations and nonsense mutations to categorize the effects of nonsynonymous mutations [36]. The distribution of synonymous mutations was consistent with their low impact on fitness relative to nonsynonymous mutations (Fig 1A, ‘Syn’). The average Wrel value, and standard deviation, for synonymous mutations were found to be 0. 82±0. 12 for G3,0. 84±0. 15 for G7, and 0. 91±0. 1 for G17 (Table 2). Given our hypothesis that other works overestimated the tolerance of mutations, we preferred to under- rather than over-estimate the fraction of deleterious mutations. Accordingly, for the assignment of a deleterious fitness effect, we chose a conservative threshold of two standard deviations under the mean of the relative fitness effect of synonymous mutations (X¯−2SD). The X¯−2SD values obtained were 0. 58,0. 55 and 0. 72 for G3, G7 and G17, respectively (Table 2) yielding an average of 0. 62 for all 3 ensembles. We therefore used Wrel ≤ 0. 6 as the threshold for indicating purging and consequently a deleterious fitness effect of a mutation. However, since, as detailed below, mutations with Wrel values of 0. 6–0. 8 were found to be systematically purged as the drift progressed, suggesting that in effect they are not neutral. This Wrel range was therefore classified as ‘nearly-neutral’ (Fig 1A). The X¯+2SD threshold was similarly applied for categorizing beneficial mutations, thus stetting the threshold for beneficial mutations as Wrel > 1. 1. Within this threshold, only 4 out of 321 synonymous mutations were defined as beneficial relative to 33 nonsynonymous mutations. Indeed, within the 0. 6–1. 1 Wrel range defined here as neutral, >93% of the synonymous mutations observed in the three selected ensembles (G3, G7, G17) were assigned as neutral (Table 2). The potential deleterious or beneficial effects of the remaining 7% were not analyzed here. The selection acting on synonymous mutations may, amongst other factors, relate to different codon usage in E. coli. Overall, given the applied thresholds, the likelihood of misassignment of neutral mutations as deleterious or beneficial was < 4% (Table 2). The Wrel threshold for defining ‘highly deleterious’ mutations was derived from the distributions of nonsense mutations that are, beyond doubt, deleterious (see also ‘Tolerance of nonsense mutations‘ below). The average Wrel value, and standard deviation, for nonsense mutations were found to be 0. 042±0. 15 for G3,0. 028±0. 13 for G7, and 0. 020±0. 11 for G17 (Table 2). Thus, a threshold of X¯+2SD, i. e. , Wrel ≤ 0. 3, was chosen for categorizing highly deleterious mutations. In summary, nonsynonymous mutations were categorized as ‘Deleterious’ if their Wrel values were ≤ 0. 6, and ‘Highly deleterious’ if Wrel ≤ 0. 3 (including eliminated mutations, Wrel = 0, i. e. , when the net frequency of a mutation was zero). Mutations were assigned as ‘Nearly-neutral’ if their frequencies in the selected populations were in the range of Wrel = 0. 6–0. 8 (X¯−2SD of the distribution of synonymous mutations) and ‘Neutral’ in the range of Wrel = 0. 8–1. 1. Finally, enrichment in the selected repertoires (Wrel >1. 1, X¯+2SD of the distribution of synonymous mutations) indicated a ‘Beneficial’ fitness effect. As can be seen in Fig 1A, the distribution of relative fitness effects of synonymous mutations centered near neutrality (Wrel ~1; see Table 2 for mean and standard deviation). In contrast, the distribution of the nonsynonymous mutations encompasses primarily deleterious mutations. Overall, ~67% out of all the possible nonsynonymous single nucleotide mutations (~1,310/1,957) were found to be deleterious, even within the conservative threshold of Wrel of ≤ 0. 6 (Table 2, ‘nonSyn’). Removal of the 42 mutations that were observed below background frequency in all repertoires, and assigned as eliminated, has almost no impact on the fraction (~1,268/1,915 = ~66%). Note that Fig 1A shows the derived fitness effects of all possible single nucleotide mutations regardless of their frequencies in the drifting populations. Indeed, a similar fraction was assigned as deleterious in all the three selected libraries (Fig 1A, ‘nonSyn’). The effect that selection had on purging deleterious mutations is clearly seen in the distribution of their frequencies (Fig 1B). This distribution shifted during the drift: from ~15% of all mutations denoted as deleterious in G3, to only ~3% in G17 (Fig 1C). Further, by our conservatively chosen threshold, mutations with assigned Wrel from 0. 6 to 0. 8 were considered as ‘Nearly-neutral’. However, mutations within this Wrel range were systematically purged throughout the drift, from ~22% in G3 to ~7% in G17 (Fig 1C), indicating small yet consistent deleterious effects. Further, as discussed below, these mutations were accepted at the background of beneficial mutations, most likely owing to their compensatory effect (as discussed in the section below). If all mutations with Wrel ≤ 0. 8 are considered, then ~83% of all possible mutations in M. HaeIII have a deleterious effect. In agreement with the reduction in the frequency of deleterious mutations, the total frequency of beneficial mutations (Wrel > 1. 1) increased consistently, from ~2% in G3 to ~10% in G17. Consistent with the distribution of Wrel values being the same along the drift (Fig 1A), we also observed that the Wrel values per given mutation remain largely the same along the drift, i. e. , when derived from the sequenced frequencies in G3, G7 or G17 (S3B and S3C Fig). This was despite the fact that the average number of mutations per gene increased form 2. 4 in G3 to 9. 6 in G17. It therefore seems that hitchhiking and/or clonal interference did not significantly bias our data, and that the derived Wrel values that average the effect of a given mutation over different genetic background are relevant. The discrepancy between our results and the results of other experimental mappings is likely due to differences between measuring the effect of a single, or at most a few mutations, to measuring the effect of accumulated mutations over a long mutational drift. To support this hypothesis, we examined the rate of accumulation of mutations along the various rounds of the drift. To this end, in addition to deep-sequencing of the gene ensembles of 3 rounds (G3, G17 and G17), we also performed conventional Sanger sequencing of the full length ORFs of randomly picked variants from each round. As expected, the total number of mutations per gene (Nt) increased linearly throughout the drift, certainly up to the 15th round, at the average of 1. 16 mutations/gene/round (Fig 2A). However, as the drift progressed, the intensity of purged nonsynonymous mutations increased as indicated by the ratio of nonsynonymous to synonymous mutations (Na/Ns; whereby Na denotes the average frequency of nonsynonymous mutations and Ns the average frequency of synonymous mutations; Fig 2B). Specifically, the first round (G1) exhibited a Na/Ns ratio of 2. 6, only mildly lower than 3. 2 –the ratio in G0, the unselected repertoire. However, by the 5th round, the Na/Ns ratio dropped to a value of ~1. By the 14th round, the accumulation of nonsynonymous mutations (Na) had slowed down, in addition to a slowdown in the total accumulation of mutations (Nt; Fig 2A). That the tolerance to mutations decreased as the drift progressed is also reflected in the continuous decline in frequency of deleterious mutations along the drift. About a third of all possible mutations were eliminated by the 3rd round (Wrel = 0, ‘Eliminated’) whereas deleterious mutations (Wrel = 0. 01–0. 6) were observed in the drifting ensembles throughout the 17 rounds. However, their frequency was small and remained constant throughout (~0. 6 deleterious mutations/gene; Fig 2C). Thus, as mutations gradually accumulated, the relative frequency of deleterious mutations became increasingly lower–in G17 their fraction became 0. 06 (0. 6 out of the ~9. 6 nonsynonymous mutations/gene) relative to 0. 26 in G3 (0. 6 out of the ~2. 4 nonsynonymous mutations/gene). In effect, the majority of mutations that did accumulate beyond G3 were neutral (Figs 1B and 2C). This finding also indicates that hitchhiking and/or clonal interference does not significantly bias our data. In addition to the accumulation of the neutral mutations (Wrel = 0. 8–1. 1) beyond G3, the later rounds were accompanied by the enrichment of beneficial mutations (Fig 1C). The beneficial mutations (Wrel > 1. 1) are likely to be compensatory mutations, increasing the global stability of M. HaeIII, or locally interacting with a specific deleterious mutation. The applied sequencing method does not reveal the specific mutational composition of individual genes, and thus, there is no way of detecting enriched, beneficial mutations that have a specific, local compensatory effect. However, as previously shown [53], mutations that were enriched in a prolonged neutral drift were experimentally confirmed to have global, stabilizing effects that compensate for a wide range of deleterious destabilizing mutations. The global compensatory effect can also be deduced from the identification of most enriched mutations as consensus mutations (S2 Table; see also Ref. [53]). Further, under selection for the acquisition of five different new DNA target specificities [51], the same mutations were rapidly fixed in all the evolved lines irrespective of which new specificity was selected (S2 Table). Compensatory mutations are essential for the acquisition of new functions because mutations that confer new functions tend to severely undermine protein stability [62–64]. By G17, each gene carried, on average, 1. 99 beneficial mutations relative to 0. 08 in G3 (Fig 2C). Conversely, the fraction of enriched, beneficial mutations in the drifting genes (out of all nonsynonymous mutations) became 0. 21 in G17 relative to 0. 03 in G3 (Fig 1B and 1C). Thus, not only was the mutational tolerance limited beyond G3, but also, the acceptance of nearly-neutral mutations (Wrel = 0. 61–0. 8, Fig 1) was dependent on the co-accumulation of compensatory mutations. The ability to tolerate highly deleterious mutations at the onset of the drift, but not once mutations further accumulate, is also vividly exemplified by the tolerance of nonsense mutations–mutations leading to stop codons (S4 Fig) or frameshifing insertions/deletions (InDels). The occurrence and tolerance of InDels in the selected G17 M. HaeIII library has been described [52], indicating that certain nonsense mutations were tolerated to some degree due to translational slippage that results in a correctly translated protein despite a frame-shifted gene. However, the levels of full length, functional proteins translated from frame-shifted genes is obviously much lower than for wild-type, in some cases as little as 1% [52]. The second form of nonsense mutations are stop codons. At the onset of the drift, at least one stop codon mutation, in position 176, was moderately tolerated (Wrel in G3 = 0. 74). Stop codons in other positions were also found, although with lower Wrel values (S4A Fig). However, once other mutations that reduce protein dose and/or function accumulated, nonsense mutations were almost entirely purged (e. g. Wrel for position 176 in G7 = 0. 47, and in G17 = 0. 24). By the later rounds, stop codons were found almost only after position 324– a region that is not under functional selection (S4 Fig). The stronger purging effect of nonsense mutations in later rounds was also manifested in an increasing fraction of nonsense mutations being assigned a Wrel value of 0 (Figs 1, ‘nonSense’, and S4B). Several laboratory mutational tolerance experiments indicated the acceptance of mutations in positions that are highly conserved in natural orthologs [21,25–27]. We therefore examined to what degree M. HaeIII' s orthologs predict acceptance in our experiment; namely, do the measured relative fitness effects of mutations (Wrel) correlate with the degree of divergence of the corresponding position in M. HaeIII’s natural orthologs? To address this question we first compared the experimental Wrel values to the natural evolutionary rates of the respective positions (Fig 3A). We used Rate4Site whereby the calculated rates relate to the degree of physico-chemical change exerted by sequence exchanges, and to the phylogenetic distances [65]. We found that positions that exhibit slow evolutionary rates (i. e. , are highly conserved in nature; log2μ ≤ −2, or μ ≤ 0. 25 [44]) show no, or low acceptance to mutations in the laboratory drift. Conversely, positions with high evolutionary rates tend to show high experimental tolerance to mutations. This trend is seen along the primary sequence (Fig 3A) as well as in the 3-dimensional structure (Fig 3B and 3C). Given that we have mapped the effects of all possible single nucleotide mutations in M. HaeIII, we further examined how well their effects could be predicted from an alignment of orthologous sequences. There are many ways of predicting the effects of mutations from multiple sequence alignments. Certain biases are inevitable; foremost, prediction is highly dependent on sequence sampling–the number and the phylogenetic distribution of available sequences that are evolutionary related to the protein in question. Other biases relate to phylogenetic relatedness of the orthologs to the reference sequence and the manner by which the degree of divergence is calculated. A meaningful measure uses profile scores (position-specific scoring matrices) that take into account not only the frequency of sequences in which a given position varies, but also the physico-chemical nature of exchanges (reviewed in [5,66]). Given the epistatic nature of sequence evolution, tolerance of mutations is largely not a matter of' if' but of' when' –namely, given enough drift, exchanges in even the most conserved sites may be tolerated [2,45,67]. PROVEAN (Protein Variation Effect Analyzer, http: //provean. jcvi. org) is a predictor that takes into account phylogenetic distances [68]. Thus, the PROVEAN score function considers the physiochemical impact of amino acid exchanges alongside the evolutionary distance between the reference protein and the homolog (s) in which a given exchange is observed. We submitted to PROVEAN the reference sequence (the wild-type M. HaeIII gene) and a multiple sequence alignment of 105 orthologs (S5 Fig). The program computed a score predicting how deleterious each possible amino acid exchange in M. HaeIII might be. The algorithm’s default thresholds are: scores ≤ -2. 5 are predicted as deleterious, and scores > -2. 5 as neutral [68]. We then compared the predicted PROVEAN score to the measured Wrel values for the 1,957 single nucleotide nonsynonymous mutations. Overall, a clear-cut trend is seen (Fig 3D) –mutations found to be deleterious in the laboratory drift (Wrel ≤ 0. 6) tend to show low PROVEAN scores (≤ -2. 5), whereas the accepted ones show high scores (> -2. 5). From PROVEAN’s point of view as a predictor of deleterious mutations, true positives occurred at a rate of 83. 3% (S3 Table). Namely, out of the 1,234 mutations that were evidently deleterious in the laboratory drift (Wrel ≤ 0. 6), 1,028 were correctly categorized by PROVEAN as deleterious (score ≤ -2. 5). True negatives–mutations predicted by PROVEAN as neutral and found to be so in the drift, occurred at a rate of 63. 5% (459 out of the 723 accepted mutations in the laboratory drift, Wrel > 0. 6, were scored with PROVEAN values of > -2. 5). When excluding mutations with borderline effects (mutations categorized as nearly-neutral, with Wrel values 0. 61–0. 8), the effects of the remaining set of mutations (1,649 out of 1,957; Wrel ≤0. 6, or > 0. 8) were, as expected, better predicted by PROVEAN (Fig 3E). Specifically, the ability to predict the effect of neutral mutations (Wrel > 0. 8) increased to 72%, (accuracy of 80. 5%, S3 Table). Notably, SIFT, a predictor similar to PROVEAN but that with no phylogenetic correction, showed lower prediction accuracy than PROVEAN (75. 3% accuracy, with 73% true positives for deleterious mutations, and 82. 2% true negatives for neutral mutations, Wrel > 0. 8; S6 Fig and S3 Table). Furthermore, the effects of mutations are best described on a continuum scale rather than a binary classification of deleterious versus neutral. The inclusion of phylogenetic distances, as in PROVEAN, also generates a continuous score that in turn seems to correlate well with the experimental Wrel values (S6C and S6D Fig). Further support to the conclusion that phylogenetic distance is a crucial factor in prediction is provided by the fact that neutral/beneficial drift mutations (Wrel > 0. 8) are decreasingly observed in orthologs as their sequences further diverge from M. HaeIII’s. In total, 39% of the single nucleotide exchanges observed in orthologs were found as neutral/beneficial in the background of M. HaeIII (Wrel > 0. 8, S7A Fig). Of these, 54% (i. e. , 21% of all single nucleotide ortholog-observed exchanges) appear in close orthologs (≤35% divergence relative to M. HaeIII, S7B Fig). Thus, sequence exchanges observed in orthologs with higher sequence identity are more indicative of a neutral fitness effect in the reference sequence (M. HaeIII in our case) than when the very same exchanges appear in diverged orthologs (>50% divergence). Finally, we examined how the experimentally measured fitness effects correlate with predicted structural and functional constrains, and to what degree these constraints apply to the sequence diversity seen in orthologs of M. HaeIII. To this end, a multiple sequence alignment was derived for M. HaeIII (S5 Fig). The alignment gave a set of 2,000 exchanges that are observed in at least one ortholog–dubbed ‘ortholog-observed’ (845, single nucleotide exchanges, and the remaining 1,155 being double and triple exchanges). The complementary set of ‘ortholog-unobserved’ was accordingly derived, and included 4,251 exchanges that are not observed in any known ortholog (of which 1,112 are single nucleotide exchanges). We then compared the experimental set of ‘neutral mutations’ (Wrel ≥ 0. 8, i. e. , excluding ‘nearly-neutral’) to ‘ortholog-observed’, and the set of ‘deleterious mutation’ Wrel ≤ 0. 6) to ‘ortholog-unobserved’. This comparison is a priori problematic. The background at which these two sets of mutations occurred differs fundamentally: the maximal divergence in the laboratory drifted G17 sequences was ~3% (an average of 9. 6 mutations per a length of 329 amino acids gene). Accordingly, whereas our experimental set comprises only single nucleotide mutations, most natural amino acid exchanges relate to two or three nucleotide exchanges within the same codon (1,155 out of 2,000 orthologs-observed, and 3,139 out of 4,251 orthologs-unobserved). Despite the above caveat, we found that drift mutations with ‘neutral/beneficial’ effect (Wrel > 0. 8), and accordingly, orthologs-observed exchanges, share the same biophysical constraints with respect to M. HaeIII’s configurational stability and enzymatic function (Fig 4). Specifically, mutations predicted using the FoldX force field to be highly destabilizing (ΔΔG ≥ 2 kcal/mol; [69,70]) were purged in the laboratory drift (‘deleterious’ Wrel ≤0. 6) and also in the natural diversity (‘orthologs-unobserved’, Fig 4A and 4B). Mutations in positions close to M. HaeIII’s active-site followed the same trend (Fig 4C and 4D). The overlap between the biophysical constrains acting in nature and in the laboratory constraints was also indicated by ‘local closeness’–a structural measure of the degree of structural connectivity of a residue to other residues [71,72] (Fig 4E and 4F). Furthermore, exchanges found in close orthologs appear to obey the above biophysical constraints to a larger extent than those in more diverged ones [41] (S7B and S7C Fig).
Our results illustrate the limitations inherent to the experimental methodologies used for measuring the fitness effects of mutations in the laboratory, and in deducing from these experiments how proteins evolve in nature. In general, the current state-of-the-art experimental mappings artificially widen the threshold for acceptance of mutations, such that the early accumulating mutations have no apparent effect on the protein’s fitness [73]. This wider experimental threshold is driven by various factors, including: (i) higher protein stability (e. g. , the stabilizing mutations we included in, or fusion tags known to increase solubility (e. g. [26]); (ii) Gene and protein copy numbers that are typically orders-of-magnitude higher than the natural ones; (iii) Growth environments that are less demanding than the natural ones; (iv) Proteins being under selection for one task out of several tasks they perform in nature under variable conditions. Such wider thresholds result in a higher, if not unrealistic tolerance of mutations relative to nature [74,75]. Once this threshold is exhausted, the loading of additional mutations results in a rapid collapse [47]. Indeed, the dynamics of our neutral drift experiment indicate that the very same deleterious mutations, including nonsense mutations, which are tolerated in the early rounds, are completely purged as the drift progressed (Figs 1C and S4; see also [76]). The continuous loading of mutations, as applied in our study, appears to portray a more realistic picture with respect to the fraction of deleterious mutations. The distribution of fitness effects of mutations (DFE) derived from this experiment is different from the distributions derived from previous experiments, namely that ~30% of mutations are deleterious, and the remaining largely neutral [25–27,33,46]. In contrast, this experiment indicates the anticipated continuum, rather than the generally assumed bimodal distribution [49,74,77] (Fig 1A). Further, even under the most conservative threshold, 67% of the mutations have evident deleterious effects (Wrel ≤ 0. 6). However, purging is also consistently seen for mutations we categorized as ‘nearly-neutral’ (Wrel = 0. 61–0. 8, Fig 1B). Individually these mutations may be close to neutrality, but collectively they impact fitness, analogously to a population’s drift load [50,78]. This is apparent by the acceptance of ‘nearly-neutral’ mutations being accompanied by the enrichment of compensatory mutations (Fig 1) [79]. If what we categorized as nearly-neutral mutations are included, ~81% of all possible amino acid mutations that derived from single nucleotide mutations are potentially deleterious. The much higher fraction of nonsynonymous mutations with deleterious effects observed in our experiment as compared to other experiments may relate to variations in the mutational tolerance of one protein vs. another. However, M. HaeIII does not seem to be a particularly slow evolving protein–the distribution of the positional evolutionary rates, and specifically the relative histogram areas of the fast versus slow evolving positions, are in agreement with a fast evolving protein [44]. Regardless, it is clear that, at present, comparing the DFEs obtained for different proteins is problematic because the experimental methodologies used to obtain these DFEs vary so much. The sensitivity of detection of fitness effects is also limited in laboratory setups by high noise levels as well as by the limited number of generations along which fitness is examined [74]. We also note that in reality, protein sequences drift in a gradual manner and via single nucleotide exchanges. Thus, the fitness effects measured for all 19 possible amino acids per position often reflect leaps in sequence space that are not taken by natural evolution. The results of our laboratory drift also support the hypothesis that natural protein drift is punctuated by deleterious and compensatory mutations. The order of their accumulation may differ, also in relation to the mutational rates. At high mutational rates, as applied here, compensatory substitutions may follow the deleterious ones [80–83]. At low mutation rates, however, mutations that initially accumulated as neutral may enable the fixation of deleterious ones [44,83]. In any case, the DFE obtained here suggests that, whereas upon drifting in nature, exchanges may be fixed by chance (the neutral theory), their fitness effects are rarely neutral–they are nearly always deleterious or compensatory [41,84] (Fig 1A). Compared to previous reports (for example see [25–27,33,36]), tolerance vs. purging of mutations in our prolonged drift shows much better correlation to the positional evolutionary rates, and to specific exchanges observed in the natural diversity, i. e. , in M. HaeIII’s orthologs (Fig 4). Such a correlation is a priori problematic. The representation of natural sequences is sporadic, especially with horizontally transferred genes that encode specialized functions such as M. HaeIII. Thus, that a certain exchange is not observed, or rarely observed in the currently known sequence does not necessarily mean it is deleterious. Nonetheless, our data seems to coincide with what had been deduced from other analyses of orthologous sequences, namely that at a given background, the vast majority of mutations are deleterious [39–41,43] (Fig 1). Our data also support the notion that the exchanges found in close orthologs are more likely to be neutral than those in more diverged ones [41] (S7 Fig). Exchanges in highly diverged orthologs are tolerated by virtue of being compensated by exchanges at other positions [41,45,67] and therefore tend to be context-specific. However, despite the above caveats, it seems that the effects of mutations can be predicted from the natural diversity of orthologs with relatively high accuracy, particularly when ‘nearly-neutral’ mutations with borderline effects are excluded (Fig 3E). A systematic exploration of the performances of various predictors and prediction parameters is beyond the scope of this work. Nonetheless, it appears that the prediction seems improved when the phylogenetic distance of orthologs is taken into account [85,86] (S6 and S7 Figs). Likewise, comparing the results of different experimental mappings of mutational effects is inherently problematic. These mapping experiments used different proteins, different mutagenesis and screening, or selection, strategies, and different ways of assigning the ‘fitness’ values to mutations. As experimental approaches of systematic mapping develop further, standard experimental and data analysis procedures may develop that will enable more meaningful comparisons. Further, the biophysical constraints acting to limit drift both in the laboratory and in nature overlap, indicating universal constraints that dictate purging of sequence exchanges [87] (Fig 4). Thus, including structural considerations, possibly ‘local closeness’ as an integrated parameter [71], may greatly improve prediction of the effect of mutations [88,89], as already shown for certain predictors [8,89–93]. Overall, the application of the experimental setup described here provides a better understanding of how protein sequences diverge in nature, as well as a new dataset that can be used for improving the prediction accuracy of the effects of protein mutations, and specifically of single nucleotide polymorphisms.
A modified M. HaeIII wild-type gene, carrying four stabilizing mutations [51], and no GGCC sites in its open reading frame, was cloned with an N-terminal His-tag into pASK-IBA3+ vector (IBA, using NcoI and NotI; the vector also carried 14 GGCC sites, 3 of which were located within the ampicillin resistance gene (See supplementary Fig 3 in [51]). Plasmids were transformed into E. coli ER2267 (EcoK r- m- McrA- McrBC-Mrr-) in which GGCC DNA methylation is not toxic [94]. Transformants were selected by growth on ampicillin. Random mutagenesis was performed as described previously [52]. Briefly, M. HaeIII' s ORF (open reading frame) was amplified by PCR with an error-prone polymerase (GeneMorphII Mutazyme, Stratagene). The mutagenic PCR was optimized to an average of 2. 2 mutations per gene. Each round of evolution, or generation (noted as' G' ), included the following steps (S1 Fig): (i) The pool of M. HaeIII genes from the previous round was randomly mutated, recloned using the NcoI and NotI sites, transformed into E. coli and plated on agar plates containing ampicillin. About 106 individual transformants were obtained in each round. (ii) Colonies grown at 37°C overnight were combined, plasmid DNA was extracted and digested with HaeIII (10–20 units, in 50 μl of NEB buffer 2, for 2 hours at 37°C), and re-purified (PCR purification kit, QIAGEN). (iii) The recovered plasmid DNA was re-transformed for another round of enrichment. Each round of drift included one cycle of mutagenesis and three cycles of enrichment (transformation, growth, plasmid extraction and digestion). The naive library, G0, relates to the transformed plasmid DNA derived from cloning of a repertoire of ~105 individual M. HaeIII genes after the first round of mutagenesis and prior to selection by HaeIII digestion. The samples of the naive (G0) and the selected libraries from Rounds 3,7 and 17 (assigned as G3, G7 and G17) were prepared as described previously [52]. Briefly, the pools of M. HaeIII' s open reading frame were PCR-amplified, purified, and concatenated by self-ligation (using XhoI restriction sites at both ends of the PCR [51]). Sequencing libraries were prepared and sequenced according to manufacturer' s protocol at the Weizmann Institute' s high throughput-sequencing core facility. The obtained sequencing reads (~40 nts) were mapped to the reference sequence of wild-type M. HaeIII with two methods: (i) Using NCBI blastn v2. 2. 20 [95] with parameters: e-value cutoff 0. 0001, word size 7, and allowing up to 6 mismatches and requiring a minimal alignment length of 24 consecutive nts, as previously described [96,97]; and (ii) Using Novoalign v2. 07. 00 with parameters: c 4 Hash step-size 6 [96]. Point mutations, insertions and deletions were assigned based on the mapping of the sequencing reads to the reference sequence as previously described [97,98]. Every mismatch or gap in the reads alignment relative to the wild-type reference was recorded per each nucleotide position, and further analyzed using custom Perl scripts (available at: https: //github. com/tawfiklab/HTS_codon_analyzer). Only codons that were intact within the 40 nt reads were included. Processing of the observed mutation counts per codon was done primarily with Excel (see S1 File). All possible single nucleotide mutations were detected in the raw data of the unselected G0 library (329X9 = 2,961 possible single nucleotide mutated codons that in turn comprise the 1,957 possible single nucleotide amino acid mutations; S2 File and S8 Fig). However, Illumina sequencing exhibits a considerable background level of mutagenesis due to PCR amplifications as well as sequencing errors. Potential sequencing artifacts, specifically mutations that were observed at the edges of reads (where sequencing errors are more frequent), were filtered out (S1 File). The background rate was determined using the region upstream of the randomly mutated open reading frame of M. HaeIII (the N-terminal fused His tag that was not subjected to mutagenesis, S2 File, S8 Fig and S4 Table). Thus, the average background frequency was subtracted from the mutational frequencies to give the net positional frequencies that were used to calculate the Wrel value of each amino acid mutation (the final analyzed data can be found in S3 File, including for double and triple mutations that were not analyzed here). Mutational frequencies were determined for every possible codon mutation (63 including single, double and triple nucleotide mutations) as the number of reads with a given mutation (s) divided by the total number of reads that mapped the corresponding position. The frequencies of all mutational events that led to the same amino acid were combined. Orthologous sequences to M. HaeIII were collected using BLASTP search within the REBASE database [99]. Within the range of 25–75% identity, 105 non-redundant family members were aligned using MUSCLE [100] (S5 Fig). The maximum likelihood phylogenetic tree was calculated using PhyML with the LG matrix [101]. The position-specific evolutionary rates (μ) were calculated by Rate4Site [65]. The positional rate is calculated that indicates how fast this site evolves relative to the average rate across all sites in the input alignment. Mutation-specific scores were calculated using PROVEAN and SIFT programs using the default parameters and M. HaeIII sequence as a reference. Both software are available on the homepage of the J. Craig Venter Institute: the SIFT tool is at http: //sift. jcvi. org [102], and the PROVEAN tool is at http: //provean. jcvi. org [68]. To ensure the consistency of this analysis, we provided the set of orthologs sequences in FASTA format (S5 Fig) rather than using the BLAST search in the PROVEAN webserver. The PROVEAN calculated scores were subsequently provided by Dr. Yongwook Choi. FoldX was used to predict the stability effects of mutations relative to wild-type M. HaeIII. The crystal structure of M. HaeIII (PDB id 1dct) [103] was first optimized using the FoldX RepairPDB function. Subsequently, all possible single mutants (19 different amino acids, at each position) were calculated by the BuildModel mutation engine, and relative stability of mutants was obtained (ΔΔG = ΔGWT-ΔGMUT). Distances of residues from the reaction center were defined as the shortest distance between the closest residue atom and either the sulfur of the catalytic cysteine or the methyl group of the SAM cofactor. These were calculated based on M. HaeIII in complex with the DNA (PDB id 1dct) [103]. The coenzyme distances were derived from a homology model based on M. HhaI in complex with SAM (PDB id 2hr1). Local closeness was calculated SPACER web server (available at http: //allostery. bii. a-star. edu. sg/) [71] using default parameters and M. HaeIII structure as a reference (PDB id 1dct) [103].
Understanding and predicting the effects of single nucleotide polymorphisms (SNPs) is of fundamental importance in many fields. Systematic experimental mappings of the effects of such mutations within a given gene/protein comprise an essential experimental tool for determining protein function and for refining models of protein evolution, as well as an important resource for improving prediction algorithms. Here, we present the results of a laboratory system that mimics the manner by which protein sequences diverge in nature: a prolonged process of gradually accumulating random mutations that retain the protein’s structure and function. The change in frequencies of mutations over generations, as obtained by deep sequencing, enabled us to assess the relative effects of all possible SNPs at the background of an accumulating number of mutations. Compared to previous reports, we found that > 80% of all possible amino acid exchanges have potential deleterious effects, with 67% being clearly deleterious. Tolerance vs. purging of mutations in our prolonged drift also showed better correlation with natural diversity. Overall, our experimental setup provides a better understanding of how protein sequences diverge in nature, plus a new basis for improving the prediction accuracy of the effects of protein mutations, and specifically of SNPs.
Systematic Mapping of Protein Mutational Space by Prolonged Drift Reveals the Deleterious Effects of Seemingly Neutral Mutations