Pervasive contingency and entrenchment in a billion years of Hsp90 evolution

The Historical Trajectory of Hsp90 Sequence Evolution.

We inferred the maximum-likelihood (ML) phylogeny of Hsp90 protein sequences from 261 species of Amorphea [the clade comprising Fungi, Metazoa, Amoebozoa, and related lineages (29)], rooted using green algae and plants as an outgroup (Fig. 1A, Fig. S1, and Datasets S1–S3). We reconstructed ancestral NTD sequences at all nodes along the trajectory from the common ancestor of Amorphea (ancAmoHsp90) to extant S. cerevisiae (ScHsp90) and identified substitutions as differences between the most probable reconstructions at successive nodes (Dataset S4).

Along this entire trajectory, substitutions occurred at 72 of the 221 sites in the NTD; because of multiple substitutions, 98 unique ancestral amino acid states existed at these sites at some point in the past and have since been replaced by the ScHsp90 state. The vast majority of these 98 ancestral states are reconstructed with high confidence (posterior probability > 0.95) in one or more ancestors along the trajectory (Fig. 1B), and every ancestral sequence has a mean posterior probability across sites of >0.95 (Fig. S2 A–C).

Entrenchment and Irreversibility.

To measure the fitness effects of ancestral amino acids when they are reintroduced into an extant Hsp90, we created a library of ScHsp90 NTD variants, each of which contains 1 of the 98 ancestral states. We determined the per-generation selection coefficient (s) of each mutation to an ancestral state relative to ScHsp90 via bulk competition monitored by deep sequencing (Dataset S5), a technique with highly reproducible results (Fig. S3). Our assay system reduces Hsp90 expression to ∼1% of the endogenous level (30), which magnifies the fitness consequences of Hsp90 mutations, enabling us to detect effects of small magnitude.

We found that the vast majority of reversions to ancestral states in ScHsp90 are deleterious (Fig. 1C). Using a mixture model to account for experimental noise in fitness measurements, we estimate that 92% of all reversions reduce the fitness of ScHsp90 (95% CI, 83−99%; Fig. 1C and Fig. S4). Three other statistical methods that differ in their assumptions yielded estimates that between 54% and 95% of reversions are deleterious (Fig. S5A). Two reversions cause very strong fitness defects (s = −0.38 and −0.54), but the typical reversion is only mildly deleterious (median s = −0.010, P = 1.2 × 10⁻¹⁶, Wilcoxon rank sum test; Fig. 1C). This conclusion is robust to excluding ancestral states that are reconstructed with any statistical ambiguity (P = 4.5 × 10⁻¹⁴). The magnitude of each mutation’s negative effect on fitness correlates with indicators of site-specific evolutionary, structural, and functional constraint, corroborating the view that they are authentically deleterious (Fig. S6).

These results do not imply that reversions can never happen. Twelve sites did undergo substitution and reversion at some point along the lineage from ancAmoHsp90 to ScHsp90. Rather, our observations indicate that, at the current moment in time, most ancestral states are selectively inaccessible, irrespective of whether they were available at some moment in the past or might become so in the future (31).

Intramolecular Versus Intermolecular Epistasis.

Reversions to ancestral states may be deleterious because the derived states were entrenched by subsequent substitutions within Hsp90 (intramolecular epistasis) (1, 4); alternatively, they may be incompatible with derived states at other loci in the S. cerevisiae genome (intermolecular epistasis), or the derived states may unconditionally increase fitness. Entrenchment because of intramolecular epistasis predicts that introducing into ScHsp90 sets of deleterious ancestral states that existed together at ancestral nodes should not reduce fitness as drastically as would be predicted from the individual mutations’ effects. To test this possibility, we reconstructed complete NTDs from two ancestral Hsp90s on the phylogenetic trajectory (Fig. 2A and Fig. S2) and assayed their relative fitness in S. cerevisiae as chimeras with ScHsp90’s other domains. This design provides a lower-bound estimate of the extent of intramolecular epistasis, because it does not eliminate interactions between substitutions in the NTD and those in other domains of Hsp90.

