Universal distribution of mutational effects on protein stability, uncoupling of protein robustness from sequence evolution and distinct evolutionary modes of prokaryotic and eukaryotic proteins

Protein-coding genes in any organismal lineage evolve at widely different rates, spanning a range of about three orders of magnitude [1, 2]. According to the molecular clock model, each gene is endowed with a characteristic evolutionary rate (ER) that remains approximately constant over long time intervals, up to the entire span of the history of life, in the case of universal genes [3]. However, subsequent measurements have shown that molecular clock is strongly over-dispersed, i.e. the ERs vary to a much greater extent than expected from sampling error alone, under the assumption of a Poisson mutational process [4–7]. The distributions of the rates across sets of orthologous genes in diverse life forms, from bacteria to mammals, which reflect the relative rates of gene evolution, show a notably greater degree of conservation than the absolute rates [1, 2]. This observation inspired the universal pacemaker, a more general model of evolution, that postulates genome-wide, synchronous changes in the ERs of genes. The universal pacemaker model yields a better fit between thousands of individual gene trees and the species tree than the strict molecular clock model [8, 9].

The conservation of the ER distribution implies simple, universal underlying factors. Traditionally, it has been assumed that ER is a multiplicative function of two terms, one of which reflects the intrinsic structural–functional constraints that affect the given protein, whereas the second one corresponds to the biological importance of the same protein [10]. Although this concept is fundamentally plausible and straightforward, for many years, it remained effectively inaccessible to empirical assessment. However, this situation has changed with the advances of functional genomics and systems biology, when prominent correlations have been shown to exist between several evolutionary and molecular phenomic variables [11–15]. For instance, the ER and propensity for gene loss are positively correlated; by contrast, each of these variables is negatively correlated with the gene expression level. Surprisingly, little if any correlation was detected between the essentiality of genes for the reproduction of organisms and the ER: at best, non-essential genes evolve slightly faster than essential genes [15–21]. Among all the detected connections, by far the most consistent and strongest one is the universal anticorrelation between the expression level of a gene and its sequence evolution rate: in all model organisms for which detailed expression data are available, highly expressed genes evolve significantly slower than lowly expressed ones [20, 22–24].

The universal link between gene expression and sequence evolution inspired the hypothesis that protein abundance or, more precisely, translation rate of protein-coding genes is the key determinant of the sequence evolution rate [23, 25, 26]. Specifically, the mistranslation-induced misfolding (MIM) hypothesis posits that the major cause of the covariation between the sequence evolution rate and expression level is the selection for robustness to protein misfolding that is increasingly important for highly expressed genes owing to the toxic effects of misfolded proteins [25–28]. Detailed computer simulations of protein evolution seem to indicate that the toxic effect of protein misfolding, indeed, could suffice to explain the observed covariation of expression level and sequence evolution rate [26]. A direct indication that translation rate substantially contributes to the rate of sequence evolution has been obtained through comparative analysis of the ERs of domains in multidomain proteins that revealed significant homogenization of the ERs compared to the same domain combinations in separate proteins [29]. Quantitative analysis of these results has suggested that the contribution of translation rate was comparable in magnitude to the contribution of structural–functional constraints. Subsequently, a more precise mathematical model has been developed to estimate the relative contributions of structural–functional constraints and translation rate to the rate of evolution of protein-coding genes through a joint comparative analysis of proteomic data on protein abundance and short term ERs in different lineages of animals [30]. The results of this analysis indicate that the contribution of structural–functional constraints is typically 3 to 5 times greater than the contribution of translation rate. Subsequent studies have refined the MIM hypothesis through the demonstration that selection for avoidance of spurious protein–protein interactions additionally contributes to the anticorrelation between expression and ER by constraining the evolution of surface residues in highly expressed proteins [31].

A parallel series of analyses has revealed a strong dependence between the protein core size and the ER of solvent-exposed amino acid residues, and hence complete proteins. It has been shown that increasing the size of the buried protein core results in increased rate of evolution of solvent-exposed residues, which typically constitute the majority of a protein molecule, whereas buried residues are less affected if at all [32–34]. Furthermore, it has been demonstrated that in yeast, expression level and core size exert largely independent effects on the rate of evolution. At the whole protein level, expression has been found to make a substantially greater contribution to the ER than core size although the contributions appeared comparable at the level of individual exposed residues [34, 35].

