1. Evolutionary Biology
Download icon

Molecular function limits divergent protein evolution on planetary timescales

  1. Mariam M Konaté
  2. Germán Plata  Is a corresponding author
  3. Jimin Park
  4. Dinara R Usmanova
  5. Harris Wang
  6. Dennis Vitkup  Is a corresponding author
  1. Columbia University, United States
  2. Division of Cancer Treatment and Diagnosis, National Cancer Institute, United States
Research Communication
  • Cited 1
  • Views 1,295
  • Annotations
Cite this article as: eLife 2019;8:e39705 doi: 10.7554/eLife.39705

Abstract

Functional conservation is known to constrain protein evolution. Nevertheless, the long-term divergence patterns of proteins maintaining the same molecular function and the possible limits of this divergence have not been explored in detail. We investigate these fundamental questions by characterizing the divergence between ancient protein orthologs with conserved molecular function. Our results demonstrate that the decline of sequence and structural similarities between such orthologs significantly slows down after ~1–2 billion years of independent evolution. As a result, the sequence and structural similarities between ancient orthologs have not substantially decreased for the past billion years. The effective divergence limit (>25% sequence identity) is not primarily due to protein sites universally conserved in all linages. Instead, less than four amino acid types are accepted, on average, per site across orthologous protein sequences. Our analysis also reveals different divergence patterns for protein sites with experimentally determined small and large fitness effects of mutations.

Editorial note: This article has been through an editorial process in which the authors decide how to respond to the issues raised during peer review. The Reviewing Editor's assessment is that all the issues have been addressed (see decision letter).

https://doi.org/10.7554/eLife.39705.001

Introduction

As proteins evolve from a common ancestor, their sequences and structures diverge from each other (Chothia and Lesk, 1986; Povolotskaya and Kondrashov, 2010). Multiple previous studies have investigated the relationship between the conservation of protein molecular function, sequence identity (Lee et al., 2007; Tian and Skolnick, 2003; Worth et al., 2009) and structural similarity (Chothia and Lesk, 1986; Wilson et al., 2000). For example, the likelihood that two proteins share the same molecular function, given their sequence (Tian and Skolnick, 2003) or structural (Wilson et al., 2000) similarity, has been used to investigate the emergence of new protein functions (Rost, 2002; Conant and Wolfe, 2008), and to perform functional annotations of protein sequences (Lee et al., 2007; Wilson et al., 2000). In this work, we focused on a different and currently unaddressed set of questions. Namely, how far can two sequences diverge while continuously maintaining the same molecular function? What are the temporal patterns of this divergence across billions of years of evolution? And how different protein sites contribute to the long-term divergence between orthologs with the same molecular function? We note that the requirement for the continuous conservation of molecular function is crucial in this context, as multiple examples of convergent evolution and protein engineering demonstrate that the same molecular function, such as catalysis of a specific chemical reaction, can in principle be accomplished by proteins with unrelated sequences and even different folds (Bork et al., 1993; Galperin et al., 1998; Omelchenko et al., 2010).

It was previously demonstrated that proteins with the same structural fold frequently diverge to very low (~10%) levels of sequence identity (Rost, 1997). These results suggest that conservation of protein fold, that is, the overall arrangement and topological connections of protein secondary structures (Murzin et al., 1995), exerts relatively minor constraints on how far protein sequences can diverge. In contrast to protein folds, it is possible that conservation of specific molecular functions will significantly limit the long-term divergence of protein orthologs. While only a relatively small fraction of protein residues (~3–5%) are usually directly involved in catalysis (Lehninger et al., 2013), recent analyses have demonstrated that even sites located far from catalytic residues may be significantly constrained in evolution. Because substitutions at these sites can have substantial effects on molecular function (Firnberg et al., 2016), it is likely that functionally-related sequence constraints extend far beyond catalytic residues.

In this study, we explored the long-term divergence patterns of protein orthologs by characterizing their pairwise sequence and structural similarity as a function of their divergence time. We used several models of molecular evolution to calculate the divergence rates, defined as the decrease in pairwise sequence identity or structural similarity per unit time, between orthologous proteins with the same molecular function. We also characterized the long-term divergence patterns at protein sites with different levels of evolutionary conservation, different locations in protein structures, and different experimentally measured fitness effects of amino acid substitutions. Finally, we explored how the limits of sequence and structural divergence after billions of years of evolution depend on the degree of functional conservation between orthologs.

Results

To study the evolution of proteins with the same molecular function, we initially focused our analysis on enzymes because their molecular function is usually well defined. The Enzyme Commission (EC) classifies enzymatic functions using a hierarchical 4-digit code (Bairoch, 1999), such that two enzymes that share all four EC digits catalyze the same biochemical reaction. We used protein sequences representing 64 EC numbers from 22 diverse model organisms across the three domains of life (Supplementary file 1). The considered activities include members of all six major enzyme classes: oxidoreductases, transferases, hydrolases, lyases, isomerases and ligases.

To investigate whether the conservation of enzymatic function limits the divergence between orthologous sequences, we first calculated global pairwise sequence identities between orthologs as a function of their divergence times (Figure 1, Figure 1—figure supplement 1). The pairwise divergence times reported in the literature (Hedges et al., 2006) between the considered 22 species (Supplementary file 1) were used as a proxy for the divergence times between corresponding orthologous proteins. For each enzymatic activity, we constructed phylogenetic trees based on the orthologous protein sequences and compared them to the corresponding species’ trees. Protein sequences with clear differences to the species'phylogenetic tree topologies, suggesting cases of horizontal gene transfer, were excluded from the analysis (see Materials and methods).

Figure 1 with 7 supplements see all
Sequence divergence of enzyme orthologs as a function of time.

The global pairwise sequence identities between pairs of orthologous enzymes are shown as a function of divergence times between the corresponding species. Three models of amino acid substitution were used to fit the divergence data. Model 1 (black lines) assumes independent and equal substitution rates across all protein sites. Model 2 (red lines) assumes, in addition to independent and equal substitution rates, that a given fraction of protein sites remains identical at large divergence distances. Model 3 (blue lines) assumes a gamma distribution of amino acid substitution rates across sites. Best fits of the models are shown for four representative EC numbers: (a) 1.5.1.3, (b) 2.7.4.3, (c) 4.2.1.2, (d) 6.3.4.2. The horizontal dashed black lines represent the average sequence identity for the global alignment of unrelated protein sequences. The data and corresponding model fits for the other EC numbers considered in the analysis are given in Figure 1—figure supplement 1 and Supplementary file 2a.

https://doi.org/10.7554/eLife.39705.002

We next considered two simple models of long-term protein evolution, one without a long-term limit of sequence divergence and the other with an explicit divergence limit. The first model corresponds to sequence divergence with equal and independent amino acid substitution rates across all proteins sites (Dickerson, 1971; Zuckerkandl and Pauling, 1965); see Equation 1, where y represents global sequence identity, t represents divergence time, and R0 represents the average substitution rate (Dickerson, 1971). In this model, back substitutions are not allowed, and sequence divergence slows down with time simply due to multiple substitutions at the same protein sites and progressively fewer non-mutated sites. The second model corresponds to sequence divergence where, in addition to sites with equal and independent substitution rates, there is a minimal fraction of identical sites at long divergence times; the fraction of identical sites is represented by Y0 in Equation 2.

(1) y=100*e-R0*t
(2) y=Y0+100-Y0*e-R0*t

We applied the two models to fit the sequence divergence for each of the considered enzymatic functions. The best model fits for four representative metabolic activities are shown in Figure 1 (black for the first model and red for the second); the fits for the remaining metabolic activities are shown in Figure 1—figure supplement 1. In 62 of the 64 cases, the second model fits the divergence data significantly better than the first model (F-test p-value<0.05, Supplementary file 2a). Moreover, in 95% of the cases (61/64) the maximum likelihood value of the parameter Y0 is significantly higher (Wald test p-value<0.05) than the average sequence identity between random protein sequences based on their optimal global alignment (~13.5%, shown in Figure 1 and Figure 1—figure supplement 1 by dashed black lines). The distribution of the fitted parameter Y0 suggests a long-term sequence identity >25% (with average ~40%) between considered orthologs (Figure 2a); this demonstrates that conservation of a specific enzymatic function significantly limits long-term protein sequence divergence. Notably, model two is mathematically equivalent (see Materials and methods) to a divergence model with equal substitution rates across sites, a limited number of amino acid types accepted per site, and allowed back substitutions (Tajima and Nei, 1984; Gilson et al., 2017; Yang, 2006). In this model, the parameter Y0 represents the inverse of the effective number of acceptable amino acid types per site during protein evolution. Our results thus suggest that, on average, only 2 to 4 amino acids are acceptable per site for proteins that strictly conserve their molecular function (Figure 2a, top blue X axis); we note that this low average number does not contradict the fact that more than four amino acid types could be observed at a given protein site at low frequencies (Breen et al., 2012).

The limit of long-term protein sequence divergence between orthologous proteins.

(a) The distribution of Y0 parameter values across 64 EC numbers for fits based on Model 2 (Equation 2). The Y0 parameter represents the minimum percentage of protein sites that remain identical at long divergence times. The parameter Y0 (considered as a fraction) can also be interpreted as the inverse of the average number of amino acid types accepted per protein site during long-term protein evolution (top blue X axis). (b) Similar to panel a, but for 29 ancient protein families annotated with non-enzymatic functions. In panels a and b, the vertical red dashed lines represent the median values of the distributions (39% and 30%, respectively).

https://doi.org/10.7554/eLife.39705.014