We found that intramolecular epistasis is the predominant cause of entrenchment. The first reconstruction, ancAscoHsp90, from the ancestor of Ascomycota fungi [estimated age, ∼450 My (32)], differs from ScHsp90 at 42 NTD sites. If the fitness effects of these ancestral states when combined were the same as when introduced individually, they would confer an expected fitness of 0.65 (95% CI, 0.61–0.69; Fig. 2B). When introduced together, however, the actual fitness is 0.99 (Fig. 2D and Fig. S7A), indicating that the current fitness deficit of ancestral states is caused primarily by deleterious epistatic interactions within the NTD.

The older ancestor, ancAmoHsp90 [estimated age, ∼1 billion years (33)], differs from ScHsp90 at 60 NTD sites. When combined, these differences would confer an expected fitness of 0.23 (95% CI, 0.21–0.26; Fig. 2B), but the actual fitness of the NTD is 0.43 (Figs. S2D and S7A), again indicating strong epistasis within the NTD. We hypothesized that the remaining fitness deficit caused by the ancAmoHsp90 NTD could be attributed to intramolecular epistasis between the NTD and substitutions in the other Hsp90 domains. We identified a candidate substitution in the protein’s middle domain that physically interacts with NTD residues to form the ATP-binding site (34); reverting this substitution to the ancAmorphea state (L378i) together with the ancAmorphea NTD increases fitness to 0.96 (Fig. 2D and Fig. S7A).

These findings indicate that virtually all of the context-dependent deleterious effects of ancestral states are caused by intramolecular interactions within the NTD and with one other site in the Hsp90 protein. Derived states that emerged along the Hsp90 trajectory have been entrenched by subsequent substitutions within the same protein, which closed the direct path back to the ancestral amino acid without causing major changes in function or fitness (1).

Contingency and Permissive Substitutions.

We next determined whether the derived states that evolved during the protein’s history were contingent on prior permissive substitutions. We constructed a library of variants of the ancAmoHsp90 NTD, the deepest ancestor of the trajectory, each of which contains 1 of the 98 forward mutations to a derived state. We cloned this library into yeast as a chimera with ScHsp90’s other domains (with site 378 in its ancestral state) and used our deep-sequencing–based bulk fitness assay to measure the selection coefficient of each mutation relative to ancAmoHsp90 (Fig. S3 D and E and Dataset S6).

We found that about half of mutations to derived states were selectively unfavorable (Fig. 3A). After accounting for experimental noise using a mixture model, we estimate that 48% of derived states reduce ancAmoHsp90 fitness (95% CI, 29–74%), and 32% are neutral (95% CI, 0–57%; Fig. S8). Three other statistical approaches gave similar results, with 48 to 66% of derived states inferred to be deleterious (Fig. S5A). Twenty percent of the derived states are beneficial in our assay (95% CI, 9–42%), which could be because they are unconditionally advantageous or because of epistatic interactions with other loci in S. cerevisiae or other regions of ScHsp90. Two derived states had very strong fitness defects, but the typical derived state is weakly deleterious (median s = −0.005, P = 5.8 × 10⁻⁴, Wilcoxon rank sum test; Fig. 3A).

As with the reversions to ancestral states, the effects of individual derived states, as measured in the ancestral background, predict fitness consequences far greater than observed when the derived states are combined in the Hsp90 genotypes that existed historically along the phylogeny (Fig. 2 C and D and Fig. S7B). Thus, many derived states would have been deleterious if they had occurred in the ancestral background, but they became accessible following subsequent permissive substitutions that occurred within Hsp90. Taken together, the data from the ancestral and derived libraries indicate that 77% of the amino acid states that occurred along this evolutionary trajectory were contingent on prior permissive substitutions, entrenched by subsequent restrictive substitutions, or both (Fig. 3B).

Specificity of Epistatic Interactions.

Epistatic effects on fitness can emerge from specific genetic interactions between substitutions that directly modify each other’s effect on some molecular property, or from nonspecific interactions between substitutions that are additive with respect to bulk molecular properties [e.g., stability (2, 35)] if those properties nonlinearly affect fitness (36–38).

