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.

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

Download asset Open asset

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 R₀ 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 Y₀ in Equation 2.

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 Y₀ 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 Y₀ 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 Y₀ 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).

Figure 2

Download asset Open asset

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).

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).

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 Y₀ 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

Download asset Open asset

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

Download asset Open asset

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.

Figure 5

Download asset Open asset

Figure 5—source data 1

RMSD versus divergence times for proteins with the same enzymatic function.

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

Download elife-39705-fig5-data1-v2.txt

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⁻²⁰) and structural similarities (2.4 vs. 1.8 Å RMSD, P-value 2 × 10⁻¹⁸) 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⁻²⁰), 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 b, Dawson et al., 2017).

Figure 6

Download asset Open asset

Figure 6—source data 1

Sequence identities between ancient orthologs sharing the same EC number and only the first three digits of their EC numbers.

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

Download elife-39705-fig6-data1-v2.txt

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.

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).

Sequence reads were deposited to the SRA database with accession number: SRP152339. Source data files have been provided for figure 1, figure 5, figure 6 and figure 1 supplements 1, 3, 5 and 6.

1. 1. Konaté MM
   2. Plata G
   3. Park J
   4. Usmanova DR
   5. Wang HH
   6. Vitkup D

   (2018) SRA

   ID SRP152339. Molecular function limits divergent protein evolution on planetary timescales.

   https://www.ncbi.nlm.nih.gov/sra/SRP152339

1. 1. Eric D Kelsic
   2. Hattie Chung
   3. Niv Cohen
   4. Jimin Park
   5. Harris H Wang
   6. Roy Kishony

   (2016) Supplementary material of: https://doi.org/10.1016/j.cels.2016.11.004

   Data S1.

   https://doi.org/10.1016/j.cels.2016.11.004

1. 1. Adkar BV
   2. Manhart M
   3. Bhattacharyya S
   4. Tian J
   5. Musharbash M
   6. Shakhnovich EI

   (2017) Optimization of lag phase shapes the evolution of a bacterial enzyme

   Nature Ecology & Evolution 1:149.

   https://doi.org/10.1038/s41559-017-0149

   • PubMed
   • Google Scholar
2. 1. Altschul SF
   2. Madden TL
   3. Schäffer AA
   4. Zhang J
   5. Zhang Z
   6. Miller W
   7. Lipman DJ

   (1997) Gapped BLAST and PSI-BLAST: a new generation of protein database search programs

   Nucleic Acids Research 25:3389–3402.

   https://doi.org/10.1093/nar/25.17.3389

   • PubMed
   • Google Scholar
3. 1. Arnold C
   2. Stadler PF

   (2010) Polynomial algorithms for the maximal pairing problem: efficient phylogenetic targeting on arbitrary trees

   Algorithms for Molecular Biology 5:25.

   https://doi.org/10.1186/1748-7188-5-25

   • PubMed
   • Google Scholar
4. 1. Bairoch A

   (1999) The ENZYME data bank in 1999

   Nucleic Acids Research 27:310–311.

   https://doi.org/10.1093/nar/27.1.310

   • PubMed
   • Google Scholar
5. 1. Barreto FS
   2. Burton RS

   (2013) Evidence for compensatory evolution of ribosomal proteins in response to rapid divergence of mitochondrial rRNA

   Molecular Biology and Evolution 30:310–314.

   https://doi.org/10.1093/molbev/mss228

   • PubMed
   • Google Scholar
6. 1. Battistuzzi FU
   2. Feijao A
   3. Hedges SB

   (2004) A genomic timescale of prokaryote evolution: insights into the origin of methanogenesis, Phototrophy, and the colonization of land

   BMC Evolutionary Biology 4:44.

   https://doi.org/10.1186/1471-2148-4-44

   • PubMed
   • Google Scholar
7. 1. Bell EA
   2. Boehnke P
   3. Harrison TM
   4. Mao WL

   (2015) Potentially biogenic carbon preserved in a 4.1 billion-year-old zircon

   PNAS 112:14518–14521.

   https://doi.org/10.1073/pnas.1517557112

   • PubMed
   • Google Scholar