Another line of investigation has demonstrated a positive correlation between protein contact density and the ER such that robust (designable) proteins with a high contact density evolve measurably faster than proteins with a low contact density [36].

Protein stability that is obviously linked to mutational robustness is considered to be one of the key factors of evolution that seems to impose limits on the ER of protein coding-genes [37–40]. The free energy of protein folding is typically in the range of 5–15 kcal mol⁻¹ [41]. Thus, proteins appear to be marginally stable so that amino acid replacements that result in free energy increase by over 3 kcal mol⁻¹ can be highly destabilizing and hence deleterious [42]. Experimental determination of protein stability remains a relatively difficult task so that genome/proteome scale analysis of the mutational effects on it is currently impractical. However, powerful and apparently accurate computational methods for estimating the free energy of protein globules and the effects of mutations on it have been developed [43–46]. Using one of these methods, FoldX [43, 44], Tokuriki et al have shown that the distributions of mutational effects on protein stability were remarkably similar for 21 diverse proteins and closely matched the available experimental data [47].

We sought to investigate the relationships between protein robustness to mutation, abundance and selection pressure on protein-coding genes (a correlate of ER) by taking the maximum advantage of the growing collection of protein structures from diverse model organisms. We come to the conclusion that variation of mutational robustness is largely neutral with respect to the selective pressure on the sequences of protein-coding genes, especially in prokaryotes. In contrast, in eukaryotes, a moderate dependence exists between the relative protein core size, which is a proxy of robustness, and selection pressure, particularly when measured separately for the surface residues. In addition, we show that proteins from a hyperthermophilic archaeon are on average significantly less robust to mutation than proteins of mesophiles, conceivably due to stronger constraints on protein structure exerted by extreme conditions.

Coverage of genomes of model organisms with protein structures and the universal dependence of dN/dS on protein abundance

Effects of mutations in the core on protein stability

For each protein, we explored the distribution of mutation effects in the core positions on the protein stability. Each core position was computationally mutated independently to each of the 19 other amino acids and the ΔΔG values (i.e. the difference between the ΔG values before and after the mutation) were estimated using the FoldX method [44] (see methods for details). The ΔΔG values were grouped into 1 kcal mol⁻¹ bins from −10 kcal mol⁻¹ to 16 kcal mol⁻¹, and the fraction of mutation resulting in the respective magnitude of ΔΔG was computed for each bin. We generated and analyzed the distributions of the fractions of mutations associated with different ΔΔG values for all possible amino acid replacements (19 mutations per position) and separately for replacements resulting from single nucleotide substitutions (9 possible mutations). In the latter distribution, two additional fractions were included, namely (i) silent mutations (mutations leading to the same amino acid) and (ii) nonsense mutations (mutations leading to a stop codon).

To our knowledge, this is the first genome-wide comparative analysis of the distributions of mutation effects on protein stability. In a previous, small scale study, similar distributions have been generated using FoldX and compared with experimental data for 21 proteins. Tokuriki et al have shown that the typical asymmetric distribution of the ΔΔG values was independent of protein sequence and fold, and that most of the mutations yielded relatively small, positive ΔΔG values (i.e. slightly destabilized the respective proteins) whereas stabilizing mutations (ΔΔG < 0) were rare [47]. The distributions for single nucleotide substitutions included fewer destabilizing mutations and more neutral mutations (mutations with a ΔΔG around 0 kcal mol⁻¹) than the complete distributions due to the well-characterized mutational robustness of the genetic code [49, 50].

The universal distribution of mutation effects across proteomes