To explore which type of epistasis predominates in the long-term evolution of Hsp90, we first investigated the two strongest cases of entrenchment, the strongly deleterious reversions V23f and E7a (uppercase letters indicate the ScHsp90 state, and lowercase, the ancestral state). We sought candidate restrictive substitutions for each of these large-effect reversions by examining patterns of phylogenetic co-occurrence. Substitution f23V occurred not only along the trajectory from ancAmoHsp90 to ScHsp90 but also in parallel on another fungal lineage; in both cases, candidate epistatic substitution i378L co-occurred on the same branch (Fig. S9 A and B). As predicted if i378L entrenched f23V, we found that introducing the ancestral state i378 in ScHsp90 relieves the deleterious effect of the ancestral state f23 (Fig. 4A). These two residues directly interact in the protein’s tertiary structure to position a key residue in the ATPase active site (Fig. S9 C and D), explaining their specific epistatic interaction.

In the case of E7a—the other reversion strongly deleterious in ScHsp90—the ancestral state was reacquired in a closely related fungal lineage. We reasoned that the substitutions that entrenched a7E on the lineage leading to ScHsp90 must have themselves reverted or been further modified on the fungal branch in which reversal E7a occurred. We identified two candidates (n13T and a151N) that met these criteria (Fig. S10 A–C). As predicted, experimentally introducing the ancestral states n13 or a151 into ScHsp90 relieves much of the fitness defect caused by the ancestral state a7, indicating that substitutions n13T and a151N entrenched a7E (Fig. 4B). These three sites are on interacting secondary structural elements that are conformationally rearranged when Hsp90 converts between ADP- and ATP-bound states (Fig. S10 D and E).

To test whether these states specifically restrict particular substitutions or are general epistatic modifiers, we asked whether the restrictive substitutions that entrenched one substitution also entrenched the other (2). As predicted if the interactions among these sets of substitutions are specific, introducing L378i does not ameliorate the fitness defect caused by E7a, and introducing T13n or N151a does not ameliorate the fitness defect caused by V23f (Fig. 4 C and D). These data indicate that specific biochemical mechanisms underlie these large-effect restrictive interactions.

Finally, we investigated whether the epistatic interactions among the set of small-effect substitutions in this trajectory are also specific or the nonspecific result of a threshold-like relationship between fitness and some bulk property such as stability (2, 35). If epistasis is mediated by a nonspecific threshold relationship, mutations that decrease fitness in one background will never be beneficial in another, although they can be neutral if buffered by the threshold (Fig. 5A) (2, 20, 25). In contrast, specific epistatic interactions can switch the sign of a mutation’s selection coefficient in different sequence contexts (Fig. 5B) (37). As predicted for specific epistasis, we found that for most differences between ancAmoHsp90 and ScHsp90 (65%), the ancestral state confers increased fitness relative to the derived state in the ancestral background but decreases it in the extant background (Fig. 5C). The selection coefficients of mutations are negatively correlated between backgrounds (P = 0.009), indicating that the substitutions that became most entrenched in the present also required the strongest permissive effect in the past. This pattern is expected if the structural constraints that determine the selective cost of having a suboptimal state at some site are conserved over time, but the specific state that is preferred depends on the residues present at other sites.

Taken together, these findings indicate that most epistasis during the long-term evolution of Hsp90 involved specific one-to-one (or few-to-few) interactions among sites, not general effects on the protein’s tolerance to mutation.

Relation to Prior Work.

We observed widespread and specific epistasis over the course of a billion years of Hsp90 evolution, during which the protein’s function, physical architecture, and fitness were conserved. The fraction of historical substitutions that were contingent on permissive substitutions, entrenched by restrictive substitutions, or both—about 80%—is considerably higher than suggested by previous experimental work (2, 22, 24) and some computational analyses (14), rivaling the highest estimates from computational studies (1, 16, 17). One explanation for the more widespread epistatic interactions in our study may be our method’s capacity to detect much smaller growth deficits than have been discernable in previous experimental studies.