The two aforementioned models simplify the process of sequence divergence by considering the same substitution rates across protein sites. A more realistic and commonly used model of protein evolution assumes a gamma distribution (Yang et al., 2000) of substitution rates across protein sites; see Equation 3; (Ota and Nei, 1994), where α represents the shape parameter of the gamma distribution. The best fits of such a variable-rate model (blue in Figure 1 and Figure 1—figure supplement 1) showed that the rates of protein sequence divergence between orthologous enzymes have decreased by more than 10 times during ~4 billion years of evolution (see Materials and methods and Supplementary file 2b). Although the third model does not explicitly consider a long-term divergence limit, the obtained model fits also show that the vast majority of orthologous enzymes with the same function will remain above 25% sequence identity on the timescales when Earth environments will be hospitable to life (1–3 billion years from the present O'Malley-James et al., 2014) (Figure 1—figure supplement 2).

(3) y=100*R0*tα+1-α

The observed divergence limit is not an artifact due to difficulty detecting remote protein homologs, as it occurs at relatively high sequence identities (Figure 1 and Figure 1—figure supplement 1), for which corresponding orthologs can be easily identified by computational sequence comparison methods. Furthermore, the results remained similar when we restricted the analysis to orthologous enzyme pairs with experimentally validated molecular functions (Figure 1—figure supplement 3), based on publications referenced in the BRENDA database (Chang et al., 2015). The results also remain robust towards the variance in the estimates of divergence times between considered species (see Materials and methods). We note that the divergence limit between orthologs with the same molecular function does not imply that the rates of molecular substitutions decrease in evolution. It is also not simply due to the curvilinear relationship between time and sequence identity caused by multiple mutations at the same sites; specifically, the observed decrease in divergence rates is substantially higher (by >10 fold) than the one expected in model one simply due to multiple substitutions at the same protein sites. Instead, the effective limit is reached when, due to a small number of amino acids types accepted per protein site and back substitutions, additional amino acid replacements do not lead to a substantial further increase in protein sequence and structural divergence (Meyer et al., 1986).

Interestingly, following the previously introduced metaphor of the expanding protein universe (Povolotskaya and Kondrashov, 2010; Dokholyan et al., 2002), we can use the third model (Equation 3) to express the divergence rate between orthologs as a function of protein distance (D = 1 – y, where y is the fractional sequence identity ranging from 0 to 1), see Equation 4. This equation, similarly to Hubble’s law of universe expansion (Hubble, 1929), describes how the divergence rate depends on the distance between protein orthologs. In contrast to the real universe, the expansion rate of the protein universe significantly decreases with divergence time and with distance between protein orthologs. For example, the divergence rate between protein orthologs drops, on average, to only ~2% sequene identity decrease per billion years when their mutual sequence identity reaches 30% (corresponding to protein distance of 70%; Figure 1—figure supplement 4).

(4) Dt= R0*(1-D)(α+1)α

The analyses described above focused on the divergence of enzymes with the same molecular function. In order to investigate whether the observed divergence patterns are not specific to enzymes, we repeated the analysis using non-enzymatic ancient orthologs (Figure 1—figure supplement 5, Supplementary file 2c). The set of analyzed 29 protein families included ribosomal proteins, heat shock proteins, membrane transporters, and electron transfer flavoproteins (Supplementary file 2d). Based on the same set of 22 species used in the analysis of enzyme families, we found that model two fitted the data significantly better than model one, and that the parameter Y0 was >25% for the majority (23/29) of the non-enzymatic protein families (Figure 2b, Supplementary file 2c). Interestingly, we also identified 19 additional orthologous groups showing two clearly different divergence patterns (Figure 1—figure supplement 6), with pairs of eukaryotic orthologs diverging substantially faster and farther than prokaryotic orthologs in the same protein family. The orthologous groups with this behavior included mitochondrial ribosomal proteins and initiation factors of mitochondrial translation (Supplementary file 2e). It has been previously postulated that mitochondrial ribosomal proteins diverged significantly faster in eukaryotes, compared to the divergence between their bacterial orthologs, due to compensatory protein substitutions following the accumulation of slightly deleterious substitutions in the mitochondrial ribosomal RNA (Barreto and Burton, 2013).

Having established, in the first half of the manuscript that conservation of molecular function significantly limits long-term sequence evolution, we investigated, in the second half, how different protein sites contribute to the observed divergence constraints. Specifically, whether the same protein sites are conserved between ancient orthologs in different phylogenetic lineages, how sites with different fitness effects of amino acid substitutions contribute to the divergence limit, and how structural locations of protein sites affect their long-term divergence patterns. We also explored how different levels of functional specificity constrain sequence and structural divergence.

To investigate whether the same protein sites are usually conserved between orthologs in different phylogenetic lineages, we aligned the sequences of ancient enzyme orthologs with the same molecular function (see Materials and methods). We then calculated how often each protein site is occupied by identical amino acids across pairs of orthologs from phylogenetically independent linages (Figure 3—figure supplement 1). Orthologous protein pairs from independent lineages were obtained using species' pairs that do not share any edges in the corresponding phylogenetic tree (Arnold and Stadler, 2010) (Figure 3a); for example, in Figure 3a the pair D-H is independent of the pair A-B but not of the pair E-F. We performed the above analysis using 30 enzymatic activities for which at least 20 independent pairs of orthologs could be identified based on annotations in the KEGG database (Kanehisa et al., 2016) (see Materials and methods). The results demonstrated that only a relatively small fraction of protein sites (10–20%) are universally conserved, that is, they are identical in a majority (>90%) of independent lineages (Figure 3b). Therefore, the observed long-term divergence limit is not primarily due to sets of universally conserved protein sites; instead, different sites contribute to the limit in independent phylogenetic lineages. By comparing the fractions of universally conserved sites to the average sequence identity between distant orthologs (~40%, Figure 2a) we found that these sites account, on average, for only ~35% of the observed sequence identity at long divergence distances. The analysis also revealed that different protein families show different probability distributions of identical sites (Figure 3—figure supplement 1). This is likely a consequence of diverse structural and functional requirements across protein families, leading to protein-family specific constraints on protein sites.

Figure 3 with 1 supplement see all
Conservation of protein sites in phylogenetically independent lineages.

To identify the fractions of protein sites that are universally conserved ― defined as sites that are identical in at least 90% of orthologs ― we considered phylogenetically independent lineages. (a) Illustration of pairs of species (e.g. A–B and D–H) representing phylogenetically independent lineages. In the figure, A-B and D-H are pairs of species that diverged within a certain time window (illustrated by the blue shaded region); these species pairs do not share more recent edges in the phylogenetic tree. (b) The distribution of the fraction of universally conserved sites across 30 enzymatic families. The analysis was performed using 30 enzymatic families for which at least 20 independent pairs of orthologs with the same molecular function could be identified based on annotations in the KEGG database (Kanehisa et al., 2016) (see Materials and methods); pairs of orthologs that diverged >2 billion years ago were selected for this analysis. Error bars represent the S.E.M. based on three replicates using different sets of orthologous pairs. The dashed red line represents the median of the distribution (~13%).

https://doi.org/10.7554/eLife.39705.015

We next investigated the long-term divergence patterns at protein sites with different fitness effects of amino acid substitutions. To that end, we experimentally measured the fitness effects of all possible single amino acid substitutions in a representative enzyme, the Escherichia coli dihydrofolate reductase (FolA, EC 1.5.1.3). We selected FolA for the experiments due to its small size (159 amino acids) and essential role in the E. coli metabolism (Benkovic et al., 1988); also, the long-term protein sequence identity between FolA orthologs (~32%, see Figure 1a) is similar to other analyzed enzymes (Figure 2a). Following a recently described strategy (Kelsic et al., 2016), we used the Multiplex Automated Genome Engineering (MAGE) approach (Wang et al., 2009) to introduce every possible amino acid substitution at each FolA site in E. coli. To evaluate the relative fitness effects of protein substitutions we measured the growth rate of strains containing each protein variant compared to the ‘wild type’ (WT) strain into which substitutions were introduced. Relative growth rates were measured in parallel by performing growth competition experiments between the pooled mutants. Amplicon sequencing of the folA gene was then used to measure the relative changes of mutant and WT abundances as a function of time (see Materials and methods, Supplementary file 3).

Using the MAGE growth measurements in E. coli, we investigated the patterns of long-term sequence divergence at protein sites with different fitness effects of amino acid substitutions. Specifically, we sorted FolA protein sites into several groups according to their experimentally measured average fitness effects (Figure 4—figure supplement 1), and explored sequence divergence for sites within each fitness group (Figure 4a, different colors). We evaluated sequence identity between FolA orthologs across divergence times using all pairwise comparisons between ~300 orthologous sequences from the COG database (Galperin et al., 2015). Although, as expected, sites with stronger fitness effects diverged more slowly, our analysis revealed interesting differences in temporal divergence patterns for sites with small and large fitness effects. For sites in the least deleterious fitness group (Figure 4a, blue) we observed, similar to the global sequence identity, a substantial decrease (~10 fold, see Equation 5 in Materials and methods) in mutual divergence rates after ~1.5 billion years of evolution. Notably, even for FolA sites with mild fitness effects, sequence identity remains above 25% at long divergence distances. In contrast to sites with mild fitness effects, sites with the most deleterious mutations (Figure 4a, black) displayed a much slower, but approximately constant average divergence rate throughout evolutionary history. This pattern suggests that, in contrast to divergence at sites with small fitness effects, the divergence at sites with large effects is not yet close to saturation.

Figure 4 with 2 supplements see all
Sequence divergence of protein sites with different fitness effects of mutations measured in E. coli.

(a) The divergence of sequence identity for FolA protein sites with different average fitness effects of mutations (measured in E. coli) is shown using different colors. The average sequence identities were calculated using bacterial FolA orthologs available in the COG database (Galperin et al., 2015); divergence times were estimated using bacterial 16S rRNA sequences (see Materials and methods). Error bars represent the S.D. of sequence identity in each bin. (b) Similar to panel (a), but for the sequence divergence at protein sites of the E. coli translation initiation factor InfA. (c) The probability that protein sites in FolA orthologs (upper panel) and InfA orthologs (lower panel) are occupied by identical amino acids as a function of the average mutant fitness (measured in E. coli) at the corresponding protein sites. The probability represents the fraction of phylogenetically independent pairs of orthologs in which sites are identical at long divergence times (2 ± 0.25 billion years for FolA, and 2.5 ± 0.25 billion years for InfA). Error bars represent the S.E.M. across sites.

https://doi.org/10.7554/eLife.39705.017

To assess the generality of the FolA results we used another dataset (Kelsic et al., 2016), obtained using MAGE, of fitness values for all possible amino acids substitutions in the E. coli translation initiation factor InfA (Figure 4b). Consistent with the relatively higher level of sequence conservation of InfA, we observed lower average mutant growth rates and lower rates of sequence divergence in each fitness group. Nevertheless, the long-term divergence patterns were qualitatively similar between the two proteins. For sites in the least deleterious InfA fitness group (Figure 4b, blue), we observed a substantial decrease in the divergence rate after ~2 billon years of evolution. In contrast, sites with strongest fitness effects (Figure 4b, pink) displayed a slower but approximately constant divergence rate.

Because the fitness effects of mutations at a protein site may change in evolution (Lunzer et al., 2010; Chan et al., 2017), it is interesting to investigate how fitness effects measured in one species, such as E. coli, correlate with the site conservation in other species at the divergence limit. To explore this question, we calculated the probability that a protein site is occupied at long evolutionary times (~2 billion years for FolA and ~2.5 billion years for InfA) by the same amino acid in phylogenetically independent lineages (Figure 3a). We then investigated how this probability changes as a function of the average fitness effects of substitutions at the site measured in E. coli (Figure 4c). For both FolA and InfA, the probability that a protein site is identical at large divergence distances first increases, approximately linearly with increasing average fitness effects, and then begins to saturate for sites with large (>30% growth rate decrease) fitness effects. Thus, the fitness effects at a protein site correlate with the site’s conservation even after billions of years of evolution.

The sequence constraints revealed by our analysis likely arise due to the conservation of corresponding protein structures required for efficient catalysis and molecular function (Wilson et al., 2000; Watson et al., 2005). Therefore, in addition to sequence divergence we also investigated the long-term structural divergence of orthologous proteins with the same function. For this analysis we used >1000 orthologous pairs of enzymes sharing all 4 EC digits with known 3D structures (Berman et al., 2000) (see Materials and methods); structures of orthologous enzymes were aligned using the TM-align algorithm (Zhang and Skolnick, 2005). This analysis demonstrated that the average root mean square deviation (RMSD) between C-alpha atoms of the orthologous enzymes increases (Spearman’s r = 0.44, p-value<1e-20) with divergence time (Figure 5a). Nevertheless, the C-alpha RMSD between orthologs rarely increases beyond 3 Å, even at long evolutionary distances. Consistent with sequence evolution (Figure 1), we also observed a substantial decrease in the rate of structural divergence after ~1.5 billion years of divergent evolution.

Long-term structural evolution of enzymes with the same molecular function.

(a) The pairwise C-alpha root mean square deviation (RMSD) as a function of the divergence time between pairs of orthologs annotated with the same EC number. RMSD values were calculated based on structural protein alignments using the TM-align algorithm (Zhang and Skolnick, 2005). Gray dots represent pairs of considered orthologs, boxes indicate the median and 25–75 RMSD percentiles for the corresponding divergence times, the vertical lines indicate the 5–95 percentiles, and the red line shows the moving average of the data. (b) Long-term divergence of sequence identity of protein sites located at different distances to enzymes’ active site residues. In this analysis we considered the same species and enzymatic activities used to explore the global sequence divergence (Figure 1 and Figure 1—figure supplement 1); the average sequence identities within each distance shell (shown using different colors) were calculated across all pairs of orthologs annotated with the same EC number (see Materials and methods). Error bars represent the S.E.M. across ortholog pairs.

https://doi.org/10.7554/eLife.39705.020

Only a small fraction of all enzyme residues forms an active site and directly participates in catalysis. Therefore, we investigated next how the sequence divergence depends on the spatial proximity of protein positions to active site residues. It was previously demonstrated that evolutionary rates of amino acid substitutions correlate with protein sites’ spatial distance to catalytic residues (Jack et al., 2016). Notably, differences in short-term evolutionary rates do not imply the existence of a divergence limit, nor do they inform how site-specific divergence patterns correlate with sites' structural locations. Thus, the main goal of our analysis was to investigate the temporal patterns of the long-term divergence, and the effective divergence limit for sites at various distances to the active site. We considered catalytic site annotations available from the Protein Data Bank (Berman et al., 2000), UniProt-KB (UniProt Consortium, 2015) and the Catalytic Site Atlas (Porter, 2004) and quantified the average sequence divergence for sites at various distances from catalytic residues (see Materials and methods, Figure 5b). We based this analysis on the same set of enzymatic activities used to study global sequence divergence (Figure 1 and Figure 1—figure supplement 1). Although, as expected, residues close to the active site were the most highly conserved (Jack et al., 2016; Halabi et al., 2009), even distant sites displayed a substantial divergence limit (~40%) at long evolutionary distances. This result suggests that the spatial constraints required for specific molecular function usually propagate throughout the entire protein structure and significantly limit the long-term divergence even at sites distant from catalytic residues.

Finally, we investigated how various degrees of functional conservation affect the long-term divergence between orthologs. To that end, we compared the long-term sequence and structural similarities of enzymes sharing their full EC classification to those sharing only the first three digits of their EC classification (Figure 6a and b); for this comparison we only used orthologs from species with divergence times > 2 billion years (see Materials and methods). In contrast to enzymes sharing all four EC digits, conservation of the first three digits indicates only a general class of substrates or cofactors (Bairoch, 1999). This analysis revealed significantly lower long-term sequence identities (27% vs. 37% identity, Mann-Whitney p-value<10−20) and structural similarities (2.4 vs. 1.8 Å RMSD, P-value 2 × 10−18) between orthologs sharing only partial EC numbers. Notably, orthologs sharing only the first three EC digits are still substantially more conserved, both in sequence and structure (p-values<10−20), than pairs of enzymes with the same structural fold, but unrelated enzyme activities (i.e. enzymes sharing no digits in the EC classification) (Figure 6a and bDawson et al., 2017).

Effect of functional specificity on long-term sequence and structural divergence between orthologs.

(a) Sequence identities between orthologous pairs of enzymes that diverged over two billion years ago. The long-term sequence divergence between pairs of orthologs sharing the same EC number (gray, n = 272) is compared to the divergence between pairs only sharing the first three digits of their EC numbers (red, n = 265), that is, enzymes conserving only a general class of substrates or cofactors. The results are based on enzyme COGs for the 22 species used to analyze global sequence divergence (Supplementary file 1). Blue points show sequence identities between pairs of proteins with the same structural fold (Dawson et al., 2017) but unrelated enzyme activities, that is, activities sharing no digits in the EC classification (n = 298, see Materials and methods). The right blue Y axis represents the average number of amino acid types accepted per protein site during long-term protein evolution. (b) Similar to panel a, but showing the corresponding C-alpha structural divergence (RMSD) between protein pairs. (c) Long-term sequence identities between orthologous enzyme pairs at the same levels of structural similarity. Results are shown for pairs of enzymes sharing their full EC classification (gray dots), or only sharing the first three digits of their EC classification (red dots). In all panels: *(p<0.05), **(p<1e-4), ***(p<1e-10) for the Mann-Whitney test.

https://doi.org/10.7554/eLife.39705.022

We also investigated the sequence constraints at the same level of protein structural divergence for proteins with different degrees of functional conservation. To that end, we calculated the sequence identity between orthologs, sharing either their full or partial EC numbers, at different bins of long-term structural similarity (Figure 6c). Interestingly, we observed that even at the same level of C-alpha RMSD divergence, orthologs sharing full EC numbers usually have higher levels of sequence identity compared to orthologous pairs with the same level of structural divergence but sharing only three EC digits. This result indicates that conservation of molecular function constrains sequence divergence even beyond the requirement to maintain a specific spatial structure.

Discussion

Our analysis demonstrates that, in contrast to proteins with the same fold (Rost, 1997), the requirement to strictly conserve the same molecular function significantly limits the long-term sequence and structural divergence of protein orthologs. Although we confirmed the result by Povolotskaya and Kondrashov (2010) that ancient protein orthologs are still diverging from each other, our study reveals that the rate of this divergence becomes increasingly slow for orthologs that strictly conserve their function. Even a slight relaxation of functional specificity, for example from full to partial EC conservation (Figure 6a and b), leads to substantially more pronounced long-term sequence and structural divergence. Similarly, a substantial sequence identity between homologous restriction endonucleases is usually limited to isoschizomers, that is, proteins specific to the same target DNA sequence (Pingoud et al., 2014).

We believe that the observed divergence patterns can be explained by the following mechanistic model. Proteins with the same molecular function usually conserve the identity of their chemical and biological substrates and interaction partners. This conservation leads to functional pressure to closely preserve the spatial positions and dynamics of key protein residues necessary for efficient catalysis and function (Lehninger et al., 2013). In turn, the requirement to continuously preserve structural properties and functional dynamics of key protein residues likely imposes a strict requirement to preserve the overall protein structure, in other words, structural optimality necessary for efficient and specific protein function. The evolutionary pressure to preserve structural optimality (to within <3 Å C-alpha RMSD) required for a given molecular function leads, in agreement with the results by Chothia and Lesk (1986) and others (Gilson et al., 2017), to substantial levels of overall sequence conservation and the observed divergence limit. Our analysis further demonstrates that less than four amino acid types are accepted, on average, per site for proteins strictly conserving their molecular function.

We note that the observed conservation may reflect both the impact of amino acid substitutions on protein activity, due to changes in equilibrium positions and dynamics of protein residues, and on protein abundance, due to changes in overall protein stability (Bershtein et al., 2015; Adkar et al., 2017). Nevertheless, direct and comprehensive biochemical experiments demonstrated that the deleterious effects of mutations primarily arise from changes in specific protein activity rather than decreases in protein stability and cellular abundance (Firnberg et al., 2016). Our results support this model, demonstrating that conservation of functional specificity imposes substantially more stringent long-term sequence constraints than simply conservation of protein folds and protein stability. Indeed, while every folded protein should maintain its stability, our results demonstrate that maintaining stable folds does not constrain long-term divergence to more than 10-15% sequence identity, which is substantially less than the requirement for maintaining specific molecular functions (~40% sequence identity). 

The presented results demonstrate that only about a third of the sequence conservation between distant orthologs with the same molecular function can be attributed to universally conserved protein sites, that is, sites occupied by identical amino acids in almost all lineages. In contrast, we found that different protein sites are usually identical between orthologs from different lineages. This result is likely due, at least in part, to the epistatic nature of protein sequence landscapes, where mutations that are neutral in one lineage are often prohibitively deleterious in another (Breen et al., 2012; Lunzer et al., 2010). In the context of the aforementioned divergence model, the evolution of mitochondrial ribosomal proteins in eukaryotes (Figure 1—figure supplement 6) provides an interesting example, suggesting that orthologs’ divergence can be substantially accelerated by co-evolution with their interaction partners or relaxation of selection pressures.

Our experimental and computational analyses also delineate two distinct stages of the long-term divergence of orthologs with the same molecular function. During the first 1–2 billion years of divergence, substitutions at protein sites with mild fitness effects lead to a substantial (40–60%) decrease in sequence identity. After the first stage, divergence at these sites effectively saturates. The saturation at sites with small fitness effects, combined with very slow divergence rate at sites with large fitness effects (Figure 4), leads to a substantially slower sequence and structural divergence during the second stage. Interestingly, as a consequence of this slowdown, for the past billion years there has not been a substantial decrease in sequence and structural similarity between ancient orthologs with the same molecular function. Further analyses of biochemical, biophysical and cellular constraints will reveal how various structural and functional properties influence proteins’ long-term evolution, and how protein functional efficiency may be compromised by deleterious mutations (Shendure and Akey, 2015).

Materials and methods

Key resources table
Reagent type
(species) or resource
DesignationSource or referenceIdentifiersAdditional
information
Strain, strain
background
(Escherichia coli EcNR2)
MG1655, bla, bio-,
λ-Red+, mutS-::cmR
PMID: 19633652Addgene #26931
Sequence-based
reagent
90 bp DNA
oligos with
phosphorothioated bases
This paperSee Supplementary file 4100 nmole DNA
Plate oligo,
Integrated DNA
Technologies
Commercial assay or kitMiseq Reagent Kit V2IlluminaMS-102–2002
Commercial assay or kitsybr greenThermoFisherS7567
Commercial assay or kitQubit HS DNA kitThermoFisherQ32854
Commercial assay or kitQ5 Hot Start High-
Fidelity Mastermix
NEBM0494S
Commercial assay or kitDNA clean and
concentration kit 5
Zymo ResearchD4013
Commercial assay or kitillustra bacteria
genomicPrep Mini
Spin kit
GE life sciences28904259
Commercial assay or kitAgilent DNA 1000 kitAgilent Genomics5067–1504
Software, algorithmSeqPrep v1.1John St. Johnhttps://github.com/jstjohn/SeqPrep
Software, algorithmBowtie2PMID: 22388286
Software, algorithmPerl scripts to count
mutant reads
This paperhttps://github.com/platyias/count-MAGE-seq
(copy archived at https://github.com/elifesciences-publications/count-MAGE-seq).
OtherTurbidostat for growth
competition assay
PMID: 23429717

Considered enzyme activities and corresponding protein orthologs

Request a detailed protocol

We selected for analysis the sequences annotated in UniProt (UniProt Consortium, 2015) with EC numbers associated with the following metabolic pathways (defined in the KEGG database) (Kanehisa et al., 2016): Glycolysis and gluconeogenesis, pentose phosphate pathway, TCA cycle, purine metabolism, pyrimidine metabolism. Using the protein sequences from 22 diverse organisms (Supplementary file 1) we constructed clusters of orthologous groups (COGs) using the EdgeSearch algorithm (Kristensen et al., 2010). Following previous studies, we considered any two proteins from different species in the same COG as orthologs (Tatusov et al., 1997). COGs were obtained using the COGsoft software (Kristensen et al., 2010), starting from an all-against-all psi-blast (Altschul et al., 1997) search, setting the database size at 108, and using a maximum considered E-value of 0.1. To obtain the largest number of likely orthologs we did not apply a filter on low complexity or composition-based statistics. Only proteins sharing the same EC number and assigned to the same COG were compared, and only COGs with sequences in 10 or more of the 22 species were used.

In order to exclude proteins clearly showing evidence of Horizontal Gene Transfer (HGT), we constructed a maximum likelihood phylogenetic tree of the 12 prokaryotes considered in our analysis using a concatenated alignment of marker genes (Wu and Scott, 2012). The species tree was then manually compared to the individual trees of the prokaryotic sequences sharing the same molecular function within each COG; COG-specific trees were built using the GAMMA model of amino-acid substitution implemented in the RAxML software (Stamatakis, 2014). Proteins that showed clear differences in tree topologies, suggesting HGT, were excluded from further analysis. Ancient gene duplications, that is, duplications occurring prior to the divergence between considered species, often lead to cases in which enzymes in the same COG but from different species have diverged for longer than the corresponding species’ divergence times; thus, we did not consider COGs with tree topologies showing evidence of ancient gene duplications. Ancient gene duplications were defined as those occurring prior to the last common ancestor of 3 or more of the 22 species considered in the analysis.

The same procedure was used to select non-enzymatic COGs for analyses (Figure 1—figure supplement 5). However, in this case we only considered COGs for which none of the proteins were annotated in UniProt with metabolic EC numbers. Naturally, UniProt functional annotations for non-enzymes vary in terms of their source and format. Therefore, it is difficult to ascertain the degree of functional specificity and conservation between non-enzymatic orthologs. To address this, we manually checked that the molecular functions associated with proteins in the same COG were related, although we could not ascertain perfect conservation of molecular function.

Models of long-term protein sequence evolution

Request a detailed protocol

Global sequence identities for pairs of proteins annotated with the same molecular function in the same COG were calculated using pairwise alignments with ClustalW2 (Larkin et al., 2007). Sequence identity was computed as the number of identical sites at aligned positions, divided by the total number of aligned sites, excluding gaps. Divergence times between organisms were obtained from the TimeTree database (Hedges et al., 2006) (November, 2015) and used as a proxy for protein divergence times; in the analysis we used the mean divergence times across studies listed in the database. Divergence times between bacteria and archaea were set to 4 billion years based on current estimates for the occurrence time of their Last Common Ancestor (Sheridan et al., 2003; Battistuzzi et al., 2004) and existing evidence of an early origin of life on Earth (Bell et al., 2015). It is likely that ancient eukaryotic genes originated through episodic endosymbiotic gene transfer events and vertical inheritance from bacterial and archaeal genomes (Ku et al., 2015; Thiergart et al., 2012). Because of the discrete nature of such transfer events, the vast majority of individual prokaryotic-eukaryotic orthologous pairs are likely to have diverged from each other long before the origin of eukaryotes (1.8 billion years ago [Parfrey et al., 2011]); specifically, because most ancient prokaryotic species would not have transferred genes to eukaryotes. Thus, based on the median divergence time between the considered prokaryotes (~4 billion years, Supplementary file 1), divergence times between eukaryotes and prokaryotes were set in our analyses at 4 billion years. The results presented in the paper remain insensitive to the exact value of this divergence estimate (within the 3–4 billion year interval). Based on the recently proposed affiliation of eukaryotes and members from the Lokiarchaeota (Spang et al., 2015), divergence times between S. solfataricus and eukaryotes were set at 2.7 billion years, that is, the estimated age of the TACK superphylum (Guy and Ettema, 2011; Betts et al., 2018).

In order to study the long-term divergence patterns of orthologs, we only used COGs containing pairs of orthologs with at least five different divergence times distributed across 4 billion years. Sequence divergence data were fitted with models 1 to 3 using the least-squares minimization algorithm implemented in the MATLAB R2017a fitnlm function (The MathWorks, Inc, Natick, MA). The best fits of model 1 and model two were compared using the F-test. To test whether the conservation of molecular function limits protein sequence divergence, the minimum sequence identity parameter in model 2 (Y0, from Equation 2) was compared, for each enzymatic activity, to the average global sequence identity between unrelated protein pairs using the Wald test.

To investigate the effect of the uncertainty of divergence times’ estimates, we repeated the analysis of the 64 enzymatic activities while randomly assigning either the maximum or minimum value of the divergence times between lineages reported in the TimeTree database. This analysis was performed for a total of 1000 independent assignment runs. Across the independent assignment runs, the expected long-term sequence identity between orthologs was higher than 25% for at least 90% of enzymes (based on model 2), and the projected sequence identity after 7.8 billion years was above 25% (based on model 3) for at least 75% of enzymes (Figure 1—figure supplement 7).

To assess the effect of computational functional annotations on the observed divergence results, we repeated the analysis using only sequences with experimentally validated molecular functions (Figure 1—figure supplement 3). To keep only sequences with validated molecular functions, we manually reviewed published references for enzyme annotations in the BRENDA database (Chang et al., 2015), and discarded any functional assignments that were based exclusively on computational or high-throughput studies. After filtering for the experimentally validated annotations, we only considered EC numbers corresponding to pairs of orthologs with at least four different divergence times distributed across 4 billion years.

Calculation of the divergence rate

Request a detailed protocol

Based on Model 3, we determined the divergence rate, that is, the rate of the decrease in sequence identity per time, at a given divergence time t by solving for the derivative of Equation 3 with respect to time:

(5) dydt=d100*R0*tα+1-αdt= -100*R0R0*tα+1-α-1

where y represents global sequence identity, t represents divergence time, R0 represents the average substitution rate, and α represents the shape parameter of the gamma distribution.

Equivalency between model two and a poisson divergence model with allowed back substitutions

Request a detailed protocol

In the Jukes-Cantor model of nucleotide divergence (Tajima and Nei, 1984; Jukes and Cantor, 1969), the expected number of substitutions per site (δ) between two sequences after a divergence time t from a common ancestor is given by:

(6) δ=-a-1aln1-aa-11-y

where y is the proportion of identical sites and a is the number of allowed nucleotide types (usually 4). The same model can be applied to the divergence of protein sequences (Yang et al., 2000; Ota and Nei, 1994), by setting a to the number of allowed amino acid types per protein site. Furthermore, δ = 2λt, where λ represents the substitution rate per site per unit time, which is assumed to be equal across all sites. Substituting δ, and solving the above equation for y yields:

(7) y= 1a+1-1aexp-2λaa-1t

which is mathematically equivalent to model 2 Equation 2, with R0=2λaa-1, and Y0 =1a. Thus, Y0 can also be interpreted as the inverse of the average number of amino acids accepted per protein site during protein evolution.

FolA competition experiment in E. coli

Request a detailed protocol

To perform competition experiments we used the EcNR2 strain derived from E. coli K12 MG1655. Mutagenesis was performed using Multiplex Automated Genomic Engineering (MAGE), as previously described (Wang et al., 2009). 90 bp DNA oligomers were designed around each folA codon using the MG1655 wild type sequence as reference (Supplementary file 4). For each codon, all possible nucleotide variants were synthesized. To avoid simultaneous mutations of multiple codons, cells were transformed targeting ten consecutive codons at a time. After four rounds of electroporation, cells were recovered and pooled together at approximately the same concentration based on cell counts. Two competition growth experiments were carried out, one for each half of the protein. For the competition experiments, cells were grown in LB media in a turbidostat while maintaining constant volume and cell density. Samples were taken every 2 hr for a period of 16 hr, spun down, washed in PBS, spun down again and stored at −20°C until all samples were collected. For each competition, the corresponding FolA region was amplified through PCR while assigning a specific DNA barcode for each time point. PCR products were then pooled and paired-end sequenced using the MiSeq Reagent Kit two from Illumina. Sequence reads were deposited to the SRA database with accession number: SRP152339.

To determine, at each time point, the abundance of each mutant relative to wild type, we joined paired-end reads using SeqPrep (v 1.1) and aligned the joined reads to the folA gene sequence using Bowtie2 (Langmead and Salzberg, 2012). We then counted the number of reads per mutant using a custom script (Plata, 2018). Reads with more than a single mutated codon were discarded. Counts were median-normalized to control for noise due to mutagenesis performed in batches of 10 codons. At each time point we calculated the ratio Rt of mutant to wild type (WT) reads. In exponential growth, the growth rate difference between a given mutant and WT was calculated based on the slope of lnRτ as a function of time:

lnRt=mi-mwt*t+ln(R0)

where mi and mwt represent the mutant and WT growth rates, respectively. Growth rate differences were calculated only for mutants with at least five time points with 20 or more reads. Relative growth rates were calculated by dividing the slopes in the equation above by the number of e-fold increases given the average dilution rate of the turbidostat (1.37/h).

To calculate a single value characterizing the effect of all possible mutations at a protein site, we first averaged the relative growth rates of mutants resulting in the same amino acid change. We then calculated the average fitness effect of mutations at each protein site by averaging across 20 possible amino acids substitutions (Supplementary file 3).

To estimate the sensitivity of our results to sequencing errors, we calculated the average fitness effect of substitutions at each FolA site using the relative growth rates of mutant strains carrying only 32 mutated codons selected at random out of 64 possible codons. We observed a high correlation (Pearson’s r: 0.95, p-value<1e-20, Figure 4—figure supplement 2) between the average growth rate effects at each site calculated using two non-overlapping subsets of 32 codons. As expected, nonsense mutations and substitutions in the folA start codon had substantially stronger average effects on growth rates compared to other substitutions (26% versus 4% slower growth than WT, respectively. Mann Whitney U, p-value<10−20). Also, the relative growth rates due to synonymous codon substitutions were usually very mild (0.2% higher growth compared to WT); 97% of synonymous substitutions had growth effects of less than 3%.

Contribution of different sites to the divergence limit

Request a detailed protocol

In order to identify phylogenetically independent pairs of species, we aligned the 16S rRNA gene sequences of bacterial species having orthologs annotated with the target 30 EC numbers (Figure 3—figure supplement 1). 16S rRNA sequences were obtained from the GreenGenes database (DeSantis et al., 2006) (October, 2016). We then built maximum likelihood phylogenetic trees based on the 16S alignments using RAxML (Stamatakis, 2014). Next, we used the Maximum Pairing Problem approach byArnold and Stadler (2010) to find the largest number of edge-disjoint pairs of species with 16S rRNA genetic distances corresponding to >2 billion years of divergence. Divergence times were estimated from the 16S genetic distances based on the linear regression of literature reported divergence times (Hedges et al., 2006) (Supplementary file 1). The F84 model of nucleotide substitution implemented in the phylip package (Felsenstein, 2005) was used to compute the genetic distances. Using the 16S alignment data, we calculated the probability that a protein site was identical across independent lineages. The probability was calculated as the fraction of orthologous pairs from phylogenetically independent species pairs with identical amino acids at the site. The amino acid identities at a given site were obtained based on the multiple sequence alignment of all orthologs associated with each EC number, obtained using ClustalW2 (Larkin et al., 2007). A similar procedure was applied to analyze FolA and InfA orthologs from the COG database (Figure 3c).

To investigate the divergence of sites with different fitness effects, we used sequences of FolA and InfA bacterial orthologs from the COG database (Galperin et al., 2015). The FolA orthologs annotated with the same EC number in UniProt (n = 311) and the InfA orthologs annotated with the same KEGG Orthology (KO) number in KEGG (n = 514) were used to build multiple sequence alignments with ClustalW2 (Larkin et al., 2007). Divergence times were estimated from the 16S genetic distances as described above. Within each divergence bin (Figure 4a,b), sequence identities of sites with different average fitness effects (represented by different colors in Figure 4a,b) were averaged across all pairs of orthologs at a given divergence time.

Analysis of global protein structural evolution

Request a detailed protocol

To study the divergence of protein structures as a function of time, we obtained PDB codes for all proteins associated with EC numbers in the BRENDA database (Chang et al., 2015). We then selected for the analysis species with experimentally solved enzyme structures for at least 10 different EC numbers. Psi-blast searches with a conservative E-value cutoff of 10−6 were used to identify orthologs (defined as bi-directional best hits) in the selected species. The 3D structures of orthologous pairs, annotated with the same EC number, were aligned using the TM-align program (Zhang and Skolnick, 2005) to obtain the C-alpha RMSD values. Pairs of proteins were not considered if more than 70% of the residues of the shortest protein could not be structurally aligned. We also removed from the analysis pairs of structures with flexibility between domains, as they could result in large RMSD values despite significant structural similarity. To identify such proteins we used the FATCAT (Ye and Godzik, 2003) software to perform flexible structural alignments of all structure pairs. We then filtered the structural pairs that were split into two or more domains by the FATCAT alignments.

Analysis of the enzyme active sites

Request a detailed protocol

To analyze divergence as a function of active site distance we used protein sequences associated with the 64 EC numbers and 22 species considered in Figure 1 and Figure 1—figure supplement 1. To that end, PDB (Berman et al., 2000) was searched for homologous sequences annotated with the same enzymatic activities and with known 3D structures. Annotations of active site residues for the corresponding structures were obtained from the Catalytic Site Atlas (Porter, 2004), PDB and UniProt-KB (UniProt Consortium, 2015). For each PDB structure with available active site information, protein sites were then stratified into different layers according to the distance between their alpha carbons and the centroid of the active site residues. Each pair of orthologs was then aligned using ClustalW2 (Larkin et al., 2007) with a homolog in PDB annotated with the same activity and with defined distance layers around the active site; the PDB sequence with the highest sequence identity to either member of the pair was used for the alignment. Sequence identities for different layers were calculated based on the structural positions in the corresponding PDB reference sequences.

Comparison of pairs of enzymes with the same structural folds

Request a detailed protocol

We used structural classifications of protein domains from the CATH database (v4.2.0) (Dawson et al., 2017). For structural comparisons, we only considered PDB structures with a single classified domain per chain. Protein pairs classified in CATH in the same homologous structural superfamily were considered as having the same fold. To obtain pairs of proteins in the same fold but with different functions, we only considered PDB structures annotated with different EC numbers in BRENDA. For this analysis we randomly selected 300 pairs of structures with the same fold that do not share any digits of their EC classification.

References

  1. 1
  2. 2
  3. 3
  4. 4
  5. 5
  6. 6
  7. 7
  8. 8
  9. 9
  10. 10
  11. 11
  12. 12
  13. 13
  14. 14
  15. 15
  16. 16
  17. 17
  18. 18
  19. 19
  20. 20
  21. 21
  22. 22
    PHYLIP (Phylogeny inference package) , version 3.6
    1. J Felsenstein
    (2005)
    University of Washington, Seattle.
  23. 23
  24. 24
  25. 25
  26. 26
  27. 27
  28. 28
  29. 29
  30. 30
  31. 31
  32. 32
    Evolution of Protein Molecules
    1. TH Jukes
    2. CR Cantor
    (1969)
    In: H. N Munro, editors. Mammalian Protein Metabolism. Academic Press. pp. 21–132.
  33. 33
  34. 34
  35. 35
  36. 36
  37. 37
  38. 38
  39. 39
  40. 40
    Lehninger Principles of Biochemistry (6th Edition)
    1. AL Lehninger
    2. DL Nelson
    3. MM Cox
    (2013)
    New York: W.H. Freeman.
  41. 41
  42. 42
  43. 43
  44. 44
  45. 45
  46. 46
  47. 47
  48. 48
  49. 49
  50. 50
  51. 51
  52. 52
  53. 53
  54. 54
  55. 55
  56. 56
  57. 57
  58. 58
  59. 59
  60. 60
  61. 61
  62. 62
  63. 63
  64. 64
  65. 65
  66. 66
  67. 67
  68. 68
    Codon-substitution models for heterogeneous selection pressure at amino acid sites
    1. Z Yang
    2. R Nielsen
    3. N Goldman
    4. AM Pedersen
    (2000)
    Genetics 155:431–449.
  69. 69
    Computational Molecular Evolution
    1. Z Yang
    (2006)
    Oxford: Oxford University Press.
  70. 70
  71. 71
  72. 72
    Evolutionary Divergence and Convergence in Proteins
    1. E Zuckerkandl
    2. L Pauling
    (1965)
    In: V Bryson, H. J Vogel, editors. Evolving Genes and Proteins. Academic Press. pp. 97–166.

Decision letter

  1. Nir Ben-Tal
    Reviewing Editor; Tel Aviv University, Israel
  2. Diethard Tautz
    Senior Editor; Max-Planck Institute for Evolutionary Biology, Germany
  3. Nir Ben-Tal
    Reviewer; Tel Aviv University, Israel

In the interests of transparency, eLife includes the editorial decision letter, peer reviews, and accompanying author responses.

[Editorial note: This article has been through an editorial process in which the authors decide how to respond to the issues raised during peer review. The Reviewing Editor's assessment is that all the issues have been addressed.]

Thank you for submitting your article "Molecular function limits divergent protein evolution on planetary timescales" for consideration by eLife. Your article has been reviewed by three peer reviewers, including Nir Ben-Tal as the Reviewing Editor and Reviewer #1, and the evaluation has been overseen by Diethard Tautz as the Senior Editor. The other reviewers remain anonymous.

The Reviewing Editor has highlighted the concerns that require revision and/or responses, and we have included the separate reviews below for your consideration. If you have any questions, please do not hesitate to contact us.

Summary:

The manuscript analyzes clusters of orthologue enzymes of the same function, defined as sharing the same 4-digits EC number, to examine the limitation of function preservation on evolutionary divergence. The results seem to indicate that function confines divergence, which is not surprising. The potential novelty here is the time dependence of the process, where phylogenetic time of species divergence is used as natural time scale "clock". The data shows that proteins diverse from each other rapidly during the first ~1.5 billion years, but saturate after that.

Opinion: The manuscript addresses an important and interesting topic, but it should be revised significantly to clarify what exactly was computed and why, what was measured and why, and how the results are related to previous studies. In view of the amount and nature of the issues raised, the authors should consider withdrawing the manuscript.

Major concerns:

(pretty much copy-pasted from the individual reports)

1) This study poses the question about the extent to which function constrains sequence divergence. The answer given is that for orthologous enzymes 25% sequence identity is the approximate limit. The major critique, is that the authors did not make an attempt to provide a clear mechanistic explanation for such phenomenon. Is it simply because a significant fraction of positions is invariant? Is it because there are various constraints on various sites? Is it because only a small subsets of amino acids are tolerated at some sites? Why not take several enzyme families with very thick alignments and analyze them position-by-position to shed light on this question? Without such deep analysis, the paper, while interesting in its general idea, reads superficial and underdeveloped, despite some experimental results.

2) The authors premise their study on the assertion made by Rost in a 1997 review that proteins of the same fold can diverge up to random sequence ID of 3-5%, i.e., that structural similarity does not impose constraint on sequence ID. As a consequence the authors attribute in their analysis all sequence identity constraints to functional similarity. However, this premise is highly problematic. While anecdotally one can find two proteins of similar fold with very low sequence ID, typically proteins with similar structure (not necessarily orthologues) have sequence ID at or above 25% – this is the essence of the famous cusp in sequence-structure relationship first reported by Lesk and Chothia in the eighties. In fact a careful analysis of structural/stability constraints on sequence divergence has been published in recent paper (Biophys J v.112, p.1350-65, 2017) where it was shown that in divergent evolution scenario protein structure is maintained up to 25-30% sequence ID and quickly deteriorates beyond that. Incidentally it is very close to 25% ID which the authors claim to stem from functional constraints. Furthermore, the above mentioned BJ paper presents a detailed analytical estimate of sequence divergence dynamics akin to exponential fits used in this PAPER. The authors should make themselves familiar with this theory and consider using it to fit their divergence curves rather than ad hoc exponential fits.

3) The high throughput experimental data is potentially interesting but its description is so cryptic and incomplete that it is very hard to assess what actually has been done there. In particular there is no assessment of noise in deep sequencing to assess fitness, of the effects of synonymous substitutions, etc. How is the data binned and why binning? The use of LB as a medium for the folA experiment is unfortunate because in LB DHFR is much less essential than in minimal media and therefore the results could be skewed by very permissive conditions. In addition the standing genetic variation is not taken into account – fitness effects could be a result of hitchhiking (i.e. originate from the variation of the background into which the mutations are introduced). These concerns aside, it is not entirely clear what does comparison with experiment tell us beyond the obvious that sites where deleterious fitness effects are greatest evolve more slowly. Was the sole purpose of the experiment to bin positions by their fitness effects? If so, what new insight is gained from slower divergence of sites with more profound fitness effects compared to faster divergence of sites that are more tolerant to substitutions in experiments. Isn't this effect expected? Maybe quantification of the effect may have some novelty, but such quantification has not been done. The obvious thing to do is to compare their alleged fitness effects with dN/dS assessment for the rates but that has not been done.