The mutation effects show highly similar, skewed bell-shaped distributions, with a heavy tail of destabilizing mutations for all analyzed proteomes (figure 1). For the purpose of this comparison, we represent the mutational effects for all proteins in a given species with a boxplot, where the box represents the distribution of fractions of mutations resulting in specific ΔΔG values for the given proteome. Figure 1 directly compares the distributions of mutation effects for five analyzed proteomes (i.e. S. cerevisiae, H. sapiens, E. coli, B. subtilis and T. gammatolerans; C. elegans and D. melanogaster proteomes were excluded from this analysis due to the small number of available structures). There was no significant difference between the eukaryotic and bacterial proteomes suggesting that the distribution of the mutation effects on protein stability is universal across the bacterial and eukaryotic domains of life. In contrast, the only archaeon and the only (hyper)thermophile in the data set, T. gammatolerans deviated from this universal pattern, with significantly more destabilizing mutations than found in other organisms (figure 1 and supplementary table S1). This feature of the T. gammatolerans proteome is compatible with the previous observations of stronger constraints on protein evolution in hyperthermophiles compared to mesophiles [51, 52].

Zoom In Zoom Out Reset image size

In the entire set of analyzed proteins, there are very few strongly stabilizing mutations (ΔΔG < − 4 kcal mol⁻¹) (figure 1). The effects of the great majority of the mutations (84%) lied between −3 and 6 kcal mol⁻¹, and 47% were between 0 and 3 kcal mol⁻¹. Given that the free energy of protein folding is distributed roughly between 5 and 15 kcal mol⁻¹ [41, 42], these findings indicate that most of the mutations are near neutral (have a negligible effect on protein stability) or slightly destabilizing. As already pointed out, the distribution is skewed toward destabilizing mutations but strongly destabilizing mutations (ΔΔG > 6 kcal mol) are rare. In T. gammatolerans, the fraction of stabilizing mutations, with effects between −3 and 0 kcal mol⁻¹, is significantly lower than in other analyzed organisms, and conversely, the fraction of strongly destabilizing mutations is significantly greater (Welch two-sample test). The distributions of the effects of single nucleotide substitutions show similar trends (compare figures 1(a) and (b)).

Although the shape of the distribution of mutation effects is universal, individual proteins in each proteome show broad variation, with some being dramatically more robust to mutations than others (figure 2). For instance, in some proteins, more than 30% percent of all mutations are slightly destabilizing, with ΔΔG between 2 and 3 kcal mol⁻¹ whereas in other proteins, only about 6% of the mutations have effects in this range. In the following sections, we investigate factors that potentially could explain this variance.

Zoom In Zoom Out Reset image size

Effects of mutations on protein stability and core size

As shown by extensive theoretical analysis, under the mean field approximation, the fraction of stabilizing amino acid replacements (i.e. replacements with ΔΔG < 0) can be used as a proxy of protein robustness to mutation [53]. This measure has the advantage of being independent of the actual ΔΔG values estimated by FoldX that can be inaccurate in some cases (Thiltgen and Goldstein 2012), but only on the fraction of these values that is below a particular threshold. In all analyzed organisms, a strong negative correlation was observed between the fraction of stabilizing mutations (pstab) and the median of the ΔΔG distribution that can be considered as a quantitative characteristic of the overall mutational effect on the stability of a protein (figure 3).

Zoom In Zoom Out Reset image size

We then examined the relationships between median ΔΔG and pstab as proxies of the mutational robustness of a protein and other structural, functional and evolutionary features. Highly significant positive correlations were observed between the relative core size, which was calculated as the proportion of buried residues in a given protein, and ΔΔG median for all 5 analyzed proteomes (figure 4). For pstab, the expected negative correlation with the relative core size was observed for four of the five analyzed proteomes; however, these correlations were weak and were significant only for S. cerevisiae and E. coli (supplementary figure S1). Furthermore, for T. gammatolerans, a significant positive correlation was observed (see discussion below). A direct proportionality between the mutational robustness of a protein and the relative core size can be expected because stable proteins with large cores can tolerate more destabilizing mutations than less stable proteins with small cores [32, 54]. In evolutionary terms, stable, large-core proteins are not subject to strong selective pressure to minimize the effect of mutations on stability, in contrast to proteins with small cores that are less robust to destabilizing mutations but are optimized by evolution such that many mutations have little or no destabilizing effect. The observation that ΔΔG median significantly and positively correlates with the relative core size, whereas pstab shows much weaker negative correlations, implies that on genome scale the former reflects protein robustness better than the latter. Accordingly, all further analysis were performed with both measures of protein robustness but we focus primarily on the results obtained with ΔΔG median.

