A protein’s genetic architecture – the set of causal rules by which its sequence determines its specific functions – also determines the functional impacts of mutations and the protein’s evolutionary potential. Prior research has proposed that proteins’ genetic architecture is very complex, with pervasive epistatic interactions that constrain evolution and make function difficult to predict from sequence. Most of this work has considered only the amino acid states present in two sequences of interest and the direct paths between them, but real proteins evolve in a multidimensional space of 20 possible amino acids per site. Moreover, almost all prior work has assayed the effect of sequence variation on a single protein function, leaving unaddressed the genetic architecture of functional specificity and its impacts on the evolution of new functions. Here we develop a new logistic regression-based method to directly characterize the global causal rules of the genetic architecture of multiple protein functions from 20-state combinatorial deep mutational scanning (DMS) experiments. We apply it to dissect the genetic architecture and evolution of a transcription factor’s specificity for DNA, using data from a combinatorial DMS of an ancient steroid hormone receptor’s capacity to activate transcription from two biologically relevant DNA elements. We show that the genetic architecture of DNA recognition and specificity consists of a dense set of main and pairwise effects that involve virtually every possible amino acid state in the protein-DNA interface, but higher-order epistasis plays only a tiny role. Pairwise interactions enlarge the set of functional sequences and are the primary determinants of specificity for different DNA elements. Epistasis also massively expands the number of opportunities for single-residue mutations to switch specificity from one DNA target to another. By bringing variants with different functions close together in sequence space, pairwise epistasis therefore facilitates rather than constrains the evolution of new functions.
This study includes important findings on protein evolution, namely that changes in function are largely attributable to pairwise rather than higher-order interactions, and that epistasis potentiates rather than constrains evolutionary paths. Compelling evidence supporting the conclusions is provided by applying a new model to a previously generated experimental dataset on deep mutational scanning of the DNA-binding domain (DBD) of steroid hormone receptor. The implications of this work are of considerable interest to protein biochemistry, evolutionary biology, and numerous other fields.
A protein’s genetic architecture – the set of causal rules by which its sequence determines its specific functions – is of fundamental importance in genetics, biochemistry, and evolution. The distribution of functions across the space of all possible protein sequences – and, in turn, the paths accessible to an evolving protein under the influence of mutation, drift, and selection – all arise from its genetic architecture. A microscopic dissection of a protein’s genetic architecture would specify how every possible amino acid at each site would affect a protein’s functions, as well as how epistatic interactions among amino acid pairs, triplets, and higher-order combinations change those effects. This kind of knowledge would allow us to develop a macroscopic overview of the most important genetic causes that generate the protein’s sequence-function relationships and its potential evolutionary trajectories.
A topic of particular interest is the extent and impact of epistatic interactions. If the effects of an amino acid state depend strongly on the states present at other sites in the protein – and especially on higher-order combinations – then the functions of a protein could become difficult to predict from its sequence or those with similar sequences. Epistasis could also constrain evolution by creating a rugged functional landscape in sequence space, in which many optimal or nearoptimal genotypes are inaccessible from others under most evolutionary scenarios(1–22). And if epistatic interactions tend to be negatively biased, making sequences less functional than anticipated from the main effects of their amino acids (23–31), then functional variants might be restricted to only those regions of sequence space where each site has the state or states with the most favorable main effects.
It is clear that epistasis is pervasive in real proteins, but its effect on evolution and the extent of higher-order interactions remains murky. Studies of single-site mutant libraries (32–42) and case studies of historical substitutions (43–49) have shown that mutations’ effects often depend on the sequence background in which they occur. Studies that measure the effects of many pairs of mutations on a single function have shown that pairwise interactions reduce the accessibility of some two-step paths from a designated wildtype protein (50–53). This work, however, has not examined higher-order epistasis or the much larger set of paths through the space of all functional proteins, thus leaving the overall complexity of genetic architecture and its evolutionary impacts largely uncharacterized.
Addressing these issues requires comprehensive analysis of higher-order combinations, which is experimentally intractable if the number of states and sites is large. Most studies have addressed this challenge by reducing the number of states and sites considered to those defined by the sequence differences between two proteins of interest (44, 46, 49, 54– 80). These binary-state studies have found that epistasis reduces the number of functional intermediates on the paths between a designated pair of starting and ending proteins (81–87). By focusing on just a few states and proteins, however, these approaches have addressed only a tiny portion of the huge ensemble of potential genotypes and the full set of evolutionary paths among them. Moreover, many protein families contain members with distinct specific functions, like the capacity to bind different substrates or other ligands, but virtually all DMS studies to date have examined only a single protein function. As a result, little is known about the genetic architecture of functional specificity and the effect of epistasis on the evolution of new functions.
A complete assessment of the genetic architecture of functional specificity requires three things: 1) an experiment that assays all combinations of 20 amino acid states at a designated set of sites for multiple functions; 2) a method to extract from these data the causal rules by which sequence determines function, and 3) a way to characterize the effect of this architecture on evolutionary trajectories through sequence space. Several studies have performed complete combinatorial mutational scans, typically for a single function and addressing three or four sites of a priori structural or functional interest (88–95). This work has established that the set of functional sequences is large and diverse, and that the direct paths between particular functional genotypes are sometimes blocked because of epistatic interactions. These studies have not, however, sought to extract from their data a characterization of the protein’s genetic architecture, how it defines which protein variants are functional in the first place, or how it affects evolutionary accessibility across sequence space as a whole.
A critical limiting factor has been the lack of tools to dissect the genetic architecture of multiple functions from 20-state combinatorial datasets. Prior work on two-state libraries has shown that modeling the rules of genetic architecture relative to a particular wild-type reference state yields globally inaccurate estimates of the effects of mutations and combinations when introduced into other backgrounds (75). Here we develop and implement a reference-free method for global dissection of 20-state sequence-function relationships. The model’s terms are encoded relative to the global functional average across all genotype rather than a particular reference sequence, and they directly express the portion of all functional variation in a protein that is attributable to any particular set of genetic determinants. This way of encoding the model also allows us to directly assess the effect of epistasis on the distribution of functions across sequence space and their accessibility under various evolutionary scenarios.
We apply this approach to a complete combinatorial DMS dataset at four sites that are critical for DNA recognition in the DNA-binding domain (DBD) of the reconstructed ancestral steroid hormone receptor (56, 93, 96). This experiment assayed the transcription factor’s ability to activate transcription from two different biologically relevant DNA response elements, enabling us to characterize the rules of genetic architecture that define functional specificity. By applying our approach to a reconstructed ancestral protein, we assessed how the protein’s historical genetic architecture shaped the maintenance of the ancestral function and the accessibility of the derived function that actually evolved during steroid hormone receptor history.
Steroid hormone receptors are a family of ligand activated transcription factors that are involved in vertebrate development, behavior, and reproduction (97). There are two major subfamilies of steroid hormone receptors, which recognize distinct palindromic DNA response elements (REs). Estrogen receptors (ERs) bind to a palindrome of the ER response element (ERE, AGGTCA), whereas the ketosteroid receptors (kSRs, including the receptors for androgens, progestogens, glucocorticoids, and mineralocorticoids) strongly prefer the steroid receptor response element (SRE, AGAACA) (98–100). Earlier experiments showed that the reconstructed DBD of the last common ancestor of the two subfamilies (called AncSR1) was ERE-specific, whereas the last common ancestor of the kSRs (AncSR2) was SRE-specific. During the interval between AncSR1 and AncSR2, three substitutions occurred in the protein’s recognition helix (RH), which inserts into the DNA major groove and makes contact with the bases that vary between ERE and SRE. These substitutions were shown experimentally to cause the shift in specificity (96).
The data that we analyze here come from a previously published deep mutational scan of AncSR1 (93). The library contained all 160,000 combinations of 20 amino acids at four sites in the RH – the three historically substituted sites plus another site that is variable in closely related nuclear receptors. The library was transformed into yeast strains carrying either an ERE-or SRE-driven GFP reporter and functionally characterized using a FACS-based Sort-seq assay.
To reduce the effect of measurement noise on the characterization of genetic architecture in this project, we transformed the continuous functional data into categorical form: each variant was categorized on each RE as a null, weak, or strong activator relative to negative control (stop-codon-containing variants) and historically-relevant positive control reference sequences (the RH sequence of AncSR1 and extant ERs on ERE, or that of AncSR2 and extant kSRs on SRE). On each RE, variants were classified as null if their mean fluorescence was indistinguishable from the negative controls, weak if they produced fluorescence between null and the historical reference state, and strong if their fluorescence was as great or greater than the reference. Across both REs, 1342 variants were classified as strong and 3166 as weak activators. Categorization substantially reduced the effect of technical noise in the assay: concordance in the activation class assigned to each variant between replicates was >97%, much better than the between-replicate correlation of mean fluorescence as a continuous variable (R2 = 0.62 for functional variants, and R2 = 0.11 when null variants are included).
Ordinal linear regression of genetic architecture
To dissect the genetic architecture of RE binding and specificity by AncSR-DBD, we developed an ordinal regression model and fit it to the experimental data. The model contains terms for the main effect of every possible amino acid at each variable site in the protein and for the epistatic effect of every pair and triplet of sites (Fig. 1A). The genetic score of each variant is defined as the sum of the main and epistatic effects of the sequence states and combinations it contains (Fig. 1B). The variant’s genetic score determines the probability that a variant is a null, weak, or strong activator through an ordinal logistic function: an increase in the genetic score causes a corresponding linear increase in the log of the odds that a variant is in a higher vs. lower activation class. The only other free parameters in the model are the average genetic score across all variants and the threshold scores that represent the boundaries between activation classes.
To incorporate functional specificity, the model contains two terms for every amino acid state or combination: one for its effect on nonspecific DNA binding (the contribution to genetic score averaged over ERE and SRE) and another that reflects its impact on specificity (the difference between its score on each RE and the average over REs). The model also contains a single term for the global effect of SRE vs. ERE, averaged across all protein variants, which represents the main effect of the RE itself. The total genetic score for any protein:RE complex is the sum of all the main and epistatic effects of its amino acids on nonspecific DNA binding, plus all the main and epistatic effects on RE specificity, plus the global effect of the RE. (The sign of the specificity and RE terms in the summation is determined by whether the complex contains SRE or ERE.)
The model is encoded to give accurate and efficient estimates of the rules of genetic architecture across the entire set of possible sequences. Rather than defining the effects of a state or combination relative to a designated wild-type sequence, we use a reference-free framework in which effects are defined relative to the global average across all variants (Fig. 1C-D; for a related approach using DNA sequences, see (101)). The main effect of any amino acid state is the mean genetic score of all variants containing that state minus the global average score. The epistatic effect of a pair of states is the mean genetic score of all variants containing that pair minus the sum of the global average and the main effects of the two states. And the third-order effect of a triplet of states is the average score of all variants containing that triplet minus the global average and the sum of the pairwise and main effects of the component states and pairs. We did not model fourth-order epistatic effects, because these cannot be distinguished from technical error in the assay measurements.
This approach has several advantages. First, the sigmoidal shape of the logistic function accounts for global non-linearity in the genotype-phenotype relationship; this is important, because failing to do so leads to spurious findings of specific epistatic interactions (5, 20, 102–106). Second, this approach reduces the effect of estimation bias and measurement error, because each effect is averaged across a large number of observations from unique genotypes at other sites (e.g. 8000 for each main effect term, 400 for each pairwise effect)(75). Finally, the model’s form allows us to use an ANOVA-like framework to directly partition the total specific genetic variance – the variance in genetic score across all protein variant-RE complexes – into the variance attributable to every causal factor in the model. The specific genetic variance attributable to any effect term in the model is the squared coefficient, weighted by the fraction of all genotypes to which it applies. The specific genetic variance explained by any set of terms – such as all the main-effect terms, or all the thirdorder epistatic terms – is simply the sum of the specific genetic variance accounted for by all the terms in that set (see Methods and Appendix 1).
We fit the model to the DMS data by least-squares ordinal logistic regression (Fig. 1E). To reduce overfitting and maximize predictive accuracy, we used L2 regression combined with cross-validation to identify the optimal regularization parameters. The model fits the experimental data well. When we used the best-fit model to predict the highest-probability activation class for each protein variant on each RE, the predicted class matched the observed class for 99.5% of all variants. The deviance of the model – an analog of the coefficient of determination used in linear regression – was 0.86 (Figure 1-Figure Supplement 1, Figure 1-Figure Supplement 2, Figure 1-Figure Supplement 3). The remaining variation in function not explained by the model is attributable to technical error in measurement and/or fourth-order epistasis. The genetic score of activators is well correlated with the apparent ΔG of dissociation from DNA, which was previously measured for a subset of variants by fluorescence anisotropy (R2 = 0.75, Figure 1-Figure Supplement 4; this correlation is slightly better than the correlation between the mean fluorescence values used to classify variants and ΔG, R2 = 0.63).
Low-order genetic architecture
We found that main effects and pairwise epistasis are important in RE recognition, but higher order epistasis is not. Of all the specific genetic variance explained by the best-fit model, more than half is attributable to the main effects of single amino acids, and virtually all of the remainder is attributable to pairwise interactions or the global effect of the RE. Main effects are the primary determinants of nonspecific DNA binding, but pairwise interactions explain the majority of RE-specificity. Third-order epistasis accounts for only 2% of total variance, with similarly small contributions to both nonspecific binding and specificity (Fig. 2A-C). About half (45%) of all the specific variance explained by epistatic amino acid terms involves sign epistasis – interactions that on average change the direction of a state’s effects in the presence of another amino acid – while the other half changes the magnitude but not the direction of effect (Fig. 2D). Although high-order interactions are unimportant, the model is very dense at low orders. Of the nearly 69,000 possible terms in the complete model, it takes >7,000 to account for 99% of specific genetic variance (“the 99% set”, Fig. 2E), including the main effect of every possible amino acid state at all four sites and >80% of all possible second-order epistatic terms. Every amino acid state participates in a minimum of 60 pairwise interactions with states at other sites (Figure 2-Figure Supplement 1). By contrast, the 99% set includes only 5% of all possible thirdorder terms.
The density of the model at low orders arises because thousands of variants activate from ERE and/or SRE, and their sequences are very diverse (Fig. 2F-G). The effects of most amino acids and pairs are therefore of small to moderate magnitude, changing the genetic score by a median of 3% and 2% of the threshold required for strong activation, respectively. Third-order effects were generally tiny, with a median absolute magnitude of just 0.02% (Fig. 2H-I, Figure 2-Figure Supplement 2, Figure 2-Figure Supplement 3). As a consequence, no single amino acid state or pair is sufficient to make a sequence into a strong activator; rather, numerous favorable states and combinations of small to moderate size at three or all four of the sites are required to generate a functional protein.
The sequence-function relationship of DNA recognition is therefore complex but not idiosyncratic. The functions of any RH sequence can be predicted with good accuracy from its main and pairwise effects alone (Figure 2-Figure Supplement 4). Higher-order epistatic terms are largely dispensable, so prediction does not require precise knowledge of the effects of the vast number of triplets and quartets (Figure 2-Figure Supplement 5). Although low-order prediction is largely sufficient, it cannot be reduced to a few simple rules, such as “must have an arginine at site 3” to bind either RE, or “must have a valine or methionine at site 4” to prefer SRE.
Rather, there are hundreds of potential “and” and “or” rules by which RH variants across the massive multidimensional space of the protein’s sequence can become a specific or nonspecific activator.
The functional effects of variation at each site can to some extent be structurally rationalized (Fig. 2J). Sites 2 and 3 have their largest effects through nonspecific binding and their side chains are positioned near the parts of the DNA that are shared between ERE and SRE; the most favorable effects come from having small-volume amino acids at site 2, particularly in combination with a basic residue at site 3, which can potentially hydrogen-bond with a guanine common to both REs. By contrast, the states at sites 1 and 4 have their largest effects on specificity, and these residues are closest to the bases in the major groove that differ between the REs; small hydrophobic residues at site 4 favor SRE binding, where they can pack against the hydrophobic methyl groups that are unique to SRE but leave ERE’s hydrogen bond (donors/acceptors) unpaired.
Epistasis generates functional activators
To determine if epistatic effects have a directional bias on functional sequences, we calculated the net contribution of all epistatic terms to the genetic score of each active variant on each RE (Fig. 3A-B). Contrary to the expectation of diminishing returns (23, 29–31), we found that epistatic terms overwhelmingly enhance function. Epistatic terms make a net positive contribution to the genetic score of 100% of strong activators. Epistatic terms contribute 64% of the threshold score for strong activation relative to the global average, and the epistatic contribution was necessary to surpass this threshold in every single case. Weak activators also benefit from epistatic interactions, which contribute positively to 98% of variants in this class, but their effect is smaller – only 39% of the strong threshold on average – explaining in large part why they are weak rather than strong. Both pairwise and thirdorder effects contribute, but second-order interactions make a far larger positive contribution than third-order effects do (Figure 3-Figure Supplement 1).
To directly test the effect of epistasis on the number, identity, and function of genotypes, we formulated alternative models of genetic architecture that exclude epistasis at various orders; we then fit these models to the experimental data and compared the sets of variants predicted to be functional under each model. Confirming the key role of epistasis in generating functional proteins, the number of strong activators in the complete third-order model is 27% greater than in a first-order only model that allows no epistasis, and weak activators grow to a similar degree (Fig. 3C). Most of this increase is caused by pairwise epistasis, because a second-order model – which contains pairwise but not third-order epistasis – contains virtually the same number of activators as the complete third-order model (98%, Figure 3-Figure Supplement 2). Epistasis also changes which sequences are functional; of variants that are strong activators in the absence of epistasis, only 64% remain strong activators when epistasis is included.
Why does epistasis increase the number of functional variants and change their identity? Even in the first-orderonly model, main effects are mostly of small magnitude (Figure 3-Figure Supplement 3): the only way to achieve a genetic score sufficient for activation is to combine the few amino acids that have the largest main effects, and the number of ways to do this is limited. With epistasis, however, positive epistatic terms can supplement main effects, improving the function of sequences that contain states that have weak main effects. Corroborating this view, states with positive main effects on binding overwhelmingly have positive pairwise and third-order epistatic interactions, enhancing function more than expected if they were independent (Fig. 3D, Figure 3-Figure Supplement 4).
The genetic architecture of mutation effects
The genetic architecture of a protein directly describes how the sequence of each protein determines its function, but evolution involves mutations – changes that turn one sequence into another. Exchanging one amino acid for another in a particular protein alters the main effect at that site and all the epistatic interactions involving that site. We used our model to characterize the impacts on the genetic score for ERE and SRE binding of all 1.2 million possible context-specific mutations, defined as one of the 1,520 mutation types (a single-residue exchange from one of the 20 possible starting residues to 19 ending states at each of four sites) introduced into any of the 8,000 possible combinations of states at the other three sites.
A huge number of mutations can transform a nonfunctional sequence into a functional one, and epistasis is critical to these transitions. There are >50,000 context-specific mutations that can directly transform a null protein into a strong activator, >17,000 that can make a weak activator strong, and >150,000 more that can turn a null activator into a weak activator (Fig. 4A). The underlying mutations are diverse, with 1160 of the 1520 unique types of amino acid exchanges moving one or more variants into a higher activation class (Fig. 4B). In nearly every case (99%), the net epistatic impact of the mutation’s effects contribute to the increase in genetic score; in 94% of cases, the epistatic contribution is necessary for the score to surpass the threshold for the next activation class, and in almost half they are sufficient (Fig. 4C-D). Sign epistasis was rampant among mutations: every single one of the 1,520 mutation types caused positive effects in some backgrounds and negative effects in others, and more than 90% of mutation types (1380 of 1520) had both positive and negative effects on more than 10% of genetic backgrounds (Fig. 4E).
Epistasis also expands the mutational opportunities to change DNA specificity. Overall, there are hundreds of context-specific mutations and mutation types that change a strong ERE-specific activator directly into a strong SRE-specific activators; these involve changes at all four sites and use dozens of different amino acid states to generate both ERE and SRE specificity (Fig. 4F). In every single case, the net change in epistatic effects caused by the mutation is necessary for the switch in specificity (the loss of ERE binding and/or the gain of SRE binding), and in almost half of cases the epistatic effects are also sufficient (Fig. 4G-H). Both pairwise and third-order epistasis contributed: in about half of cases, third-order epistasis is necessary, but it is rarely sufficient (Figure 4-Figure Supplement 1).
To directly test the role of epistasis in creating mutational opportunities to switch DNA specificity, we compared the mutations available under the best-fit third-order model to those under the best-fit first-order model. Without epistasis, the total number of specificity-switching mutations is reduced by >80%, and the number of mutation types declines from 85 to just three (Fig. 4F). Under a model that incorporates pairwise but not third-order epistasis, the number of RE-switching mutations is intermediate, but closer to that in the third-order model (Figure 4-Figure Supplement 2). Thus, epistasis – mainly second order – is the primary factor in the opportunity for mutations to change RE-specificity with a single amino acid change.
Evolutionary paths in sequence space
We next assessed trajectories of functional evolution in sequence space, and the effect of epistasis on those trajectories. The RH sequence space consists of all 160,000 possible sequences, each representing one protein variant, with edges connecting nodes that can be directly interconverted by a single nucleotide change given the standard genetic code. We assigned to each node the function(s) predicted by its genetic score under each model. Nonfunctional nodes – those that are null or weak on both REs – were excluded from the network, because purifying selection will remove these genotypes from an evolving population under most circumstances 107. We focused our analysis on evolutionary paths from ERE-specific starting points to SRE-specific endpoints, because this is the functional transition that occurred during history.
Sign epistasis is generally thought to constrain epistasis by creating a fragmented sequence space in which local optima are separated by impassable valleys. While sign epistasis is rampant in our datasets (Fig. 2D, Fig. 4E), we found that given the best-fit third-order model, 99.8% of functional sequences are connected to each other in a single mutually accessible network (Fig. 5A).
Epistasis also enlarges the network of connected functional sequences, consistent with the function-enhancing effect of epistatic terms on functional variants. The thirdorder model contains 1603 functional nodes (1600 of which are connected in a single network), whereas the first-order epistasis-free model contains only 1080 functional nodes. The second-order model is almost as large as the complete third-order model, indicating that pairwise interactions account for most of the network-enlarging effect of epistasis (Figure 5-Figure Supplement 1). In multidimensional sequence space, epistasis therefore does not isolate functional variants: rather, it expands the network of functional sequences, virtually any of which can reach any other through a series of single-nucleotide substitutions.
The probability of reaching any variant from a given starting point depends on the number of substitutions required to reach it. We characterized the evolutionary accessibility of SRE specificity from ERE-specific starting points by simulating evolution from each ERE-specific genotype and measuring the length of the paths required for SRE specificity to evolve. We used two evolutionary scenarios for these simulations. The first represents purifying selection and neutral drift: from any ERE-specific node, every step to any functional neighbor has equal probability, irrespective of which REs it activates, and the path ends as soon as an SRE-specific node is reached (Fig. 5B). We found that in this scenario epistasis shortens the average path length to SRE-specificity for 80% of ERE-specific nodes, with an average reduction of 25% (Fig. 5C). Accessibility under the second-order model was very similar to the complete model, indicating that the increased accessibility of SRE-specific nodes was primarily caused by pairwise epistasis (Figure 5-Figure Supplement 2).
The second evolutionary scenario incorporates positive selection for SRE-specificity by making the probability of visiting each neighbor determined by a selection coefficient that is proportional to the node’s preference for SRE over ERE. Compared to the complete third-order model, eliminating epistasis in the first-order model again lengthened the average path taken to SRE-specificity, although this difference shrank as selection became stronger (Fig. 5D). This difference, too, was attributable primarily to pairwise epistasis (Figure 5-Figure Supplement 3).
Epistasis therefore shortens rather than lengthens the paths taken during functional evolution. To understand why this is so, we examined the distribution of functions across sequence space and its effect on the paths that connect ERE- and SRE-specific sequences. We found several contributing factors. First, the minimum distance from ERE-specific starting points to the nearest SRE-specific variant is on average 10% shorter in the third-order model than when epistasis is eliminated in the first-order model (Fig. 5E). Second, there are far more direct single-step mutations from ERE-to SRE-specific nodes: the average number of SRE-specific neighbors per ERE-specific node is 2.5 times higher in the third-order model, and each SRE-specific variant is also more likely to have at least one ERE-specific single-step neighbor (Fig. 5F). Third, epistasis reduces the average number of ERE-specific neighbors around ERE-specific nodes, further increasing the probability that a step to SRE specificity will be taken if one is available (Fig. 5F). Finally, epistasis increases the number of paths by which SRE specificity can evolve: in the third-order model, each ERE-specific node can reach its closest SRE-specific nodes by an average of 65 different paths; excluding epistasis shrinks this number by about 20%, and these paths are less distinct (Fig. 5G-H, Figure 5- Figure Supplement 4).
Taken together, these data indicate that epistasis facilitates the evolution of new protein functions, rather than constraining it. Epistasis does constrain evolution on a fine scale: by making some of the paths between designated pairs of sequences inaccessible, epistasis can and does make the evolution of particular outcomes contingent on intermediate steps. But epistasis also shapes which pairs exist in the network of functional sequences in the first place and how many steps apart those sequences are. When the genetic architecture includes epistasis, there are more functional sequences; further, similar sequences are more likely to have different functions, because a single amino acid difference can combine with the particular residues at the other sites, sometimes dramatically improving one function and impairing the other. By contrast, if the genetic architecture consists entirely of main effects, variants that share functions will tend to have sequences that are more similar to each other and more different from those with different functions: in this case, each function will be distributed more smoothly across sequence space, and neighboring sequences are less likely to have distinct functions than when epistasis is present.
Genetic architecture of a historical change in function
Finally, we sought to understand how the DBD’s genetic architecture shaped the historical transition from ERE-to SRE-specificity that occurred during the phylogenetic interval between AncSR1 and AncSR2. The three amino acid changes from the ancestral RH sequence egKa to the derived GSKV (ancestral in lower case, derived in upper case) can each be conferred by a single nucleotide change. The functions of the intermediates along the direct path from EGKA to GSKV suggest the order in which the substitutions are likely to have occurred: of the three single-mutant variants, only one (GGKA) is functional on either RE – indicating that this was probably the first step – and the same is true of the possible second steps (GGKV being the only accessible functional node, Fig. 6A).
The ancestral and derived RH variants that occurred during history – and persist to the present – were accessible only because of pairwise epistasis. Under the main-effects only model, neither EGKA nor GSKV is a strong activator on either RE. Including the terms from the pairwise or third-order models restores the function of these variants, because interactions make large and necessary contributions to their functions. For example, the glutamate at site 1 in the ancestral variant has virtually no main effect on ERE activation, but it makes a major contribution via pairwise interactions with sites 2 and 3 (Fig. 6B). Similarly, the serine at site 2 in the derived variant has no main effect on SRE specificity, but its interactions with sites 3 and 4 are necessary to make this protein an SRE activator. These RH sequences are conserved today in all known ERs and kSRs: therefore the functionality of both historical and present-day steroid receptors depends critically on pairwise epistasis.
The historical transition from ERE-to SRE-specificity also altered the architecture by which REs are recognized, changing the sites and pairs that confer DNA binding rather than just exchanging amino acid states within an otherwise fixed architecture (Fig. 6B). Moreover, this architecture continued to change once the RH became SRE-specific: the final substitution in the trajectory (g2S, which had no net effect on specificity) replaced the initial historical determinants of SRE specificity with a different set of favorable interactions. As a result, the final RH sequence – which is now found in all extant kSRs – has determinants of function that are entirely different from those in the ancestral RH and even those that first mediated the acquisition of SRE specificity during history.
Many of our findings differ from other studies of the genetic architecture of protein function and evolution. These differences appear to arise primarily because most prior work has 1) taken a local approach and focused on the effects of mutating a designated reference protein to a particular derived state at each variable site, 2) worked within a two-state framework, which considers a tiny fraction of the 20N proteins that comprise sequence space when all possible amino acids are considered, and 3) examined a single protein function rather than functional specificity.
A global view of genetic architecture
Most prior studies have attempted to decompose protein genetic architecture by estimating the effects of mutations (and combinations) on a designated reference sequence (27, 81, 84, 86, 87, 104, 107). This strategy provides an accurate local picture of mutation effects when they are introduced into the reference or nearby proteins, but accuracy is poor across other backgrounds and is very sensitive to experimental error (75). One method reduces this problem by averaging estimates of mutation effects across sets of variants containing a state or combination of interest, but it has been implemented only for binary data, and it defines mutation effects relative to a particular reference state, which may make it sensitive to error in the measurement of variants containing the reference state (75). In contrast, our approach provides a globally optimal dissection of genetic architecture by minimizing the total deviation between predicted and observed functions across the entire ensemble of variants, regardless of the number of states allowed at each site.
The shift from a reference-based account to a global model of genetic architecture reframes questions about the relationship of sequence to function and its impacts on evolution. In the local model, the wild-type sequence is taken as a given; the goal is to understand the effect of mutations as the protein moves along paths away from that starting point (usually to a defined endpoint). The starting point itself has no genetic determinants at all; a pair of mutations, for example, changes the function of the original protein via two main effects plus an epistatic deviation.
A reference-free model of genetic architecture, in contrast, asks how sequence states per se – not just a change from a particular beginning state to a particular end state – determine function. The function of every variant in the entire ensemble of sequences is caused by all the states it contains, not by the differences that separate it from a wild-type sequence. Epistasis not only affects the available paths through sequence space from the reference to other nodes; it also determines which nodes are in the space in the first place – the entire set of possible starting, ending, and intermediate genotypes. The effect of mutations is also conceived differently: in the global model, a pair of mutations replaces the main effect of the ancestral states with two new derived main effects, and replaces the original epistatic interactions with every site in the protein – including with those sites at which the state does not change. A mutation does not simply replace one brick in a static genetic architecture; it also changes many of the second-order interactions that determined the protein’s functions but that were invisible if the starting point is viewed as an already functional black-box. As a protein evolves, each substitution changes the effect on function of the state at many other sites, and it modifies the effect that scores of other substitutions may have in the future.
Low-order, dense genetic architecture
Like previous work we found pervasive epistasis (81, 84, 87, 104, 107). Unlike most previous studies, we found that higher-order epistasis plays only a tiny role (but see (87)). We did not measure fourth-order epistasis, but it is unlikely that it will account for much global variation: we obtained good estimates of function from main and pairwise effects alone, and for fourthorder terms to matter, they would have to be considerably larger on average than third-order and second-order terms, but we and others have observed the opposite pattern – effects of shrinking magnitude as epistatic order increases (87).
Why do our findings concerning higher-order epistasis differ from previous studies? A likely cause is that prior work used methods that rely on a reference sequence or a reference state, and/or they did not adequately account for global non-linearities; both of these factors are known to cause spurious inference of higher-order epistasis (75, 104, 107).
Another difference from previous work is that we found a dense genetic architecture at first and second orders, such that main and pairwise effects involving virtually all possible states at all sites make substantial contributions to genetic architecture. Some other studies have found a very sparse architecture, with just a few states or combinations explaining the distribution of functions across variants (75). One cause of this difference may be that we analyzed all 20 amino acid states rather than just two, so the number of possible terms in our model is far larger than studies that analyze only pairs of states present in two high-functioning proteins. Consistent with this idea, combinatorial DMS studies have found that there are often numerous functional sequences that share few to no sequence states at the mutated sites, which necessitates a denser genetic architecture to account for the entire set of functional sequences compared to a set of sequences where only two amino acid states are allowed at each site. Another factor may be that all the sites we examined are in the protein’s “active site” – in this case, making direct contact with the DNA to which the protein binds.
Epistasis and evolution
Prior work has viewed epistasis as an evolutionary constraint, but we found that the primary effect of epistasis is to facilitate functional evolution. Virtually all of those studies, however, considered the effect of epistasis on direct paths through a binary sequence space between two functional sequences designated a priori as starting and ending points, and they have focused on optimization of a single function. In that context, the major effect of epistasis is to create idiosyncrasies in the functions of the intermediates between two high-functioning proteins, which by necessity can only constrain evolution along those paths (108, 109) 108,109. But this framework never considers the effect of epistasis in making the designated sequences functional in the first place or its effect on opportunities for functional innovation, and it excludes from view the vast number of additional potential starting, ending, and intermediate nodes in a 20-state combinatorial sequence space.
Our global analysis reveals the facilitating effect of epistasis on evolution. The global genetic architecture of function generates a huge network of densely connected functional RH sequences that evolution can navigate (110, 111). Epistatic interactions are critical in making protein variants functional, so it dramatically increase the number of functional nodes that can be visited during evolution. Moreover, by determining the protein’s DNA specificity, epistasis brings variants with different specificities closer together in sequence space; this makes regions of sequence space with new functions far easier to access from those with ancestral functions. By increasing the number of functional nodes and shortening the paths between those with different functions, epistasis dramatically increases the probability that evolution will arrive at a new function from anywhere in the network of potential starting points.
Overall, epistasis makes the topography of sequence space idiosyncratic enough that new functions are easily accessible from nodes with different functions – but not so idiosyncratic that functional nodes become inaccessible from each other. Although epistatic interactions do constrain trajectories when we limit our view to arrival at a particular destination from a particular origin, interactions open opportunities to evolve new functions that are often just a step away from the trajectories that were realized during evolution.
Limitations and future work
We are aware of several limitations in our study. First, we used an ordinal model to analyze the genetic architecture of functional data transformed into categorical form. Compared to analysis of continuous data, categorical analysis might reduce precision in estimating the terms of genetic architecture or dampen the epistatic effects of combinations of states that at lower orders are all strongly beneficial or deleterious. Exploring this possibility will require applying reference-free global analysis to continuous data with an appropriate nonlinear correction and then using it to analyze DMS data that are less noisy than those we studied here.
Second, we studied just four residues that directly contact the DNA’s major groove in a single protein. It is possible that our results will not generalize to other proteins or even other sites in the same protein. For example, genetic architecture away from the protein’s active site could be sparser than that in the RH, where every possible state at all sites affected function through main and pairwise effects. However, complete single-mutant scans of all sites in other proteins (and the few studies of all double mutants) shows that most mutations at most sites affect function, consistent with a generally dense genetic architecture – although perhaps less extremely so than in the RH sites that we studied here (33, 51– 53, 112–126). Higher-order epistasis might be more important in other parts of the protein than in the active site, but we know of no structural rationale that might support this expectation: the RH sites are all packed closely in space and bind adjacent nucleotides in the response element, a physical architecture that might be expected to enrich for rather than deplete higher-order effects. We know of no structural or functional reason to expect that high-order epistasis will be more important in other proteins; application of our approach to other proteins can resolve this question. It is possible that third-(and higher-) order interactions could become more important when larger numbers of sites are considered, because the number of third-order terms that affect each variant will increase as the number of variable sites grows (87); this possibility could be evaluated by using reference-free analysis to analyze third-order combinations at a larger number of sites.
Perhaps the greatest challenge to a broader understanding of protein genetic architecture is the immense size of sequence space when more than just a few sites are considered. Covering all combinations of states using DMS is currently practical for at most five variable sites at a time, so a complete characterization of genetic architecture at all orders is possible for only tiny protein fragments. Our finding that high-order epistasis is relatively unimportant implies that modeling and-first second-order terms should be sufficient to account for most genetic variation in a protein’s functions. If generalizable, the focus could shift from complete sampling of all high-order combinations to covering the much smaller set of pairs. It would still be necessary, however, to analyze the effects of pairs across sufficiently diverse backgrounds at other sites to allow effective global averaging of secondorder terms. Evaluating strategies to efficiently accomplish this end without collapsing back into the local tunnel-vision of reference-based, binary state analyses is a key goal for the future.
To what extent are our findings – particularly the lack of high-order epistasis – generalizable from proteins to the genetic architecture of other biological and biosocial traits? Our complete and precise analysis of genetic architecture was possible only because we reduced the problem to the capacity of a few sites in a single protein domain to carry out a simple specific function that can be measured in the laboratory using a high-throughput assay. Our strategy’s biggest virtue is therefore also its greatest liability. Complex phenotypes that are generated by interactions among multiple molecules or loci may have far more complex architectures (127). Supporting this possibility, although third-order amino acid interactions are rare in the DBD, we found that most functional specificity – preferential interactions with one RE or another – is attributable to pairwise interactions between amino acids: this phenomenon represents third-(or potentially fourth-) order intermolecular epistatic interactions between amino acid pairs and the variable nucleotides in the RE. Moreover, our assay was carried out in a controlled and constant environment, and our library was combinatorically complete; this allowed us to dissect genetic architecture free of any environmental variance, population structure, gene-environment interactions, or gene-environment correlations – all factors that complicate the genotype-phenotype relationship in natural populations. Within the DBD recognition helix, the relationship between genetic score and function is complex but potentially tractable, given further advances in experimental and analytical techniques. This lends hope to biochemistry but less to genetics, and none to those who would predict the biological future – or justify the present – of any complex organism by reference to its genotype.
The experimental data analyzed here were first reported in Starr et al. 2017 (93). Libraries of all combinations of all 20 amino acid states at four variable sites in the SR DBD recognition helix (RH, 160,000 total protein variants) were prepared in two different reconstructed ancestral backgrounds: AncSR1 and AncSR1+11P, which contains 11 historical permissive substitutions, which nonspecifically increase DNA affinity without strongly affecting specificity (96). The libraries were transformed into yeast strains in which a fluorescent reporter is driven by the ancestral estrogen receptor response element (ERE) or the derived ketosteroid receptor response element (SRE). Activation by each variant on each RE was measured using a FACS-assisted sort-seq assay and then discretized into three ordered activation categories on each RE: Null if their mean fluorescence was not significantly greater than stop-codon-containing variants, Strong if their mean fluorescence was significantly greater than stop-codon variants and not less than 80% of the historical reference for that RE (AncSR1 on ERE, AncSR1+11P+GSKV on SRE, where GSKV is the sequence of extant ketosteroid receptors in the RH), or Weak if they were neither Null nor Strong. Further details of the categorization can be found in Starr et al. 2017 (93). For model fitting, we only used data from the combinatorial library created on the AncSR1+11P background as it contained the most balanced number of variants across the three classes.
We formulated an ordinal logistic regression model to dissect the effects of amino acid states and combinations on the activation class of RH variants on the two REs. The model is a forward-cumulative model with a proportional log-odds assumption – a standard approach to modeling data discretized into ordered classes of an otherwise continuous variable (in this case, mean fluorescence). The model imposes a linear relationship between the predictors (the genetic states/combinations) and the log-odds that a variant is in a set of higher classes vs. a set of lower classes; log-odds is defined as the logarithm of the ratio of the probabilities of the higher vs. the lower classifications. For three activation classes, the relationship between the predicted class Y (g) of a genetic variant with sequence g is:
where tNW and tW S are thresholds that delimit the boundaries between the Null/Weak and Weak/Strong binding classes, X is a matrix of indicator variables that represent all the possible sequence states and combinations, and β (the variables to be estimated) is a set of coefficients that describe the effect of each state or combination on the log-odds of classification. We evaluated the proportional odds assumption by fitting two simple logistic regression models to the observed activation class data, with variants reclassified as null vs. weak-or-strong, and again as null-or-weak vs. strong. Estimated effects across the two thresholds were similar, as required by the proportional odds assumption. (Figure 1-Figure Supplement 3). However, this assumption can be relaxed if there is evidence that the effects of amino acid states or combinations is dependent on the binding class.
We used a reference-free encoding of indicator variables and effect coefficients. First, we describe the form of the model that would be used if only a single phenotype were measured; we then expand the description to incorporate two different phenotypes. For the single-phenotype model, the indicator matrix X describes the amino acid states and combinations contained by each of the 160,000 protein variants. Each row in X specifies a protein variant with sequence q1q2q3q4 at the four variable sites, and each column represents a particular state at a site (e.g., q1 = A) or a particular combination of states at a combination of sites (e.g., q1q2 = AA, or q1q2q3 = AAA). Each cell is assigned value 1 or 0 to specify whether that state or combination is present or absent in the variant. Each coefficient represents the effect of having one of these states or combinations on the log-odds of being in an activation class, relative to the average across all variants. The complete model therefore consists of 80 first-order effects (20 amino acid states at each of the four sites), 2400 second-order effects (400 pairs of states at each of the six pairs of sites), and 32,000 third-order effects (8000 triplets of states at each of the four triplets of sites). A single global intercept, β0 describes the average across all variants (and is associated with indicator variable with value 1 in all variants). In principle, the model could contain all fourth-order combinations and effects, but we did not include them, because each fourth-order coefficient would be estimated from the observed activation class of a single variant and so would be indistinguishable from measurement error.
To incorporate the effects of genetic states on two different functions, the indicator matrix is expanded to uniquely identify all 320,000 protein-RE complexes based on the RE it contains and the amino acid states/combinations in the protein variant. Each amino acid state/combination is now associated with two coefficients – one for its average effect across the two REs (β which represents the RE-nonspecific effect on activation) and the other for the difference in effect on ERE vs. SRE (σ representing the effect on specificity). Indicator variables for the RE-nonspecific effects are assigned value 1 if an amino acid state/combination is present in a complex and 0 if absent; indicators for specificity effects are assigned value 1 if the complex contains ERE (and the amino acid state/combination is present), or -1 if it contains SRE. The net effect of any amino acid state or combination on ERE activation is therefore its binding coefficient plus its specificity coefficient; its effect on SRE is the binding coefficient minus the specificity coefficient. A global RE coefficient σ0 represents the main effect of having ERE vs SRE in the complex, averaged over all protein variants (with indicator 1 or -1 for complexes containing ERE or SRE, respectively). As in the single-phenotype model, a global intercept coefficient represents the global average across all complexes (β0, with indicator 1 in all complexes).
Although the number of total coefficients in the model is large, each of the 320,000 complexes is described by just 30 non-zero indicator variables: one for the global intercept, one for the main RE effect, and a binding and specificity coefficient for each of the four amino acid states, six pairs, and four triplets that it contains. For a particular variant with sequence q1q2q3q4 the complete model is specified as follows, after coefficients are multiplied by the relevant indicator variables:
where | indicates the RE in the complex, and i, j, and k index the four sites in the RH sequence.
Specific vs non-specific epistasis and the genetic score
In general, epistatic interactions among states can be classified as specific or non-specific (20, 103). Non-specific epistasis occurs if a global non-linear relationship between an underlying physical property and the measured phenotype affects all combinations of states in the same way. Specific epitasis occurs when combinations of states have distinct non-additive effects on the underlying physical properties. Failure to account for non-specific epistasis when modeling epistatic interactions can artificially inflate inferences of specific epistasis (104, 107). Nonspecific epistasis can arise from intrinsically nonlinear relationships between physical properties and measured phenotypes (such as the energy of binding and occupancy of the bound state), limited dynamic range in an assay or biological phenotype, non-linear scaling within the dynamic range, or the imposition of thresholds within the dynamic range for classifying the functionality of variants.
Our model estimates specific and non-specific epistatic effects in a single fitting procedure by including both forms of epistasis in the model. We do so by defining the genetic score y(g) of any variant g with sequence q1q2q3q4 in complex with ERE or SRE as the sum of all the specific genetic effects of the states in g and the complex:
The specific epistatic interactions are estimated as the second- and third-order epistatic terms of the model, which have an additive effect on the genetic score (together with all the non-epistatic first- and zero-order effects). Nonspecific epistasis is then incorporated through the logit function, which makes classification a nonlinear outcome of the genetic score.
Epistatic coding and interpreting coefficient estimates
Most efforts to model sequence-function relationships have encoded genetic effects as the effects of mutations (singly or in combination) relative to a designated reference or “wildtype” sequence. Others have used a Fourier transformation to recode genotypes and estimate the effects of the transformed states relative to the average across all genotypes (76, 101, 128). Here we develop a reference-free approach that directly estimates the effect of any amino acid state or combination on the genetic score relative to the global average across all variants. In our model, the zero-order term is the global average:
where y(g) is the genetic score of a particular genotype g and G is the set of all possible protein sequence/RE element complexes; the factor before the sum averages over all protein genotypes times 2 to reflect the number of REs. The main effect of a particular amino acid state q at a particular site i on the log-odds of activation (not specific to an RE) is the average difference between the genetic score of all variants with that state and the global average:
where Gqi is the set of genotypes with amino acid state q at site i. Similarly, a pairwise epistatic effect is the difference between the average genetic score of Gqiqj – the set of variants containing states qi and qj at a particular pair of sites i and j – and the expected score if there was no pairwise interaction:
Third-order interactions are the difference between the expected genetic score based on the relevant lower-order effects:
The specificity coefficients are encoded similarly. The global RE coefficient reflects the average difference between the genetic score of all complexes containing ERE and the global average (which is of equal magnitude but opposite sign as the average difference of all complexes containing SRE from the global average). This equals half the difference between the average genetic scores for all complexes on ERE and all complexes on SRE:
where y(gE) is the genetic score when a protein variant with sequence g is bound to ERE and y(gS) is the genetic score when a genotype g is bound to SRE. The remaining specificity coefficients follow the same pattern, reflecting differences in specificity beyond what is expected from lowerorder effects. For the main specificity coefficients:
For pairwise specificity coefficients:
For third order specificity coefficients:
Model selection, regularization, and cross-validation
We fit the model to the DMS data using regularized logistic regression in R (129). To avoid overfitting model parameters to measurement noise, we used regularization via ridge regression (L2 norm). To identify the optimal regularization penalty lambda, we used iterated 10-fold cross-validation. Specifically, we fit a series of models with 900 lambda values of decreasing magnitude from 10−1 to 10−8 on a log10 scale. The data were randomly divided into 10 equally sized blocks, and the model at each lambda was fit to 90% of the data, leaving out 10% of the data as a test set, and the activation class of each variant in the test set was predicted on each RE. We compared these predictions to the observed class for each variant, counting misclassification into an adjacent category (weak-null or weak-strong) as a single error and misclassification between null and strong as two errors. We repeated this procedure, using each of the ten blocks as the test set in turn and calculated the mean misclassification rate. We then repeated this entire procedure ten times, dividing the data into different blocks each time, and using the variation in the estimated misclassification ratio to determine the standard error in our estimate. We applied this procedure to all values of lambda and selected the value with the lowest average misclassification rate.
To implement the cumulative odds model with three ordered classes and the proportional odds assumption in a framework that allowed for regularization and cross-validation, we made several modifications to the glmnetcr function in the glmnetcr package (130).
Proportional Odds assumption via logistic regression
In general, the cumulative odds model can be viewed as generating a set of binary variables around each threshold, each of which distinguishes between the set of classes lower and the set of classes higher than the threshold. For our three activation classes, we recoded the categories using two binary variables – one that distinguishes the Null class from Weak-or-Strong and another that distinguishes Null-or-Weak from Strong. Null corresponds to the 00 class, weak to the 10 class, and strong to the 11 class. Logistic regression can then be used to estimate the best-fit model parameters that predict these recoded classes. This produces equivalent estimates of model coefficients as ordinal regression.
Sparse Matrix Representation
Manipulating the large model matrix imposed a large computational burden. Since a large portion of the matrix consists of zeros (each row contains only 30 non-zero values), we used the MatrixModels R package (131) and modification of functions in the glmnetcr package to implement the handling of sparse matrices.
Long vector accommodation
The version of R we used fit logistic regression models with regularization using Fortran. To accommodate the large number of observations and covariates, we modified the lognet function from the glmnet R package to handle long vectors and passed this to the underlying Fortan code using the dotCall64 R package (132).
Zero-sum constraint on coefficients
The reference-free framework requires the sum of model coefficients for a given site (or combination of sites) to sum to zero. All marginal subsets of pairwise or third-order coefficients – defined as the subset of coefficients for a given pair (or triplet) of sites that share a particular state (or pair of states) – must also sum to zero. For example, of the 400 pairwise interactions between sites 1 and 2, the entire set must sum to zero; of these, 20 have an A at site 1, and this marginal set must also sum to zero. This constraint guarantees that all coefficients represent the average effect of having some sequence state or combination relative to prediction from applicable lower-order terms and the global average.
Imposing this constraint during regularized regression is computationally expensive. We therefore performed unconstrained regularized regression and then adjusted the unconstrained coefficients post-hoc to impose this constraint. Specifically, we calculated the average of each set of unconstrained coefficients and each marginal set of unconstrained coefficients. We then subtracted the relevant averages from each unconstrained coefficient, thus recentering each set and marginal subset around an average of zero. For third-order interactions, the constrained estimates are:
Where and are the constrained coefficients, β and σ the unconstrained coefficients, and the averages are over the possible states at a site (or combination) indicated by the *. Pairwise coefficients follow the same pattern, with the total marginal effect of each combination of amino acids to be zero:
Similarly for the main effect terms:
The intercepts then capture the average deviations from zero:
This procedure does not change the genetic score of any variant, and it does not change the likelihood of any model. Constraining the model this way merely centers and scales coefficients so that each represents the average effect of an amino acid or interaction relative to the average of the entire data set (and to the average predicted by lower order terms). In practice, we found that L2 regularization brings the average of most sets of coefficients close to zero anyway, so the post-hoc constraint changes most terms by less than 1% of the unconstrained estimates (Figure 1-Figure Supplement 5).
Alternative models without epistasis
To understand the effects of specific epistasis on protein function and evolution, we compared the third-order model described above to models that do not incorporate any epistasis (first-order model) or no higher-order epistasis (second-order model). To implement these lower-order models, we modified the matrix of indicator variables by removing all indicators for combinations at the orders to be excluded. Model-fitting and selection procedures were otherwise identical.
Comparison to observed DMS data
To estimate the quality of the model fit, we compared the observed class of each protein variant to the observed classification from the empirical data. In each case, we calculated the number of true positives, true negatives, false positives, and false positives for each of the three activation categories, as well as a “non-strong” category (union of null and weak) and an “activators” category (union of weak and strong).
Comparison to measured ΔG
We compared the estimated genetic score for a subset of protein variants to their energy of binding to ERE and SRE, previously measured using fluorescence anisotropy (133). We used simple least-squares linear regression to estimate the relationship between the genetic score and the ΔG of binding. Because both genetic score and ΔG accrue additively, the relationship between these quantities should be linear. The slope and intercept of the relationship between genetic score and ΔG was then used it to estimate the ΔG of binding for all variants in the data set from their genetic score.
Model variance explained by individual terms
The reference-free coding of model coefficients leads to a simple relationship between the effect size of a term and the fraction of specific genetic variance in phenotype that it explains. A detailed proof is provided in Appendix 1 and summarized here.
The specific genetic variance is defined as the variance in phenotype attributable to all the effects of individual genetic states and combinations; it excludes genetic variance caused by nonspecific epistasis (which is captured by the logit function in our model). The total specific genetic variance is the average squared deviation of each variant’s genetic score from the global mean:
where G is the set of all protein variants on either RE, y(g) is the genetic score of a particular protein-RE complex, and 2 · 204 is the total number of complexes.
The total variance explained by a model can be partitioned into the partial variances explained by each causal factor. In the case of our model, the total specific genetic variance can be partitioned into the partial variances attributable to the possible states at each site or (combination of states at multiple sites). Organizing these partial variances by the epistatic order and effect type (nonspecific or RE-specific):
where i, j, and k index the four variable sites, the variances are over all the possible states (or combinations), and the summations are over all sites (or combination).
The model coefficients have a straightforward relationship to this partitioned variance. Main-effect coefficients are defined as the average deviation of a subset of variants containing a state or combination relative to the global mean, and higher-order coefficients are defined as deviations from the sum of lower-order terms. The variance attributable to any state (or combination) is therefore simply the square of its coefficients, and the variance attributable to any set of states is the average of the squared coefficients (see proof in Appendix 1). Thus, for any given site i (or site combination i,j or i,j,k), the variance attributable to the set of model terms applicable to that site can be written as follows:
The leading factors in equations 25-31 can can all be expressed as a function of O, the epistatic order of the effect they represent, where O = 0 for the global average effect of the RE, and 1, 2, or 3 for the main, pairwise, and thirdorder epistatic effects of each amino acid combination. Substituting into equation 24, the total specific genetic variance is therefore the sum of the squared effects of every model term at every order (except for the global average), each divided by the number of model terms at that order:
We can use this same approach to calculate F (V ar(θ)), the fraction of the total specific genetic variance explained by any coefficient θ (a particular β or σ of any order, except the global average):
The fraction of specific genetic variance explained by any set of coefficients is simply the sum of F (V ar(θ)) over all the coefficients in the set.
This formalizes the intuition that within an epistatic order, model terms of large magnitude explain more variation than smaller terms; if two terms at different order have the same magnitude, the lower-order term explains more variation than the higher-order term, because the former affects more genotypes than the latter. The leading terms effectively weight each set of coefficients in a set by the number of genotypes to which they apply.
We used this simple relationship between the magnitude of a model coefficient, its order of epistatic effect, and the fraction of specific genetic variance explained by that coefficient to directly calculate the percent of variance explained by every model term. We then summed the percent variance explained by individual terms to determine the percent of variance explained by sets of terms.
Fraction of genetic score attributable to epistasis or specificity
To find the relative impact of epistasis terms on each variant’s genetic score, we calculated the sum of the absolute value of all pairwise and third-order model terms for that genotype and divided it by the sum of the absolute value of all model terms for that genotype. For specificity, we summed the absolute value of only the specificity terms divided by the total absolute sum of all terms.
Types of epistasis
Epistasis can be divided into several subtypes depending on the signs and magnitude of mutational effects on different genetic backgrounds. These distinctions are most commonly made for pairwise epistatic interactions. For example, sign epistasis arises when the direction of effect of a mutation depends on the genetic background. Reciprocal sign epistasis occurs when the direction of effect of both mutations depends on the genetic background. All other instances of epistasis are classified as magnitude epistasis. Within magnitude epistasis, there are two broad forms, diminishing returns and diminishing costs epistasis. The former, diminishing returns, also includes common forms of epistasis such as amplifying costs, differing only on which genotype is considered ‘wild-type’ or the starting point of a mutant cycle. Likewise, diminishing costs epistasis also covers amplifying returns epistasis, only differing in which genotype is considered the beginning of a mutant cycle.
To classify the types and frequencies of different forms of epistasis, we compared the estimated effect of each pairwise interaction to the estimated effect of each single amino acid involved in that pairwise interaction. The fraction of specific genetic variance explained by the particular pairwise effect was then added to the appropriate category, either sign, reciprocal sign, diminishing returns, or diminishing costs epistasis. To extend this approach to third order interactions, we sequentially fixed each amino acid in the third order interaction, taking that genotype as the starting genotype and then varied the other two amino acids in a manner analogous to the pairwise epistasis calculation. The fraction of specific genetic variance explained by a particular third order interaction was then divided by three (as there are three states that can be fixed) and added to the appropriate category.
Direction of the effect of mutations’ effects
The genetic score of a specific genotype is the summation of the effects in that sequence, i.e. the single global term, four single site effects, six pairwise interactions, and four three-way interactions for both binding and specificity. The model thus gives the effect of each amino acid or interaction relative to the average of all genotypes, not the effect of individual mutations from one state to another. To calculate the effect of a mutation in a given genetic background, we calculated the difference in the genetic score between two genotypes with a single amino acid difference. This entails a change in one main effect, three pairwise interactions, and three third-order interactions. In this way, the model directly captures the fact that even a single amino acid change has the potential to alter numerous epistatic interactions. As a result, to determine the effect of a mutation requires measuring its impact across numerous genetic backgrounds. To determine the range of possible effects of a single mutation, we compared the genetic score between two genotypes differing by a single amino acid on each of the 8000 genetic backgrounds on which the particular change from one amino acid to another could occur and determined the frequency in which that particular mutation increased or decreased binding on a particular DNA element. To minimize the effect of small changes in genetic score upon mutation being counted, we only considered changes in genetic score greater than 5% of the difference between the average genotype and the threshold needed to be classified as a strong activator. We also used the fact that the model identifies which epistatic factors (main, pairwise, third-order) are responsible for the change in genetic score between variants, as well as how this effect is partitioned between binding and specificity, to identify when changes in epistasis were necessary, sufficient, or both for the change in genetic score.
Sequence space analysis
To determine how epistasis affected the evolution of DBD binding specificity, we generated connected networks for each model of epistasis using the R package igraph (134). In each case, genotypes were connected if they differed by a single amino acid (hamming distance) or if the single amino acid difference could occur via a single nucleotide mutation (genetic code distance). In the latter case, we accounted for the disjoint nature of serine codons in the genetic code, treating the two groups as non-interchangeable by coding all serine amino acids in both an S and a Z form. Thus, for a sequence containing n serine amino acid, there are 2n versions of that amino acid sequence with the same phenotype that inhabit different locations in sequence space with different connectivity to neighboring sequences. We created networks for the full model, the model containing up to pairwise epistasis, and the model containing no specific epistasis. We used only those genotypes that were functional in a particular model, i.e. were either weak or strong on either ERE, SRE, or both. In addition, we created a network of all 160,000 genotypes. This latter network contained no information about the function of each genotype and was used to model evolution over sequence space unconstrained by function.
Network distances were calculated as the shortest path distance between two genotypes, constrained to only take steps to neighboring genotypes in a network. Because only functional genotypes are present in a network, the distance between genotypes in a network can exceed four, even when connections are made by hamming distance, if there are no direct paths between the genotypes. We also calculated the distance from each ERE-specific genotype to the nearest SRE-specific genotype, of which there may be several, by finding the minimum distance for each ERE-specific genotype.
Shared vs unique network genotypes
One of the main effects of epistasis is to influence whether a genotype is functional or not. Thus, many genotypes are only functional in one or a subset of networks. However, even if the same genotypes are functional with and without epistasis, epistasis can influence the paths between them that evolution is likely to take. To distinguish between these effects, we compared the effects of epistasis on path lengths for all functional genotypes and for the set of genotypes that are functional in all three networks.
One-step functional transitions and epistatic necessity/sufficiency
To determine how epistasis influenced changing binding classes or specificity, we identified all instances in which two neighboring genotypes in a network differ in binding class or specificity. These represent portions of sequence space where a single amino acid mutation can move a genotype from one binding class to another or alter specificity from ERE-specific to SRE-specific. For changes between binding classes, we identified all such mutations on the full network with third order epistasis. For changes in specificity, we identified the mutation responsible for all such instances in the full network, as well as the networks that had no third-order epistasis or had no epistatic effects.
We calculated whether epistasis was necessary, sufficient, or both for transitions in binding classes and specificity. For binding classes, epistasis was categorized as necessary if the main effect of the amino acid change was not large enough to move the genotype into the new binding class without contribution from the epistatic effects. Similarly, epistasis was categorized as sufficient if epistatic changes were large enough to move the genotype into the new binding class without contribution from the main effects. For single step transitions in specificity, epistasis was necessary if the main, non-epistatic, effect of an amino acid change was not large enough to either reduce ERE binding or increase SRE binding enough to switch the binding class from ERE-specific to SRE-specific. Epistasis was sufficient if the epistatic effects alone were large enough to both reduce ERE binding and increase SRE binding enough to switch the binding class from ERE-specific to SRE-specific. In the full model, a similar distinction was made for third order epistasis relative to main effects and pairwise epistasis and calculated in an analogous manner.
Number of paths and path distinctiveness
The number of paths between two genotypes was calculated as the number of distinct shortest paths between them. Two paths were considered distinct if at least one genotype along the paths were different. To account for the fact that many paths traverse similar sets of genotypes, we also developed a metric for distinctiveness of a set of paths between two genotypes. To calculate this metric, we extracted from the larger network the set of genotypes along any shortest path between the starting and ending genotype. Intuitively, for a path with S genotypes, the set of paths between the starting and ending genotypes should be considered more distinct if they contain more unique genotypes or fewer shared edges between intermediate genotypes. To capture these aspects, we calculated two metrics. The first metric relates the effective number of paths (Pg) to the number of unique genotypes (G) and the path length (S - 1). The second metric relates a different measure of the effective number of paths (Pe) to the number of unique edges (E) and the number of genotypes along a single path (S).
The first calculation provides the maximum number of completely distinct paths with no shared genotypes that could be created given the path length and the number of genotypes present. It thus provides an upper bound on the number of distinct paths possible. If a set of paths are completely distinct, i.e. sharing no genotypes, then Pg will equal Pe. However, if the paths are not completely distinct, then there will be more edges among the set of paths than if all paths were distinct. Thus, for a given number of genotypes and path length, the greater the number of edges connecting those genotypes, the less distinct the set of paths are. By comparing the difference in these metrics to the maximum possible value, we can estimate the effective number of distinct paths (D) in a set.
We used this equation to calculate the distinctiveness of paths from each ERE-specific genotype to each SRE-specific genotype over the network.
To determine whether epistasis increased or decreased the ability of evolution to find a protein variant with altered specificity, we performed forward evolutionary simulations over the previously described genotype networks. These simulations captured neutral wandering among a set of functional variants and are similar to the verbal model of protein evolution described by Maynard Smith (19). In each simulation, evolution was constrained to move among only functional variants, i.e. those classified by the particular model as either ERE-specific, SRE-specific or promiscuous. Genotypes were classified as either ERE-specific if they bound ERE and not SRE, SRE-specific if they bound SRE and not ERE, and promiscuous if they bound both ERE and SRE. Genotypes were considered neighbors if a single nucleotide change could mutate one variant into the other (genetic code) or if they differed by a single amino acid (hamming distance). To estimate the global effect of epistasis, evolution was initiated from each ERE-specific genotype and allowed to take a random step to a neighboring genotype, regardless of DNA binding specificity, until an SRE-specific genotype was reached. This was repeated 100 times for each starting genotype. For the subset of ERE-specific genotypes found in each of the three models with varying amounts of epistasis (Full model, pairwise only, no epistasis), an additional 1000 simulations were performed to estimate how epistasis altered particular evolutionary paths. For each simulation, we recorded the number of steps, ending genotype identity, and the identity of each intermediate genotype during the evolutionary walk.
Model with selection for increased specificity
To estimate how epistasis might interact with selection for increased specificity, we performed additional forward evolutionary simulations. As before, evolution was initiated from ERE-specific genotypes and terminated when an SRE-specific genotype was first encountered. Unlike the neutral simulations, however, we incorporated selection for increased specificity. At each step, the specificity of each neighboring genotype was calculated from the underlying genetic scores as the difference in binding to SRE and ERE. We then used Kimura’s formula for fixation probability to calculate the relative probability of fixation among all possible neighbors (135). The effects of specificity were scaled such that the difference between a non-binder and a strong binder corresponded to a scaled selection coefficient (i.e. Ns) of 1 in a population of size 100. We initiated evolution with a range of population sizes between 1 and 10,000 to capture the transition between neutrality and selection driven changes, performing 100 simulations for each population size on each network using the genetic code to connect genotypes. We noticed during these simulations, that for a small subset of starting genotypes, evolution would often become ‘stuck’ on local optimum in the pairwise model, continuously cycling between two genotypes. This was more common for larger population sizes and had the effect of drastically increasing the average evolutionary path length away from the median, or typical, path length. To account for this effect, we also estimated evolutionary path lengths after removing all evolutionary simulations in which more than 20 steps were taken (Figure 5-Figure Supplement 3).
The coding scripts for running all analyses are located on Github (https://github.com/JoeThorntonLab/DBD.GeneticArchitecture).
We thank members of the Thornton lab for helpful comments on the manuscript and the University of Chicago’s Research Computing Center for computational resources.
Supplementary Note 1: Supplementary Figures
Supplementary Note 2: Appendix 1
Partitioning of Variance
In our modeling formalism, the importance of each model term – the fraction of total specific genetic variance that it accounts for – can be calculated directly from its magnitude and order. In addition, these microscopic contributions can be summed over any set of terms to yield the fraction of variance accounted for by the set. Specifically, the fraction of variance accounted for by any term is equal to the square of the term, weighted by the number of genotypes to which it applies, which is a function of the term’s epistatic order (main, pairwise, third-order, etc.). This property flows ultimately from the constraint we placed on the model coefficients to sum to zero at each site and combination of sites. Here we show why this relationship holds.
We begin with the definition of specific genetic variance, which is the variance in the genetic score among all genotype:
Where g is a particular genotype on a particular DNA element, and y(g) the genetic score (phenotype) of that genotype, and G the set of all possible protein sequences/RE element complexes. The variance is calculated over all genotypes for both ERE and SRE. The phenotype of a particular genotype y(g) can be expressed as the sum of the model coefficients for which that genotype has a particular amino acid state or combination of states. We retain a form of y(g) with both binding and specificity coefficients, although this distinction is not necessary for the proof. The binding coefficients (β s) are the average effect on the genetic score of an amino acid or combination of amino acids, regardless of which DNA element is being considered, while specificity coefficients (σ s) are one half of the difference in the genetic score between binding to ERE and SRE. In the convention we use, binding to ERE adds the specificity effects, while binding to SRE subtracts the specificity terms. Thus, for a particular genotype g with sequence q1q2q3q4 we have:
Here gE and gS are genotypes bound to ERE or SRE respectively, i,j, and k index the available sites, and βq /σq are the corresponding model terms based on the state of genotype g at site i.
The global binding intercept β0 is defined as the average phenotype of all genotypes on the two DNA elements:
Similarly, the global specificity term σ0 reflects the difference in the average phenotype of genotypes on one element relative to the global average, or one half of the difference in average phenotypic binding between the two DNA elements:
The remaining model terms capture deviations from the expected phenotype in the absence of progressively higher order interactions:
Where Gqi is the set of genotypes with state q at site i. Because these terms capture deviations from an average, the set of terms at a site or combination of sites sum to zero by definition:
Where q2 and q3 are the set of all possible pairs or triplets of amino acid states respectively. Using the definition of β0 and inserting equations (2) and (3) into equation (1), we can rewrite the specific genetic variance as:
As expected, the specific genetic variance is independent of the global binding coefficient β0. However, it is not independent of the global specificity coefficient σ0. We can rewrite (20) using a more compact notation as:
Where a is the set of sites or combinations of sites for all orders other than the intercepts. Thus, a includes all of the qi first order terms, the qiqj second order terms, and the qiqjqk third order terms. We next expand the squared term, rewriting the sum as the sum of each squared coefficient plus the sum of the product of all pairwise combinations of coefficients:
Where b also indexes the sites/combinations of sites among the different epistatic orders. Distributing the outer sum over these terms gives:
To simplify this expression, we note that all of the cross-product terms in (23) turn out to be zero. To demonstrate this, we note that the set of all genotypes can be represented by conditioning on a particular amino acid state at a particular site and then taking the union over all possible amino acid states for that site:
We can thus reorganize a sum of cross product terms over all genotypes in (10) as a sum of cross product terms conditioned on a particular set of amino acid states and then sum over all possible amino acid states. For example, the sum of the crossproducts of the binding coefficients with a = q1 and b = q2 can be rewritten by conditioning on the amino acid states at sites 1 and 2, and then summing over all possible combinations of amino acids for those sites:
Where is the set of genotypes with state q1 at site 1 and state q2 at site 2 bound to ERE. From here, we note that once conditioned on a particular set of amino acid states at a set of sites, that combination of amino acid states will appear on all possible combinations of amino acids at the remaining sites. We can thus replace a sum over all genotypes with number of such combinations and then factor the result:
Since each of these latter sums is over the entire set of coefficients at a site, then by construction they sum to zero, and their product is thus also zero:
A similar logic holds for higher order terms, for specificity terms, and for the combination of binding and specificity terms; as long as the conditional sums are over the entire set of amino acid combinations for those sites then we can replace summing over all genotypes with the number of genetic backgrounds a particular combination of states will appear on. For example, for the sum of cross products between binding and specificity terms with a = q1q2 and b =
Similar logic can be applied to all terms in (23) that are the sum of the product of coefficients at other orders: the terms for any pair or triplet of sites always sum to zero, so the sum of the product of those terms (times any other set of terms) also sums to zero.
Combining equation (23) with this result, we can rewrite the specific genetic variance as:
Where the second equality follows because the sets of squared coefficients are the same for binding ERE and SRE. To simplify further, we again use the fact that the set of all genotypes can be partitioned into smaller sets based on the amino acid state at a site. For example, for the squared binding coefficients at site 1:
Similarly, when coefficients involve multiple sites:
And the global specificity term:
Which is the sum of the squared coefficients, each normalized by the number of terms at a particular epistatic order. We can rewrite this more compactly as:
where β and σ represent the individual binding and specificity terms, and O is the epistatic order of each term, with O = 0 for global (intercept) terms, 1 for main effects, 2 for pairwise interactions, and 3 for third order interactions. Said another way, the total variance is the sum of each term’s squared effect, weighted by the number of genotypes that the term affects.
One consequence of the relationship between the effect size of a term, its order of epistasis, and the total specific genetic variance explained, is that the fraction of specific genetic variance explained by any single term is readily calculable. For example, the fraction of specific genetic variance explained by each first order binding term at a particular site is:
In addition, since the total specific genetic variance is the simple weighted sum of squared terms, the fraction of specific genetic variance explained by any set of terms is also easy to calculate. For example, the fraction of specific genetic variance explained by all first order binding terms at a particular site is:
Thus, knowing only the value of a model coefficient and the epistatic order it pertains to, we can readily calculate the fraction of the total specific genetic variance contributed by an individual term and thus any set of terms in the model, including across sites, epistatic orders, mechanisms of action, and the particular form of epistasis.
- 1.Kondrashov. Epistasis as the primary factor in molecular evolutionNature 490:535–538https://doi.org/10.1038/nature11510
- 2.Adaptive landscapes and protein evolutionProceedings of the National Academy of Sciences 107:1747–1751https://doi.org/10.1073/pnas.0906192106
- 3.A multilevel neuromolecular architecture that uses the extradimensional bypass principle to facilitate evolutionary learningPhysica D: Nonlinear Phenomena 75:417–437https://doi.org/10.1016/0167-2789(94)90295-X
- 4.Towards High Evolvability Dynamics IntroductionEvolutionary Systems :33–43
- 5.The Causes and Consequences of Genetic Interactions (Epistasis)Annual Review of Genomics and Human Genetics 20:433–460https://doi.org/10.1146/annurev-genom-083118-014857
- 6.Evolutionary Accessibility of Mutational PathwaysPLoS Computational Biology 7https://doi.org/10.1371/journal.pcbi.1002134
- 7.Evolution and speciation on holey adaptive landscapesTrends in ecology & evolution 12:307–312https://doi.org/10.1016/S0169-5347(97)01098-7
- 8.Percolation on the fitness hypercube and the evolution of reproductive isolationJournal of Theoretical Biology 184:51–64https://doi.org/10.1006/jtbi.1996.0242
- 9.Sequence Entropy and the Absolute Rate of Amino Acid Substitutionsbiorxiv
- 10.Towards a General Theory of Adaptive Walks on Rugged LandscapesJournal of Theo 128:11–45
- 11.The NK model of rugged fitness landscapes and its application to maturation of the immune responseJournal of Theoretical Biology 141:211–245
- 12.Topological features of rugged fitness landscapes in sequence spaceTrends in Genetics 31:24–33https://doi.org/10.1016/j.tig.2014.09.009
- 13.The role of epistasis in protein evolutionNature 497:7–9https://doi.org/10.1038/nature12219
- 14.Epistasis and intramolecular networks in protein evolutionCurrent Opinion in Structural Biology 69:160–168https://doi.org/10.1016/j.sbi.2021.04.007
- 15.The causes of evolvability and their evolutionNature Reviews Genetics 20:24–38https://doi.org/10.1038/s41576-018-0069-z
- 16.Amino acid coevolution induces an evolutionary Stokes shiftProceedings of the National Academy of Sciences 109:E1352–E1359https://doi.org/10.1073/pnas.1120084109
- 17.Exploring protein fitness landscapes by directed evolutionNature reviews. Molecular cell biology 10:866–876https://doi.org/10.1038/nrm2805
- 18.Contingency and entrenchment in protein evolution under purifying selectionProceedings of the National Academy of Sciences 112:E3226–E3235https://doi.org/10.1073/pnas.1412933112
- 19.Natural selection and the concept of a protein spaceNature 225:563–564https://doi.org/10.1038/225563a0
- 20.Epistasis in protein evolutionProtein Science 25:1204–1218https://doi.org/10.1002/pro.2897
- 21.A Model of Substitution Trajectories in Sequence Space and Long-Term Protein EvolutionMolecular biology and evolution 32:542–554https://doi.org/10.1093/molbev/msu318
- 22.Multiple Fitness Peaks and EpistasisAnnual Review of Ecology and Systematics 26:601–629https://doi.org/10.1146/annurev.es.26.110195.003125
- 23.Diminishing returns epistasis among beneficial mutations decelerates adaptationScience 332:1190–2https://doi.org/10.1126/science.1203799
- 24.Experimental Studies of Evolutionary Dynamics in MicrobesTrends in Genetics 34:693–703https://doi.org/10.1016/j.tig.2018.06.004
- 25.Higher-fitness yeast genotypes are less robust to deleterious mutationsScience 366:490–493https://doi.org/10.1126/science.aay4199
- 26.Global epistasis makes adaptation predictable despite sequence-level stochasticityScience 344:1519–1522https://doi.org/10.1126/science.1250939
- 27.Idiosyncratic epistasis creates universals in mutational effects and evolutionary trajectoriesNature Ecology and Evolution https://doi.org/10.1038/s41559-020-01286-y
- 28.Global epistasis emerges from a generic model of a complex traiteLife 10https://doi.org/10.7554/ELIFE.64740
- 29.Diminishing returns and tradeoffs constrain the laboratory optimization of an enzymeNature Communications 3https://doi.org/10.1038/ncomms2246
- 30.Patterns and Mechanisms of Diminishing Returns from Beneficial MutationsMolecular Biology and Evolution 36:1008–1021https://doi.org/10.1093/molbev/msz035
- 31.Diminishing-returns epistasis decreases adaptability along an evolutionary trajectoryNature Ecology and Evolution 1https://doi.org/10.1038/s41559-016-0061
- 32.Mutational effects on stability are largely conserved during protein evolutionProceedings of the National Academy of Sciences 110:21071–21076https://doi.org/10.1073/pnas.1314781111
- 33.A Systematic Survey of an Intragenic Epistatic LandscapeMolecular Biology and Evolution 32:229–238https://doi.org/10.1093/molbev/msu301
- 34.Co-evolution of interacting proteins through non-contacting and non-specific mutationsNature Ecology and Evolution 6:590–603https://doi.org/10.1038/s41559-022-01688-0
- 35.Impact of in Vivo Protein Folding Probability on Local Fitness LandscapesMolecular Biology and Evolution 36:2764–2777https://doi.org/10.1093/molbev/msz184
- 36.Stability-mediated epistasis constrains the evolution of an influenza proteineLife 2https://doi.org/10.7554/eLife.00631
- 37.Epistatically Interacting Substitutions Are Enriched during Adaptive Protein EvolutionPLoS genetics 10https://doi.org/10.1371/journal.pgen.1004328
- 38.Pervasive cryptic epistasis in molecular evolutionPLoS Genetics 6https://doi.org/10.1371/journal.pgen.1001162
- 39.Epistatic drift causes gradual decay of predictability in protein evolutionScience 376:823–830
- 40.An experimental assay of the interactions of amino acids from orthologous sequences shaping a complex fitness landscapePLoS Genetics 15https://doi.org/10.1371/journal.pgen.1008079
- 41.Pervasive contingency and entrenchment in a billion years of Hsp90 evolutionProceedings of the National Academy of Sciences 115:4453–4458https://doi.org/10.1073/pnas.1718133115
- 42.Statistical analysis reveals co-expression patterns of many pairs of genes in yeast are jointly regulated by interacting lociPLoS genetics 9https://doi.org/10.1371/journal.pgen.1003414
- 43.Thornton. An epistatic ratchet constrains the direction of glucocorticoid receptor evolutionNature 461:515–519https://doi.org/10.1038/nature08249
- 44.Evolution of hormone-receptor complexity by molecular exploitationScience 312:97–101https://doi.org/10.1126/science.1123348
- 45.Contingency between Historical Substitutions in the Acetylcholine Receptor PoreACS chemical neuroscience 11:2861–2868https://doi.org/10.1021/ACSCHEMNEURO.0C00410
- 46.Stability-Mediated Epistasis Restricts Accessible Mutational Pathways in the Functional Evolution of Avian HemoglobinMolecular Biology and Evolution 34:1240–1251https://doi.org/10.1093/molbev/msx085
- 47.Crystal structure of an ancient protein: evolution by conformational epistasisScience 317:1544–1548https://doi.org/10.1126/science.1142819
- 48.Shifting mutational constraints in the SARS-CoV-2 receptor-binding domain during viral evolutionScience 377:420–424
- 49.Epistasis Constrains Mutational Pathways of Hemoglobin Adaptation in High-Altitude PikasMolecular biology and evolution 32:287–298https://doi.org/10.1093/molbev/msu311
- 50.Epistasis in a fitness landscape defined by antibody-antigen binding free energyCell systems 8https://doi.org/10.1016/J.CELS.2018.12.004
- 51.The genetic landscape of a physical interactioneLife 7:1–31https://doi.org/10.7554/eLife.32472
- 52.A comprehensive biophysical description of pairwise epistasis throughout an entire protein domainCurrent Biology 24:2643–2651https://doi.org/10.1016/j.cub.2014.09.072
- 53.Coevolution-based inference of amino acid interactions underlying protein functioneLife 7https://doi.org/10.7554/eLife.34300
- 54.Evolving New Protein-Protein Interaction Specificity through Promiscuous IntermediatesCell 163:1–13https://doi.org/10.1016/j.cell.2015.09.055
- 55.The adaptive landscape of a metallo-enzyme is shaped by environment-dependent epistasisNature Communications 12https://doi.org/10.1038/s41467-021-23943-x
- 56.Intermolecular epistasis shaped the function and evolution of an ancient transcription factor and its DNA binding siteseLife 4https://doi.org/10.7554/eLife.07864
- 57.Compensatory mutations restore fitness during the evolution of dihydrofolate reduc-taseMolecular Biology and Evolution 27:2682–2690https://doi.org/10.1093/molbev/msq160
- 58.Exploring the effect of sex on empirical fitness landscapesAmerican Naturalist 174:S15–S30https://doi.org/10.1086/599081
- 59.Determinants of BH3 Binding Specificity for Mcl-1 versus Bcl-xLJournal of Molecular Biology 398:747–762https://doi.org/10.1016/j.jmb.2010.03.058
- 60.Distinct modes of regulation by chromatin encoded through nucleosome positioning signalsPLoS computational biology 4https://doi.org/10.1371/journal.pcbi.1000216
- 61.Epistatic mutations in PUMA BH3 drive an alternate binding mode to potently and selectively inhibit anti-apoptotic Bfl-1eLife 6:1–23https://doi.org/10.7554/eLife.25541
- 62.RNA-guided editing of bacterial genomes using CRISPR-Cas systemsNat Biotechnol 31:233–239https://doi.org/10.1038/nbt.2508
- 63.Negative epistasis between beneficial mutations in an evolving bacterial populationScience 332:1193–6https://doi.org/10.1126/science.1203801
- 64.Equally parsimonious pathways through an RNA sequence space are not equally likelyJournal of Molecular Evolution 45:278–284https://doi.org/10.1007/PL00006231
- 65.Stepwise acquisition of pyrimethamine resistance in the malaria parasiteProceedings of the National Academy of Sciences of the United States of America 106:12025–12030https://doi.org/10.1073/pnas.0905922106
- 66.The biochemical architecture of an ancient adaptive landscapeScience 310:499–501https://doi.org/10.1126/science.1115649
- 67.Ancestral lysozymes reconstructed, neutrality tested, and thermostability linked to hydro-carbon packingNature 345:86–89https://doi.org/10.1038/345086a0
- 68.Quantitative Description of a Protein Fitness Landscape based on Molecular FeaturesMolecular Biology and Evolution 32:1774–1787https://doi.org/10.1093/molbev/msv059
- 69.How mutational epistasis impairs predictability in protein evolution and designProtein Science 0:1–13https://doi.org/10.1002/pro.2876
- 70.Step-wise enhancement of catalytic performance of haloalkane dehalogenase LinB towards β-hexachlorocyclohexaneAMB Express 4https://doi.org/10.1186/S13568-014-0072-5
- 71.Epistasis Among Adaptive Mutations in Deer Mouse HemoglobinScience 340:1324–1327https://doi.org/10.1126/science.1236862
- 72.Intramolecular epistasis and the evolution of a new enzymatic functionPLoS ONE 7https://doi.org/10.1371/journal.pone.0039822
- 73.Quan-titative exploration of the catalytic landscape separating divergent plant sesquiterpene synthasesNature chemical biology 4:617–23https://doi.org/10.1038/nchembio.113
- 74.Delayed commitment to evolutionary fate in antibiotic resistance fitness landscapesNature Communications 6https://doi.org/10.1038/ncomms8385
- 75.Learning the pattern of epistasis linking genotype and phenotype in a proteinNature Communications 10https://doi.org/10.1038/s41467-019-12130-8
- 76.The Context-Dependence of Mutations: A Linkage of FormalismsPLOS Computational Biology 12https://doi.org/10.1371/journal.pcbi.1004771
- 77.Constructing and analyzing the fitness landscape of an experimental evolutionary processChemBioChem 9:2260–2267https://doi.org/10.1002/cbic.200800371
- 78.Darwinian Evolution Can Follow Only Very Few Mutational Paths to Fitter ProteinsScience 312:111–114https://doi.org/10.1126/science.1123539
- 79.Higher-order epistasis shapes the fitness landscape of a xenobiotic-degrading enzymeNature Chemical Biology 15:1120–1128https://doi.org/10.1038/s41589-019-0386-3
- 80.Fisher’s geometrical model of fitness landscape and variance in fitness within a changing environmentEvolution 66https://doi.org/10.1111/j.1558-5646.2012.01610.x
- 81.Higher-order epistasis creates idiosyncrasy, confounding predictions in protein evolutionbioRxiv https://doi.org/10.1101/2022.09.07.505194
- 82.Empirical fitness landscapes and the predictability of evolutionNature Reviews Genetics 15:480–490https://doi.org/10.1038/nrg3744
- 83.Inferring a complete genotype-phenotype map from a small number of measured phenotypesPLoS Computational Biology 16https://doi.org/10.1371/journal.pcbi.1008243
- 84.High-order epistasis shapes evolutionary tra-jectoriesPLoS Computational Biology 13https://doi.org/10.1371/journal.pcbi.1005541
- 85.Quantitative analyses of empirical fitness landscapesJournal of Statistical Mechanics: Theory and Experiment 2013https://doi.org/10.1088/1742-5468/2013/01/P01005
- 86.FISHER’S GEOMETRIC MODEL OF ADAP-TATION MEETS THE FUNCTIONAL SYNTHESIS: DATA ON PAIRWISE EPISTASIS FOR FITNESS YIELDS INSIGHTS INTO THE SHAPE AND SIZE OF PHENOTYPE SPACEEvolution 67:2957–2972https://doi.org/10.1111/evo.12156
- 87.The Influence of Higher-Order Epistasis on Biological Fitness Landscape TopographyJournal of Statistical Physics 172:208–225https://doi.org/10.1007/s10955-018-1975-3
- 88.Diversification of DNA-Binding Specificity by Permissive and Specificity-Switching Mutations in the ParB/Noc Protein FamilyCell Reports 32https://doi.org/10.1016/j.celrep.2020.107928
- 89.Uncovering the basis of protein-protein interaction specificity with a combinatorially complete libraryeLife 9https://doi.org/10.7554/eLife.60924
- 90.Engineering orthogonal signalling pathways reveals the sparse occupancy of sequence spaceNature 574:702–706https://doi.org/10.1038/s41586-019-1639-8
- 91.Podgornaia and Michael T Laub. Pervasive degeneracy and epistasis in a protein-protein interfaceScience 347:673–677https://doi.org/10.1126/science.1257360
- 92.Origins of Allostery and Evolv-ability in Proteins: A Case StudyCell 166:468–480https://doi.org/10.1016/j.cell.2016.05.047
- 93.Exploring protein sequence–function landscapesNature Biotechnology 35:125–126https://doi.org/10.1038/nbt.3786
- 94.Translation dynamics of single mRNAs in live cells and neuronsScience https://doi.org/10.1126/science.aaf1084
- 95.Bridging non-overlapping reads illuminates high-order epistasis between distal protein sites in a GPCRNature Communications 11https://doi.org/10.1038/s41467-020-14495-7
- 96.Evolution of DNA Specificity in a Transcription Factor Family Produced a New Gene Regulatory ModuleCell 159:58–68https://doi.org/10.1016/j.cell.2014.09.003
- 97.Comparative Vertebrate Endocrinology
- 98.Analysis of estrogen response element binding by genetically selected steroid receptor DNA binding domain mutants exhibiting altered specificity and enhanced affinityJournal of Biological Chemistry 274:23591–23598https://doi.org/10.1074/jbc.274.33.23591
- 99.Determinants of cell- and gene-specific transcriptional regulation by the glu-cocorticoid receptorPLoS Genetics 3:927–938https://doi.org/10.1371/journal.pgen.0030094
- 100.Genomic actions of estrogen receptor alpha: What are the targets and how are they regulated?Endocrine-Related Cancer 16:1073–1089https://doi.org/10.1677/ERC-09-0086
- 101.Maximally efficient modeling of DNA sequence motifs at all levels of complexityGenetics 187:1219–1224https://doi.org/10.1534/genetics.110.126052
- 102.Missense meanderings in sequence space: A biophysical view of protein evolutionNature Reviews Genetics 6:678–687https://doi.org/10.1038/nrg1672
- 103.Biophysical mechanisms for large-effect mutations in the evolution of steroid hormone receptorsProceedings of the National Academy of Sciences 110:11475–11480https://doi.org/10.1073/pnas.1303930110
- 104.Inferring the shape of global epistasisProceedings of the National Academy of Sciences 115:E7550–E7558https://doi.org/10.1073/pnas.1804015115
- 105.Epistasis–the essential role of gene interactions in the structure and evolution of genetic systemsNature reviews Genetics 9:855–67https://doi.org/10.1038/nrg2452
- 106.Molecular ensembles make evolution unpredictableProceedings of the National Academy of Sciences :1–34https://doi.org/10.1073/pnas.1711927114
- 107.Detecting high-order epistasis in nonlinear genotype-phenotype mapsGenetics 205:1079–1088https://doi.org/10.1534/genetics.116.195214
- 108.Empirical fitness landscapes reveal accessible evolutionary pathsNature 445:383–386https://doi.org/10.1038/nature05451
- 109.Perspective: Sign epistasis and genetic constraint on evolutionary trajectoriesEvolution 59:1165–1174
- 110.Genotype network intersections promote evolutionary innovationPLOS Biology 17https://doi.org/10.1371/JOURNAL.PBIO.3000300
- 111.Latent phenotypes pervade gene regulatory circuitsBMC systems biology 8https://doi.org/10.1186/1752-0509-8-64
- 112.A fundamental protein property, thermodynamic stability, revealed solely from large-scale measurements of protein functionProceedings of the National Academy of Sciences 109:16858–16863https://doi.org/10.1073/pnas.1209751109
- 113.An Experimentally Determined Evolutionary Model Dramatically Improves Phylogenetic FitMolecular biology and evolution 31:1956–1978https://doi.org/10.1093/molbev/msu173
- 114.Ultra-low input single tube linked-read library method enables short-read second-generation sequencing systems to generate highly accurate and economical long-range sequencing information routinelyGenome Research https://doi.org/10.1101/gr.260380.119
- 115.Predictive Bcl-2 family binding models rooted in experiment or structureJournal of Molecular Biology 422:124–144https://doi.org/10.1016/j.jmb.2012.05.022
- 116.A Comprehensive, High-Resolution Map of a Gene’s Fitness LandscapeMolecular biology and evolution 31:1581–1592https://doi.org/10.1093/molbev/msu081
- 117.High-resolution mapping of protein sequence-function relationshipsNature methods 7:741–746https://doi.org/10.1038/nmeth.1492
- 118.The fitness landscape of the codon space across environmentsbioRxiv https://doi.org/10.1101/252395
- 119.Experimental illumination of a fitness landscapeProceedings of the National Academy of Sciences 108:7896–901https://doi.org/10.1073/pnas.1016024108
- 120.The spatial architecture of protein function and adaptationNature 490:138–142https://doi.org/10.1038/nature11500
- 121.Analyses of the effects of all ubiquitin point mutants on yeast growth rateJournal of Molecular Biology 425:1363–1377https://doi.org/10.1016/j.jmb.2013.01.032
- 122.Local fitness landscape of the green fluorescent proteinNature 533:397–401https://doi.org/10.1038/nature17995
- 123.Activity-enhancing mutations in an E3 ubiquitin ligase identified by high-throughput mutagenesisProceedings of the National Academy of Sciences of the United States of America 110:1263–1272https://doi.org/10.1073/pnas.1303309110
- 124.Deep Mutational Scanning of SARS-CoV-2 Receptor Binding Domain Reveals Constraints on Folding and ACE2 BindingCell 182:1295–1310https://doi.org/10.1016/j.cell.2020.08.012
- 125.The inherent mutational tolerance and antigenic evolvability of influenza hemagglutinineLife 3https://doi.org/10.7554/eLife.03300
- 126.Optimization of affinity, specificity and function of designed influenza inhibitors using deep sequencingNature biotechnology 30:543–8https://doi.org/10.1038/nbt.2214
- 127.Idiosyncratic epistasis leads to global fitness–correlated trendsScience 376:630–635https://doi.org/10.1126/SCIENCE.ABM4774
- 128.On the sparsity of fitness functions and implications for learningProceedings of the National Academy of Sciences of the United States of America 119
- 129.R: A language and environment for statistical computing
- 130.glmnetcr: An R Package for Ordinal Response Prediction in High-dimensional Data Settings
- 131.MatrixModels: Modeling with Sparse and Dense Matrices
- 132.dotCall64: An R package providing an efficient interface to compiled C, C++, and Fortran code supporting long vectorsSoftwareX 7:217–221https://doi.org/10.1016/J.SOFTX.2018.06.002
- 133.Strong selection genome-wide enhances fitness trade-offs across environments and episodes of selectionEvolution 68:16–31https://doi.org/10.1111/evo.12259
- 134.The igraph software package for complex network researchInterJournal
- 135.ON THE PROBABILITY OF FIXATION OF MUTANT GENES IN A POPULATIONGenetics