8. 1. Benkovic SJ
   2. Fierke CA
   3. Naylor AM

   (1988) Insights into enzyme function from studies on mutants of dihydrofolate reductase

   Science 239:1105–1110.

   https://doi.org/10.1126/science.3125607

   • PubMed
   • Google Scholar
9. 1. Berman HM
   2. Westbrook J
   3. Feng Z
   4. Gilliland G
   5. Bhat TN
   6. Weissig H
   7. Shindyalov IN
   8. Bourne PE

   (2000) The protein data bank

   Nucleic Acids Research 28:235–242.

   https://doi.org/10.1093/nar/28.1.235

   • PubMed
   • Google Scholar
10. 1. Bershtein S
    2. Serohijos AW
    3. Bhattacharyya S
    4. Manhart M
    5. Choi JM
    6. Mu W
    7. Zhou J
    8. Shakhnovich EI

    (2015) Protein homeostasis imposes a barrier on functional integration of horizontally transferred genes in Bacteria

    PLOS Genetics 11:e1005612.

    https://doi.org/10.1371/journal.pgen.1005612

    • PubMed
    • Google Scholar
11. 1. Betts HC
    2. Puttick MN
    3. Clark JW
    4. Williams TA
    5. Donoghue PCJ
    6. Pisani D

    (2018) Integrated genomic and fossil evidence illuminates life's early evolution and eukaryote origin

    Nature Ecology & Evolution 2:1556–1562.

    https://doi.org/10.1038/s41559-018-0644-x

    • PubMed
    • Google Scholar
12. 1. Bork P
    2. Sander C
    3. Valencia A

    (1993) Convergent evolution of similar enzymatic function on different protein folds: the hexokinase, ribokinase, and galactokinase families of sugar kinases

    Protein Science 2:31–40.

    https://doi.org/10.1002/pro.5560020104

    • PubMed
    • Google Scholar
13. 1. Breen MS
    2. Kemena C
    3. Vlasov PK
    4. Notredame C
    5. Kondrashov FA

    (2012) Epistasis as the primary factor in molecular evolution

    Nature 490:535–538.

    https://doi.org/10.1038/nature11510

    • PubMed
    • Google Scholar
14. 1. Chan YH
    2. Venev SV
    3. Zeldovich KB
    4. Matthews CR

    (2017) Correlation of fitness landscapes from three orthologous TIM barrels originates from sequence and structure constraints

    Nature Communications 8:14614.

    https://doi.org/10.1038/ncomms14614

    • PubMed
    • Google Scholar
15. 1. Chang A
    2. Schomburg I
    3. Placzek S
    4. Jeske L
    5. Ulbrich M
    6. Xiao M
    7. Sensen CW
    8. Schomburg D

    (2015) BRENDA in 2015: exciting developments in its 25th year of existence

    Nucleic Acids Research 43:D439–D446.

    https://doi.org/10.1093/nar/gku1068

    • PubMed
    • Google Scholar
16. 1. Chothia C
    2. Lesk AM

    (1986) The relation between the divergence of sequence and structure in proteins

    The EMBO Journal 5:823–826.

    https://doi.org/10.1002/j.1460-2075.1986.tb04288.x

    • Google Scholar
17. 1. Conant GC
    2. Wolfe KH

    (2008) Turning a hobby into a job: how duplicated genes find new functions

    Nature Reviews Genetics 9:938–950.

    https://doi.org/10.1038/nrg2482

    • PubMed
    • Google Scholar
18. 1. Dawson NL
    2. Lewis TE
    3. Das S
    4. Lees JG
    5. Lee D
    6. Ashford P
    7. Orengo CA
    8. Sillitoe I

    (2017) CATH: an expanded resource to predict protein function through structure and sequence

    Nucleic Acids Research 45:D289–D295.

    https://doi.org/10.1093/nar/gkw1098

    • Google Scholar
19. 1. DeSantis TZ
    2. Hugenholtz P
    3. Larsen N
    4. Rojas M
    5. Brodie EL
    6. Keller K
    7. Huber T
    8. Dalevi D
    9. Hu P
    10. Andersen GL

    (2006) Greengenes, a chimera-checked 16S rRNA gene database and workbench compatible with ARB

    Applied and Environmental Microbiology 72:5069–5072.

    https://doi.org/10.1128/AEM.03006-05

    • PubMed
    • Google Scholar