Another difference from previous research is that we primarily observed specific epistasis, whereas several studies have found a dominant role for nonspecific stability-mediated epistasis, particularly during the short-term evolution of viruses (2, 11, 20). This disparity could be attributable to a difference in selective regime or in timescale: the epistatic constraints caused by specific interactions are expected to be maintained over far longer periods of time than those caused by nonspecific interactions, which are easily replaced by other substitutions because of the many-to-many relationship between permissive and permitted amino acid states (2, 20, 37). The prevalence and type of epistasis may also vary because of differences in proteins’ physical architectures. Additional case studies will be necessary to evaluate the causal role of these and other factors in determining the nature of epistatic interactions during evolution.

Limitations.

Our strategy has some known limitations, but none is likely to change our major conclusions. For example, we assayed the effect of long-past substitutions in the context of extant yeast cells. But our experiments indicate that there is only very weak epistasis for fitness between historical substitutions within Hsp90 and those at other loci, because the reconstructed ancestral Hsp90 chimeras cause a fitness deficit of only 0.01–0.04 when introduced into S. cerevisiae cells—much smaller than the sum of intramolecular incompatibilities revealed by introducing the ancestral states individually. Our finding of widespread contingency and entrenchment is therefore not an artifact of incompatibilities between ancestral Hsp90 states and the genotype of present-day S. cerevisiae at other loci.

A second potential limitation is that the ancestral states we tested were reconstructed phylogenetically, not known empirically. But the vast majority of states were inferred with high statistical confidence, because the Hsp90 NTD is well conserved and we used a densely sampled alignment. The ambiguity that was present primarily concerned the specific ancestral node at which an inferred ancestral state was present, not whether or not it was ancestral somewhere along the trajectory, which is the key inference for our purposes. Even when all states with any degree of statistical uncertainty in the ancestral reconstruction were excluded from the analysis, the remaining data strongly supported our conclusions concerning contingency and entrenchment.

Finally, we measured fitness under a particular set of experimental conditions. Our assay system reduces Hsp90 expression to ∼1% of the endogenous level (30). Based on previous work quantifying the relationship between Hsp90 function, expression, and growth rate (30), we estimate that the average selection coefficient of −0.01 we observed among contingent or entrenched substitutions corresponds to a fitness deficit of approximately s = −5 × 10⁻⁶ under native-like expression levels. Mutations with selection coefficients in this range would likely be subject to purifying selection in large microbial populations (39–41). Our assay also tests fitness under log-phase growth conditions in rich media. A more heterogeneous or demanding environment would likely increase the magnitude of selective effects of Hsp90 mutations, because stress should amplify the fitness consequences of mutations in the proteostasis machinery.

Implications.

Our observation that contingency and entrenchment affected the majority of historical substitutions suggests a daisy-chain model by which genetic interactions structured long-term Hsp90 evolution (Fig. 5D). A permissive mutation becomes entrenched and irreversible once a substitution contingent upon it occurs; if the contingent substitution subsequently permits a third substitution, it too becomes entrenched (1, 16).

Most of the substitutions along the trajectory from ancAmoHsp90 to ScHsp90 were both contingent and entrenched, suggesting that they occupy an internal position in this daisy chain. Each of these changes closed reverse paths at some sites and opened forward paths at others, which, if taken, would then entrench the previous step. Evolving this way over long periods of time, proteins come to appear exquisitely well-adapted to the conditions of their existence, with most present states superior to past ones. But the conditions that make today’s states so fit include—or are even dominated by—the transient internal organization of the protein itself.

Phylogenetic Inference and Ancestral Reconstruction.

We inferred the ML phylogeny for 261 Hsp90 protein sequences from the Amorphea clade, with Viridiplantae as an outgroup, under the LG+Γ+F model in RAxML, version 8.1.17 (42). Most probable ancestral NTD sequences were reconstructed on this ML phylogeny using the AAML module of PAML, version 4.4 (43). The trajectory of sequence change was enumerated from the amino acid sequence differences between successive ancestral nodes on the lineage from the common ancestor of Amorphea (ancAmoHsp90) to S. cerevisiae Hsp82 (ScHsp90). Ancestral states are defined as amino acid states not present in ScHsp90 that occurred in at least one ancestral node on the lineage from ancAmoHsp90 to ScHsp90. Derived states are defined as amino acid states not present in the reconstructed ancAmoHsp90 sequence that occurred in at least one descendent node on the lineage to ScHsp90.