4) The word "protein" should be replaced with "enzyme" in the title. "We focused on enzymes because their molecular function is usually well defined," doesn't sound quite true. Reading the manuscript, it seems that the goal of this work is to investigate divergence of enzymes. The divergence of non-enzymes will most likely follow a different pattern and will generally be much lower than 25% identity. If the authors want to generalize to non-enzymes, they need to perform appropriate analyses. Currently, only one experiment is dedicated to a non-enzyme. It is not fair to the readers to generalize beyond the scope of analysis that was carried out.

5) The sentence in the Abstract "divergence rates of orthologous enzymes decrease substantially after ~1-2 billion years of independent evolution" is confusing. "Rate" means number of amino acid substitutions per site per unit of time. Moreover, the authors write in the Introduction "a divergence limit does not imply that the rate of amino acid substitutions slows down in evolution," thus contradicting the Abstract. The Abstract should be edited and the misleading sentence corrected. From the context of the paper is seems like by "divergence rate" the author mean the rate of decrease in pairwise sequence identity with time. This definition poses additional question. This "divergence rate" will always decrease due to multiple and back substitutions and it is trivial. The authors need to make it clear that what they mean goes beyond a trivial curvilinear relationship between time and sequence identity.

6) Further on this, evolutionary rate is much more than simply a count of changes in the identity of the amino acids. To measure it properly one needs to rely on the phylogeny. Indeed, since the study is based on several clusters of orthologous proteins, why not do just that? For example, why not use the dN/dS formalism?