Zoom In Zoom Out Reset image size

Mutational effects on stability, protein abundance and purifying selection

A comparison of the ΔΔG median to the protein abundance showed no significant correlation for S. cerevisiae, H. sapiens and T. gammatolerans. In contrast, a significant negative correlation was observed in bacteria (figure 4) where highly abundant proteins appear to have evolved to allow fewer destabilizing mutations (i.e. are less mutationally robust) than low abundance proteins.

The ΔΔG median for each protein was further compared to the dN/dS value, i.e. the strength of purifying selection, for the respective gene. All proteomes showed some positive correlation that, however, was found to be significant only in H. sapiens (figure 4). Thus, slow-evolving proteins tolerate fewer destabilizing mutations, i.e. are somewhat less mutationally robust, than faster-evolving proteins. The correlation increased when only surface residues were included in the comparison, and in this case, the correlation became significant also for E. coli (figure 4). Apparently, mutations in the core have a greater effect on the stability and evolution of the surface areas of proteins that on the core itself, in accord with previous observations [34, 35].

No significant correlations with either protein abundance or dN/dS were observed for pstab (supplementary figure S1). Taken together, these observations unexpectedly indicate that robustness of proteins to mutation inferred from the ΔΔG distribution obtained using FoldX is linked to abundance in bacteria but not in eukaryotes and is largely decoupled from sequence evolution in all tested organisms. As shown previously, the link between ER and robustness could be stronger at the per-residue level than it is on the whole protein level [35, 55, 56].

Relative protein core size, evolution and abundance of proteins

The relative core size (proportion of amino acid residues in the core) is a complementary measure of the mutational robustness of a protein that has been reported to be linked to the evolution of protein-coding genes [51, 52]. For each structure, we identified solvent-inaccessible positions, i.e. the core positions that primarily contribute to the protein fold maintenance [57] (see methods). A moderate but significant positive correlation between dN/dS value and proportion of core positions (relative core size) (methods) was observed in S. cerevisiae, D. melanogaster and H. sapiens but not in prokaryotes (table 2 and figure 5). This positive correlation was observed also after correction for protein size using partial correlation analysis although in this case, it remained significant only for S. cerevisiae. Using the dN/dS values estimated for the core positions only, we found a significant correlation in eukaryotes (except for C. elegans and D. melanogaster) but again no correlation in prokaryotes. Finally, we estimated the dN/dS only from the surface (exposed) positions and detected a significant positive correlation with the relative core size in all species except for T. gammatolerans; the correlation coefficients for the surface residues were greater than those for core residues or complete proteins. Notably, however, the correlation coefficients were on average higher in eukaryotes (i.e. 0.24 for H. sapiens, 0.26 for C. elegans, 0.38 for D. melanogaster and 0.32 for S. cerevisiae) than in bacteria (i.e. 0.24 for B. subtilis and 0.12 for E. coli). Using more restrictive thresholds for solvent accessibility (5% or 0) yielded similar but somewhat weaker trends than those detected with the 25% (data not shown). These findings imply that increase in the relative core size results in an increased robustness of the respective protein to mutations, and as a consequence, lowers the pressure of purifying selection in eukaryotic proteins. This effect of the core size is most pronounced in the solvent-exposed, surface parts of proteins, and for prokaryotic proteins, this is the only portion for which the (in this case, weak) connection between the core size and selective pressure is detectable at all.

Zoom In Zoom Out Reset image size