20. 1. Dickerson RE

    (1971) The structures of cytochrome c and the rates of molecular evolution

    Journal of Molecular Evolution 1:26–45.

    https://doi.org/10.1007/BF01659392

    • PubMed
    • Google Scholar
21. 1. Dokholyan NV
    2. Shakhnovich B
    3. Shakhnovich EI

    (2002) Expanding protein universe and its origin from the biological big bang

    PNAS 99:14132–14136.

    https://doi.org/10.1073/pnas.202497999

    • PubMed
    • Google Scholar
22. Software

    1. Felsenstein J

    (2005) PHYLIP (Phylogeny inference package) , version 3.6

    University of Washington, Seattle.

    http://evolution.genetics.washington.edu/phylip.html
23. 1. Firnberg E
    2. Labonte JW
    3. Gray JJ
    4. Ostermeier M

    (2016) A comprehensive, High-Resolution map of a gene's Fitness Landscape

    Molecular Biology and Evolution 33:1378.

    https://doi.org/10.1093/molbev/msw021

    • PubMed
    • Google Scholar
24. 1. Galperin MY
    2. Walker DR
    3. Koonin EV

    (1998) Analogous enzymes: independent inventions in enzyme evolution

    Genome Research 8:779–790.

    https://doi.org/10.1101/gr.8.8.779

    • PubMed
    • Google Scholar
25. 1. Galperin MY
    2. Makarova KS
    3. Wolf YI
    4. Koonin EV

    (2015) Expanded microbial genome coverage and improved protein family annotation in the COG database

    Nucleic Acids Research 43:D261–D269.

    https://doi.org/10.1093/nar/gku1223

    • PubMed
    • Google Scholar
26. 1. Gilson AI
    2. Marshall-Christensen A
    3. Choi JM
    4. Shakhnovich EI

    (2017) The role of evolutionary selection in the dynamics of protein structure evolution

    Biophysical Journal 112:1350–1365.

    https://doi.org/10.1016/j.bpj.2017.02.029

    • PubMed
    • Google Scholar
27. 1. Guy L
    2. Ettema TJ

    (2011) The archaeal 'TACK' superphylum and the origin of eukaryotes

    Trends in Microbiology 19:580–587.

    https://doi.org/10.1016/j.tim.2011.09.002

    • PubMed
    • Google Scholar
28. 1. Halabi N
    2. Rivoire O
    3. Leibler S
    4. Ranganathan R

    (2009) Protein sectors: evolutionary units of three-dimensional structure

    Cell 138:774–786.

    https://doi.org/10.1016/j.cell.2009.07.038

    • PubMed
    • Google Scholar
29. 1. Hedges SB
    2. Dudley J
    3. Kumar S

    (2006) TimeTree: a public knowledge-base of divergence times among organisms

    Bioinformatics 22:2971–2972.

    https://doi.org/10.1093/bioinformatics/btl505

    • PubMed
    • Google Scholar
30. 1. Hubble E

    (1929) A relation between distance and radial velocity among extra-galactic nebulae

    PNAS 15:168–173.

    https://doi.org/10.1073/pnas.15.3.168

    • PubMed
    • Google Scholar
31. 1. Jack BR
    2. Meyer AG
    3. Echave J
    4. Wilke CO

    (2016) Functional sites induce Long-Range evolutionary constraints in enzymes

    PLOS Biology 14:e1002452.

    https://doi.org/10.1371/journal.pbio.1002452

    • PubMed
    • Google Scholar
32. Book

    1. Jukes TH
    2. Cantor CR

    (1969)

    Evolution of Protein Molecules

    In: Munro H. N, editors. Mammalian Protein Metabolism. Academic Press. pp. 21–132.

    • Google Scholar
33. 1. Kanehisa M
    2. Sato Y
    3. Kawashima M
    4. Furumichi M
    5. Tanabe M

    (2016) KEGG as a reference resource for gene and protein annotation

    Nucleic Acids Research 44:D457–D462.

    https://doi.org/10.1093/nar/gkv1070

    • PubMed
    • Google Scholar