7) "billion years of evolution, we observed a significant decrease in mutual divergence rates," see item 5 above. Also, shouldn't "significant" be quantified is some way?

8) "Species' divergence times were used to estimate the times of divergence between corresponding orthologous proteins": How reliable is this time estimate?

9) Discussion "The presented results demonstrate that, in contrast to proteins with the same fold [Rost, 1997], the requirement to maintain the same molecular function significantly constrains the long-term divergence of protein orthologs": This difference is surprising. Maybe it has to do with how changes were estimated in the two studies? To state it with confidence direct comparison is needed, where the exact same methodology is used.

10) Finally, the manuscript should be edited for clarity by a professional scientific editor.

Separate reviews (please respond to each point):

Reviewer #1:

Summary of paper:

The manuscript analyzes clusters of orthologue enzymes of the same function, defined as sharing the same 4-digits EC number, to examine the limitation of function preservation on evolutionary divergence. The results seem to indicate that function confines divergence, which is anticipated- perhaps also shown before. The potential novelty here is the time dependence of the process. The data shows (to the best of my understanding) that proteins diverse from each other rapidly during the first ~1.5 billion years, but saturate after that. The data further suggests (again, assuming that I understood correctly) that the rapid evolutionary phase mostly reflects random drift, i.e., amino acid sites that can freely mutate, and that the saturation value reflects their relative fraction of the entire enzyme sequence.

Opinion: The manuscript could be important but it should be revised significantly to clarify what exactly was computed and why, what was measured and why, and how the results are related to previous studies. I also believe the manuscript should be edited for clarity by a professional scientific editor.

Major issues:

1) Results "In the first model, all protein sites have equal and independent substitution rates; see Equation 1": Equation 1 obviously represents an unrealistic evolutionary model because it is well known that the evolutionary rates are not equal at all amino acid sites. For example, catalytic residues hardly change.

2) Where does Equation 2 come from? It also doesn't look particularly realistic.

3) Equation 3 is the commonly used evolutionary model I think. So why not start with it?

4) "Species' divergence times were used to estimate the times of divergence

between corresponding orthologous proteins": How reliable is this time estimate?

5) Figure 2: How is "fitness cost" defined?

6) "phylogenetically independent pairs of species": What does it mean and how is "independent" defined? Is the pair A-B in Figure 3A phylogenetically independent? And the pair D-H?

7) "Notably, the probability that a protein site is identical, and thus contributes to the divergence limit, first increases linearly with increasing fitness effects at the site, and then begins to saturate for sites with high (>30%) fitness effects (Figure 3B)." This is perhaps true for FolA but I do not see saturation for InfA.

8) The distributions of the probability of identical sites in FolA and InfA (Figures 4A and B) are very different, but they are discussed together without reference to the difference.

9) Discussion "The presented results demonstrate that, in contrast to proteins with the same fold [Rost, 1997], the requirement to maintain the same molecular function significantly constrains the long-term divergence of protein orthologs": I am surprised about this difference. Maybe it has to do with how changes were estimated in the two studies? To state it with confidence I would have liked to see direct comparison where the exact same methodology is used.

10) Evolutionary rate is much more than simply a count of changes in the identity of the amino acids. To measure it properly one needs to rely on the phylogeny. Indeed, since the study is based on several clusters of orthologous proteins, why not do just that?

Reviewer #2:

This study poses the question about the extent to which function constrains sequence divergence. The answer given is that for orthologous enzymes 25% sequence identity is the approximate limit. The major critique, is that the authors did not make an attempt to provide a clear mechanistic explanation for such phenomenon. Is it simply because a significant fraction of positions is invariant? Is it because there are various constraints on various sites? Is it because only a small subsets of amino acids are tolerated at some sites? Why not take several enzyme families with very thick alignments and analyze them position-by-position to shed light on this question? Without such deep analysis, the paper, while interesting in its general idea, reads superficial and underdeveloped, despite some experimental results.

1) I suggest replacing the word "protein" with "enzyme" in the title. "We focused on enzymes because their molecular function is usually well defined," doesn't sound quite true. Reading the manuscript, it seems that the goal of this work is to investigate divergence of enzymes. I am positive that the divergence of non-enzymes will follow a different pattern and will generally be much lower than 25% identity. If the authors want to generalize to non-enzymes, they need to perform appropriate analyses. Currently, only one experiment is dedicated to a non-enzyme. It is not fair to the readers to generalize beyond the scope of analysis that was carried out.

2) I think that the sentence in the Abstract "divergence rates of orthologous enzymes decrease substantially after ~1-2 billion years of independent evolution" is confusing. "Rate" usually means number of amino acid substitutions per site per unit of time. Moreover, the authors write in the Introduction "a divergence limit does not imply that the rate of amino acid substitutions slows down in evolution," thus contradicting the Abstract. I suggest to edit the Abstract and correct the misleading sentence. From the context of the paper is seems like by "divergence rate" the author mean the rate of decrease in pairwise sequence identity with time. This definition poses additional question. This "divergence rate" will always decrease due to multiple and back substitutions and it is trivial. The authors need to make it clear that what they mean goes beyond a +trivial curvilinear relationship between time and sequence identity.

3) "billion years of evolution, we observed a significant decrease in mutual divergence rates," see item 2 above. Also, shouldn't "significant" be quantified is some way?

4) It was not made clear why these experiments were performed and how they integrate with the rest of the study. Addition of these experiments seems preliminary and no confident conclusions follow. Was the sole purpose of the experiment to bin positions by their fitness effects? If so, what new insight is gained from slower divergence of sites with more profound fitness effects compared to faster divergence of sites that are more tolerant to substitutions in experiments. Isn't this effect expected? Maybe quantification of the effect may have some novelty, but such quantification has not been done.

Reviewer #3:

In this paper the authors aim to assess how does the requirement of conserved function (i.e. orthology) constraints sequence evolution. The author premise their study on the assertion made by Rost in 1997 review that proteins of the same fold can diverge up to random sequence ID of 3-5% i.e. that structural similarity does not impose constraint on sequence ID. As a consequence the authors attribute in their analysis all sequence identity constraints to functional similarity. To that end they analyze sequence divergence of orthologs using phylogenetic time of species divergence as natural time scale "clock". They observed saturating time dependencies of sequence divergence with time and conclude that the rate of sequence divergence decreases with decreasing divergence time. Further, they carried out a high throughput mutational experiment on fitness effects of substitutions in 2 genes – folA and InfA – in E. coli and find, perhaps unsurprisingly, that sites where mutational effects on fitness are stronger are less constrained in evolution.

While the paper contains some interesting bioinformatics observations and analysis I have serious concerns about its premise and interpretation of the results. Apparently some issues stem from author's apparent gaps in knowledge of modern literature on biophysical determinants of protein evolution.

1) The premise that structural similarity does not constraint sequence identity is highly problematic. While anecdotally one can find two proteins of similar fold with very small sequence ID typically proteins with similar structure (not necessarily orthologues) have sequence ID at or above 25% – this is the essence of the famous cusp in sequence-structure relationship first reported by Lesk and Chothia in the eighties. In fact a careful analysis of structural/stability constraints on sequence divergence has been published in recent paper (Biophys J v.112, p.1350-65 (2017) where it was shown that in divergent evolution scenario protein structure is maintained up to 25-30% sequence ID and quickly deteriorates beyond that. Incidentally it is very close to 25% ID which the authors claim to stem from functional constraints Furthermore, the above mentioned BJ paper presents a detailed analytical estimate of sequence divergence dynamics akin to exponential fits used in this PAPER. The authors should make themselves familiar with this theory and use it to for their divergence curves rather ad hoc exponential fits.

2) The high throughput experimental data is potentially interesting but its description is so cryptic and incomplete that it is very hard to assess what actually has been done there. In particular there is no assessment of noise in deep sequencing to assess fitness., of the effects of synonymous substitutions etc. How is the data binned and why binning? The use of LB as a medium for the folA experiment is unfortunate because in LB DHFR is much less essential than in minimal media and therefore the results could be skewed by very permissive conditions. In addition the standing genetic variation is not taken into account – fitness effects could be a result of hitchhiking (i.e. originate from the variation of the background into which the mutations are introduced).These concerns aside, it is not entirely clear what does comparison with experiment tell us beyond the obvious that sites where deleterious fitness effects are greatest evolve more slowly. The obvious thing to do is to compare their alleged fitness effects with dN/dS assessment of the rates but that has not been done.