It seems likely that the difference in the effects of the relative core size on protein evolution in prokaryotes and eukaryotes are due to the greater overall compactness of prokaryotic proteins. Indeed, a comparison of the distributions of core sizes in eukaryotes and prokaryotes showed that the mode of the distributions was shifted toward larger relative core size in prokaryotes whereas the tail of small cores was much heavier in eukaryotes (figure 6(a)). Statistical analysis indicates that the difference between the distributions of the relative protein core sizes in prokaryotes and eukaryotes is highly significant (supplementary table S2) whereas between different eukaryotes and prokaryotes the distributions are statistically indistinguishable (figure 6(b)). Thus, it appears likely that the much stronger connection between the relative core size that is observed in eukaryotes compared to prokaryotes results from the 'fragility' of eukaryotic proteins with small cores that slows down their evolution, in agreement with recent observations [58].

Zoom In Zoom Out Reset image size

Notably, there was a significant negative correlation between the relative core size and protein abundance that was substantially more pronounced in bacteria (but not in the archaeon T. gammatolerans) than it was in eukaryotes, and in the former case, survived the correction for the universal correlation between abundance and dN/dS (table 2). Conceivably, protein abundance is associated with an increased number of protein–protein interactions which require a large exposed area to be available in a protein globule, and this effect is particularly pronounced in compact bacterial proteins. For eukaryotic proteins, the selective pressure correlated with the relative core size much stronger than with the median ΔΔG. This difference potentially could be attributed to low accuracy of ΔΔG estimation by FoldX. However, given that such a difference was not observed for prokaryotic proteins, it appears more likely that the relative core size reflects aspects of protein robustness that are distinct from those reflected in the median ΔΔG.

Structural and evolutionary trends in individual sets of orthologous proteins

What determines the rate of protein evolution (or more precisely, selection on protein sequences), is a storied, complex and obviously difficult problem. Clearly, protein abundance (or translation rate) is linked to the sequence evolution rate but can hardly be considered the single or even the most important determinant [15, 59]. Efforts on quantitative dissection of the effects of the rate of translation and domain-specific structural and functional constraints suggest that the contributions of both factors are important but the latter is measurably greater [30]. What is the nature of these constraints on protein evolution?

Here, we took advantage of the steadily increasing database of protein structures from diverse organisms to explore some plausible candidates, in particular mutational robustness of protein domains. If protein misfolding is a key component of the fitness cost of mutations, a direct connection should be expected to exist between the robustness of protein folds to mutation and protein stability, on the one hand, and abundance, evolution rate and evolvability of proteins, on the other hand [26, 28, 30, 39, 40]. Despite the striking progress in experimental methods for sequence manipulation, protein expression and especially experimental measurement of protein stability still lag far behind, making genome-wide analysis of protein robustness impractical for any foreseeable future. Accordingly, at present, this problem can only be addressed using computational simulation of protein folding. FoldX employs an empirical force field that was developed to rapidly calculate the free energy (ΔG) of proteins and evaluate the effect of amino acid replacements on ΔG and accordingly the stability of proteins; this effect is measured as ΔΔG, i.e. the difference between the free energies of the native and mutated states [43, 44]. The limited available comparisons with experimental data seem to attest to a reasonable accuracy of the computational estimates [47], and the universality of the ΔΔG distributions observed here is compatible with this assessment. Therefore, in the absence of practicable direct approaches, FoldX appears to be the method of choice for a genome-wide analysis of protein robustness and stability and their evolutionary correlates.

The first result of this analysis to be emphasized is the significant positive correlation between the median ΔΔG of a protein with its relative core size. This correlation seems to indicate that mutational robustness of proteins increases with the enlargement of the core. In other words, proteins with a large proportion of amino acid residues in the core are not 'fragile' and thus have not evolved to minimize the effect of mutations. Given that the core size and ΔΔG are determined from the protein structures independently, this finding appears indicative of the structural–functional relevance of the ΔΔG estimates. Larger protein cores can accommodate more destabilizing mutations than smaller cores, i.e. are associated with a greater mutational robustness of proteins. This definition of robustness might seem somewhat counter-intuitive given the greater effect of destabilizing mutations on protein stability but the strong correlation between the median ΔΔG and relative core size implies that it reflects the realities of protein evolution.