34. 1. Kelsic ED
    2. Chung H
    3. Cohen N
    4. Park J
    5. Wang HH
    6. Kishony R

    (2016) RNA structural determinants of optimal codons revealed by MAGE-Seq

    Cell Systems 3:563–571.

    https://doi.org/10.1016/j.cels.2016.11.004

    • PubMed
    • Google Scholar
35. 1. Kristensen DM
    2. Kannan L
    3. Coleman MK
    4. Wolf YI
    5. Sorokin A
    6. Koonin EV
    7. Mushegian A

    (2010) A low-polynomial algorithm for assembling clusters of orthologous groups from intergenomic symmetric best matches

    Bioinformatics 26:1481–1487.

    https://doi.org/10.1093/bioinformatics/btq229

    • PubMed
    • Google Scholar
36. 1. Ku C
    2. Nelson-Sathi S
    3. Roettger M
    4. Sousa FL
    5. Lockhart PJ
    6. Bryant D
    7. Hazkani-Covo E
    8. McInerney JO
    9. Landan G
    10. Martin WF

    (2015) Endosymbiotic origin and differential loss of eukaryotic genes

    Nature 524:427–432.

    https://doi.org/10.1038/nature14963

    • PubMed
    • Google Scholar
37. 1. Langmead B
    2. Salzberg SL

    (2012) Fast gapped-read alignment with bowtie 2

    Nature Methods 9:357–359.

    https://doi.org/10.1038/nmeth.1923

    • PubMed
    • Google Scholar
38. 1. Larkin MA
    2. Blackshields G
    3. Brown NP
    4. Chenna R
    5. McGettigan PA
    6. McWilliam H
    7. Valentin F
    8. Wallace IM
    9. Wilm A
    10. Lopez R
    11. Thompson JD
    12. Gibson TJ
    13. Higgins DG

    (2007) Clustal W and clustal X version 2.0

    Bioinformatics 23:2947–2948.

    https://doi.org/10.1093/bioinformatics/btm404

    • PubMed
    • Google Scholar
39. 1. Lee D
    2. Redfern O
    3. Orengo C

    (2007) Predicting protein function from sequence and structure

    Nature Reviews Molecular Cell Biology 8:995–1005.

    https://doi.org/10.1038/nrm2281

    • PubMed
    • Google Scholar
40. Book

    1. Lehninger AL
    2. Nelson DL
    3. Cox MM

    (2013)

    Lehninger Principles of Biochemistry (6th Edition)

    New York: W.H. Freeman.

    • Google Scholar
41. 1. Lunzer M
    2. Golding GB
    3. Dean AM

    (2010) Pervasive cryptic epistasis in molecular evolution

    PLOS Genetics 6:e1001162.

    https://doi.org/10.1371/journal.pgen.1001162

    • PubMed
    • Google Scholar
42. 1. Meyer TE
    2. Cusanovich MA
    3. Kamen MD

    (1986) Evidence against use of bacterial amino acid sequence data for construction of all-inclusive phylogenetic trees

    PNAS 83:217–220.

    https://doi.org/10.1073/pnas.83.2.217

    • Google Scholar
43. 1. Murzin AG
    2. Brenner SE
    3. Hubbard T
    4. Chothia C

    (1995) SCOP: a structural classification of proteins database for the investigation of sequences and structures

    Journal of Molecular Biology 247:536–540.

    https://doi.org/10.1016/S0022-2836(05)80134-2

    • PubMed
    • Google Scholar
44. 1. O'Malley-James JT
    2. Cockell CS
    3. Greaves JS
    4. Raven JA

    (2014) Swansong biospheres II: the final signs of life on terrestrial planets near the end of their habitable lifetimes

    International Journal of Astrobiology 13:229–243.

    https://doi.org/10.1017/S1473550413000426

    • Google Scholar
45. 1. Omelchenko MV
    2. Galperin MY
    3. Wolf YI
    4. Koonin EV

    (2010) Non-homologous isofunctional enzymes: A systematic analysis of alternative solutions in enzyme evolution

    Biology Direct 5:31.

    https://doi.org/10.1186/1745-6150-5-31

    • Google Scholar