3) The interpretation of the experimental fitness effects in terms of function is also questionable. The authors are apparently unaware of series of experimental works from Shakhnovich lab where determinants of fitness effects of mutations are addressed. In particular it has been shown using both point mutations and orthologous chromosomal replacements for folA gene (PLOS Genetics 2015 DOI:10.1371/journal.pgen.1005612) and adk gene (Nature Ecology Evolution, 2017 http://dx.doi.org/10.1038/s41559-017-0149) that fitness is determined by product of folded protein abundance A and activity kcat/KM. Mutations may affect stability and through that parameter A (by changing the balance between protein production and degradation, see Mol Cell v.49, pp133-44 (2013). Therefore interpretation of the experimental trends entirely in functional terms is not warranted.

A minor comment: The concept and metaphor of expanding protein universe has been introduced 10 years before Kondrashov's work in the paper "Expanding protein universe and its origin from biological big bang" PNAS 2002, v.99 pp. 14132-6.

[Editors' note: further revisions were suggested before publication, as described below.]

Thank you for resubmitting your work entitled "Molecular function limits divergent protein evolution on planetary timescales" for further consideration at eLife. Your revised article has been favorably evaluated by Diethard Tautz as the Senior Editor, a Reviewing Editor, and two reviewers.

As you can see from the reports below, the reviewers appreciated the revisions. However, there are still major outstanding issues. While some of these can be resolved by changes in the presentation, others are fundamental. We would strongly encourage you to address all of these prior to publication.

Reviewer #2:

I find that the authors did a thorough revision of the manuscript. At least now I think I understood the main conclusion of the paper. In enzymes and other proteins with very strong functional conservation, the number of different amino acids acceptable at a position is about 4. It is not because many sites are invariant and some are variable (not a dirty trick of average temperature in a hospital), but because most sites (except the invariant ones) are constrained to use a library of 3 to 5 different amino acids, not more than that. If this is not the bottom line, then the authors need to do better job at crystallizing their main claim and result.

If it is, I think it is a meaningful finding that could be explained better to the readers. I guess the second claim is that 3-5 amino acid limit is universal to all enzymes and (conserved!) non-enzymes. I do doubt (as in the original review) the validity of such a strong claim, which could be a result of the authors' bias in selecting families for their analysis. At least for non-enzymes, they selected most conserved proteins (like ribosomal proteins), so of course such selection is biased to get proteins that saturate easily in evolution. The authors try to justify this biased selection suggesting that it is difficult to find orthologs. But that statement by itself totally discredits this study. Why? Because if you cannot find your orthologs, wouldn't it mean that they already diverged beyond (lower than) your claimed 25% identity as the universal limit, and the author's conclusions do not apply to such families? Maybe for the enzymes too, the 25% limit is simply a reflection of the search methods the authors used to find orthologs, that fail to find more distant ones.

To improve the presentation and make this paper quite interesting (well, reviewers are not supposed to direct the study, but the authors seem passionate about their work and also seem rather inexperienced in both logical thinking and putting a paper together, so maybe this recommendation could be helpful), I would suggest to base it on two plots.

1) Between-protein variability.

The first plot is the histogram of the average estimated number of amino acids allowable per site (without invariant sites, or with) for enzyme families and other protein families. Well, the authors would have to try harder to find orthologs for proteins that manage to diverge below 25%, even my rotation students can do that. This histogram would be expected to have a maximum around 3 to 5, and for some variable proteins it could be around 10 or 15? Right? Then protein families with very low and very high numbers could be discussed with an attempt to give explanations about their uniqueness.

2) Within-protein variability.

The second plot is the histogram of estimated number of amino acids allowable per site within a protein family. The authors could either normalize to the average per protein, or select a bin with the maximum count from previous plot (let's say 4). These histograms can be averaged for all families with mean and SD showing for each bin. I would assume there will be a large count for invariant sites for enzymes (at 1), counts for 2 amino acids used, 3, 4, etc. Will this histogram have a single mode? Maybe around 4? Several modes? Will there be sites using more than 10 amino acids? How are these distributed in spatial structure and relative to the active site? Discussion of these details could be quite insightful and interesting.

If these authors do not wish to make these plots, since this review will be published, maybe someone else will, and we will learn something interesting as a result.

Currently, there are so many plots in this paper and many of them are not particularly helpful to get the point across quickly. Also I still find the usage of non-standard terms (like divergence rate) confusing, not necessary, and not insightful about the mechanism. Yes, the terms are defined, but they are just masking the reality, which is simpler: constrained usage of amino acids in positions of enzyme molecules. Not all 20, not even 10, but 4! Why not state and illustrate this clearly? Don't you agree that the impact of the paper will increase because it will be easier to understand?

The Abstract is still very poor and misleading. For instance, the statement "The effective divergence limit (>25% sequence identity) is not primarily due to multiple substitutions at the same sites" is completely wrong. According to the authors' explanations, this "effective divergence limit" is exactly due to multiple substitutions at the same site! If you have only 4 amino acids acceptable at a site, due to multiple substitutions sequence identity will saturate are 25%. Like in DNA. Why not write a precise and clear Abstract about what this work presents? Due to so many logical flaws in the authors' thinking and presentation (as illustrated above), I would be scared to publish this paper without a careful read. One statement at a time. And trying to figure out why the statement is wrong. If clearly not wrong, then move on to the next one.

And, finally, why not compare your results with this paper more thoroughly PMID: 27138088? Isn't it a bit similar? I guess I don't understand the meaning behind "to investigate the temporal patterns of the long-term divergence."

Do the authors still claim that the sites are saturated at the usage of 3 to 5 amino acids because the time passed was not sufficient to gain more changes? Or the time was enough and protein of the same function simply cannot tolerate additional amino acids well and still keep the function? Which one is right? I got an impression it was the latter. And then what I said at the beginning of this review holds. If it is the former, it needs to be convincingly justified.

If the authors disagree with my review, then I did not understand the paper, which is possible, and still suggests that the authors need to improve the presentation.

Reviewer #3:

This version of the manuscript is a significant improvement over the previous version in terms of added details (e.g. experimental procedure of folA mutagenesis re now described in sufficient detail to be understandable and/or potentially reproducible).

The authors attempted to address my (and other reviewers) main concern by presenting the analysis in new Figure 5 that shows that full orthologues (all EC numbers coincide) are more sequence-constrained than partial orthologues (3 EC numbers coincide). However this analysis fails to control for different structural divergence between full and partial orthologues.

Therefore, the most problematic aspect of the analysis – that authors attribute observed sequence conservation in diverging clades to conservation of function between orthologues has not been fully addressed. A clear alternative explanation that such conservation is explained by maintenance of structure and stability regardless of function still stands.

The authors used a truism by pointing out that proteins of the same FOLD (i.e. topologically similar arrangement of elements of secondary structure) can diverge to 10% or less sequence ID again citing an old work with data collected on very limited number of structures available at that time.

However, their full functional orthologues are much more similar structurally than proteins sharing the same fold. The correct control which was suggested in my initial review has not been done satisfactorily. Specifically, the authors should compare their sets of proteins with sets that have similar degree of structural similarity (measured as distribution of TM-scores) but different function. There are such examples where structures are quite similar (TM-score-wise) but functions differ significantly. Good examples of that kind ate TIM-barrels which are almost exclusively enzymes with wide variety of specificity and Igfold proteins again with very conserved structures but broadly diverged functional annotations.

If the authors find that conservation of structure in functionally divergent proteins imposes less sequence divergence constraints than same degree of structural conservation in orthologues – that will be a clear demonstration of additional constraints imposed by functional conservation which is the main message of this work. In the absence of such analysis the current data does not justify the conclusion.

Minor point: The authors severely misquote Firnberg et al., 2016. They say: "Nevertheless, direct and comprehensive biochemical experiments demonstrated that the deleterious effects of protein mutations primarily result from changes in specific protein activity rather than decreases in protein stability and cellular abundance [Firnberg et al., 2016]". In fact, the authors of Firnberg et al., 2016, say directly the opposite: '… These DFEs provide insight into the inherent benefits of the genetic code's architecture, support for the hypothesis that mRNA stability dictates codon usage at the beginning of genes, an extensive framework for understanding protein mutational tolerance, and evidence that mutational effects on protein thermodynamic stability shape the DFE..…" (cited from the Abstract of Firnberg et al., 2016). The authors misrepresent the main result of Firnberg et al., 2016: According to Figure 5B of Firnberg et al., 2016, the product of abundance and catalytic activity shapes fitness effects in TEM1, not abundance alone. This is exactly what is established in Bershtein et al., 2015 and Adkar et al., 2017, and shows that abundance (which is a function of stability) enters fitness landscape on equal footing with catalytic activity, i.e. there is as much selection for abundance (i.e. stability) as it is for kcat/KM and related measures of activity. This misinterpretation somewhat undermines the major premise of the present paper that there is separate selection for activity/function and stability/structure. In reality one cannot disentangle the two because fitness landscape depends on the function (i.e. the product) of the two factors.

https://doi.org/10.7554/eLife.39705.034

Author response

The manuscript analyzes clusters of orthologue enzymes of the same function, defined as sharing the same 4-digits EC number, to examine the limitation of function preservation on evolutionary divergence. The results seem to indicate that function confines divergence, which is not surprising. The potential novelty here is the time dependence of the process, where phylogenetic time of species divergence is used as natural time scale "clock". The data shows that proteins diverse from each other rapidly during the first ~1.5 billion years, but saturate after that.

Opinion: The manuscript addresses an important and interesting topic, but it should be revised significantly to clarify what exactly was computed and why, what was measured and why, and how the results are related to previous studies. In view of the amount and nature of the issues raised, the authors should consider withdrawing the manuscript.

We sincerely thank the reviewing editor and reviewers for their thoughtful suggestions and comments. Following the reviewers’ suggestions, we have substantially revised the manuscript to clarify the main results, their novelty and importance. Although the fact that functional conservation constrains protein divergence is indeed not surprising, our manuscript describes interesting long-term temporal patterns and the limits of this divergence for proteins with the same molecular function. We welcome the opportunity provided by this innovative peer review model to share with the readers our replies to the reviewers’ comments. We believe that this scientific dialogue will be an important contribution to understanding of long-term evolution of protein sequences and structures.

Major concerns:

(pretty much copy-pasted from the individual reports)

1) This study poses the question about the extent to which function constrains sequence divergence. The answer given is that for orthologous enzymes 25% sequence identity is the approximate limit. The major critique, is that the authors did not make an attempt to provide a clear mechanistic explanation for such phenomenon. Is it simply because a significant fraction of positions is invariant? Is it because there are various constraints on various sites? Is it because only a small subsets of amino acids are tolerated at some sites? Why not take several enzyme families with very thick alignments and analyze them position-by-position to shed light on this question? Without such deep analysis, the paper, while interesting in its general idea, reads superficial and underdeveloped, despite some experimental results.

We apologize for the possible confusion. The mechanistic explanation of our results is the following. We demonstrate in the manuscript that the long-term evolution of proteins that continuously maintain the same molecular functions generally requires a strict conservation of protein structures (<3Å C-α root mean square deviation). The structural conservation is likely necessary for proper structural positioning of active site residues required for efficient function and catalysis. This level of structural conservation, in agreement with the Chothia and Lesk analysis (Chothia and Lesk, 1986), requires a substantial conservation of protein sequences, which leads to the observed sequence divergence limit. Furthermore, our analysis shows that only about 35% of the sequence conservation between distant orthologs can be attributed to universally conserved protein sites, i.e. sites occupied by identical amino acids in almost all lineages. In contrast, the observed divergence limit primarily originates from a small number of acceptable amino acids at each protein site under the continuous constraint to strictly conserve protein structure and function. Following the reviewers’ suggestion, we now present an analysis of sequence conservation across protein sites using multiple sequence alignments for 30 diverse enzyme families (Figure 2B, Figure 2—figure supplement 1).

2) The authors premise their study on the assertion made by Rost in a 1997 review that proteins of the same fold can diverge up to random sequence ID of 3-5%, i.e., that structural similarity does not impose constraint on sequence ID. As a consequence the authors attribute in their analysis all sequence identity constraints to functional similarity. However, this premise is highly problematic. While anecdotally one can find two proteins of similar fold with very low sequence ID, typically proteins with similar structure (not necessarily orthologues) have sequence ID at or above 25% – this is the essence of the famous cusp in sequence-structure relationship first reported by Lesk and Chothia in the eighties. In fact a careful analysis of structural/stability constraints on sequence divergence has been published in recent paper (Biophys J v.112, p.1350-65, 2017) where it was shown that in divergent evolution scenario protein structure is maintained up to 25-30% sequence ID and quickly deteriorates beyond that. Incidentally it is very close to 25% ID which the authors claim to stem from functional constraints. Furthermore, the above mentioned BJ paper presents a detailed analytical estimate of sequence divergence dynamics akin to exponential fits used in this PAPER. The authors should make themselves familiar with this theory and consider using it to fit their divergence curves rather than ad hoc exponential fits.

We believe that there may be a potential misunderstanding here in terms of the definition of protein folds and structures, and in terms of the nature of evolutionary constraints. A protein fold is defined as a particular set of protein secondary structures, their arrangement and topological connections (Murzin et al., 1995). Both the Rost 1997 paper and the Biophysical J. paper (cited by the reviewer) actually show that a low sequence identity (<25%) between proteins with the same fold is very common. Thus, the sequence constraints necessary to simply maintain a protein fold cannot be responsible for the divergence limit we observe. We also would like to emphasize that selection does not specifically act to maintain structural similarity; instead, structural and sequence conservation is a consequence of the selection for a molecular function necessary to maintain species fitness. In contrast to protein folds, our analysis demonstrates that the continuous conservation of specific molecular functions usually requires high structural similarity (<3Å C-α root mean square deviation), and that the requirement for the structural conservation likely constrains the divergence of sequence identity (>25%).

In terms of the analytical model presented in the Biophysical J. 2017 paper, we note that the exponential formulas used in our fits are not ad hoc. They directly correspond to a divergence model with equal probability of mutations across protein sites, as in the Biophysical J. model. The BJ model is numerically identical to our second model. Specifically, the fraction of identical sites (y) after a certain divergence time (t) in the BJ model is given by:

y=1-((a-1)/a)(1-(1-a/l(a-1)μt)

according to Equation 7 in the BJ paper, where the parameter l is the length of the protein, a is the number of acceptable amino acid types per protein site, and μ is the per-protein substitution rate. After re-arranging, and given that the substitution rate per site μ0 =μ × l-1:

y=1a+(1-1/a)(1+((-a/(a-1))/l))0t

Given the well-known limit: exp(x)=limn→∞⁡(1+x/n)n, and because protein length l is substantially larger than a/(a-1):

y=1/a+(1-1/a)exp(-R0t); with R0=(a/(a-1))μ0

This equation is numerically identical to our model 2. Furthermore, this demonstrates that the parameter Y0 in model 2 can be also interpreted as the inverse of the effective number of amino acid types acceptable per protein site when molecular function is strictly conserved. Our results thus suggest that, on average, only 2 to 4 amino acids types are usually accepted per protein site as long as protein molecular function is conserved.

3) The high throughput experimental data is potentially interesting but its description is so cryptic and incomplete that it is very hard to assess what actually has been done there. In particular there is no assessment of noise in deep sequencing to assess fitness, of the effects of synonymous substitutions, etc. How is the data binned and why binning? The use of LB as a medium for the folA experiment is unfortunate because in LB DHFR is much less essential than in minimal media and therefore the results could be skewed by very permissive conditions. In addition the standing genetic variation is not taken into account – fitness effects could be a result of hitchhiking (i.e. originate from the variation of the background into which the mutations are introduced). These concerns aside, it is not entirely clear what does comparison with experiment tell us beyond the obvious that sites where deleterious fitness effects are greatest evolve more slowly. Was the sole purpose of the experiment to bin positions by their fitness effects? If so, what new insight is gained from slower divergence of sites with more profound fitness effects compared to faster divergence of sites that are more tolerant to substitutions in experiments. Isn't this effect expected? Maybe quantification of the effect may have some novelty, but such quantification has not been done. The obvious thing to do is to compare their alleged fitness effects with dN/dS assessment for the rates but that has not been done.

