{'fraction_non_alphanumeric': 0.02930027204209963, 'fraction_numerical': 0.023993221246042012, 'mean_word_length': 5.011662198391421, 'pattern_counts': {'":': 0, '<': 13, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Variation in Mutational Robustness between Different Proteins and the Predictability of Fitness Effects

Random mutations in genes from disparate protein classes may have different distributions of ﬁtness effects (DFEs) depending on different structural, functional, and evolutionary constraints. We measured the ﬁtness effects of 156 single mutations in the genes encoding AraC (transcription factor), AraD (enzyme), and AraE (transporter) used for bacterial growth on L -arabinose. Despite their different molecular functions these genes all had bimodal DFEs with most mutations either being neutral or strongly deleterious, providing a general expectation for the DFE. This contrasts with the unimodal DFEs previously obtained for ribosomal protein genes where most mutations were slightly deleterious. Based on theoretical considerations, we suggest that the 33-fold higher average mutational robustness of ribosomal proteins is due to stronger selection for reduced costs of translational and transcriptional errors. Whereas the large majority of synonymous mutations were deleterious for ribosomal proteins genes, no ﬁtness effects could be detected for the AraCDE genes. Four mutations in AraC and AraE increased ﬁtness, suggesting that slightly advantageous mutations make up a signiﬁcant fraction of the DFE, but that they often escape detection due to the limited sensitivity of commonly used ﬁtness assays. We show that the ﬁtness effects of amino acid substitutions can be predicted based on evolutionary conservation, but those weakly deleterious mutations are less reliably detected. This suggests that large-effect mutations and the fraction of highly deleterious mutations can be computationally predicted, but that experiments are required to characterize the DFE close to neutrality, where many mutations ultimately ﬁxed in a population will occur.


Introduction
The distribution of fitness effects (DFE) of random mutations for a whole genome or individual genes is a key factor in determining how evolution progresses in response to mutation, selection, and drift. Fitness constraints can be present at many levels, including protein and mRNA expression levels, protein stability, and structure as well as interactions with other cellular components, which makes prediction of the DFE extremely complex (Travisano and Shaw 2013;de Visser and Krug 2014). Several methods have been used to gain insight into the parameters of the DFE, including mutation accumulation experiments, fitness measurements of single induced mutations, and DNA sequence data (Eyre-Walker and Keightley 2007;Bataillon and Bailey 2014). Distributions inferred from sequence analysis concern mutational effects small enough that SNPs can be detected in populations whereas experimental approaches can only detect mutations with relatively large effects, with the limit set by the sensitivity of the fitness assay used. Whereas many studies focus on the structure-function relationship and phenotypic effects of mutations (Suckow et al. 1996;Jacquier et al. 2013;Jimenez et al. 2013;Firnberg et al. 2014;Sarkisyan et al. 2016), sometimes providing data for thousands of mutations, the evolutionary relevance of such mutational landscapes is uncertain. For example, Kondrashov and co-workers (Sarkisyan et al. 2016) measured the effects of mutations in Green Fluorescent Protein (GFP). In this case, the gene was moved from its original organism, the jellyfish Aequorea victoria, into E. coli and mutated, codon optimized, fused to another protein, and expressed from a plasmid at high level. "Gene fitness" was then measured with low sensitivity as fluorescence, with excitation at a wavelength not primarily selected for in its original organism, in an external environment and temperature very different from where the gene had evolved, and with no connection between the function of the gene and cell fitness. Whereas this and similar studies certainly contribute to understanding constraints on the molecular level, their relevance for evolution or for organismal fitness landscapes is debatable.
Relatively few studies have directly measured the fitness effects of random single mutations in protein-coding genes in their native context with high sensitivity (Lind et al. 2010;Sanjuan 2010;Hietpas et al. 2011) and provided explanations for which factors influence the shape of DFEs. We previously reported that two ribosomal protein (RP) genes (Lind et al. 2010) had unimodal DFEs peaked at slightly deleterious mutations, with an average fitness effect of s ¼ À0.012, no complete loss of function, and very few neutral mutations. This finding contrasts with the distributions reported for viruses, b-lactamase, and GFP where a high fraction of complete loss-of-function mutations was observed (Sanjuan 2010;Hietpas et al. 2011;Jacquier et al. 2013;Firnberg et al. 2014;Sarkisyan et al. 2016). To address these divergent results and their underlying evolutionary rationale, we used highly sensitive competition experiments to measure the DFE of random mutations for three genes from different protein classes in their native context to investigate if the high robustness of RPs is an exception and if any general properties of the DFE can be established. Specifically we wanted to determine the fraction of deleterious, neutral and advantageous mutations, investigate the molecular basis for mutational effects and the potential for in silico prediction, and examine if there is variation depending on differences in protein function.

Experimental System
Three genes in Salmonella typhimurium LT2, required for rapid growth on L-arabinose, were chosen for random mutagenesis based on their different molecular functions, conditional dispensability, and large fitness costs in a selective environment, i.e., growth on arabinose as the only carbon and energy source, when disrupted ( fig. 1A). We introduced 156 single base pair substitutions, excluding stop codons, into araC (n ¼ 55), araD (n ¼ 47), and araE (n ¼ 54) in their native chromosomal position in S. typhimurium and estimated their fitness (s) in one-to-one competition with wild-type controls ( fig. 1B) in the presence of arabinose minimal medium. This set-up allowed detection of fitness differences of jsj > 0.004. For araD the mutations were introduced by mutagenic primers where the position and type of mutation were randomized with uniform probabilities. For araC and araD, errorprone PCR was used to generate mutations and the mutational spectrum is biased towards transitions (supplementary table S4, Supplementary Material online). Mutants with stop codons were excluded for both methods.
AraD is a homotetrameric epimerase that catalyzes the final step in the degradation pathway for L-arabinose to Dxylulose-5-phosphate using a Zn 2þ as cofactor. AraC is a transcription factor that binds L-arabinose and activates transcription of the araBAD operon and araE in the absence of glucose and represses its own transcription and that of araBAD in the absence of arabinose (Schleif 2010). AraC belongs to the AraC/XylS family of transcriptional regulators that are widespread among bacteria and are involved in regulation of carbon metabolism, stress response, pathogenesis, and antibiotic resistance (Gallegos et al. 1997). AraE is an arabinose/proton symporter that transports L-arabinose across the inner membrane (Daruwalla et al. 1981). It belongs to the major facilitator (MFS) superfamily of importers and exporters, which are present in all domains of life and have been implicated in antimicrobial resistance and in human diseases including cancer and type 2 diabetes (Quistgaard et al. 2016).
The AraCDE proteins represent three major functional classes, transcription factor, enzyme, and transporter, which are expected to have different biophysical constraints. Although all these proteins need to recognize and interact with a small molecule to perform their function there are FIG. 1. Experimental system and methods. (A) Three genes (araC, araD and araE) encoding proteins required for rapid growth on L-arabinose were chosen for random mutagenesis. AraE transports L-arabinose into the cell, AraC binds L-arabinose and activates transcription of araE, araC, and araBAD. AraD catalyzes the final step in the degradation pathway for L-arabinose to D-xylulose-5-phosphate, which then enters the pentose phosphate pathway. (B) The change in ratio between mutant and wild-type populations is determined after growth competition under a serial transfer protocol. The fitness of each mutant relative to the wild-type is calculated as a selection coefficient (s) per generation that is 0 for neutral, negative for deleterious and positive for advantageous mutations under growth in minimal medium with arabinose.
Variation in Mutational Robustness . doi:10.1093/molbev/msw239 MBE additional constraints in interactions with membrane components for AraE, DNA binding and RNA polymerase interaction for AraC and oligomerization for AraD and AraC. There are also differences in system level properties, where AraC regulates transcription of the entire system, the level of AraE determines the flow of L-arabinose across the membrane and the production of D-xylulose-5-phosphate by AraD is the final biochemical determinant for flow into the pentose phosphate pathway ( fig. 1A).
The fitness effects of a complete loss-of-function mutation in each of the three genes were determined by examining growth on arabinose when the respective gene was deleted. Loss of AraD is lethal when grown on L-arabinose due to accumulation of a toxic intermediate (Englesberg et al. 1962) and s for DaraD was set to À1. Deletion of AraC results in loss of activation of transcription for araBAD and araE (s ¼ À0.36). As AraE is not the only way for the cells to take up L-arabinose, deletion of this gene had less dramatic consequences on fitness (s ¼ À0.25).
Synonymous Mutations in araC, araD, and araE Do Not Reduce Fitness Synonymous substitutions accounted for 47 of the 156 single base pair substitutions and all of them were classified as "neutral", which here means they cannot be reliably distinguished from wild-type controls (jsj < 0.004, in the shaded area in fig.  2A-C, full details in supplementary table S1, Supplementary Material online). The mean cost for all synonymous substitutions was s ¼ 2 Â 10 À5 (SD ¼ 0.0014), suggesting that their effects are very close to neutrality on average and the frequency of non-neutral synonymous mutations is likely to be below 6% (P ¼ 0.05 hypergeometric test, assuming random sampling and correct classification). In contrast, for the RP genes only two out of 38 synonymous substitutions were "neutral" and the average cost was s ¼ À0.0096 (SD ¼ 0. 0064). This difference could be explained by higher expression levels of RPs (Bakshi et al. 2012), more optimized codon usage and that a minor reduction in the amount of an RP reduces growth to a significant degree (Knöppel et al. 2016). Firnberg et al. (2014) also reported significant effects of single synonymous mutations on TEM-1 activity when expressed from a tac promoter, suggesting that these effects might vary substantially depending on gene, expression level, and selective conditions. The AraCDE Proteins Are Less Robust than Ribosomal Proteins, but Have a Higher Fraction of Advantageous Mutations For the amino acid substitution mutants, the fractions of complete loss-of-function mutants were similar for AraD (7/36) and AraC (9/38) whereas only one of the AraE mutants was close in fitness to the deletion mutant ( fig. 2D). Deleterious mutants that still maintained some function made up a similar fraction as the complete loss-of-function mutants for AraD (9/36) and AraC (15/38), whereas they were more common for AraE (20/38) ( fig. 2D).
Amino acid substitution mutants that were indistinguishable from wild-type fitness (jsj < 0.004) made up 33-56% of all mutations ("neutral" in fig. 2D). This provides an upper limit for the fraction of mutants that are nearly neutral in a population genetics sense with fitness effects smaller than the inverse of the effective population size, which for many bacteria would be in the order of jsj < 10 À5 to 10 À9 (Berg 1996;  Charlesworth 2009). We also found mutants with significantly higher fitness, one in AraC (M164T) and three in AraE (T338I, V31I, R257H) (two-sided t-tests P < 10 À4 ) ( fig. 2A and C; sup plementary table S2, Supplementary Material online). This suggest that several percent of random mutations might be advantageous when the mutated genes are under strong selection, but the few mutations found here makes estimates of frequency uncertain (f ¼ 2.6%, f > 0.9%, P ¼ 0.05, hypergeometric test) assuming that our sampling is random and that mutations are classified correctly as advantageous or not). This frequency is similar to that reported for mutations in TEM-1 b-lactamase that increased resistance to cefotaxime and increased fitness under selection for function (3.4%) (Schenk et al. 2012). The fitness increases for the mutations in araC and araE were relatively small (s ¼ 0.0061-0.019) and would not have been detected using less sensitive fitness assays, but as most mutations are expected from theory to have small effects (Fisher 1930) this again highlights the importance of investigating the DFE close to neutrality. Similarly high rates of advantageous mutations with small effects have previously been reported (Desai et al. 2007;Perfeito et al. 2007), but they are often not detected in experimental evolution studies where they are lost through population bottlenecks and competition with rare high fitness mutants (Sniegowski and Gerrish 2010).

MBE
Whereas 95% of random substitutions in the genes of RP S20 and RP L1 were weakly deleterious this is expected to be largely due to constraints on the mRNA level given that synonymous substitutions had effects of similar magnitude (Lind et al. 2010). If the average effect of the synonymous mutations in RPs is subtracted from the nonsynonymous, the average cost is only s ¼ À0.0037 Nonsynonymous mutations in ribosomal proteins never reduced fitness more than s ¼ À0.1, even though deletion of the genes resulted in much larger fitness costs (estimated maximum frequency for all nonsynonymous mutations 3.3% P ¼ 0.05, hypergeometric test) (Lind et al. 2010). For the AraCDE protein, these highly deleterious mutations represented 27% of amino acid substitution mutants (95% confidence interval for all nonsynonymous mutations 23-31%, assuming random sampling, and correct classification) ( fig. 3A and B).

Selection to Reduce Costs of Misreading Leads to High Robustness for Ribosomal Proteins
What could be the evolutionary reason for the approximately 30-fold higher robustness of ribosomal proteins as compared with the AraCDE proteins? One possibility is that genes central to growth, such as RPs, might evolve robustness to reduce the costs of deleterious mutations, but based on the calculation outlined below we argue that this is unlikely. The fitness of a population under the pressure of random mutations [rate u % 10 À10 (Foster et al. 2015)] is given by e -uM , where M is the number of positions in the genome where a mutation leads to a reduction in fitness (Kimura and Maruyama 1966). The increase in population fitness if one of the deleterious positions were transformed into a neutral one, would be Ds ¼ ue -uM % u. This is an extremely small selective advantage that could be effective only in a population that has a very large effective size, i.e., for N e > 10 10 or when the mutation rate is very high. The effective population size for many bacteria is usually estimated to be smaller than the inverse of the mutation rate (N e 10 9 ) (Berg 1996;Charlesworth 2009), which means that mutation pressure alone would have an extremely small effect on the selection for an increased mutational robustness at each individual site. If we allow the removal of one deleterious position to also reduce the cost of mutations in i other positions, the fitness gain for the population is Ds ¼ iue -uM % iu, which is still very small unless the number of other stabilized sites are exceedingly large. Thus, although the effect of a mutation in an FIG. 3. Distribution of fitness effects for amino acid substitutions mutations in (A) ribosomal proteins S20 and L1 is unimodal and (B) AraC, AraD, and AraE is bimodal. The average cost of a synonymous mutation was subtracted from selection coefficients previously obtained for RP S20 and L1 (Lind et al. 2010).
Variation in Mutational Robustness . doi:10.1093/molbev/msw239 MBE individual cell can be large as indicated by the DFE, the effect of deleterious mutations on the fitness of the whole population is exceedingly small and selection on the mutation pressure is unlikely to shape the DFE.
In contrast, transcription/translation errors will occur at a much higher rate p ($10 À3 to 10 À4 per amino acid; Drummond and Wilke 2009;Meyerovich et al. 2010;Imashimizu et al. 2013) in all cells and the effect on the population fitness can be large. We can estimate the effects of misreading errors on fitness in the following simplified model. Consider a protein of size n amino acids that is present in m copies in the cell. When all copies of the protein carry the same mistake (like a mutation), the fitness effect is s as measured in the DFE. Assume that the effects are additive so that when only one copy of the protein carries the error, the effect on growth is s/m. Assume further that the average effect of all possible missense errors in the protein is the same as the average of the DFE, i.e., s=m per copy, which would be the case if most translation errors are due to single base-pair misreading. The expected number of errors at all codons and in all copies of the protein is n Â p Â m. Thus, the overall average effect on fitness from missense errors in this protein is Ds ¼ n Â p Â s. For a ribosome with n % 7500 (Nikolay et al. 2001), s¼ À0.0037, and assuming p ¼ 2 Â 10 À 4 , this gives Ds ¼ À0.0056. The fitness cost of misreading errors is similar for a protein of average size, with mutational robustness like the AraCDE proteins [Ds ¼ À0.0072 assuming n % 300 (Brocchieri and Karlin 2005) and s ¼ À0:12]. However, if the early primitive ribosomes had the same DFE as an average protein under strong selection (exemplified by the AraCDE proteins with s¼ À0.12), then fitness would be greatly reduced with Ds ¼ À0.18. Thus, for large protein complexes (like ribosomes), the selective advantage of shifting the average of the fitness spectrum to smaller values, thereby increasing robustness, would be very large and could explain the high robustness of modern ribosomes.

Inability of Mutations in Ribosomal Proteins to Disrupt Intermolecular Interactions Provides a Mechanism for Mutational Robustness
What could be the mechanistic basis for the experimentally observed difference in robustness between AraCDE versus RP S20 and L1? A structure-fitness analysis of the AraCDE proteins using homology models (supplementary fig. S1, Supplementary Materials online) suggests that mutations with large deleterious effects appear to disrupt specific functions (supplementary table S2, Supplementary Material online) rather than being generally destabilizing. For the AraD protein, conditional lethal mutations are in the direct proximity of the active site (L16R, G172C, and H118Y), including Zn 2þ binding residues, or in the C-terminal end (Y229N and G230R) that is directed across the active site of the adjacent subunit (Luo et al. 2001;Schleif 2010). Although AraC and AraE lack enzymatic active sites, a similar pattern emerges where highly deleterious mutations in AraE appear mainly in the transport channel (A46V, D39H, L155P, T158R, and Y131C), disrupt salt bridges (P202S), or are close to residues involved in substrate recognition (G273V, S400L, and W405C) (Sun et al. 2012;Iancu et al. 2013). In AraC, highly deleterious mutations are near residues implicated in arabinose binding (I57V, T50I, G55D, and R84Q), DNA binding (S207L and R178C), interaction with RNA polymerase (G250E, A191S, and D256V), or dimerization (L160S, L108P, L161F, Y103H, and P39L) (supplementary table S2, Supplementary Material online) (Soisson et al. 1997;Schleif 2010).
In contrast to the AraCDE proteins many RPs primarily serve a structural function, which could allow them to be stabilized by numerous interactions with rRNA. Specific functions are difficult to assign to individual ribosomal proteins, but RP S20 and L1 do not interact directly with mRNA or tRNA or form part of the peptidyl transferase center and their roles might be described as mainly structural. S20 is a primary rRNA binding protein that initiates proper folding of the 5 0 domain of the 23S rRNA and is crucial for binding of ribosomal protein S16. S20 is crucial for anchoring rRNA helix 44, which is part of the A-and P sites (Schluenzen et al. 2000) and forms bridging contacts with the 50S subunit (Selmer et al. 2006). Deletion of RP S20 leads to reduced mRNA binding and inhibits association with the 50S subunit (Tobin et al. 2010). RP L1 binds directly to 23S rRNA and is not known to influence the binding of other ribosomal proteins. It is located on a flexible stalk region that is believed to be involved in structural rearrangements that is coupled to movement of deacetylated tRNA in and out of the E-site (Ning et al. 2014). Deletion of L1 leads to aberrant ribosome biogenesis and L1 is required for balanced formation of subunits (Kåhrström et al. 2012). Mutational robustness of RPs might be increased by introduction of specific amino acids at the protein-RNA interface (Jones et al. 2001) so that function may not be significantly disrupted by single amino acid substitutions. In addition, strong selection for proper folding and ribosome assembly of proteins containing amino acid substitutions could potentially lead to the recruitment of molecular chaperones (Fares 2015) and assembly factors (Shajani et al. 2011).
An alternative explanation for the increased robustness of RPs is to evolve a codon usage that avoids single mutations that lead to amino acid substitutions with large effects on protein structure. Our data does not support this hypothesis based on similar distribution of values from substitution matrices for the mutations in RPs compared with the AraCDE proteins (supplementary fig. S2A and B, Supplementary Material online). The predicted effects of mutations on changes in thermodynamic stability between the two classes are also similar (supplementary fig. S2C and D, Supplementary Material online). However, it is possible that ribosomal proteins have evolved robustness by increasing their stability threshold by selection for generally stabilizing mutations (Kather et al. 2008).

Predicting Mutational Effects of Amino Acid Substitutions on Fitness
Prediction of the DFE from sequence or structural data could play a major part in evolutionary forecasting efforts and complement experimental studies. However, it is not clear to what degree fitness effects are predictable, which methods would Lind et al. . doi:10.1093/molbev/msw239 MBE be most successful and how this depends on the functional class the protein belongs to. One fast and readily available method is to use amino acid sequence conservation in general, in the form of substitutions matrices such as PAM and BLOSUM, which were reported to be correlated with fitness for TEM-1 b-lactamase (r ¼ 0.4) (Firnberg et al. 2014). For AraCDE, RP S20 and RP L1, amino acid substitution matrices did not show a significant predictive ability ( fig. 4A, r ¼ 0.04-0.39, P > 0.01) with the exception of PAM10 and AraD (r ¼ 0.50, P ¼ 0.0015) (supplementary fig. S3A and B and table S3, Supplementary Material online). Possibly this is due to the fact that AraD is the only enzyme among the five proteins and that other classes of proteins are underrepresented in the sequences used to construct amino acid substitution matrices.
Methods based on multiple alignments of homologous sequences proved more successful and for AraCDE these have a significant predictive ability (Consurf r ¼ 0.47-0.58, P < 0.004, Provean r ¼ 0.47-0.61, P < 0.002), whereas for the ribosomal proteins evolutionary conservation did not correlate significantly with fitness cost (r < 0.26, P > 0. An alternative approach is to use structural information from experimental structures or homology models to make computational predictions of changes in thermodynamic protein stability, which has been recognized as an important determinant of the fitness costs of mutations (DePristo et al. 2005;Wylie and Shakhnovich 2011). In contrast to the alignment-based methods, information about active site residues or vital intermolecular interactions is not included, which is expected to reduce correlation with fitness because many complete loss-of-function mutations will not be predicted. The simplest of structure-based parameters is the relative solvent accessibility (RSA) of the mutated amino acid. A low RSA represents a buried amino acid and that mutations in that position are more likely to disrupt protein structure. Given the indication that many mutations classified as neutral were at the surface of the proteins, we calculated the exposed surface area of each amino acid in homology models of both ribosomal proteins and AraCDE and found a significant correlation between fitness cost and RSA ( There are a large variety of publically accessible programs to estimate the effect of mutation on protein stability [reviewed in Kucukkal et al. (2014)]. A comprehensive evaluation of these methods for predicting fitness effects was not performed because several of the studied proteins did not have experimental structures for S. typhimurium and initial efforts using three different methods generally performed worse than for RSA ( fig. 4C; supplementary table S3 and fig.  S3F and G, Supplementary Material online). In addition, structure-based methods are expected to provide very low correlations for the ribosomal proteins given that a major part of the fitness constraints is on the mRNA level.
We used several publically available methods that predict phenotypic effects of mutations by combining evolutionary conservation and structural effects as well as a range of other parameters using neural network classifiers (SNAP2; supple  . Whereas these methods were not primarily designed to predict fitness effects, but rather the relevance of mutations for human disease, they could have substantial predictive ability also for nonhuman proteins. It is however reasonable to expect these methods to be better suited for predicting large effect mutations that are have been shown to be strongly associated with disease, have large effects on protein function or are in highly conserved positions in sequence alignments. For AraD, AraE and RPS20 there were minor improvements from the methods using only conservation and SNPs&GO using homology models of the proteins produced the highest correlations (AraD r ¼ 0.71, P < 10 À4 ; AraE r ¼ 0. 68, P < 10 À4 ; and S20 r ¼ 0. Clearly some of these methods can produce highly significant correlations between score and fitness effect and this is largely due to their ability to identify mutations with very large effects, which does not necessarily mean that they are useful for distinguishing between deleterious and neutral mutations. The most successful methods for classifying mutations in the five proteins as either deleterious or neutral were Provean (72% correct) and SNAP2 (73% correct) and several methods produced higher than 60% correct classifications (supplementary table S3, Supplementary Material online). Up to 89% of the highly deleterious (s < À0.03) mutations were predicted as costly whereas the proportions were lower for slightly deleterious (% 70%) and about 60% of mutations predicted as neutral had no significant fitness effects (supple mentary table S3, Supplementary Material online). It should be noted that the proportion of truly neutral mutations is expected to be substantially lower than what we were able to determine as not different from wild-type, due to the limited sensitivity of the competition experiments (i.e., mutants with jsj < 0.004 are not reliably distinguishable from wild-type).
In summary, none of the methods provide meaningful predictions for mutations in ribosomal proteins ( fig. 4), which can be explained by the strong constraints on the mRNA level and high robustness of ribosomal proteins. Whereas several methods can, at first glance, predict the fraction of deleterious mutations in ribosomal proteins, this is likely to be due to the very low frequency of apparently neutral mutations (5%), which means that the method that predicts the highest fraction of deleterious mutations will perform best. For the AraCDE proteins, predictions were more successful and quite equally so despite the different functions of the proteins. Methods based on evolutionary conservation by multiple alignments performed best. There appears to be little to gain by also including additional information, as done by the hybrid/machine learning methods (SNAP2, Suspect, MutPred and SNP&GO) compared with Provean and Consurf that only use conservation ( fig. 4B vs. fig. 4D). This is true both when comparing correlation between prediction score and size of fitness effect (average correlation 0.51 for hybrid methods vs. 0.53 for conservation methods) and mutations classification as neutral or deleterious (average correctly classified for hybrid methods 0.65 vs. 0.71 for conservation methods). Thus, the choice of method among those based on evolutionary conservation might be guided mainly by the usability of user interfaces and speed.

Conclusions
An essential part of understanding evolutionary processes and the ability to forecast mutational trajectories is a deeper knowledge of DFEs and their underlying mechanistic explanations. If each gene has a unique DFE there is little hope to arrive at any general principles for forecasting, but this seems unlikely given the high functional and structural conservation observed between biological systems. The similarities between the DFEs of AraC, AraD, and AraE, despite their diverse molecular functions, suggest that a bimodal distribution is widespread among different classes of proteins and that the fractions of deleterious mutations are similar (% 50%). This provides a general expectation for the shape of the DFE for individual genes under strong selection, but the results for the RPs demonstrate that some classes of genes have completely different DFEs.
We suggest that selection to reduce the costs of misreading errors is the main evolutionary force behind the high robustness of ribosomal proteins and possibly other proteins in large complexes. This hypothesis is based on wellestablished population genetics theory and experimental measurements of misreading rates and provides a general expectation for the robustness of proteins across functional classes. The importance of misreading errors and evolution towards increased mutational robustness has previously been demonstrated experimentally in mutants with increased error rates (Goldsmith and Tawfik 2009;Bratulic et al. 2015) under stringent selection for function. An additional factor that could work with the misreading hypothesis to drive RPs towards higher robustness than the AraCDE proteins is a difference in the strength of selection. Thus, in natural settings, the arabinose proteins are most likely under strong selection only during short periods when arabinose is limiting for growth (as is the case in the present experiments). Under conditions where arabinose is not limiting, most mutations in the ara genes would be (near-) neutral and their DFE would look robust with a much smaller value of s than that measured here. In contrast, ribosomal proteins are in constant use and would be under selection under all growth conditions. In general, the shape of the DFE and its robustness must depend on the growth conditions under which it is measured.
The frequency of nonsynonymous mutations that were classified as deleterious was 56% for the AraCDE proteins and 73% if also including the ribosomal proteins. This is similar to the fraction strongly deleterious effect in mammalian species of about 70% (Eyre-Walker et al. 2002;Yampolsky et al. 2005). We classified 40% of the amino acid substitutions in AraCDE as "neutral", but the truly neutral mutations with 1 < N e Â jsj is likely to represent a minor part of this given that we can only detect jsj > 0.004 whereas fitness costs 10 4 times smaller could be subject to natural selection in many bacterial species. In humans, with an effective population size of about 10 4 (Charlesworth 2009), the fraction of slightly deleterious (s > À10 À3 ) and neutral mutations has been estimated to 44-57% (Bataillon and Bailey 2014), which suggest that the major differences in the fractions of neutral and deleterious substitutions between species is largely an effect of variation in population sizes.
Based on the pattern of highly deleterious mutations in AraCDE we suggest that these strong fitness effects, and the resulting bimodal fitness distributions, are not primarily due to reduced thermodynamic stability. This been suggested as the main cause of highly deleterious mutational effects in TEM-1 beta lactamase and GFP, but it is possible that disruption the buried functional sites in these proteins also contributes to the correlation (Firnberg et al. 2014;Sarkisyan et al. 2016). Instead our hypothesis is that the amino acid Lind et al. . doi:10.1093/molbev/msw239 MBE substitutions disrupt intermolecular interactions and directly reduce function through loss of proper binding affinity to a small molecule substrate (AraC, AraD, and AraE), loss of DNA or protein recognition (AraC) and disruption of oligomerization (AraC and AraD). If this interpretation would prove correct, studies of monomeric enzymes or fluorescent proteins might not be representative for the majority of proteins. The mechanistic reason for the lack of highly deleterious mutations in ribosomal proteins could then be a result of their binding to ribosomal RNA through numerous amino acid residues and this interaction cannot be completely disrupted by single substitutions. Although rates of molecular evolution are largely determined by biophysical constraints on thermodynamic stability and protein abundance (DePristo et al. 2005;Drummond et al. 2005;Serohijos et al. 2013), the shapes of DFEs are, under this hypothesis, mainly determined by the potential for loss of intermolecular interactions. This idea is similar to the concept of functional density, which was introduced by Zuckerkandl (1976) as an hypothesis to explain rates of molecular evolution. However, we suggest that the proportion of molecular sites involved in specific functions as compared with general functions (i.e., functional density) determines the fraction of highly deleterious mutations and forms the basis for bimodal DFEs, but it does not determine evolutionary rate at the molecular level. Our results also show that very highly conserved proteins, like ribosomal proteins, are not necessarily fragile as previously suggested (Assis and Kondrashov 2014), and that fitness constraints on the mRNA level might also reduce evolutionary rates.
The effects of individual nonsynonymous mutations on fitness can be predicted with some success using evolutionary conservation or machine learning methods, achieving a correlation of up to 0.71 and correctly classifying up to 73% of mutations as deleterious or neutral, but with highly varying results for different genes depending on methods used (sup plementary table S3, Supplementary Material online). The difficulty of correctly predicting mutations with small effects of s ¼ À0.03 to À0.004, a range that is not possible to study using less sensitive fitness assays, suggests that a substantial fraction of deleterious mutations are not correctly classified using these methods. This could have relevance for the genetic basis of complex human disease where alleles with small deleterious effects are not easily detected, but could make a significant contribution to the etiology of the disease when combined (Gibson 2011). If the DFEs of different protein classes are conserved within and between species then prediction of mutational effects would be greatly simplified, which could inform evolutionary forecasting by reducing the region of sequence space that must be experimentally explored to estimate the robustness of a gene or functional network.

Strains and Media
All strains used are derived from S. enterica subsp. enterica serovar Typhimurium str. LT2 (designated S. typhimurium in the text). Luria-Bertani agar (LA) was used as solid media, supplemented with sucrose 6%, kanamycin 30 mg/l, chloramphenicol 12.5 mg/l, or tetracycline 7.5 mg/l for selection following transformation or maintenance of the pSIM5 k Red plasmid. LA plates without added sodium chloride were used for selection against sacB. Lysogeny broth (LB), SOC, or M9 minimal media supplemented with 0.4% (w/v) L-arabinose were used as liquid media.

Strain Construction
Mutations in araD (STM0101), araC (STM0104), and araE (STM3016) were introduced using k Red recombineering (Datsenko and Wanner 2000). For araD, mutations were selected by computer randomization with uniform probability for mutations location and mutation type, but excluding stop codons. Mutagenic primers containing the desired mutation 6 bp from the 3 0 -end and with 40 bp of homology 5 0 to the mutation were used to PCR amplify part of the araD open reading frame and an adjacent kanamycin resistance marker using Phusion high fidelity DNA polymerase (Thermo Scientific). For araC and araE we used the Genemorph II random mutagenesis kit (Agilent Technologies) to introduce mutations in the open reading frames using pZE21 plasmids with araC or araE in the multiple cloning site as template (0.5-0.8 lg per 20 ll reaction for 21-25 cycles).
Mutated DNA was electroporated into k Red competent cells and for araD kanamycin resistance was used to select for recombinants. For araC and araE, we constructed strains with a sacB marker, conferring sucrose sensitivity, inserted into the chromosome replacing the araC and araE open reading frames. This allowed us to create marker-less mutations by introducing the mutated DNA into k Red competent cells and then select for sucrose resistance and verify inserts with PCR and sequencing. This was the main reason for changing from the more labor-intensive mutagenic primers method used for araD, which had the advantage of introducing a less biased spectrum of mutations (supplementary table S4, Supplementary Materials online). All mutations were moved using P22 transduction into strains with fluorescent protein marker cassettes (amp-PLlacO-mTagBFP2/sYFP2, Genbank: KM018302/KM018301 (Gullberg et al. 2014)) inserted at a neutral position (cbiA, STM2035; Knöppel et al. 2014) for the competition experiments. For the araD mutants kanamycin resistance was used to select for successful transductions and for the araC and araE mutants the strains with fluorescent proteins had sacB cassettes inserted into the araC or araE genes which allowed selection for sucrose resistance after P22 transduction. This protocol reduces the impact of secondary mutations/rearrangements introduced during k Red recombineering and adaptive mutations during laboratory manipulation and creates two independent strains to ascertain that the fitness effects are dependent on the mutation introduced. It should be noted that mutations in araCDE causing extreme toxicity leading to lethality when grown on LB agar cannot be detected in our experiments, but there is little reason to expect these kinds of mutations to exist given the low expression level in the absence of arabinose and the molecular functions of the genes.
Variation in Mutational Robustness . doi:10.1093/molbev/msw239 MBE Independent control strains were isolated from the same populations as the mutants with the only difference being the lack of mutations (i.e., failed mutagenesis) when sequenced. All colonies were selected at random, only dependent on their position on the agar plate, after transformation, transductions, and for competitions. All strains were stored in 10% dimethyl sulfoxide and 90% LB at À80 C and all cultures were frozen immediately after completion of transformation or transduction protocols.

Fitness Assay
Strains were taken from the freezer and incubated overnight at 37 C on LA plates. Single colonies were used to inoculate 180 ll LB cultures in 96-well round bottom plates (Thermo Scientific #262162) and grown for 6 h with 900 rpm shaking (Grant PHMP-4) at 37 C. Mutant and control strains containing opposite fluorescent markers were mixed 1:1 (5 ll each) and grown overnight under the same conditions. After 16 h, 1 ll of each culture was transferred to 180 ll M9 arabinose and grown for 24 h (900 rpm, 37 C). At the same time, each culture was diluted 200-fold in phosphate buffered saline (PBS) for flow cytometry. In order to estimate the ratio of wild-type to mutant for each culture we counted 100,000 cells and determined the fraction of BFP and YFP positive cells using flow cytometry (MACSQuant VYB, Miltenyi Biotec). After 24 h growth the protocol was repeated to obtain the ratio after approximately seven generations of growth. We maintained the cultures for three additional serial transfers of 1.4 ll culture to 180 ll M9 arabinose before measuring final ratios after 28 generations of growth.
The selection coefficients were determined using the regression model s ¼ [ln(R(t)/R(0))]/[t], as previously described (Dykhuizen 1990) where R is the ratio of mutant to wild-type. For mutants with very high fitness costs (s < À0.15) the ratios in the final count were too low and selection coefficients were calculated using the ratios from the first flow cytometry measurements (after growth in LB) and after 24 h growth in M9 arabinose, representing approximately seven generations. For the majority of mutants the selection coefficients were calculated from the second flow cytometry assay (after 24 h growth in M9 arabinose) and after three serial transfers, approximately 21 generations of growth.
The fitness of each mutant was measured in triplicate for both BFP and YFP markers, where each competition was inoculated with independent single colonies of mutant and control strains. The differences in fitness cost of the fluorescent markers were measured using control competitions included for each competition plate (minimum of six cultures) and used to adjust the selection coefficients. Three (araC, araE) or four (araD) wild-type control strains that were independently constructed on different days, were used to measure the standard variation between genetically identical strains. The average cost of the BFP marker relative to the YFP marker were s ¼ À0.0018 (SD ¼ 0.0021, n ¼ 24) for araD ¼ À0.0028 (SD ¼ 0.0020, n ¼ 17) for araC, and s ¼ À0.0015 (SD ¼ 0.0019, n ¼ 18) for araE (supplementary fig.  S5, Supplementary Material online). Based on these results, we decided that a mean selection coefficient between À0.004 and 0.004 (approximately 62 SD) for six cultures is designated as "neutral" in the sense that they cannot be reliably distinguished from the wild-type. After adjustment for the differential cost of the marker, the fitness effects in the BFP and YFP strains were highly correlated (r ¼ 0.979, median difference sBFP-sYFP ¼ 4 Â 10 À5 indicating no systematical differences between strains compared with the wild-type reference).
Adaptive or compensatory mutations can arise during competition experiments and to reduce the influence of such events on estimated fitness we removed a maximum of one extreme outlier if it had a selection coefficient more than two standard deviations from the average of all replicates. For highly deleterious mutations that could clearly be classified as costly, outliers were not removed. In cases where more than one replicate diverged from the others or there was a large difference between YFP and BFP strains, the competition was repeated. If differences between YFP and BFP remained, the strains were independently reconstructed using k Red recombineering. This was also done for synonymous mutations with non-neutral effects or for any mutations with significantly higher fitness than the wild-type control.

Supplementary Material
Supplementary tables S1-S4 and figures S1-S5 are available at Molecular Biology and Evolution online.