46. 1. Ota T
    2. Nei M

    (1994) Estimation of the number of amino acid substitutions per site when the substitution rate varies among sites

    Journal of Molecular Evolution 38:642–643.

    https://doi.org/10.1007/BF00175885

    • Google Scholar
47. 1. Parfrey LW
    2. Lahr DJ
    3. Knoll AH
    4. Katz LA

    (2011) Estimating the timing of early eukaryotic diversification with multigene molecular clocks

    PNAS 108:13624–13629.

    https://doi.org/10.1073/pnas.1110633108

    • PubMed
    • Google Scholar
48. 1. Pingoud A
    2. Wilson GG
    3. Wende W

    (2014) Type II restriction endonucleases--a historical perspective and more

    Nucleic Acids Research 42:7489–7527.

    https://doi.org/10.1093/nar/gku447

    • PubMed
    • Google Scholar
49. Software

    1. Plata G

    (2018) Count-MAGE-Seq, version 4f827f6

    GitHub.

    https://github.com/platyias/count-MAGE-seq
50. 1. Porter CT

    (2004) The catalytic site atlas: a resource of catalytic sites and residues identified in enzymes using structural data

    Nucleic Acids Research 32:129D–133.

    https://doi.org/10.1093/nar/gkh028

    • Google Scholar
51. 1. Povolotskaya IS
    2. Kondrashov FA

    (2010) Sequence space and the ongoing expansion of the protein universe

    Nature 465:922–926.

    https://doi.org/10.1038/nature09105

    • PubMed
    • Google Scholar
52. 1. Rost B

    (1997) Protein structures sustain evolutionary drift

    Folding and Design 2:S19–S24.

    https://doi.org/10.1016/S1359-0278(97)00059-X

    • PubMed
    • Google Scholar
53. 1. Rost B

    (2002) Enzyme function less conserved than anticipated

    Journal of Molecular Biology 318:595–608.

    https://doi.org/10.1016/S0022-2836(02)00016-5

    • PubMed
    • Google Scholar
54. 1. Shendure J
    2. Akey JM

    (2015) The origins, determinants, and consequences of human mutations

    Science 349:1478–1483.

    https://doi.org/10.1126/science.aaa9119

    • PubMed
    • Google Scholar
55. 1. Sheridan PP
    2. Freeman KH
    3. Brenchley JE

    (2003) Estimated minimal divergence times of the major bacterial and archaeal phyla

    Geomicrobiology Journal 20:1–14.

    https://doi.org/10.1080/01490450303891

    • Google Scholar
56. 1. Spang A
    2. Saw JH
    3. Jørgensen SL
    4. Zaremba-Niedzwiedzka K
    5. Martijn J
    6. Lind AE
    7. van Eijk R
    8. Schleper C
    9. Guy L
    10. Ettema TJG

    (2015) Complex archaea that bridge the gap between prokaryotes and eukaryotes

    Nature 521:173–179.

    https://doi.org/10.1038/nature14447

    • PubMed
    • Google Scholar
57. 1. Stamatakis A

    (2014) RAxML version 8: a tool for phylogenetic analysis and post-analysis of large phylogenies

    Bioinformatics 30:1312–1313.

    https://doi.org/10.1093/bioinformatics/btu033

    • PubMed
    • Google Scholar
58. 1. Tajima F
    2. Nei M

    (1984) Estimation of evolutionary distance between nucleotide sequences

    Molecular Biology and Evolution 1:269–285.

    https://doi.org/10.1093/oxfordjournals.molbev.a040317

    • PubMed
    • Google Scholar
59. 1. Tatusov RL
    2. Koonin EV
    3. Lipman DJ

    (1997) A genomic perspective on protein families

    Science 278:631–637.

    https://doi.org/10.1126/science.278.5338.631

    • PubMed
    • Google Scholar
60. 1. Thiergart T
    2. Landan G
    3. Schenk M
    4. Dagan T
    5. Martin WF

    (2012) An evolutionary network of genes present in the eukaryote common ancestor polls genomes on eukaryotic and mitochondrial origin

    Genome Biology and Evolution 4:466–485.

    https://doi.org/10.1093/gbe/evs018

    • PubMed
    • Google Scholar