We agree with the reviewers that the correlation between the fitness effects of mutations and evolutionary conservation is indeed expected. But this particular result was neither a goal nor a major conclusion of our analysis. Our primary goals were to characterize the long-term temporal divergence patterns at sites with different fitness effects, to explore the limits of the sequence divergence at different sites, and to investigate the contribution of sites with various fitness effects to the divergence limit. For example, our analysis demonstrated (Figure 3), for the first time to our knowledge, that the divergence at sites with relatively small fitness effects usually saturates after 1-2 billion years of evolution, while the divergence at sites with higher fitness effects still continues, although at a substantially smaller rate.

In our analysis, we binned the sites (based on their average fitness effects) in order to explore the sequence divergence for sites within each fitness category. We used logarithmic bins in our analysis due to the skewed distribution of fitness effects in the experimental data (Figure 3—figure supplement 2).

Following the reviewers’ suggestion, we now use several approaches to demonstrate the reliability of our experimental measurements: 1.) we show that mutants targeting start and stop codons have, on average, significantly larger growth effects compared to amino acid replacements; 2.) we show that synonymous substitutions yield very small fitness effects; and 3.) by calculating average growth defects per site (based on strains with non-overlapping sets of introduced codons), we demonstrate a high reproducibility (r = 0.95) of our experimental findings. Overall, these results demonstrate that neither sequencing noise nor background mutations have a significant effect on the experimental measurements.

4) The word "protein" should be replaced with "enzyme" in the title. "We focused on enzymes because their molecular function is usually well defined," doesn't sound quite true. Reading the manuscript, it seems that the goal of this work is to investigate divergence of enzymes. The divergence of non-enzymes will most likely follow a different pattern and will generally be much lower than 25% identity. If the authors want to generalize to non-enzymes, they need to perform appropriate analyses. Currently, only one experiment is dedicated to a non-enzyme. It is not fair to the readers to generalize beyond the scope of analysis that was carried out.

We indeed believe that for enzymes, specific molecular functions are usually well-defined. Following the reviewers’ suggestion, we now also present the results of similar sequence analyses for orthologous proteins that are not enzymes. In order to explore the long-term divergence limit, protein orthologs need to be present across very long evolutionary distances. Therefore, the analyzed proteins mostly represent ribosomal proteins, chaperones and electron transport flavoproteins. Notably, the results obtained for these protein families are consistent with the divergence limit observed for enzymes (Figure 1—figure supplement 5). We also now describe an interesting and specific case of several mitochondrial ribosomal proteins which display much faster divergence rates compared to their prokaryote orthologs; a likely consequence of their co-evolution with the fast-evolving mitochondrial-encoded rRNA.

5) The sentence in the Abstract "divergence rates of orthologous enzymes decrease substantially after ~1-2 billion years of independent evolution" is confusing. "Rate" means number of amino acid substitutions per site per unit of time. Moreover, the authors write in the Introduction "a divergence limit does not imply that the rate of amino acid substitutions slows down in evolution," thus contradicting the Abstract. The Abstract should be edited and the misleading sentence corrected. From the context of the paper is seems like by "divergence rate" the author mean the rate of decrease in pairwise sequence identity with time. This definition poses additional question. This "divergence rate" will always decrease due to multiple and back substitutions and it is trivial. The authors need to make it clear that what they mean goes beyond a trivial curvilinear relationship between time and sequence identity.

The reviewer is correct, and we apologize for being unclear. We specifically define the sequence divergence rate (now in the third paragraph of the Introduction) as the decrease in the pairwise sequence identity per unit time. Notably, this divergence rate is different from the substitution rate, which corresponds to the number of amino acid substitutions per site per unit time. We have also edited the Abstract and the main text to further clarify this distinction. As we discuss in the paper (fifth paragraph of the Results) the observed divergence limit is not a trivial consequence of the curvilinear relationship between divergence time and sequence identity. For example, according to the data fits using model 1 and model 3, the decrease in the sequence divergence rate after ~1-2 billion years of evolution is more than an order of magnitude higher than expected simply due to multiple substitutions at the same sites. We show in the paper that the observed sequence divergence limit is primarily due to a small number of acceptable amino acids at each protein site under the continuous and strict conservation of protein function.

6) Further on this, evolutionary rate is much more than simply a count of changes in the identity of the amino acids. To measure it properly one needs to rely on the phylogeny. Indeed, since the study is based on several clusters of orthologous proteins, why not do just that? For example, why not use the dN/dS formalism?

We would like to clarify that the primary goal of our manuscript was not to determine the short-term evolutionary rates across different sites, which can be indeed quantified using the dN/dS formalism. Our main goals were to characterize the patterns of long-term sequence divergence across protein sites and the contribution of sites to the divergence limit for proteins with the same molecular function. Unfortunately, short-term substitution rates (based on the dN/dS formalism) are not informative across long evolutionary distances (billions of years) due to saturation of substitutions at synonymous sites. We rely in our analyses on phylogenetic information by using the consensus species' divergence times (Supplementary file 1).

7) "billion years of evolution, we observed a significant decrease in mutual divergence rates," see item 5 above. Also, shouldn't "significant" be quantified is some way?

Following the reviewer’s comment, we have re-worded this sentence to clarify that we are specifically referring to the decline of sequence identity per unit time; we also now report that, based on model 3 and Equation 5, the sequence divergence rate for proteins with the same molecular function decreases more than ten times after ~1.5 billion years of evolution.

8) "Species' divergence times were used to estimate the times of divergence between corresponding orthologous proteins": How reliable is this time estimate?

We thank the reviewer for this question. We now present in Supplementary file 1 data describing the variance in the estimated divergence times based on multiple independent publications. Overall, most estimates of species divergence differ by tens to several hundred million years. To assess the sensitivity of our results to this variation, we repeated the analysis using ~1000 independent runs by randomly sampling either the maximum or minimum reported divergence times between each pair of lineages. The fits of the divergence model 2 based these simulations demonstrate that for ~90% of proteins the independent assignments of the divergence times support a long-term sequence identity limit higher than 25% (Figure 1—figure supplement 7).

9) Discussion "The presented results demonstrate that, in contrast to proteins with the same fold [Rost, 1997], the requirement to maintain the same molecular function significantly constrains the long-term divergence of protein orthologs": This difference is surprising. Maybe it has to do with how changes were estimated in the two studies? To state it with confidence direct comparison is needed, where the exact same methodology is used.

We agree with the reviewer that this result is both interesting and surprising. This is one of the key results of our manuscript. This result is especially interesting given that proteins with the same fold often have very low levels of sequence similarity. Our analysis demonstrates that, in contrast to protein that simply share the same fold, orthologs that continuously conserve specific molecular functions usually display high degree of structural and sequence similarity even after billions of years of divergent evolution. We also show (Figure 5) that even a small relaxation of functional conservation (for example, from sharing all 4 digits of the EC number to only sharing the first 3 digits of the EC number) results in a significant increase of long-term sequence divergence.

10) Finally, the manuscript should be edited for clarity by a professional scientific editor.

Following the reviewer’s suggestion, we have extensively revised the text to improve clarity, and to highlight the novelty and importance of our findings.

Separate reviews (please respond to each point):

Reviewer #1:

Summary of paper:

The manuscript analyzes clusters of orthologue enzymes of the same function, defined as sharing the same 4-digits EC number, to examine the limitation of function preservation on evolutionary divergence. The results seem to indicate that function confines divergence, which is anticipated- perhaps also shown before. The potential novelty here is the time dependence of the process. The data shows (to the best of my understanding) that proteins diverse from each other rapidly during the first ~1.5 billion years, but saturate after that. The data further suggests (again, assuming that I understood correctly) that the rapid evolutionary phase mostly reflects random drift, i.e., amino acid sites that can freely mutate, and that the saturation value reflects their relative fraction of the entire enzyme sequence.

Opinion: The manuscript could be important but it should be revised significantly to clarify what exactly was computed and why, what was measured and why, and how the results are related to previous studies. I also believe the manuscript should be edited for clarity by a professional scientific editor.

Major issues:

1) Results "In the first model, all protein sites have equal and independent substitution rates; see Equation 1": Equation 1 obviously represents an unrealistic evolutionary model because it is well known that the evolutionary rates are not equal at all amino acid sites. For example, catalytic residues hardly change.

2) Where does Equation 2 come from? It also doesn't look particularly realistic.

3) Equation 3 is the commonly used evolutionary model I think. So why not start with it?

All scientific models are ultimately approximations to reality. They are essentially useful tools for answering specific questions and testing hypotheses. For example, despite being clearly an approximation to reality, model 2 fits the long-term sequence divergence data pretty well (Figure 1). Notably, increasing the complexity of the models (e.g. model 2 compared to model 3) does not significantly improve the fits to the data (see R2 values in Supplementary file 2A). Model 2 describes sequence divergence with independent and equal probability of amino acid substitutions across protein sites and a divergence limit at long evolutionary distances.

Our main reason for first considering models 1 and 2 is to test whether the observed decline in protein sequence identity (i.e. the sequence divergence rate) is simply due to multiple substitutions at the same sites (a trivial explanation for the curvilinear relationship between time and sequence identity), or whether it is necessary to assume additional constraints, such as the long-term divergence limit. Testing the goodness of fits between these two simple models allows us to answer this specific question. We now describe in the text that model 2 represents a simple modification to model 1 in which the exponential decay of sequence identity is bounded by a minimum identity at long evolutionary distances.

4) "Species' divergence times were used to estimate the times of divergence

between corresponding orthologous proteins": How reliable is this time estimate?

We thank the reviewer for this question. We now present in Supplementary file 1 data describing the variance in the estimated divergence times based on multiple independent publications. Overall, most estimates of species divergence differ by tens to several hundred million years. To assess the sensitivity of our results to this variation, we repeated the analysis using ~1000 independent runs by randomly sampling either the maximum or minimum reported divergence times between each pair of lineages. The fits of the divergence model 2 based these simulations demonstrate that for ~90% of proteins the independent assignments of the divergence times support a long-term sequence identity limit higher than 25% (Figure 1—figure supplement 7).

5) Figure 2: How is "fitness cost" defined?

Fitness cost was defined as the percent decrease in bacterial growth rate relative to the original strain, i.e. the “wild type”, in which mutations were introduced. We have now clarified this definition in the manuscript.

6) "phylogenetically independent pairs of species": What does it mean and how is "independent" defined? Is the pair A-B in Figure 3A phylogenetically independent? And the pair D-H?

Phylogenetically independent pairs of species are defined as pairs of species that do not share any edges in the phylogenetic tree with each other since each pair diverged from a common ancestor. In Figure 2A of the updated manuscript, the pair A-B is independent from the pair D-H because they do not share edges in the phylogenetic tree since divergence from a common ancestor. By contrast, the pair D-H is not independent from the pair E-F. We now provide these examples in the main text to clarify these points.

7) "Notably, the probability that a protein site is identical, and thus contributes to the divergence limit, first increases linearly with increasing fitness effects at the site, and then begins to saturate for sites with high (>30%) fitness effects (Figure 3B)." This is perhaps true for FolA but I do not see saturation for InfA.

We agree with the reviewer that the saturation is more prominent for FolA. However, we also observe it in InfA; for example, for InfA there is a larger (>2-fold) change in the probability of site identity in the 10-40% interval of the relative growth rate decrease (x-axis), compared to the 40-70% interval.

8) The distributions of the probability of identical sites in FolA and InfA (Figures 4A and B) are very different, but they are discussed together without reference to the difference.

We apologize for being unclear. Likely due to family-specific constraints, the distributions of the probability of identical sites are indeed quite different between FolA and InfA. But our main point in this analysis is that in both cases the fraction of universally conserved sites, i.e. sites conserved in almost all lineages, is relatively small. We now present in Figure 2B and Figure 2—figure supplement 1 a more comprehensive analysis based on 30 different enzyme families. These results demonstrate that sites conserved in the majority of the linages are responsible for only about a third of the overall sequence identity between distant orthologs.

9) Discussion "The presented results demonstrate that, in contrast to proteins with the same fold [1997], the requirement to maintain the same molecular function significantly constrains the long-term divergence of protein orthologs": I am surprised about this difference. Maybe it has to do with how changes were estimated in the two studies? To state it with confidence I would have liked to see direct comparison where the exact same methodology is used.

We agree with the reviewer that this result is interesting and surprising. This is one of the key results of our manuscript. Previous studies demonstrated (for example, Rost, 1997), that proteins sharing the same fold often have very small degrees of sequence and structural similarity. Our analysis shows that, in contrast to proteins that simply share the same fold, orthologs that continuously conserve specific molecular functions usually display high degree of structural and sequence similarity even after billions of years of divergent evolution. We also show (Figure 5) that even a small relaxation of functional conservation (for example, from sharing all 4 digits of the EC number to only sharing the first 3 digits of the EC number) results in a significant increase of long-term sequence divergence.

10) Evolutionary rate is much more than simply a count of changes in the identity of the amino acids. To measure it properly one needs to rely on the phylogeny. Indeed, since the study is based on several clusters of orthologous proteins, why not do just that?

We agree with the reviewer that evolutionary rate is more than a count of amino acid changes, which is why we use multiple divergence models to fit the data. We note that we indeed use phylogenetic information by considering published divergence times between different lineages. However, we do not fit an evolutionary model to the clusters of orthologous proteins because the short-term evolutionary rates derived in this way would not be informative about divergence patterns over long geological timescales. Our study is focused on understanding the patterns of long-term divergence of sequence and structural identity for proteins conserving molecular function, rather than estimating short-term evolutionary rates across protein sites.

Reviewer #2:

This study poses the question about the extent to which function constrains sequence divergence. The answer given is that for orthologous enzymes 25% sequence identity is the approximate limit. The major critique, is that the authors did not make an attempt to provide a clear mechanistic explanation for such phenomenon. Is it simply because a significant fraction of positions is invariant? Is it because there are various constraints on various sites? Is it because only a small subsets of amino acids are tolerated at some sites? Why not take several enzyme families with very thick alignments and analyze them position-by-position to shed light on this question? Without such deep analysis, the paper, while interesting in its general idea, reads superficial and underdeveloped, despite some experimental results.

We thank the reviewer for suggestions and apologize for the potential confusion. The mechanistic explanation of our results is the following. We demonstrate in the manuscript that the long-term evolution of proteins that continuously maintain the same molecular functions generally requires a strict conservation of protein structures (<3Å C-α root mean square deviation). The structural conservation is likely necessary for proper structural positioning of active site residues required for efficient function and catalysis. This level of structural conservation, in agreement with the Chothia and Lesk analysis (Chothia and Lesk, 1986), requires a substantial conservation of protein sequences, which leads to the observed sequence divergence limit. Furthermore, our analysis shows that only about 35% of the sequence conservation between distant orthologs can be attributed to universally conserved protein sites, i.e. sites occupied by identical amino acids in almost all lineages. In contrast, the observed divergence limit primarily originates from a small number of acceptable amino acids at each protein site under the continuous constraint to strictly conserve protein structure and function. Following the reviewers’ suggestion, we now present an analysis of sequence conservation across protein sites using multiple sequence alignments for 30 diverse enzyme families (Figure 2B, Figure 2—figure supplement 1).

1) I suggest replacing the word "protein" with "enzyme" in the title. "We focused on enzymes because their molecular function is usually well defined," doesn't sound quite true. Reading the manuscript, it seems that the goal of this work is to investigate divergence of enzymes. I am positive that the divergence of non-enzymes will follow a different pattern and will generally be much lower than 25% identity. If the authors want to generalize to non-enzymes, they need to perform appropriate analyses. Currently, only one experiment is dedicated to a non-enzyme. It is not fair to the readers to generalize beyond the scope of analysis that was carried out.

We indeed believe that for enzymes, specific molecular functions are usually well-defined. Following the reviewers’ suggestion, we now also present the results of similar sequence analyses for orthologous proteins that are not enzymes. In order to explore the long-term divergence limit, protein orthologs need to be present across very long evolutionary distances. Therefore, the analyzed proteins mostly represent ribosomal proteins, chaperones and electron transport flavoproteins. Notably, the results obtained for these protein families are consistent with the divergence limit observed for enzymes (Figure 1—figure supplement 5). We also now describe an interesting and specific case of several mitochondrial ribosomal proteins which display much faster divergence rates compared to their prokaryote orthologs; a likely consequence of their co-evolution with the fast-evolving mitochondrial-encoded rRNA.

