Introduction
============

"Prediction is very difficult, especially about the future. Niels Bohr"

While commonly attributed to Niels Bohr, it is difficult to determine with confidence the primary source of the above quote, demonstrating that the statement is also true for predictions involving the past.  Predicting the future is difficult because, in lieu of direct observations, we must extrapolate based on present-day information.  For the same reason, it is also difficult to reconstruct past events that occurred sufficiently long ago that little or no direct record of them remains.  This is often the case in evolutionary biology, which studies the past to understand the present.  As the past cannot be directly observed, we must rely on methods, such as phylogenetics, that make inferences about the past to describe the patterns and comprehend the processes that have shaped biodiversity. 

Reconstructing past evolution is hard because it is difficult to disentangle signal from noise, and because our understanding of the biological process is imperfect.  Moreover, inferences about the past can almost never be validated experimentally.  The validation of phylogenetic methods and evolutionary models depends almost exclusively on simulations ([@msu305-B3]).  Such in silico experiments can readily produce simulated data based on models of evolutionary processes.  Phylogenetic methods or models can then be tested and/or compared in their ability to accurately reconstruct evolutionary events or estimate parameters of the evolutionary process that generated the simulated sequences.  However, our models of evolutionary process are overly simplistic and by extension limited in their ability to reproduce the emergent properties of complex systems ([@msu305-B42]; [@msu305-B2]). 

For instance, phylogenetic methods and evolutionary models can be used to infer the ancestral molecular sequences of extant protein-coding genes (at the DNA or protein level) ([@msu305-B40]; [@msu305-B55]; [@msu305-B24]).  Although the performance of these methods and models in terms of ancestral sequence reconstruction (ASR) can be evaluated through simulation experiments, current models cannot anticipate the emergent properties of protein ancestors in their native state (correct folding, 3D structure, enzymatic characteristics, etc).  It is only when these ancestors are resurrected in vitro or in vivo and their functionality is verified that one can make an evaluation of the performance of evolutionary models.  In this article, we experimentally validate the computational predictions of the relative performance of evolutionary models in terms of ASR accuracy in order to answer two questions: 1) are more complex evolutionary models able to infer more accurate historical trajectories of proteins and 2) if so, do these improved genotypes translate into more accurate phenotypes? 

Ancestral protein resurrection holds great potential for understanding how evolutionary processes and biochemical properties interplay to produce the structures and functions of extant proteins ([@msu305-B6]; [@msu305-B23], [@msu305-B24]).  Fifty years ago, [@msu305-B40] proposed that the resurrection of ancestral sequences inferred in silico could open the possibility of experimentally studying the ancestors of modern proteins.  This is possible because, given a set of homologous sequences, a corresponding phylogenetic tree, and a model of sequence evolution, one can infer ancestral sequences for any node of the phylogeny.  These putative ancestral sequences can then be "resurrected" in the laboratory using standard molecular biology techniques, giving access to extinct proteins and their phenotypes.  Since the work of [@msu305-B37] and [@msu305-B46], who first implemented this idea in practice, numerous studies combining ASR with experimental resurrection have investigated diverse biological questions, ranging from ancient adaptations to temperature ([@msu305-B17], [@msu305-B16]; [@msu305-B25]), to ancestral ecological adaptations ([@msu305-B7]; [@msu305-B38]), the emergence of protein function ([@msu305-B4]; [@msu305-B39]), the influence of gene duplication on functional divergence ([@msu305-B51]), the evolution of molecular complexes ([@msu305-B15]), and industrial, technological, or biomedical applications of ancestral proteins ([@msu305-B28]; [@msu305-B8]; [@msu305-B9]). 

With the increase in popularity of the ASR approach, several methodological improvements have been proposed ([@msu305-B55]; [@msu305-B29]; [@msu305-B44]; [@msu305-B52]; [@msu305-B43]).  Using maximum likelihood (ML), [@msu305-B55] proposed the marginal reconstruction algorithm that we have employed in this study and which is used in almost all modern ASR studies.  With this approach, at a given site in the sequence alignment and at a given internal node, posterior probabilities (PPs) for all possible states are computed.  The state having the highest PP is considered as the ancestral state.  It is worth noting that PPs provide confidence in the reconstruction inference ([@msu305-B55]).  Despite the flexibility afforded by such a probabilistic approach, and the correspondingly wide range of available tree reconstruction algorithms, few studies ([@msu305-B22]) have focused on the effect of the phylogenetic tree on ASR. 

