Should evolutionary geneticists worry about higher-order epistasis?

Abstract

Natural selection drives evolving populations up the fitness landscape, the projection from nucleotide sequence space to organismal reproductive success. While it has long been appreciated that topographic complexities on fitness landscapes can arise only as a consequence of epistatic interactions between mutations, evolutionary genetics has mainly focused on epistasis between pairs of mutations. Here we propose a generalization to the classical population genetic treatment of pairwise epistasis that yields expressions for epistasis among arbitrary subsets of mutations of all orders (pairwise, three-way, etc.). Our approach reveals substantial higher-order epistasis in almost every published fitness landscape. Furthermore we demonstrate that higher-order epistasis is critically important in two systems we know best. We conclude that higher-order epistasis deserves empirical and theoretical attention from evolutionary geneticists.

Figure 1

Epistatic coefficients as a function of order. Because epistatic coefficients may be positive or negative (Box 1), mean squared values are shown. The zeroth-order epistatic coefficient is the mean fitness across all genotypes (see Box 1); in each case here, fitness values were normalized to make this quantity equal to 1.0. First-order and second-order coefficients are analogous to classical selection coefficients and classical pairwise epistasis terms, respectively (see Box 1). Error bars represent standard deviation among coefficients of given order; those that extend to the x-axis overlap 0. (a) Theoretical fitness landscapes. Additive (circle) genotypic fitness values are the sum of the fitness effects of constituent mutations, which in turn were drawn uniformly on the interval [0, 1]. Here, all pairwise and higher epistatic coefficients are zero. I.I.D. (squares) genotypic fitness values were drawn independently and identically from a uniform distribution over [0, 1]. Here, magnitude of mean squared epistatic coefficient increases exponentially with order. We expect empirical results to lie in between these two extremes. Enzyme folding stability model (triangles) considers the fitness landscape defined by 1/(1 + e^ΔG/k_bT) over 7 missense mutations with identical and additive ΔΔG = 1 kcal/mol [44••]. Here ΔG is the free energy of folding, k_b is Boltzman’s constant and T is temperature. See text for further details. (b–o) Empirical fitness landscapes in Table 1; citations given in square brackets. Growth rate (panels c, e, h, l, n and o) and drug resistance (f, i and k) values were log-transformed before epistatic coefficients were computed. In cases where more than one combinatorially complete subset of mutations was identified (panels d, g, k, m) results for a randomly selected subset is shown.

Figure 2. The Walsh transform of the fitness landscape W⃑ into Walsh coefficients E⃑

Here we consider the fitness landscape defined by all combinations of L = 3 mutations in the avian lysozyme characterized by Malcolm et al. [[25]; melting temperature is used as a proxy for fitness]. Each row is ordered by a binary string whose bits left-to-right correspond to the T40S, I55V and S91T mutations. (Here, T, S, I and V stand for threonine, serine, isoleucine and valine, respectively, and the number is the mutated residue in the enzyme.) In the case of the fitness landscape vector W⃑, each ‘1’ in the string signals a contribution from that mutation to the corresponding fitness value. In the Walsh coefficients vector W⃑, each ‘1’ in the string signals a contribution from that mutation to the corresponding interaction coefficient. Thus for example we observe that the Walsh coefficient corresponding to the S91T mutation is equal to −1.53°C (second line) and that Walsh coefficient corresponding to the I55V and S91T mutations (fourth line) is equal to zero. ψ can be written for arbitrary L, as for example with the hadamard() function in the software package Matlab (Mathworks, Natick, MA).