2) I think that the sentence in the Abstract "divergence rates of orthologous enzymes decrease substantially after ~1-2 billion years of independent evolution" is confusing. "Rate" usually means number of amino acid substitutions per site per unit of time. Moreover, the authors write in the Introduction "a divergence limit does not imply that the rate of amino acid substitutions slows down in evolution," thus contradicting the Abstract. I suggest to edit the Abstract and correct the misleading sentence. From the context of the paper is seems like by "divergence rate" the author mean the rate of decrease in pairwise sequence identity with time. This definition poses additional question. This "divergence rate" will always decrease due to multiple and back substitutions and it is trivial. The authors need to make it clear that what they mean goes beyond a +trivial curvilinear relationship between time and sequence identity.

The reviewer is correct, and we apologize for being unclear. We specifically define the sequence divergence rate (now in the third paragraph of the Introduction) as the decrease in the pairwise sequence identity per unit time. Notably, this divergence rate is different from the substitution rate, which corresponds to the number of amino acid substitutions per site per unit time. We have also edited the Abstract and the main text to further clarify this distinction. As we discuss in the paper (fifth paragraph of the Results) the observed divergence limit is not a trivial consequence of the curvilinear relationship between divergence time and sequence identity. For example, according to the data fits using model 1 and model 3, the decrease in the sequence divergence rate after ~1-2 billion years of evolution is more than an order of magnitude higher than expected simply due to multiple substitutions at the same sites. We show in the paper that the observed sequence divergence limit is primarily due to a small number of acceptable amino acids at each protein site under the continuous and strict conservation of protein function.

3) "billion years of evolution, we observed a significant decrease in mutual divergence rates," see item 2 above. Also, shouldn't "significant" be quantified is some way?

Following the reviewer’s comment, we have re-worded this sentence to clarify that we are specifically referring to the decline of sequence identity per unit time; we also now report that, based on model 3 and Equation 5, the sequence divergence rate for proteins with the same molecular function decreases about ten times after ~1.5 billion years of evolution.

4) It was not made clear why these experiments were performed and how they integrate with the rest of the study. Addition of these experiments seems preliminary and no confident conclusions follow. Was the sole purpose of the experiment to bin positions by their fitness effects? If so, what new insight is gained from slower divergence of sites with more profound fitness effects compared to faster divergence of sites that are more tolerant to substitutions in experiments. Isn't this effect expected? Maybe quantification of the effect may have some novelty, but such quantification has not been done.

We agree with the reviewer that the correlation between the fitness effects of mutations and evolutionary conservation is expected. But this particular result was neither a goal nor a major conclusion of our analysis. Our primary goals were to characterize the long-term temporal divergence patterns at sites with different fitness effects, to explore the limits of the sequence divergence across sites, and to investigate the contribution of sites with various fitness effects to the divergence limit. For example, our analysis demonstrated (Figure 3), for the first time to our knowledge, that the divergence at sites with relatively small fitness effects usually saturates after 1-2 billion years of evolution, while the divergence at sites with higher fitness effects still continues, although at a substantially smaller rate.

Reviewer #3:

[…] While the paper contains some interesting bioinformatics observations and analysis I have serious concerns about its premise and interpretation of the results. Apparently some issues stem from author's apparent gaps in knowledge of modern literature on biophysical determinants of protein evolution.

1) The premise that structural similarity does not constraint sequence identity is highly problematic. While anecdotally one can find two proteins of similar fold with very small sequence ID typically proteins with similar structure (not necessarily orthologues) have sequence ID at or above 25% – this is the essence of the famous cusp in sequence-structure relationship first reported by Lesk and Chothia in the eighties. In fact a careful analysis of structural/stability constraints on sequence divergence has been published in recent paper (Biophys J v.112, p.1350-65 (2017) where it was shown that in divergent evolution scenario protein structure is maintained up to 25-30% sequence ID and quickly deteriorates beyond that. Incidentally it is very close to 25% ID which the authors claim to stem from functional constraints Furthermore, the above mentioned BJ paper presents a detailed analytical estimate of sequence divergence dynamics akin to exponential fits used in this PAPER. The authors should make themselves familiar with this theory and use it to for their divergence curves rather ad hoc exponential fits.

We believe that there may be a potential misunderstanding here in terms of the definition of protein folds and structures, and in terms of the nature of evolutionary constraints. A protein fold is defined as a particular set of protein secondary structures, their arrangement and topological connections (Murzin et al., 1995). Both the Rost 1997 paper and the Biophysical J. paper (cited by the reviewer) actually show that a low sequence identity (<25%) between proteins with the same fold is very common. Thus, the sequences constraints necessary to simply maintain a protein fold cannot be responsible for the divergence limit we observe. We also would like to emphasize that selection does not specifically act to maintain structural similarity; instead, structural and sequence conservation is a consequence of the selection for a molecular function necessary to maintain species fitness. In contrast to protein folds, our analysis demonstrates that the continuous conservation of specific molecular functions usually requires high structural similarity (<3Å C-α root mean square deviation), and that the requirement for the structural conservation likely constrains the divergence of sequence identity (>25%).

In terms of the analytical model presented in the Biophysical J. 2017 paper, we note that the exponential formulas used in our fits are not ad hoc. They directly correspond to a divergence model with equal probability of mutations across protein sites, as in the Biophysical J. model. The BJ model is numerically identical to our second model. Specifically, the fraction of identical sites (y) after a certain divergence time (t) in the BJ model is given by:

y=1-((a-1)/a)(1-(1-a/l(a-1))μt)

according to Equation 7 in the BJ paper, where the parameter l is the length of the protein, a is the number of acceptable amino acid types per protein site, and μ is the per-protein substitution rate. After re-arranging, and given that the substitution rate per site μ0 = μ × l-1:

y=1/a+(1-1/a)(1+(-a/(a-1))/l)0t

Given the well-known limit: exp(x)=limn→∞⁡(1+x/n)n, and because protein length l is substantially larger than a/(a-1):

y=1/a+(1-1/a)exp(-R0t); with R0=(a/(a-1))μ0

This equation is numerically identical to our model 2. Furthermore, this demonstrates that the parameter Y0 in model 2 can be also interpreted as the inverse of the effective number of amino acid types acceptable per protein site when molecular function is strictly conserved. Our results thus suggest that, on average, only 2 to 4 amino acids types are usually accepted per protein site as long as protein molecular function is conserved.

2) The high throughput experimental data is potentially interesting but its description is so cryptic and incomplete that it is very hard to assess what actually has been done there. In particular there is no assessment of noise in deep sequencing to assess fitness., of the effects of synonymous substitutions etc. How is the data binned and why binning? The use of LB as a medium for the folA experiment is unfortunate because in LB DHFR is much less essential than in minimal media and therefore the results could be skewed by very permissive conditions. In addition the standing genetic variation is not taken into account – fitness effects could be a result of hitchhiking (i.e. originate from the variation of the background into which the mutations are introduced).These concerns aside, it is not entirely clear what does comparison with experiment tell us beyond the obvious that sites where deleterious fitness effects are greatest evolve more slowly. The obvious thing to do is to compare their alleged fitness effects with dN/dS assessment of the rates but that has not been done.

We agree with the reviewers that the correlation between the fitness effects of mutations and evolutionary conservation is indeed expected. But this particular result was neither a goal nor a major conclusion of our analysis. Our primary goals were to characterize the long-term temporal divergence patterns at sites with different fitness effects, to explore the limits of the sequence divergence at different sites, and to investigate the contribution of sites with various fitness effects to the divergence limit. For example, our analysis demonstrated (Figure 3), for the first time to our knowledge, that the divergence at sites with relatively small fitness effects usually saturates after 1-2 billion years of evolution, while the divergence at sites with higher fitness effects still continues, although at a substantially smaller rate.

We would also like to clarify that the primary goal of our manuscript was not to determine the short-term evolutionary rates across different sites, which can be indeed quantified using the dN/dS formalism. Our main goals were to characterize the patterns of long-term sequence divergence for proteins with the same molecular function. Short-term substitutions rates (calculated using the dN/dS formalism) cannot reveal the patterns of sequence divergence across billions of years of evolution due to saturation of substitutions at synonymous sites.

In our analysis, we binned the sites (based on their average fitness effects) in order to explore the sequence divergence for sites within each fitness category. We used logarithmic bins in our analysis due to the skewed distribution of fitness effects in the experimental data (Figure 3—figure supplement 2).

Following the reviewers’ suggestion, we now use several approaches to demonstrate the reliability of our experimental measurements: 1.) we show that mutants targeting start and stop codons have, on average, significantly larger growth effects compared to amino acid replacements; 2.) we show that synonymous substitutions yield very small fitness effects; and 3.) by calculating average growth defects per site (based on strains with non-overlapping sets of introduced codons), we demonstrate a high reproducibility (r = 0.95) of our experimental findings. Overall, these results demonstrate that neither sequencing noise nor background mutations have a significant effect on the experimental measurements.

3) The interpretation of the experimental fitness effects in terms of function is also questionable. The authors are apparently unaware of series of experimental works from Shakhnovich lab where determinants of fitness effects of mutations are addressed. In particular it has been shown using both point mutations and orthologous chromosomal replacements for folA gene (PLOS Genetics 2015 DOI:10.1371/journal.pgen.1005612) and adk gene (Nature Ecology Evolution, 2017 http://dx.doi.org/10.1038/s41559-017-0149) that fitness is determined by product of folded protein abundance A and activity kcat/KM. Mutations may affect stability and through that parameter A (by changing the balance between protein production and degradation, see Mol Cell v.49, pp133-44 (2013). Therefore interpretation of the experimental trends entirely in functional terms is not warranted.

We are aware of the experimental work from the Shakhnovich lab, and agree with the reviewer that the effects of mutations on fitness could arise either due to effects on protein activity or protein stability, which affects the product of protein abundance and activity. We note however that direct and comprehensive biochemical experiments demonstrate that the deleterious effects of protein mutations primarily arise from changes in specific protein activity rather than decreases in stability or protein cellular abundance; see, for example, Firnberg et al., 2014. We now clarify this point in the Discussion.

A minor comment: The concept and metaphor of expanding protein universe has been introduced 10 years before Kondrashov's work in the paper "Expanding protein universe and its origin from biological big bang" PNAS 2002, v.99 pp. 14132-6.

We thank the reviewer for bringing this manuscript to our attention, and we now cite this paper in the main text.

[Editors' note: further revisions were suggested before publication, as described below.]

As you can see from the reports below, the reviewers appreciated the revisions. However, there are still major outstanding issues. While some of these can be resolved by changes in the presentation, others are fundamental. We would strongly encourage you to address all of these prior to publication.

We again thank the editors and reviewers for their remarks and suggestions. We have updated the manuscript to fully address the comments and remaining concerns. Following the recommendation of the second reviewer, we highlighted the low average number of acceptable amino acids per site associated with the conservation of protein molecular function. Following the suggestions of the third reviewer, we further clarified the effects of functional constraints on the divergence limit both for protein sequences and structures. We also demonstrated, as the third reviewer requested, that even at the same level of structural divergence, orthologs sharing more specific molecular functions display substantially higher levels of long-term sequence identity.

Reviewer #2:

I find that the authors did a thorough revision of the manuscript. At least now I think I understood the main conclusion of the paper. In enzymes and other proteins with very strong functional conservation, the number of different amino acids acceptable at a position is about 4. It is not because many sites are invariant and some are variable (not a dirty trick of average temperature in a hospital), but because most sites (except the invariant ones) are constrained to use a library of 3 to 5 different amino acids, not more than that. If this is not the bottom line, then the authors need to do better job at crystallizing their main claim and result.

We thank the reviewer for the comments. Indeed, the observation that only a few (less than 4) amino acids are accepted, on average, per site for proteins with conserved molecular function is one of our main results, which we now highlight in the Abstract. However, as we explained in our previous reply, this is not the only main result or conclusion of our work.

We apologize for repeating our statements from a previous reply. In the paper we asked several important and interrelated questions: 1) How far can two sequences diverge while continuously maintaining the same specific molecular function? 2) What are the temporal patterns of this divergence across billions of years of evolution? and 3) How do different protein sites contribute to the long-term divergence between orthologs with the same molecular function? Our paper provides corresponding main conclusions: 1) Orthologs that share the same molecular function usually do not diverge beyond ~3.5 Å RMSD and ~25% sequence identity (this translates to 3-4 amino acids, on average, per site). 2) The decrease of sequence identity between diverging orthologs becomes increasingly slow past 1-2 billion years of divergent evolution, such that ancient orthologs have retained approximately the same level of similarity for the past billion years. And 3) the observed long-term sequence similarity is not primarily due to universally conserved sites across all orthologs. Instead, as the reviewer correctly comments, most sites are on average constrained to a library of less than 4 amino acids. We also show that long-term sequence constraints can be partially explained by a site’s fitness effects of mutations and its distance to catalytic residues.

If it is, I think it is a meaningful finding that could be explained better to the readers. I guess the second claim is that 3-5 amino acid limit is universal to all enzymes and (conserved!) non-enzymes. I do doubt (as in the original review) the validity of such a strong claim, which could be a result of the authors' bias in selecting families for their analysis. At least for non-enzymes, they selected most conserved proteins (like ribosomal proteins), so of course such selection is biased to get proteins that saturate easily in evolution. The authors try to justify this biased selection suggesting that it is difficult to find orthologs. But that statement by itself totally discredits this study. Why? Because if you cannot find your orthologs, wouldn't it mean that they already diverged beyond (lower than) your claimed 25% identity as the universal limit, and the author's conclusions do not apply to such families? Maybe for the enzymes too, the 25% limit is simply a reflection of the search methods the authors used to find orthologs, that fail to find more distant ones.

We apologize for the possible confusion. But we discussed this particular concern in the manuscript, and are afraid that the related discussion escaped the attention of the reviewer. First, for almost all analyzed enzymatic protein families, sequence divergence saturates at a relatively high sequence identity (25%-40%); at these sequence identities, protein orthologs can be easily identified by computational sequence comparison methods. Second, when we relax the stringency of the molecular function conservation (e.g. from sharing all digits of an EC number (EC4) to sharing only the first three digits (EC3)), we immediately observe orthologous protein pairs with significantly lower levels of sequence identity (Figure 6). For non-enzymatic protein families, we need to investigate families that conserved specific molecular functions across billions of years of evolution. That is why we are limited in the number of protein families we can consider, not because we cannot find orthologs, but because there are not many such non-enzymatic protein families (beyond known examples such as ribosomes, flavoproteins or chaperones) that continuously conserved specific molecular functions across billions of years of evolution.

To improve the presentation and make this paper quite interesting (well, reviewers are not supposed to direct the study, but the authors seem passionate about their work and also seem rather inexperienced in both logical thinking and putting a paper together, so maybe this recommendation could be helpful), I would suggest to base it on two plots.

1) Between-protein variability.

The first plot is the histogram of the average estimated number of amino acids allowable per site (without invariant sites, or with) for enzyme families and other protein families. Well, the authors would have to try harder to find orthologs for proteins that manage to diverge below 25%, even my rotation students can do that. This histogram would be expected to have a maximum around 3 to 5, and for some variable proteins it could be around 10 or 15? Right? Then protein families with very low and very high numbers could be discussed with an attempt to give explanations about their uniqueness.

2) Within-protein variability.