In most, if not all, previous studies where ASR and protein resurrection have been performed, ancestral sequences were inferred using a gene tree reconstructed using only the multiple sequence alignment of existing sequences ([@msu305-B23]); we refer to such gene trees as *Species-tree-unaware trees*, thereafter named *S-unaware trees*.  Individual sequences alone contain limited signal, and as a result phylogenetic reconstruction almost always involves choosing between statistically equivalent or weakly distinguishable relationships.  Furthermore, while each set of homologous genes has its own unique story, they are all related by a shared species history, which could be helpful for gene tree inference.  To exploit this possibility, genome evolutionary processes such as duplication, horizontal transfer, and loss must be modeled to reconcile the gene tree with the species tree ([@msu305-B47]).  The advantage of such "species tree aware" methods is that they allow the detection and the correction of tree reconstruction errors resulting from the finite size of alignments or the inadequacy of the substitution model employed, while at the same time retaining bona fide phylogenetic discord produced by genome evolutionary processes ([fig.  1](#msu305-F1){ref-type="fig"}).  In many simulation studies, methods that combine the substitution model with models of genome evolution to reconstruct *Species-tree-aware trees*, thereafter named *S-aware trees*, have been proved to increase the accuracy of gene trees ([@msu305-B1]; [@msu305-B45]; [@msu305-B5]; [@msu305-B48]; [@msu305-B53]).  F[ig]{.smallcaps}.  1.Schematic illustration of the impact of gene tree/species tree reconciliation on the topology of the gene tree.  The gene family under consideration evolves along the species tree shown in (*a*).  In this example, two copies of a gene are present in each of the species, genes in species A are denoted *a*~1~ and *a*~2~, genes in *B* by *b*~1~ and *b*~2~ etc.  The gene tree reconstructed using the traditional S-unaware method that optimizes a sequence evolution-based score is shown in (*b*).  It contains poorly supported phylogenetic relationships, including branches in red that conflict with the species tree.  On the left part of the tree in (*b*), the conflict is strongly supported by the sequence information and conserved in the reconciled tree shown in (*c*), suggesting that a horizontal gene transfer likely occurred between species B and C. In the middle part of the tree in (*b*), the support for a possible transfer is low (30/100), indicating that a gene tree with a higher joint likelihood probably exists.  In the reconciled S-aware tree (*c*) that optimizes a joint sequence evolution and gene family evolution score (here the joint likelihood that considers both the sequence substitutions and the duplication, transfer, and loss of genes), some of the conflicts are resolved (blue branches), because, while the sequence evolution component of joint likelihood is slightly lower, the gene family evolution component is significantly improved.  Reconciling the gene tree (*c*) with the species tree (*a*) requires a horizontal transfer and a gene duplication predating the divergence of species A, B, and C. In contrast, the reconciliation of the gene tree (*b*) requires, beside the statistically supported gene transfer, at least two losses, another transfer, and a duplication. 

The purpose of this study is to investigate to what extent both ASR and protein resurrection can benefit from the use of such biologically realistic models of tree reconstruction. 

Results
=======

Impact of the Phylogenetic Tree on ASR
--------------------------------------

We first investigated the influence of the phylogenetic tree reconstruction method on ASR accuracy through simulation experiments.  We evaluated the impact of using S-aware trees in comparison with S-unaware trees on ASR accuracy.  To do so, we made use of the data set of [@msu305-B48], comprising 1,099 gene families from 36 cyanobacterial genomes.  For each of these biological gene families, a reconciled tree was computed in their original study ([@msu305-B48]).  In this work, we randomly chose 100 families out the 1,099 and we simulated sequences along these reconciled tree topologies, thereafter considered as "true" gene trees.  To measure reconstruction accuracy, we considered both the raw and Grantham ([@msu305-B18]) distances when comparing inferred ancestral sequences to true sequences recorded during simulations (see Material and Methods).  As the patterns between the two distance metrics were highly similar, only raw distances are discussed in this article. 

With the 100 simulated alignments, the corresponding S-unaware trees were reconstructed either with PhyML ([@msu305-B21]) and the site-homogeneous LG model ([@msu305-B32]) or with PhyML-CAT ([@msu305-B34]) and the site-heterogeneous C60 model (originally used to simulate sequences, see Material and Methods).  To compute the S-aware trees, that is, reconciled gene trees that maximize the joint sequence-reconciliation likelihood, the amalgamated likelihood estimation (ALE) program ([@msu305-B48]) was used, with the cyanobacterial species tree computed by [@msu305-B47].  Ancestral sequences were then inferred along these reconstructed S-unaware or reconstructed S-aware trees, as well as along the "true trees."  For nodes defining similar monophyletic clades between the S-unaware or S-aware tree and the true tree, these ancestral sequences were compared with the true ancestral sequences recorded during the simulation. 

[Figure 2](#msu305-F2){ref-type="fig"}*a* shows that, on average, the S-unaware trees reconstructed either with LG or with C60 contain significantly more topological errors than the S-aware trees, in comparison with the true trees.  These results confirm the findings of [@msu305-B48], showing that S-aware trees are more accurate than S-unaware trees, even when they are reconstructed with the complex model used to simulate the sequences (C60).  Furthermore, this has a direct impact on the ASR accuracy: When ancestral sequences are reconstructed along the S-aware trees, the accuracy is greatly and significantly improved ([fig.  2](#msu305-F2){ref-type="fig"}*b*) and is close to the accuracy obtained with the true trees.  F[ig]{.smallcaps}.  2.Impact of the phylogenetic tree on ASR. (* a*) Phylogenetic reconstruction accuracy.  Robinson--Foulds distances were computed between S-unaware trees (LG or C60) or S-aware trees and the "true" tree.  The exODT model is the reconciliation model described in [@msu305-B49] (*b*) ASR accuracy depending on the phylogenetic tree.  Distances between inferred and true ancestral sequences were computed for nodes defining similar monophyletic clades between the S-unaware or S-aware tree and the true tree.  \*\*\**P* value \< 0.001; NS, nonsignificant. 

We then investigated the patterns of incorrectly inferred sites.  We only report results with the LG S-unaware trees and ALE S-aware trees, as results obtained with trees reconstructed with the C60 model are highly similar to those obtained with LG S-unaware trees.  We used the Grantham matrix to measure the biochemical properties of the differences between inferred and true amino acids.  The average Grantham scores of amino acid differences are 66.0 and 64.8 for the S-unaware and S-aware tree, respectively, and the overall distributions of Grantham Scores between S-unaware and S-aware trees are very close to each other ([supplementary fig.  S1](http://mbe.oxfordjournals.org/lookup/suppl/doi:10.1093/molbev/msu305/-/DC1), [Supplementary Material](http://mbe.oxfordjournals.org/lookup/suppl/doi:10.1093/molbev/msu305/-/DC1) online).  Nonetheless, this difference is statistically significant (Wilcoxon test, *P* value \< 0.001), indicating that S-unaware trees tend to lead to inference errors with more important biochemical consequences than S-aware trees.  The average score (65) corresponds to pairs of amino acids that have either a similar polarity and different molecular volumes (e.g., M--W) or the opposite case (e.g., S--D).  As expected, [supplementary figure S1](http://mbe.oxfordjournals.org/lookup/suppl/doi:10.1093/molbev/msu305/-/DC1) and [table S1](http://mbe.oxfordjournals.org/lookup/suppl/doi:10.1093/molbev/msu305/-/DC1), [Supplementary Material](http://mbe.oxfordjournals.org/lookup/suppl/doi:10.1093/molbev/msu305/-/DC1) online, show that a large proportions of amino acid differences concern amino acids with very similar properties (e.g., L--I or F--Y).  However, they also show that many reconstruction errors involve pairs of amino acids that are biochemically dissimilar, for example, L--H or A--Q. We observed that biochemical differences increase with the height of the internal node.  For instance, the average Grantham score par quartile of node height is 64.5, 65.0, 66.2, and 68.1 for the S-unaware trees.  All these characteristics of the distributions of Grantham scores highlight the impact that inference errors may have on the biochemical properties of resurrected proteins and the importance of favoring methods that increase the accuracy of the reconstruction. 

Reconstruction errors were also investigated in light of site-specific evolutionary rates.  With both S-unaware and S-aware trees, reconstruction errors occur more frequently in fast-evolving sites (Correlation test, $r^{2} = 0.2$, *P* value \< 0.001).  However, [figure 3](#msu305-F3){ref-type="fig"} clearly shows that the excess of errors obtained with S-unaware trees is not distributed uniformly with evolutionary rates.  Reconstruction errors tend to occur more frequently at slow-evolving sites (average posterior rate of 1.3) with the S-unaware trees in comparison with the S-aware trees (average posterior rate of 1.46, *P* value \< 0.001).  This demonstrates how topological errors can have a profound impact on ASR, as even conserved sites can be subject to wrong ancestral amino acid inferences.  F[ig]{.smallcaps}.  3.Reconstruction errors also affect slow-evolving sites with S-unaware trees.  For a given gene tree reconstruction method, the pattern of reconstruction errors for all internal nodes over the 100 simulations was analyzed in the light of site-specific (posterior) evolutionary rates.  Rates were computed a posteriori with the use of the Gamma distribution of site rates used during the ML reconstruction of ancestral sequences.  Orange: Species tree-unaware gene trees.  Blue: Species tree-aware gene trees.  The two data series are overlapped. 

We finally examined the PP for residues inferred differently with the S-unaware trees and with the S-aware trees.  The average PP reaches 0.82 and 0.81 for the S-unaware trees and S-aware trees, respectively.  This shows that the difference in inferences can involve residues that are unambiguously reconstructed with the S-unaware trees, and that the use of S-aware trees can radically change ancestral predictions. 

Resurrection and Experimental validation
----------------------------------------

We previously used the biochemical and biophysical properties of reconstructed ancestral LeuB enzymes to investigate thermal adaptation in *Bacillus* ([@msu305-B25]).  Furthermore, we used the biochemical and biophysical properties of the resurrected enzymes as a measure of their accuracy (e.g., a high Michaelis--Menten constant suggests a biologically unrealistic, and therefore inaccurate, ancestral enzyme).  Here, we have used the same approach to compare two versions of the same ancestral LeuB enzyme from the last common ancestor of the Firmicutes, the bacterial phylum to which *Bacillus* belongs.  These enzymes were inferred and resurrected to investigate the influence of the phylogenetic tree on potential biological conclusions regarding protein phenotypes.  The two enzymes were reconstructed either with the LeuB S-aware tree or with the LeuB S-unaware tree and are named LeuB*~S~*~--aw~ and LeuB*~S~*~--unaw~, respectively.  The ALE program, which was used to reconcile sequence and species information, detected 0 duplications, 14 lateral gene transfers, and 15 losses.  The S-aware tree has a Robinson--Foulds distance with the S-unaware tree equal to 32, which is very high.  The LeuB*~S~*~--aw~ and LeuB*~S~*~--unaw~ sequences differ by approximately 10% (36 amino acids).  Note that LeuB is the only enzyme on which we performed resurrections. 

The Michaelis--Menten constant (*K~M~*) for the substrate isopropylmalate (IPM) with LeuB*~S~*~--aw~ is similar to those measured for other thermophilic LeuB enzymes, such as the contemporary BCVX enzyme and the previously reconstructed thermophilic ancestors ANC1, ANC3, and ANC4 ([table 1](#msu305-T1){ref-type="table"}).  In contrast, the *K~M~*(IPM) for LeuB*~S~*~--unaw~ is about 4-fold higher, indicating its poorer affinity for this substrate ([table 1](#msu305-T1){ref-type="table"}).  Replicate *K~M~* determinations could not be performed for LeuB*~S~*~--aw~ and LeuB*~S~*~--unaw~ due to the expense of the substrate and the relatively high *K~M~* (IPM) of LeuB*~S~*~--unaw~, therefore we are unable to say whether the difference in *K~M~*(IPM) is statistically significant; however, a comparison of the Michaelis--Menten plots for these two enzymes ([supplementary fig.  S2](http://mbe.oxfordjournals.org/lookup/suppl/doi:10.1093/molbev/msu305/-/DC1), [Supplementary Material](http://mbe.oxfordjournals.org/lookup/suppl/doi:10.1093/molbev/msu305/-/DC1) online) illustrates that the difference in substrate affinity is considerable.  Furthermore, the *K~M~*(IPM) of LeuB*~S~*~--unaw~ is substantially higher than the highest *K~M~*(IPM) value in the BRENDA enzyme database ([www.brenda-enzyme.org](http://www.brenda-enzyme.org)).  In terms of turnover rate (*k*~cat~), LeuB*~S~*~--unaw~ exhibits a greater than 2-fold higher *k*~cat~ than LeuB*~S~*~--aw~. Although LeuB*~S~*~--unaw~ exhibits a high turnover rate, its high *K~M~* for IPM suggests that the substrate would have to be present at a very high concentration inside the cell for binding to actually occur.  Table 1.Kinetic Constants, Thermoactivity, and Biophysical Parameters for the Ancestral LeuB Enzyme from the Firmicutes Ancestor. Enzyme$K_{M}^{(\text{IPM})}$ (mM)$K_{M}^{(\text{NAD})}$ (mM)*k*~cat~ $(s^{- 1})$*T*~opt~ (°C)$\Delta G_{N - U}^{\ddagger}$ $(\text{kJ\ mol}^{- 1})$BPSYC0.20.66.54794.9 ± 0.2BSUB0.78.148.75395.9 ± 0.5BCVX1.10.853.869100.7 ± 0.2ANC11.30.5141.873100.9 ± 0.5ANC21.00.941.74991.1 ± 0.4ANC32.71.0102.36095.6 ± 0.1ANC41.71.0362.270110.8 ± 0.4*LeuB~S--aw~1.66.5181.285110.9* ± 1.6*LeuB*~S~*~--unaw~6.85.5441.27891.4* ± 0.6[^5]

The thermoactivity profiles of the two resurrected enzymes reveal that they are highly thermophilic with *T*~opt~ values greater than 75 °C ([table 1](#msu305-T1){ref-type="table"} and [fig.  4](#msu305-F4){ref-type="fig"}*a*).  We also determined the $\Delta G_{N - U}^{\ddagger}$ values for these enzymes, as we have previously found this parameter to be a useful measure of a biologically realistic enzyme ([@msu305-B25]).  $\Delta G_{N - U}^{\ddagger}$ indicates the conformational stability of a protein between its native (folded) and unfolded states and can be calculated from the measured unfolding rates of a protein in different concentrations of the denaturant urea.  Both LeuB*~S~*~--aw~ and LeuB*~S~*~--unaw~ are highly thermophilic, therefore they should exhibit some resistance to unfolding and have relatively high $\Delta G_{N - U}^{\ddagger}$ values ([supplementary fig.  S3](http://mbe.oxfordjournals.org/lookup/suppl/doi:10.1093/molbev/msu305/-/DC1), [Supplementary Material](http://mbe.oxfordjournals.org/lookup/suppl/doi:10.1093/molbev/msu305/-/DC1) online).  In accordance with its high *T*~opt~ value, LeuB*~S~*~--aw~ is very kinetically stable with a $\Delta G_{N - U}^{\ddagger}$ value of 110.9 kJ mol^−1^. In contrast, LeuB*~S~*~--unaw~ is thermophilic but unfolds rapidly in comparison with LeuB*~S~*~--aw~ ([fig.  4](#msu305-F4){ref-type="fig"}*b*) and is consequently kinetically unstable; its $\Delta G_{N - U}^{\ddagger}$ value of 91.4 kJ mol^−1^ is lower than that of contemporary and ancestral psychrophilic and mesophilic LeuB enzymes ([table 1](#msu305-T1){ref-type="table"}) and greater than 13 kJ mol^−1^ lower than would be predicted from its *T*~opt~ ([supplementary fig.  S3](http://mbe.oxfordjournals.org/lookup/suppl/doi:10.1093/molbev/msu305/-/DC1), [Supplementary Material](http://mbe.oxfordjournals.org/lookup/suppl/doi:10.1093/molbev/msu305/-/DC1) online).  As $\Delta G_{N - U}^{\ddagger}$ is related to the unfolding rate of a protein via an exponential function, a difference of 1 or 10 kJ mol^−1^ in $\Delta G_{N - U}^{\ddagger}$ equates to a 1.5-fold or 48-fold difference in the rate of unfolding, respectively.  The low kinetic stability of LeuB*~S~*~--unaw~ suggests that, while it is adapted to function at high temperatures, it would unfold rapidly in a thermophilic environment.  The structural/molecular reason(s) for the differences in $\Delta G_{N - U}^{\ddagger}$ and *K~M~*(IPM) between LeuB*~S~*~--unaw~ and LeuB*~S~*~--aw~ remains to be elucidated ([supplementary material](http://mbe.oxfordjournals.org/lookup/suppl/doi:10.1093/molbev/msu305/-/DC1) and [fig.  S4](http://mbe.oxfordjournals.org/lookup/suppl/doi:10.1093/molbev/msu305/-/DC1), [Supplementary Material](http://mbe.oxfordjournals.org/lookup/suppl/doi:10.1093/molbev/msu305/-/DC1) online).  Nonetheless, the kinetic instability of LeuB*~S~*~--unaw~ which is not concordant with its thermophilic adaptation, combined with its impaired *K~M~* for IPM, suggests that this enzyme is not biologically realistic and implies that its inferred sequence contains errors.  F[ig]{.smallcaps}.  4.Resurrection of ancestral LeuBs and impact of the phylogenetic tree on biological interpretations. (* a*) Thermoactivity profiles for ancestral LeuB enzymes.  Blue curve: LeuB*~S~*~--unaw~. Green curve: LeuB*~S~*~--aw~. (*b*) Unfolding rates of ancestral LeuB enzymes.  Unfolding rates are shown in 8 M urea as a decrease in intrinsic protein fluorescence with time.  Colors are the same as in (*a*). 

Discussion
==========

Our in silico investigations support that the use of an S-aware gene tree can have a profound impact on the inference of ancestral sequences.  This phylogenetic prediction is congruent with the conclusions obtained with our resurrection experiment, which suggest the need for reconciled gene trees (maximizing the joint sequence-reconciliation likelihood) to provide accurate substitution trajectories and ancestral protein phenotypes.  When the gene family under study has experienced a complex evolutionary history involving gene duplications, lateral transfers, and losses (such as L*euB*), it becomes necessary to account for these genomic events to reconstruct the tree along which ASR is performed.  Numerous methods that implement models of duplication, transfer, and loss of genes are now available to reconcile an S-unaware tree with a species tree ([@msu305-B1]; [@msu305-B11]; [@msu305-B12]; [@msu305-B45]; [@msu305-B48]; [@msu305-B53]).  Here, we demonstrate that the resulting gene tree is considerably more accurate than the original S-unaware tree and allows us to infer more accurately the history of protein evolution. 

Although the present results highlight how more complex evolutionary models improve ASR, potential limitations remain regarding hypotheses made by some methods that we used: Ancestral sequences were reconstructed in ML, with the marginal ASR approach ([@msu305-B55]).  With this approach, at a given position and at a given internal node, the state (amino acid in our case) having the highest PP is chosen as the ancestral state.  A well-known bias exists with this approach.  ML tends to assign to ancestral residues the state having the highest frequency at a given site ([@msu305-B54]).  With a simple contact potential used to calculate the free energy of protein ancestors of the purple acid phosphatase, [@msu305-B52] highlighted with simulations that because of this bias, ML may infer ancestral sequences that are biased toward thermostability.  Considering a Bayesian sampling approach, consisting of randomly drawing ancestral amino acids in the posterior distribution (instead of selecting the amino acid with the maximum probability), may be an ideal approach to bypass this bias.  Even though this result would need to be confirmed with further experiments (i.e., with a model allowing the 3D structure to change overtime or with the use of a more sophisticated energy potential capturing more appropriately the complexity of protein folding), and while this bias regarding thermostability has not been observed in our previous resurrection study ([@msu305-B25]), we cannot exclude that our results are not affected by a similar bias.  However, we do not anticipate any reason for which this bias would lead us to question our predictions, as we expect that the increase in ASR accuracy due to the use of S-aware trees would also apply with another ASR method, such as Bayesian Inference. To reconcile species and gene information, we used the ALE program ([@msu305-B48]).  For the moment, different uncertainties are not accounted for in the reconstruction of ancestral sequences along S-aware trees.  For instance, to what extent species tree reconstruction or incomplete lineage sorting impact ASR in our new methodology is currently unknown.  Plus, ALE makes use of a time-calibrated species tree to compute the probabilities of horizontal gene transfers during the reconstruction of S-aware trees.  The uncertainty in the estimation of species divergence times should also be accounted for in the reconstruction of ancestral sequences along S-aware trees.  Finally, despite the substantial increase in accuracy in gene tree reconstruction offered by ALE, we previously observed cases where S-unaware trees were more accurate than S-aware trees, in part due to overfitting of the species tree ([@msu305-B48]).  Although this occurred in a minority of cases (S-unaware trees were more accurate in 22.9% of cases), it could potentially impact ASR.  Additional experiments are required to investigate these specific issues. 

To date, protein resurrection studies have used species tree unaware methods of phylogenetic reconstruction methods---often producing unreliable gene trees.  Although these shortcomings did not necessarily prevent the resurrection of functional ancestors, this study suggests the potential dependence of biological conclusions regarding the phenotype of protein ancestors on the accuracy of the reconstructed phylogeny.  Further in vitro and/or in vivo investigations are needed to statistically confirm our preliminary results on additional proteins.  However, our study advocates the use of information on species history, in combination with state-of-the-art sequence evolution models ([@msu305-B19]) to accurately predict ancestral protein function and structure. 

Materials and Methods
=====================

Data Used for In Silico Experiments and Substitution Models
-----------------------------------------------------------

To perform in silico experiments to investigate the influence of the phylogenetic tree on ASR, we used the data set from [@msu305-B47].  This data set comprises 1,099 gene families from 36 cyanobacterial genomes available in the HOGENOM database ([@msu305-B41]).  The phylogenomic species tree of these 36 species that [@msu305-B47] reconstructed was also used in the present study.  With this species topology and a newly described model of gene tree/species tree reconciliation, [@msu305-B49] computed the reconciled S-aware trees for the 1,099 families.  Here, we randomly chose 100 families out the 1,099 and we simulated sequences along their corresponding S-aware tree topologies, that we considered as true gene trees.  On average, 2.17 duplications, 3.37 transfers, and 6.39 losses occurred along these true trees.  We added an outgroup species to both the species tree and true topologies.  The branch length leading to the outgroup species was set to one-half of the S-aware tree height. 

All models employed in this study are empirical Markovian substitution models and were all used in combination with a discrete Γ distribution to model the site-specific rate variation, with four categories. 

Simulations
-----------

Available substitution models may contain several parameters aiming at capturing molecular footprints left by biological processes during evolution.  Even so, they are too simplistic in comparison with the complexity of processes acting on biological data.  To mimic this gap between simplicity of substitution models and complexity of biological data, we used a relatively complex model to simulate sequences along the 100 true S-aware gene trees, and reconstructed phylogenetic trees and ancestral sequences with simpler models, described by a fewer number of parameters and constructed along different mathematical settings.  The site-heterogeneous C60 model ([@msu305-B34]), which is the most complex empirical substitution model currently available in the literature, was used to simulate data.  This model is a mixture of profiles, with a single Poisson exchangeability matrix that is assigned to all components (profiles) of the mixture (see [supplementary material](http://mbe.oxfordjournals.org/lookup/suppl/doi:10.1093/molbev/msu305/-/DC1), [Supplementary Material](http://mbe.oxfordjournals.org/lookup/suppl/doi:10.1093/molbev/msu305/-/DC1) online).  Consequently, this model contains 60 × 19 + 59 = 1,199 empirical parameters to describe the substitution process.  Alignments were simulated using the original alignment sizes of the 100 cyanobacterial families.  Simulations were performed with our own C++ program depending on Bio++ libraries ([@msu305-B20]).  For a given alignment, because sites are supposed to evolve independently, all 60 components of the mixture were used to simulate subalignments with a number of sites proportional to their empirical weight, with all subalignments being subsequently concatenated to produce the final alignment. 

Ancestral Sequence Reconstruction
---------------------------------

With the simulated data, ASR was performed with the UL3 mixture model [@msu305-B35], which is a mixture model containing fewer empirical parameters.  Indeed, contrarily to C60, which is a mixture of profiles, UL3 is a mixture of matrices (see [supplementary material](http://mbe.oxfordjournals.org/lookup/suppl/doi:10.1093/molbev/msu305/-/DC1), [Supplementary Material](http://mbe.oxfordjournals.org/lookup/suppl/doi:10.1093/molbev/msu305/-/DC1) online).  Each component of the UL3 mixture possesses its own exchangeability matrix and its own set of equilibrium frequencies.  This model possesses 3 × ($\frac{19 \times 20}{2} - 1$) + 3 × 19 + 2 = 626 empirical parameters, which is far less than the number of empirical parameters describing the C60 model. 

For both simulated and LeuB data, ML estimates of branch lengths and parameters of the substitution model were inferred with bppML, which belongs to the bppSuite of programs ([@msu305-B13]) and depends on Bio++ libraries ([@msu305-B20]).  The weight of each component of the mixture model was optimized by ML.  With all these ML estimates, ancestral sequences were then inferred with bppAncestor ([@msu305-B13]) using the marginal ASR approach ([@msu305-B55]).  A posteriori weight values of the mixture are used to perform ASR.  For a given site at a given internal node of the tree, the state having the maximum PP was inferred as the putative ancestral state. 

ASR Accuracy Measurement
------------------------

Inferred ancestral sequences were compared to true internal sequences by computing two distances: 1) the raw distance, which is simply the number of amino acid differences divided by the length of the sequence and 2) the Grantham distance ([@msu305-B18]), defined as the amino acid pair distance computed with the Grantham distance matrix, which takes into account biochemical similarities between amino acids in terms of polarity and volume. 

Gene Tree/Species Tree Reconciliations
--------------------------------------

[@msu305-B49] recently described a probabilistic reconciliation model that accounts for the duplication, transfer, and loss of genes along a species tree.  Given a fixed species tree, the model allows exploring possible paths along which a gene tree may have been generated by a series of speciations, duplications, transfers, and losses.  To efficiently explore the space of all reconciled trees according to the joint sequence-reconciliation likelihood that combines sequence information and information on the species phylogeny, [@msu305-B48] proposed the ALE algorithm.  ALE makes use of a sample of S-unaware gene trees (for instance, a sample of posterior trees produced by a Bayesian program such as PhyloBayes ([@msu305-B30])) to compute conditional clade probabilities ([@msu305-B26]), which are used to approximate the PP of all gene trees that can be amalgamated from clades present in the sample. 

ALE was used to perform all S-unaware gene tree/species tree reconciliations for both simulated and biological (see below) data sets.  For each simulated alignment, PhyloBayes (version 3.3f) was run to obtain an MCMC sample of trees using a simple F81 (Poisson) substitution model.  Two chains were run in parallel to check for convergence, with a burn-in of 1,000 samples followed by at least 10,000 samples.  These MCMC samples were then used by ALE to explore the space of reconciled trees in combination with the ML estimation of duplication, transfer, and loss rates, to eventually propose the *S-aware tree*---the reconciled gene tree that maximizes the joint sequence-reconciliation likelihood.  ALE calculations were performed with the calibrated species tree initially used to compute the true gene trees (see above). 

Experimental Resurrection of LeuB Enzymes
-----------------------------------------

### Firmicutes Species Tree and LeuB S-Unaware Tree Reconstructions

Firmicutes genomic sequences were downloaded from the NCBI, as of April 2012.  Orthologous gene families corresponding to all 53 bacterial ribosomal proteins were constructed with BLAST.  Each individual gene was aligned with Mafft ([@msu305-B27]) and ambiguous sites were trimmed by BMGE ([@msu305-B10]), using the BLOSUM30 matrix.  Only 46 out of the 53 ribosomal gene alignments were then concatenated.  The remaining seven genes (L25, L30, L32, L33, S4, S14, S21) were discarded owing to either the presence of paralogs or a patchy distribution over Firmicutes species.  To root both the species tree and the LeuB tree, we incorporated two outgroup *LeuB* sequences from two Actinobacteria species, *Corynebacterium glutamicum*, and *Streptomyces coelicolor*.  The final alignment contains 68 Firmicutes species, and the species tree ([supplementary fig.  S5](http://mbe.oxfordjournals.org/lookup/suppl/doi:10.1093/molbev/msu305/-/DC1), [Supplementary Material](http://mbe.oxfordjournals.org/lookup/suppl/doi:10.1093/molbev/msu305/-/DC1) online) was computed with PhyloBayes ([@msu305-B30]) using the CAT model ([@msu305-B31]).  Two independent chains were run in parallel to check for convergence.  The model of [@msu305-B49] used by ALE ([@msu305-B48]) to search for the S-aware gene tree needs divergence times between speciation nodes to compute the probabilities of gene transfers between branches.  Therefore, the species tree was calibrated with relative times using PhyloBayes and an arbitrary calibration of 1,000 time unit at the root.  The Log-normal autocorrelated relaxed clock model ([@msu305-B50]) was chosen to allow substitution rates to vary in time. 

The gene family corresponding to the 71 *LeuB* sequences found in the 68 species was reconstructed and a preliminary alignment was inferred using Muscle ([@msu305-B14]) and used to build a preliminary S-unaware phylogenetic tree using PhyML ([@msu305-B21]) with the LG model and a Γ distribution for rate variation.  This preliminary S-unaware tree was used as a guide tree in Prank ([@msu305-B36]) to realign LeuB sequences.  The final LeuB S-unaware tree along which ancestral sequences were reconstructed was computed with PhyloBayes, using the LG+Γ(4) model, and rooted on the branch between the Firmicutes and outgroup LeuBs ([supplementary fig.  S6](http://mbe.oxfordjournals.org/lookup/suppl/doi:10.1093/molbev/msu305/-/DC1), [Supplementary Material](http://mbe.oxfordjournals.org/lookup/suppl/doi:10.1093/molbev/msu305/-/DC1) online).  Three chains were run in parallel to ensure that convergence of the MCMC was reached. 

### LeuB S-Aware Gene Tree Reconstruction

We used the model described in [@msu305-B49] and implemented in the ALE program ([@msu305-B48]) to search for the ML S-aware reconciled tree, that is, the reconciled gene tree that maximizes the joint sequence-reconciliation likelihood ([supplementary fig.  S7](http://mbe.oxfordjournals.org/lookup/suppl/doi:10.1093/molbev/msu305/-/DC1), [Supplementary Material](http://mbe.oxfordjournals.org/lookup/suppl/doi:10.1093/molbev/msu305/-/DC1) online).  ALE used the sample of S-unaware trees produced by PhyloBayes (see above) and the calibrated species tree to compute the S-aware tree along which ASR was performed.  The S-aware tree was used as a guide tree in Prank to compute the final alignment. 

### Model Selection, Fit to the LeuB Data, and ASR

ASR of LeuB was performed with the site-heterogeneous EX_EHO mixture substitution model ([@msu305-B33]).  EX_EHO was deemed to be the best site-heterogeneous model at fitting the LeuB data according to the AIC criterion, in comparison with all other site-heterogeneous mixture models currently available in the literature ([@msu305-B34]; [@msu305-B35]; [@msu305-B33]) and implemented in the Bio++ libraries ([@msu305-B20]).  See [supplementary materials](http://mbe.oxfordjournals.org/lookup/suppl/doi:10.1093/molbev/msu305/-/DC1), [Supplementary Material](http://mbe.oxfordjournals.org/lookup/suppl/doi:10.1093/molbev/msu305/-/DC1) online, for information on the different site-homogeneous and site-heterogeneous mixture models.  As with simulations, ancestral sequences were inferred with bppAncestor ([@msu305-B13]).  When Prank was used to compute the final LeuB alignment, we used the "-anc" option to jointly infer ancestral gaps, which were subsequently incorporated into ancestral sequences inferred by bppAncestor.  This two-step approach mimics the one proposed in a previous publication ([@msu305-B15]), which makes use of the Fitch algorithm to a priori infer ancestral gap positions and then incorporates these gaps into ancestral sequences. 

### Protein Expression and Purification

Gene sequences for the two inferred versions of the ancestral Firmicutes LeuB were codon optimized for expression in *Escherichia coli* and chemically synthesised by Geneart (Life Technologies) with a 5′-*NcoI* site and a 3′-*PstI* site.  Following ligation of the genes into the protein expression vector pPROEX HTb, recombinant proteins were expressed in *E. coli* DH5*α* with 1 mM IPTG induction at 37 °C for 24 h. Proteins were purified to ≥95% purity by nickel affinity chromatography, and subsequent size-exclusion chromatography using the buffers detailed in [@msu305-B25].  Protein concentrations were determined using a NanoDrop 2000 (Thermo Scientific) and extinction coefficients calculated using ProtParam on the ExPASy server (web.expasy.org/protparam/). 

### LeuB Enzyme Characterization

LeuB activity was measured by following the reduction of NAD at 340 nm as described in [@msu305-B25].  The V~max~ and Michaelis--Menten constants for both substrates (IPM and NAD) were found using the Michaelis--Menten nonlinear fitting function in Graphpad Prism 6.  Thermoactivity profiles were determined by measuring the initial rate of activity at 1--5 °C intervals over a 20--30 °C temperature range in triplicate.  Thermoactivity profile reactions contained 15 mM IPM, 50 mM NAD, and 10-50 μM LeuB enzyme.  The free energy of unfolding, $\Delta G_{N - U}^{\ddagger}$, for each enzyme was determined from urea unfolding rates as described in [@msu305-B25]. 

Supplementary Material
======================

[Supplementary material](http://mbe.oxfordjournals.org/lookup/suppl/doi:10.1093/molbev/msu305/-/DC1), [table S1](http://mbe.oxfordjournals.org/lookup/suppl/doi:10.1093/molbev/msu305/-/DC1) and [figures S1--S7](http://mbe.oxfordjournals.org/lookup/suppl/doi:10.1093/molbev/msu305/-/DC1) are available at *Molecular Biology and Evolution* online (<http://www.mbe.oxfordjournals.org/>). 

###### Supplementary Data

The authors are grateful to Bastien Boussau, Laurent Duret, Vincent Daubin, Eric Tannier, and Nicolas Lartillot for fruitful comments and suggestions.  This work was supported by the French Agence Nationale de la Recherche (ANR) and is a contribution to the Ancestrome project (ANR-10-BINF-01-01).  G.J.S. was supported by the Marie Curie CIG 618438 "Genestory" and the Albert Szent-Györgyi Call-Home Researcher Scholarship A1-SZGYA-FOK-13-0005 supported by the European Union and the State of Hungary, cofinanced by the European Social Fund in the framework of TÁMOP 4.2.4.  A/1-11-1-2012-0001 "National Excellence Program." 

[^1]: ^‡^Present address: Department of Biological Engineering, Massachusetts Institute of Technology, Cambridge, MA

[^2]: ^§^Present address: Department of Biochemistry and Microbiology, University of Victoria, BC, Canada

[^3]: ^†^These authors equally contributed to this work. 

[^4]: **Associate editor:** Tal Pupko

[^5]: N[ote]{.smallcaps}.---Values obtained in this study for the ancestor of the Firmicutes (italic characters) were inferred using either the LeuB S-unaware tree or the LeuB S-aware reconciled tree, with the site-heterogeneous EX_EHO model.  Data for contemporary (first three lines) and other ancestral LeuBs for *Bacillus* (ANC1-4) characterized in [@msu305-B25] are shown for comparison.  Errors for $\Delta G_{N - U}^{\ddagger}$ are the standard error in the calculation of $\Delta G_{N - U}^{\ddagger}$ from ≥5 measurements of the unfolding rate in different urea concentrations.