The distribution of the ΔΔG values for protein cores was found to be universal among eukaryotes and bacteria, with the mode and median at slightly destabilizing amino acid replacements and a heavy tail of more strongly destabilizing mutations. These findings indicate that the shape of the free energy landscape of the core is a fundamental feature of globular proteins, with similar levels of mutational robustness attained by protein evolution in a broad variety of life forms. This distribution joins the growing list of universal patterns identified by quantitative comparative genomics [60]. By implication, this distribution, with the characteristic level of optimization, can be inferred to have evolved at an early stage of the evolution of each widespread protein fold, probably before the divergence of the three domains of cellular life. This conclusion is compatible with the observations that radiation of the common protein folds into major families antedates the last common ancestor of the extant cellular life forms [61]. However, the universal ΔΔG distribution appears to be significantly shifted towards destabilizing mutations in hyperthermophiles (pending validation by analysis of additional proteomes from hyperthermopihlic organisms), resulting in a greater characteristic robustness that seems to reflect the stronger stability-determined constraints on proteins of thermophiles compared to those of mesophiles [51, 52].

Unexpectedly, in most cases, no significant connection was found to exist between the median ΔΔG (a proxy of mutational robustness) and the pressure of purifying selection (dN/dS) affecting a protein sequence (and accordingly, the rate of protein evolution) although the sign of the observed weak correlation was as expected under the hypothesis that more robust (less fragile) proteins are subject to weaker constraints than less robust ones. The connection between mutational robustness and protein abundance is also weak and is significant only in bacteria. Conversely, the relative core size, which is directly linked to robustness, showed significant correlation with dN/dS only in eukaryotes although a stronger dependence was detected to exist for surface amino acid residues in all organisms. For eukaryotic proteins the relative core size was found to be a much better predictor of selective pressure (dN/dS) than it was for prokaryotic proteins. In part, this difference could be attributed to the inaccuracy of FoldX but the lack of a similar difference in prokaryotes implies that the relative core size and median ΔΔG reflect distinct aspects of protein robustness. The nature of this difference remains to be elucidated.

Taken together, these findings seem to indicate, perhaps counter-intuitively, that robustness to destabilizing effects of mutation is not a major, defining factor of protein evolution, at least on the relatively short scale that is suitable for the dN/dS measurements and on the level of whole proteins. Conceivably, as suggested above, for the common protein folds, the degree of robustness of the cores that is essential for protein stability and function was attained at early stages of evolution whereas much of the subsequent variation appears to be effectively neutral. This conclusion holds particularly well for the compact, highly optimized protein of prokaryotes that generally are subject to strong purifying selection linked to the characteristic large effective population size [62] and show some dependence on core robustness only for surface residues. In the less compact, typically more fragile eukaryotic proteins, the connection between robustness and sequence evolution is significantly stronger.

In this work, we took advantage of the growing collection of protein structures from diverse organisms to examine, on the whole genome scale, the connections between biophysical characteristics of proteins and their evolution. We find that the distribution of the mutational effects on protein stability is a universal characteristic of cellular life forms, with significant deviation detected only in a hyperthermophilic archaeon. The results of this analysis imply that a protein's robustness to mutation, at least on the whole protein level, is, to a large extent, uncoupled from its abundance and sequence evolution. Such uncoupling appears to be particularly pronounced in the compact, apparently highly optimized proteins of prokaryotes. The less compact, more fragile eukaryotic proteins show a notably stronger connection between the relative protein core size, taken as a measure of mutational robustness, and sequence evolution rate, particularly at the protein surface. The general conclusion from this analysis is that the minimum required or optimal mutational robustness was reached at early stages of protein fold evolution whereas much of subsequent variation was neutral.

Data sets

Evolutionary analysis

Identification of orthologs between distant species

Protein structures

Core positions

ΔΔG calculation

Distributions of mutational effects on protein stability

The authors thank Yuri I Wolf for critical reading of the manuscript and numerous helpful suggestions, and the Koonin group members for useful discussions. The authors' research is supported by intramural funds of the US Department of Health and Human Services (to the National Library of Medicine).

Author contributions

GF performed data analysis, GF and EVK analyzed the results, EVK wrote the manuscript that was read and approved by both authors.