The second plot is the histogram of estimated number of amino acids allowable per site within a protein family. The authors could either normalize to the average per protein, or select a bin with the maximum count from previous plot (let's say 4). These histograms can be averaged for all families with mean and SD showing for each bin. I would assume there will be a large count for invariant sites for enzymes (at 1), counts for 2 amino acids used, 3, 4, etc. Will this histogram have a single mode? Maybe around 4? Several modes? Will there be sites using more than 10 amino acids? How are these distributed in spatial structure and relative to the active site? Discussion of these details could be quite insightful and interesting.

If these authors do not wish to make these plots, since this review will be published, maybe someone else will, and we will learn something interesting as a result.

We thank the reviewer for the suggestions, but would like to first make an editorial comment. We do not think that condescending personal comments, i.e. calling into question the authors’ abilities for logical thinking (see more about it below), are either appropriate or helpful. We do not take these comments personally, but anonymously throwing personal comments is not what science or publishing should be about. We hope that the reviewer agrees. We also believe that publication of peer-reviews by eLife and other prominent journals will raise the overall level of scientific discussion, and make reviewers feel responsible for their comments.

Back to science. The comments above are especially surprising because the suggested analyses were mostly already present in the manuscript. We again think they may have simply escaped the reviewer’s attention. We have now made several further changes to highlight these results:

1) The between-protein variability in long-term sequence identity (previously Figure 1—figure supplement 2A) is now presented as a main figure in the manuscript (Figure 2A). We have added an additional x-axis (at the top of the figure) to show the corresponding average number of accepted amino acids per site; we also now present the results of an equivalent analysis for non-enzymatic protein families (Figure 2B).

2) In regards to the comparison between enzymes that conserve their function and proteins that “manage to diverge below 25% identity”, we have now modified Figure 6A to highlight the higher average number of acceptable amino acids for ancient orthologs sharing only three digits of the EC annotation. We also show the sequence identity of proteins with different molecular functions but the same structural fold. Additionally, we compared those differences to the corresponding long-term structural similarities between orthologs (Figure 6B).

3) In terms of within protein variability, we now comment on the different shapes of the distributions shown in Figure 3—figure supplement 1. We decided not to average these distributions across enzymes, as this would hide an interesting underlying variability in the number of acceptable amino acids across protein families.

As we describe in the manuscript (see discussion of Figure 6), the aforementioned results demonstrate that the average number of acceptable amino acids per site increases from the case when a specific molecular function is conserved (~3 amino acids for EC4), to conservation of less specific functions (~4 amino acids for EC3), to simply preserving the same protein fold (~10 amino acids). Naturally, a higher long-term sequence identity limit usually corresponds to a smaller average number of acceptable amino acids per site.

Currently, there are so many plots in this paper and many of them are not particularly helpful to get the point across quickly. Also I still find the usage of non-standard terms (like divergence rate) confusing, not necessary, and not insightful about the mechanism. Yes, the terms are defined, but they are just masking the reality, which is simpler: constrained usage of amino acids in positions of enzyme molecules. Not all 20, not even 10, but 4! Why not state and illustrate this clearly? Don't you agree that the impact of the paper will increase because it will be easier to understand?

We appreciate the reviewer’s opinion. As we described above, we now further highlight the aforementioned result about the number of acceptable amino acids. Nevertheless, as we stated previously, this is not the only main result in the paper. There are many other interesting results, analyses, and experiments described in the manuscript. We understand that the reviewer particularly appreciates this result. We do too! But, based on the comments of other reviewers and editors, and discussions with our colleagues, many researchers are likely to be interested in other results as well. Such as the existence of the sequence and structural divergence limits, temporal profiles of long-terms divergence, various structural and functional analyses, correlations with experimentally measured fitness effects, etc. Therefore, we cannot reduce the paper to one particular result or figure. We currently have 6 main figures in the manuscript, which is arguably not many, and, we believe, they can be easily understood by most readers interested in the topic. We hope that all of our analyses are useful and provide different and complementary insights to different researchers.

In terms of the “divergence rate”, we do not understand this comment. The terms “divergence rate” or “protein divergence rate” are in fact quite standard and popular in molecular evolution. Searches for these terms in Google or PubMed return many hundreds of usages in the same or similar context and in multiple related papers. We now also define it in the text, to prevent any possibility for misinterpretation or confusion.

The Abstract is still very poor and misleading. For instance, the statement "The effective divergence limit (>25% sequence identity) is not primarily due to multiple substitutions at the same sites" is completely wrong. According to the authors' explanations, this "effective divergence limit" is exactly due to multiple substitutions at the same site! If you have only 4 amino acids acceptable at a site, due to multiple substitutions sequence identity will saturate are 25%. Like in DNA. Why not write a precise and clear Abstract about what this work presents? Due to so many logical flaws in the authors' thinking and presentation (as illustrated above), I would be scared to publish this paper without a careful read. One statement at a time. And trying to figure out why the statement is wrong. If clearly not wrong, then move on to the next one.

We apologize for the possible confusion, but we do not believe that our statements are either wrong or illogical. Let’s consider the example given by the reviewer. It is true that DNA sequence identity saturates at 25% under the assumption that all nucleotide substitutions are equally accepted at DNA sites. In this case the limit is indeed directly driven by random and neutral substitutions. But in the case of proteins, this scenario, i.e. acceptance of all amino acids with the same probability, would lead to saturation at 5% sequence identity, and not at 25-40%, which is what we observe for proteins with conserved molecular functions. That is why we said that “the limit is not primarily driven by multiple substitutions at the same sites”. Obviously, substitutions at the same sites occur all the time, but the small number of acceptable substitutions is ultimately not the cause but the consequence of the functional constraints and the divergence limit. In any case, to eliminate any possibility of confusion here, following the reviewer’s suggestion, we rewrote the Abstract accordingly.

And, finally, why not compare your results with this paper more thoroughly PMID: 27138088? Isn't it a bit similar? I guess I don't understand the meaning behind "to investigate the temporal patterns of the long-term divergence."

Do the authors still claim that the sites are saturated at the usage of 3 to 5 amino acids because the time passed was not sufficient to gain more changes? Or the time was enough and protein of the same function simply cannot tolerate additional amino acids well and still keep the function? Which one is right? I got an impression it was the latter. And then what I said at the beginning of this review holds. If it is the former, it needs to be convincingly justified.

We apologize for possible confusion. By “the temporal patterns of the long-term divergence” we mean the changes in sequence and structural similarities between orthologs as a function of time. As we discuss in the paper (see trajectories in Figure 1 and Figure 5A), our results clearly demonstrate that most ancient orthologs are close to saturation, and are unlikely to diverge, at least on planetary timescales, beyond ~25% sequence identity.

In terms of the paper by Jack et al., as we state in the paper, Jack et al. explored the difference in short-term substitution rates as a function of sites’ distance to catalytic residues. In contrast, our study is focused on the global long-term sequence and structural divergence across billions of years of evolution, the observed divergence limit, its origins, and the related analyses, such as contribution of various sites. Importantly, as we already discussed in our previous reply, the differences in short-term evolutionary rates do not imply the existence of a divergence limit, do not inform about the timing at which the divergence of orthologs saturates, or explain the functional origins of this saturation.

If the authors disagree with my review, then I did not understand the paper, which is possible, and still suggests that the authors need to improve the presentation.

We believe that the reviewer understood the paper, and thank the reviewer again for the comments and suggestions.

Reviewer #3:

This version of the manuscript is a significant improvement over the previous version in terms of added details (e.g. experimental procedure of folA mutagenesis re now described in sufficient detail to be understandable and/or potentially reproducible).

The authors attempted to address my (and other reviewers) main concern by presenting the analysis in new Figure 5 that shows that full orthologues (all EC numbers coincide) are more sequence-constrained than partial orthologues (3 EC numbers coincide). However this analysis fails to control for different structural divergence between full and partial orthologues.

Therefore, the most problematic aspect of the analysis – that authors attribute observed sequence conservation in diverging clades to conservation of function between orthologues has not been fully addressed. A clear alternative explanation that such conservation is explained by maintenance of structure and stability regardless of function still stands.

The authors used a truism by pointing out that proteins of the same FOLD (i.e. topologically similar arrangement of elements of secondary structure) can diverge to 10% or less sequence ID again citing an old work with data collected on very limited number of structures available at that time.

However, their full functional orthologues are much more similar structurally than proteins sharing the same fold. The correct control which was suggested in my initial review has not been done satisfactorily. Specifically, the authors should compare their sets of proteins with sets that have similar degree of structural similarity (measured as distribution of TM-scores) but different function. There are such examples where structures are quite similar (TM-score-wise) but functions differ significantly. Good examples of that kind ate TIM-barrels which are almost exclusively enzymes with wide variety of specificity and Igfold proteins again with very conserved structures but broadly diverged functional annotations.

If the authors find that conservation of structure in functionally divergent proteins imposes less sequence divergence constraints than same degree of structural conservation in orthologues – that will be a clear demonstration of additional constraints imposed by functional conservation which is the main message of this work. In the absence of such analysis the current data does not justify the conclusion.

Following the reviewer’s comments, we now present the suggested analysis (Figure 6C) testing the contribution of function to sequence constrains independent of structural similarity. The results clearly demonstrate that, even at the same level of structural divergence (RMSD), there is significantly higher sequence similarity for pairs of orthologs sharing the full EC annotation, i.e. a higher degree of functional similarity. We have updated the Discussion to reflect this new result.

Importantly, we would like to emphasize that the main message of our work is not Figure 6C, but the fact that the conservation of specific molecular function limits how far orthologs can diverge both in sequence and structure (i.e. in terms of the RMSD). As we comment in the Discussion, this limit is likely due to the conservation of structural optimality (but not necessarily protein stability, see below), i.e. of the specific structural arrangements and dynamics of protein residues required for efficient catalysis and function. To further emphasize this point, we now show side by side (Figure 6A, 6B) the long-term sequence and structural divergence of orthologs sharing 1.) their full EC numbers, 2.) orthologs sharing only the first three digits of their EC numbers, and 3.) enzymes classified in the same structural fold, but with different molecular functions. These results clearly demonstrate that orthologs maintaining a less stringent level of functional conservation diverge significantly further in both sequence and structure.

Minor point: The authors severely misquote Firnberg et al., 2016. They say: "Nevertheless, direct and comprehensive biochemical experiments demonstrated that the deleterious effects of protein mutations primarily result from changes in specific protein activity rather than decreases in protein stability and cellular abundance [Firnberg et al., 2016]". In fact, the authors of Firnberg et al., 2016, say directly the opposite: '… These DFEs provide insight into the inherent benefits of the genetic code's architecture, support for the hypothesis that mRNA stability dictates codon usage at the beginning of genes, an extensive framework for understanding protein mutational tolerance, and evidence that mutational effects on protein thermodynamic stability shape the DFE..…" (cited from the Abstract of Firnberg et al., 2016). The authors misrepresent the main result of Firnberg et al., 2016: According to Figure 5B of [Firnberg et al., 2016] the product of abundance and catalytic activity shapes fitness effects in TEM1, not abundance alone. This is exactly what is established in Bershtein et al., 2015 and Adkar et al., 2017, and shows that abundance (which is a function of stability) enters fitness landscape on equal footing with catalytic activity, i.e. there is as much selection for abundance (i.e. stability) as it is for kcat/KM and related measures of activity. This misinterpretation somewhat undermines the major premise of the present paper that there is separate selection for activity/function and stability/structure. In reality one cannot disentangle the two because fitness landscape depends on the function (i.e. the product) of the two factors.

Frankly, we are quite surprised by the Reviewer’s comment. The reviewer suggests that we misquote the Firnberg et al. study and cites several sentences from their abstract. We are puzzled by this comment because immediately after the cited text, in fact in the next sentence (!) of the abstract, Fimberg et al. directly and clearly state: “Contrary to prevailing expectations, we find that deleterious effects of mutations primarily arise from a decrease in specific protein activity and not cellular protein levels”, which is exactly what we are saying in our paper.

Of course we do not dispute the fact that Fitness ~ Enzyme_Rate= kcat*[E] (under ligand saturation), which is textbook biochemical knowledge. Also we do not dispute that stability effects do contribute to protein evolution, we specifically said that in the text and previous reply. However, it is simply not true that because Fitness(x)=kcat(x)*[E(x)] (where “x” a specific mutation), it is not possible to understand the average relative contribution of random mutations to changes in fitness due to functional (kcat) or stability/abundance ([E]) effects. In direct mathematical terms: if P(x)=A(x)*B(x) it is clearly possible to disentangle the relative contribution of A and B to the variance in P due to changes in x. In fact, the experiments in Figure 5 of Firnberg et al. were specifically designed to do that, and they directly support the primary contributions of effects related to protein activity. The different, and in fact sometimes opposite, contributions of stability and activity to protein evolution are also illustrated by several examples of enzymatic optimization where protein stability is traded for an increase in catalytic efficiency (see for example PMID: 26940154, PMID: 20036254).

In our view, the discussion above is not a minor point. These results support our conceptual model of the observed limit. We also think the Reviewer maybe is confused that mutations can affect protein function (i.e. kcat) only through purely sequence effect. In fact, sequence substitutions affect the average protein structure and dynamics (but again not necessarily protein stability!) and make the structure not optimal for catalysis and function. Furthermore, we would like to emphasize again that there is no specific selection in evolution for protein stability, i.e. evolution does not just select proteins to be stable per se. Obviously, proteins need to have a certain level of stability for functioning, but, based on the Fimberg et al. results, the selection for protein stability is not the only and likely not the main target of functional selection.

https://doi.org/10.7554/eLife.39705.035

Article and author information

Author details

  1. Mariam M Konaté

    1. Department of Systems Biology, Columbia University, New York, United States
    2. Division of Cancer Treatment and Diagnosis, National Cancer Institute, Bethesda, United States
    Contribution
    Data curation, Software, Formal analysis, Investigation, Methodology, Writing—original draft
    Contributed equally with
    Germán Plata
    Competing interests
    No competing interests declared
  2. Germán Plata

    Department of Systems Biology, Columbia University, New York, United States
    Contribution
    Data curation, Software, Formal analysis, Investigation, Methodology, Writing—original draft, Writing—review and editing
    Contributed equally with
    Mariam M Konaté
    For correspondence
    gap2118@cumc.columbia.edu
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-6470-7748
  3. Jimin Park

    1. Department of Systems Biology, Columbia University, New York, United States
    2. Department of Pathology and Cell Biology, Columbia University, New York, United States
    Contribution
    Investigation, Methodology, Writing—review and editing
    Competing interests
    No competing interests declared
  4. Dinara R Usmanova

    Department of Systems Biology, Columbia University, New York, United States
    Contribution
    Formal analysis, Investigation, Methodology, Writing—review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-5031-0013
  5. Harris Wang

    1. Department of Systems Biology, Columbia University, New York, United States
    2. Department of Pathology and Cell Biology, Columbia University, New York, United States
    Contribution
    Supervision, Methodology, Project administration, Writing—review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-2164-4318
  6. Dennis Vitkup

    1. Department of Systems Biology, Columbia University, New York, United States
    2. Department of Biomedical Informatics, Columbia University, New York, United States
    Contribution
    Conceptualization, Supervision, Funding acquisition, Methodology, Writing—original draft, Project administration, Writing—review and editing
    For correspondence
    dv2121@cumc.columbia.edu
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-4259-8162

Funding

National Institute of General Medical Sciences (R01)

  • Dennis Vitkup

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

We sincerely thank Dan Tawfik, Eugene Koonin, and Fyodor Kondrashov for very helpful discussions. This work was supported in part by the National Institute of General Medical Sciences grants GM079759 and R35GM13188 to DV.

Senior Editor

  1. Diethard Tautz, Max-Planck Institute for Evolutionary Biology, Germany

Reviewing Editor

  1. Nir Ben-Tal, Tel Aviv University, Israel

Reviewer

  1. Nir Ben-Tal, Tel Aviv University, Israel

Publication history

  1. Received: June 29, 2018
  2. Accepted: August 7, 2019
  3. Version of Record published: September 18, 2019 (version 1)
  4. Version of Record updated: September 19, 2019 (version 2)

Copyright

© 2019, Konaté et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 1,295
    Page views
  • 110
    Downloads
  • 1
    Citations

Article citation count generated by polling the highest count across the following sources: PubMed Central, Crossref, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Download citations (links to download the citations from this article in formats compatible with various reference manager tools)

Open citations (links to open the citations from this article in various online reference manager services)

Further reading

    1. Computational and Systems Biology
    2. Evolutionary Biology
    Roman Sloutsky, Kristen M Naegle
    Tools and Resources
    1. Chromosomes and Gene Expression
    2. Evolutionary Biology
    Ryan Bracewell et al.
    Research Article Updated