61. 1. Tian W
    2. Skolnick J

    (2003) How well is enzyme function conserved as a function of pairwise sequence identity?

    Journal of Molecular Biology 333:863–882.

    https://doi.org/10.1016/j.jmb.2003.08.057

    • PubMed
    • Google Scholar
62. 1. UniProt Consortium

    (2015) UniProt: a hub for protein information

    Nucleic Acids Research 43:D204–D212.

    https://doi.org/10.1093/nar/gku989

    • PubMed
    • Google Scholar
63. 1. Wang HH
    2. Isaacs FJ
    3. Carr PA
    4. Sun ZZ
    5. Xu G
    6. Forest CR
    7. Church GM

    (2009) Programming cells by multiplex genome engineering and accelerated evolution

    Nature 460:894–898.

    https://doi.org/10.1038/nature08187

    • PubMed
    • Google Scholar
64. 1. Watson JD
    2. Laskowski RA
    3. Thornton JM

    (2005) Predicting protein function from sequence and structural data

    Current Opinion in Structural Biology 15:275–284.

    https://doi.org/10.1016/j.sbi.2005.04.003

    • PubMed
    • Google Scholar
65. 1. Wilson CA
    2. Kreychman J
    3. Gerstein M

    (2000) Assessing annotation transfer for genomics: quantifying the relations between protein sequence, structure and function through traditional and probabilistic scores

    Journal of Molecular Biology 297:233–249.

    https://doi.org/10.1006/jmbi.2000.3550

    • PubMed
    • Google Scholar
66. 1. Worth CL
    2. Gong S
    3. Blundell TL

    (2009) Structural and functional constraints in the evolution of protein families

    Nature Reviews Molecular Cell Biology 10:709–720.

    https://doi.org/10.1038/nrm2762

    • PubMed
    • Google Scholar
67. 1. Wu M
    2. Scott AJ

    (2012) Phylogenomic analysis of bacterial and archaeal sequences with AMPHORA2

    Bioinformatics 28:1033–1034.

    https://doi.org/10.1093/bioinformatics/bts079

    • Google Scholar
68. 1. Yang Z
    2. Nielsen R
    3. Goldman N
    4. Pedersen AM

    (2000)

    Codon-substitution models for heterogeneous selection pressure at amino acid sites

    Genetics 155:431–449.

    • PubMed
    • Google Scholar
69. Book

    1. Yang Z

    (2006)

    Computational Molecular Evolution

    Oxford: Oxford University Press.

    • Google Scholar
70. 1. Ye Y
    2. Godzik A

    (2003) Flexible structure alignment by chaining aligned fragment pairs allowing twists

    Bioinformatics 19 Suppl 2:ii246–ii255.

    https://doi.org/10.1093/bioinformatics/btg1086

    • PubMed
    • Google Scholar
71. 1. Zhang Y
    2. Skolnick J

    (2005) TM-align: a protein structure alignment algorithm based on the TM-score

    Nucleic Acids Research 33:2302–2309.

    https://doi.org/10.1093/nar/gki524

    • PubMed
    • Google Scholar
72. Book

    1. Zuckerkandl E
    2. Pauling L

    (1965)

    Evolutionary Divergence and Convergence in Proteins

    In: Bryson V, Vogel H. J, editors. Evolving Genes and Proteins. Academic Press. pp. 97–166.

    • Google Scholar

[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.

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 μ₀ =μ × l^-1:

y=1a+(1-1/a)(1+((-a/(a-1))/l))^lμ₀t

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

y=1/a+(1-1/a)exp(-R₀t); with R0=(a/(a-1))μ0

This equation is numerically identical to our model 2. Furthermore, this demonstrates that the parameter Y₀ 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 R² 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 μ₀ = μ × l^-1:

y=1/a+(1-1/a)(1+(-a/(a-1))/l)^lμ₀t

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

y=1/a+(1-1/a)exp(-R₀t); with R0=(a/(a-1))μ₀

This equation is numerically identical to our model 2. Furthermore, this demonstrates that the parameter Y₀ 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.

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

   "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

   "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

   "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

   "This ORCID iD identifies the author of this article:" 0000-0003-4259-8162

• 6,513

  Page views

• 444

  Downloads

• 15

  Citations

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