Bulk Growth Competitions.

The ScHsp90 and ancAmoHsp90 NTD (+L378i; SI Methods and Fig. S2D) protein-coding sequences were expressed together with the other domains from ScHsp90 from the p414ADHΔTer plasmid (30). Individual mutations in each variant library were introduced via PCR. Library genotypes were tagged with short barcodes to simplify sequencing steps during the bulk competition (44); barcodes were associated with variant genotypes via paired-end sequencing on an Illumina MiSeq instrument. Variant libraries were transformed into the DBY288 Hsp90 S. cerevisiae shutoff strain (45, 46) and grown for 48 h (ScHsp90 library) or 31 h (ancAmoHsp90 library) under selective conditions. Cultures were maintained in log phase by regular dilution with fresh media, maintaining a population size of 10⁹ or greater, and samples of ∼10⁸ cells were collected at regular time points over the course of the bulk competition. Plasmid DNA was isolated from each time point (47), and the frequency of each library genotype at each time point was determined via Illumina sequencing. Bulk competitions were performed in duplicate.

Determination of Selection Coefficients.

The ratio of the frequency of each variant in the library relative to wild type (ancAmoHsp90 or ScHsp90) was determined from the number of sequence reads at each time point, and the slope of the logarithm of this ratio versus time (in number of generations) was determined as the raw per-generation selection coefficient (s) (48):

where n_m and n_wt are the number of sequence reads of mutant and wild type, respectively, and time is measured in number of wild-type generations. For genotypes that were not assayed in the bulk competition, selection coefficients were determined from monoculture growth rates (48), measured over 30 h of growth with periodic dilution to maintain log-phase growth (30). Selection coefficients for both types of measurement were scaled in relative fitness space (w = e^s) such that an Hsp90 null allele, which is lethal, has a relative fitness of 0 (s = −∞).

Estimating the Fraction of Deleterious Mutations.

We estimated the fraction of mutations in each library that are deleterious by fitting the distribution of mutant selection coefficients to a mixture model of underlying Gaussian distributions. One of the underlying mixture components was required to have the mean and SD estimated from replicate measurements of wild-type sequences present in the library but represented by independent barcodes; remaining mixture components had a freely fit mean and SD, and all components had a freely fit mixture proportion. The optimal number of mixture components was determined via Akaike information criterion. The mixture component derived from the wild-type sampling distribution was taken to represent neutral mutations in the library, and the other mixture components were taken to reflect nonneutral mutations. We estimated the fraction of neutral mutations as the mixture proportion of the neutral mixture component. We then estimated the proportion of deleterious (or beneficial) mutations as the sum of the weighted cumulative density of all nonneutral components that was below (or above) zero.

To determine the robustness of the mixture model’s estimates of the fraction of deleterious (or beneficial) mutations, we used three other approaches to analyze the distribution of fitness measurements. The simplest, a nonparametric approach, estimates the proportion of deleterious (or beneficial) mutations as the proportion of mutations with observed selection coefficients less (or greater) than 0. The second approach is an empirical Bayes approach that calculates the posterior probability that each mutation is deleterious (or beneficial), given the experimentally observed selection coefficients and measurement error for wild-type and mutant genotypes; these posterior probabilities are summed to yield the estimated proportion of mutations in each fitness category. Third, we constructed 95% CIs for each mutation’s selection coefficient given its mean over two replicates and SE estimated over all mutations and counted the number of mutations with selection coefficients less (or greater) than zero whose CIs do not overlap zero.

Data and Code Availability.

Processed sequencing data and scripts to reproduce all analyses are available at https://github.com/JoeThorntonLab/Hsp90_contingency-entrenchment. Tables listing mutants and their selection coefficients are included as Datasets S5 and S6. Raw sequencing data from each bulk competition have been deposited in the National Center for Biotechnology Information Sequence Read Archive under accession number SRP126524. Tables linking barcode variants to their associated Hsp90 genotype are included as Datasets S7 and S8.