Epistasis and Adaptation on Fitness Landscapes

Probabilistic fitness landscape models allow for the definition of arbitrarily epistatic fitness landscapes with few parameters and therefore lend themselves to mathematical analysis of the consequences of epistasis (Krug, 2021). They are agnostic to any mechanistic or molecular reasons for epistasis and are only based on probability distributions. The RMF model (Aita et al., 2001, 2000) is highlighted here as representative of this class of models; see Hwang et al. (2018), de Visser & Krug (2014), Ferretti et al. (2016) for descriptions of the various other members of this class, including the prominent NK model (Kauffman & Weinberger, 1989). In the RMF model, each locus has an additive effect, either taken to be a constant or drawn from a probability distribution, and each genotype is assigned an epistatic deviation from its additive expectation, drawn from a normal distribution with mean 0.

The RMF model interpolates between a fully additive fitness landscape, in which each substitution contributes a fixed amount to its carrier’s fitness independent of the rest of the genome (i.e., there is no epistasis), and the so-called House-of-Cards fitness landscape (Kauffman & Levin, 1987), which maximizes epistasis by assigning a random fitness to each genotype (Figure 1). The difference between the two limiting cases highlights an essential consequence of epistasis: it makes fitness of a genotype less predictable. For example, for a fitness landscape with $L=10$ loci and two alleles per locus, resulting in $2^{L}=1024$ genotypes, only 10 measurements are required to correctly predict the fitness of all genotypes in an additive fitness landscape, whereas 1024 measurements are necessary to fully determine the House-of-Cards landscape.

Probabilistic fitness landscapes such as the RMF model are an ideal testing ground to study the evolutionary consequences of epistasis when the landscape is tuned from additive to epistatic (Neidhart et al., 2014). In the smoothest, non-epistatic landscape, there is only a single, global fitness peak and any population will eventually climb this peak (Figure 1A,F). At the same time, the evolutionary history of independently evolved populations will tend to be different since any order of substitutions is possible. With maximum epistasis, the RMF landscape is equivalent to the House-of-Cards landscape and displays a large number of fitness peaks, which is on average $2^{L}/(L+1)$ (Kauffman & Levin, 1987; Figure 1E,H). Overall, epistasis introduces the potential for parallel fitness evolution towards different adapted genotypes (local fitness peaks). However, an increasing number of peaks in the landscape decreases the length of adaptive trajectories and less of the full landscape is explored by evolution. Consequently, well-adapted genotypes are most divergent for intermediate ruggedness of the landscape (Figure S2). Moreover, epistasis tends to create historical contingencies of adaptation. It limits the possible orders of substitutions and thus the number of evolutionary trajectories towards a fitness peak, making evolution from the same starting point more repeatable and seemingly more predictable. (Figure 1, Figure S2). Notably, the regularity of the RMF landscape (and other probabilistic fitness landscape models, Hwang et al., 2018), in which every genotype has the same statistical properties, is an unrealistic assumption for real interaction networks between mutations or proteins. Nevertheless, the RMF model was shown to recover features of various experimental fitness landscapes (Szendro et al., 2013, de Visser & Krug, 2014, Bank et al., 2016).

Phenotype-fitness models such as Fisher’s Geometric Model (FGM; Fisher, 1930, reviewed in Tenaillon, 2014) provide a perspective on epistasis complementary to that of probabilistic fitness landscape models. Importantly, under FGM, mutations that are additive on the phenotypic level can be epistatic on the fitness level, because phenotype and fitness are connected by a nonlinear fitness function (Figure 2A; Martin et al., 2007, Gros et al., 2009). In FGM, fitness is determined by the distance of a phenotype from a single phenotypic optimum in a multidimensional trait space, usually under the assumption of a Gaussian fitness function (i.e., there is stabilizing selection for an optimum phenotype). New mutations appear as jumps in a random direction originating in the current phenotype. A mutation is considered beneficial when it decreases the distance from the phenotypic optimum. Epistasis occurs, for example, when the same mutation appears at different distances from the optimum (Figure 2B).

FGM turns into a genotype-fitness landscape by introducing an explicit underlying genotype-phenotype map, allowing for the comparison of FGM-based and experimental fitness landscapes (Blanquart et al., 2014, Hwang et al., 2017). Interestingly, the single-peaked phenotype-fitness landscape of FGM can give rise to an underlying genotype-fitness landscape with multiple peaks (Hwang et al., 2017, Park et al., 2020, Srivastava & Payne, 2022). The amount of epistasis for fitness at the genotype level depends on the distance from the fitness peak: far away from the optimum, the same mutation tested on nearby genetic backgrounds tends to be non-epistatic, whereas strong epistasis can occur when the same mutation finds itself on nearby genetic backgrounds close to the phenotypic optimum. Interestingly and importantly for experimental design, the sampling choice of mutations from the greater fitness landscape influences the shape of the combinatorily complete subset of genotypes built by those mutations: individually beneficial substitutions, measured in a reference genetic background, will tend to result in a more epistatic fitness landscape than substitutions that have occurred subsequently in an adaptive walk (Blanquart et al., 2014).

The last subsections showed that in FGM, epistasis arises from a non-linear phenotype-fitness map, whereas RMF and other probabilistic fitness landscape models directly map genotype to fitness without any assumptions about the mapping of genotype to phenotype(s) to fitness. Importantly, epistasis for fitness can arise both through an epistatic genotype-phenotype map and a non-linear phenotype-fitness map, and these levels may interact to yield a (potentially epistatic) genotype-fitness map (Figure 2C). This leads to an open question: how much epistasis for fitness can be attributed to an epistatic genotype-phenotype map, and how much is due to non-linear phenotype-fitness relationships?

Ultimately, the genotype is connected to fitness by multiple layers of multiple phenotypes, and thus the picture discussed here and sketched in Figure 2C is greatly oversimplified. In an elegant theoretical study, the propagation of epistasis across hierarchical levels in a metabolic network was recently studied by Kryazhimskiy (2021). They found that negative epistasis cannot be converted into positive epistasis at higher levels, whereas the opposite is possible. Interestingly, their work suggests that epistasis does not get diluted at higher levels but tends to become stronger.

Obtaining general mathematical results about evolution on theoretical fitness landscapes requires additional assumptions about the evolutionary dynamics. Most genotype-fitness landscape models address haploid, non-recombining organisms in the limit of strong selection and weak mutation (i.e., neglecting neutral substitutions and assuming that mutations appear and fix sequentially Gillespie, 1984; see Fragata et al., 2019 for a review considering fitness landscape models with neutral substitutions and the consequences of genetic drift). For example, when recombination is included in the study of adaptive walks on fitness landscapes, the connectedness of genotypes in the fitness landscape changes because recombination produces new combinations of (multiple) segregating alleles (Moradigaravand & Engelstädter, 2012, Nowak et al., 2014). For rugged RMF landscapes, recombination creates an initial advantage to recombining populations that is later lost during adaptation because recombination dilutes genotypes necessary for crossing fitness valleys. Thus, theory predicts that recombining populations adapting on an epistatic fitness landscape can get stalled at local optima even more severely than non-recombining populations (Nowak et al., 2014).

Studying diverse and recombining populations on fitness landscapes creates conceptual and technical questions. In a diploid fitness landscape, dominance has to be accounted for with respect to the direct effect of individual mutations and their epistasis because it affects both the duration for which epistatic alleles segregate and their ultimate fate in a population (see, e.g. Turelli & Orr, 2000, Blanckaert & Bank, 2018). For example, there is currently no commonly used biologically reasonable and parameter-light definition of the dominance of epistasis, i.e., of the effect of combinations of one or two copies of sets of epistatic alleles.

Finally, population dynamics considering multiple coexisting genotypes are difficult to characterize mathematically and to simulate computationally. If a starting population with standing genetic variation is to be considered, how should this population be distributed in the fitness landscape? Finding simple, computationally feasible, and yet biologically relevant assumptions is a major challenge underlying theoretical studies of fitness landscapes.

As an intermediate between theoretical fitness landscapes and data from nature, experimental fitness landscapes are appealing for several reasons. Firstly, depending on the design of the mutant library they may contain all genotypes of interest, not only those which are segregating naturally; thus, deleterious or weakly beneficial genotypes can be measured, and inference of the selection coefficient (i.e., the relative competitive advantage of a genotype) can be performed for each genotype, resulting in a measurement of the (potentially combinatorially complete) fitness landscape. Secondly, the lab environment is controlled such that all genotypes grow or evolve in a constant, known, and similar environment. Finally, samples can be resequenced to check for secondary mutations, which are anyways minimized in the setup of deep mutational scanning, where populations are observed for only a few generations. Together, these features result in elegant and reproducible fitness data.

Many fitness landscape studies come from systems and molecular biology, where epistasis is considered a natural consequence of interactions within and between genes (reviewed in Domingo et al., 2019, Phillips, 2008, Costanzo et al., 2019, Athavale et al., 2014). Because intra- and intergenic interactions intuitively arise from the physical organisation of proteins and the connection of genes in biological pathways, the question is not whether there is epistasis but how much, and which models of epistasis are most suitable. Among a multitude of studies that screened hundreds or thousands of engineered genotypes for (usually pairwise) epistasis, most have focused on measurements of phenotypes that may be fitness-related but do not directly translate into competitive fitness. Across these studies, pairwise epistasis was found to be common and there is increasing evidence of higher-order epistasis (Domingo et al., 2019). Moreover, an increasing number of studies have followed evolutionary trajectories in the laboratory and dissected epistatic effects along those trajectories (e.g. Kryazhimskiy et al., 2014, Johnson et al., 2019). This review focuses on engineered slices of fitness landscapes built from combinations of more than two mutational steps and with an emphasis on studies measuring competitive fitness.

Although the extent of inferred epistasis in experimental fitness landscapes differed between study system, measure of fitness, and environment, many studies shared the observation of negative (diminishing-returns) epistasis, where the effect of the same beneficial mutation was smaller when it was measured in a higher-fitness genetic background (see de Visser & Krug, 2014 for an earlier review). Diminishing-returns epistasis results in a gradual slow-down of fitness evolution and is predicted across various fitness landscape models when the population approaches a fitness peak (Martin et al., 2007, Lyons et al., 2020, Draghi & Plotkin, 2013, Greene & Crona, 2014). Diminishing-returns epistasis is predicted generally during adaptation on a fitness landscape when there is widespread “idiosyncratic” epistasis, corresponding to epistasis on the genotype-phenotype level in Figure 2C; Lyons et al. (2020), Reddy & Desai (2021)). Models of idiosyncratic epistasis can, for example, explain he fitness trajectories obtained from a long-term evolution experiment in E. coli, in which even after 50000 generations of evolution no fitness peak was reached (Wiser et al., 2013). Interestingly, the epistasis that occurs when an additive trait is translated into fitness through a non-linear phenotype-fitness map (also termed “nonspecific epistasis”; Domingo et al., 2019; see for example the genotype-phenotype-fitness map of FGM in Figure 2A) is not sufficient to guarantee persistent diminishing-returns epistasis (Reddy & Desai, 2021). Thus, theory indicates that (potentially ubiquitous) epistasis on the genotype-phenotype level may contribute to observed evolutionary trajectories.

Retracing adaptive steps during evolution in the wild is becoming feasible with the advance of experimental methods, especially in organisms that lend themselves to experimental manipulation in the laboratory. However, few studies to date have explicitly screened evolutionary trajectories for the presence of epistasis for fitness (see Storz, 2018 for a review of studies of genotype-phenotype epistasis). Pokusaeva et al. (2019), introduced above, built intragenic fitness landscapes including extant variation across various yeast species. They found ubiquitus and often higher-order epistasis. The epistasis resulted in fitness ridges such that fit genotypes tended to be located near extant variation in the sequence space. Thus, most likely epistasis played a role during evolution by defining the accessible sequence space, including potentially neutrally accumulated divergence.

In a study of a protein under strong selection, Gong et al. (2013), Gong & Bloom (2014) dissected the epistatic evolutionary trajectory of the human influenza nucleoprotein (NP) between 1968 and 2007. They used reverse genetics approaches to create genotypes carrying individual substitutions that had occurred during the evolutionary trajectory on the original genetic background from 1968, and used the growth rate of these genotypes in cell cultures as measure of their fitness. Gong et al. (2013) found that several of the substitutions that had fixed during the evolutionary history of NP were individually deleterious in the background of the original sequence, but enabled the virus to subsequently accumulate beneficial immune-system evading substitutions. Interestingly, no epistatic trajectories were observed in NP in swine influenza (Gong & Bloom, 2014). Whereas humans tend to get infected by influenza multiple times, swine are rarely infected more than once, which leads to much stronger selection for immune invasion in human influenza than in swine influenza. The results from Gong et al. (2013), Gong & Bloom (2014) indicate that strong selection was a prerequisite for observing epistasis in the evolutionary trajectory. Whereas laboratory studies in bacteria support that more strongly selected mutations have stronger interactions (e.g. Schenk et al., 2013), it remains unknown whether epistatic trajectories in nature are enriched for individually strongly selected substitutions. Whether this pattern is predicted by theory depends on the considered model scenario. For example, FGM predicts strong epistasis between strongly selected substitutions when the optimum is close, but not when it is far away.

One example of the dissection of an evolutionary trajectory in complex organisms comes from Karageorgi et al. (2019), who studied a convergent evolutionary trajectory of the adaptation of multiple insects including monarch butterflies to cardiac glycoside toxins. They overcame the challenges of retracing evolution by using CRISPR-Cas9 based genome editing, which allowed them to create and test combinations of the target mutations introduced into Drosophila melanogaster. They found strong positive epistasis for toxin tolerance between the three substitutions that had been repeatedly observed during evolution of toxin resistance. In addition, the order at which substitutions could accumulate was found to be constrained by epistasis (Karageorgi et al., 2019). Positive epistasis at fitness peaks is a general theoretical prediction that holds independent of the fitness landscape model (Greene & Crona, 2014).

In the research area of speciation, epistasis plays a major role when it manifests itself as hybrid incompatibility (reviewed in Seehausen et al., 2014, Coughlan & Matute, 2020). The classical model that explains speciation via the evolution of postzygotic hybrid incompatibilities in geographically isolated populations is the (Bateson-)Dobzhansky-Muller model. In this model, allopatrically diverging populations accumulate substitutions over time. Any combinations of substitutions not tested by natural selection because they never segregated in the same population have the potential to be incompatible, leading to so-called Dobzhanzky-Muller incompatibilities (DMIs). As pointed out by Gavrilets (1997), the resulting fitness landscape contains ridges of connected high-fitness genotypes and valleys of incompatible genotypes reflecting negative epistasis between the outermost genotypes on the ridge. Thus, a defining element of a classical DMI is that only one of the recombinant offspring of hybrids suffers from low fitness, whereas the broad definition of hybrid incompatibility includes any combination of alleles reducing fitness in hybrids (Seehausen et al., 2014).

Hybrid incompatibilities and DMIs between species have been mapped in laboratory crosses and more recently using large genomic screens, often revealing hundreds of pairwise incompatibilities between species pairs. For example, a genetic map of incompatibilities between the two fly species D. melanogaster and D. simulans resulted in an estimate of $>$100 hybrid-lethal incompatibility regions (Presgraves, 2003). Similarly, between two hybridizing swordtail fish in Mexico, genomic analyses of hybrid genomes indicated $>$100 incompatibility loci (Schumer et al., 2014). Also within species, segregating epistatic interactions may be common: Corbett-Detig et al. (2013) found various pairwise epistatic interactions when screening the global genetic variation within D. melanogaster.

Measuring a complex fitness landscape in the wild that encompasses more than pairwise interactions is difficult for various reasons. For example, diploidy creates a challenge for defining fitness, it is often difficult to obtain a sufficiently large sample to have a good representation of genotypes, fitness measurements are even less straightforward than in the lab, background variation might affect measurements, or experimental manipulation and creation of fitness variation can be impossible. Hybrid populations circumvent several of these challenges by being natural sources of genetic and fitness variation. Few studies to date have attempted to leverage this opportunity to map multi-locus fitness landscapes in the wild. In a landmark study, Martin & Wainwright (2013) mapped the phenotype-fitness landscape of hybrids of three recently diverged pupfish ecotypes in field enclosures. Based on measures of growth and survival, they found a multi-peaked phenotype-fitness landscape isolating molluscivores from generalists with no evidence of frequency-dependent selection (Martin & Wainwright, 2013, Martin & Gould, 2020, Martin, 2016). In a follow-up study, Patton et al. (2021) recently studied the genotype-fitness landscape underlying the observed fitness differences between specialist and generalist phenotypes.

Epistasis enables multiple genetic solutions to adaptation, represented by multiple peaks in the fitness landscape. Consequently, epistasis is thought to favor allopatric divergence. However, it is an open question how much epistatic alleles contribute to speciation early during divergence or in the presence of gene flow, and how (much) epistatic variation is accumulated and maintained within populations. Barton & Bengtsson (1986) argued that hybrid incompatibilities tend to be weak barriers to gene flow. Evolutionary theory has predicted that, because of the specific shape of their fitness landscape, epistatic alleles in a DMI cannot be maintained polymorphic in a population and that their accumulation in the presence of gene flow can only be a by-product, rather than a driver, of evolution (Gavrilets, 1997, Bank et al., 2012). However, more complex epistatic interactions might behave differently. Blanckaert et al. (2020) recently showed that adding a third locus to the fitness landscape of a hybrid incompatibility allows for the evolution of complete reproductive isolation in the presence of gene flow (see Paixão et al., 2015 for a network-based approach). In a series of theoretical works, Fraïsse et al. (2016), Simon et al. (2018), Schneemann et al. (2020) have studied the evolution of hybrid incompatibility using FGM, recovering many of the predictions of previous theoretical models and patterns observed in empirical data. Furthermore, they developed indices characterizing the respective influence of the environment and genetics on hybrid fitness (Schneemann et al., 2020).

Taking a fitness-landscape view of speciation models in which higher-order interactions are possible is challenging due to the complexity of the models (and the potential genomic inference of such interactions), but may be an opportunity to resolve the dichotomy between the frequency at which hybrid incompatibilities are observed and the theoretically ascribed minor role of epistatic loci early during divergence (Seehausen et al., 2014). Given the potential overrepresentation of epistatic loci at short physical distances in the genome (e.g., within proteins, or regulatory interactions), it will be important to address the role of genetic architecture and recombination distances in determining the evolutionary dynamics of epistatic loci (see, e.g., Blanckaert & Bank, 2018). Here, probabilistic fitness landscape models can be a tool to quantify how the amount of epistasis in a fitness landscape affects speciation. For example, with increasing ruggedness, the number of fitness peaks increases and with it the potential for populations to climb different fitness peaks. At the same time these fitness peaks will be located close to each other in sequence space, which should result in a smaller barrier to gene flow when the diverged populations hybridize (see also Figure S2). Thus, an intuitive prediction would be that intermediate epistasis maximizes the potential for stable divergence.

Despite great advances in the study of both theoretical and empirical fitness landscapes and epistasis, little is known about how fitness landscapes and epistasis change across environments. Whereas experimental screens have reported several cases of ubiquitous genotype-by-environment interactions for fitness (GxE; also referred to as antagonistic pleiotropy and costs of adaptation; e.g., Bank et al., 2014, Kinsler et al., 2020, Flynn et al., 2020), the presence and magnitude of environment-dependent epistasis (or genotype-by-genotype-by-environment interactions; GxGxE) for fitness has received little attention to date. The few existing examples have reported striking variation in epistasis between environments. When Hall et al. (2019) competed the 32 E.coli genotypes from Khan et al. (2011) across eight different environments, the resulting fitness landscapes looked vastly different between environments; epistasis was shown to be both strong and environment-specific (Figure 3; see also Flynn et al., 2013). In an experimental study of yeast that had evolved in the presence of a fungicide, epistatic interactions between two point mutations could change greatly dependent on the concentration of the drug (Ono et al., 2017). Gorter et al. (2018) showed that the evolution of yeast populations in the presence of increasing concentrations of nickel followed the evolutionary trajectory predicted by the changes in the fitness landscape along the dosage gradient. Lindsey et al. (2013) showed that environment-dependent epistasis prohibited evolutionary rescue of E. coli populations under high rates of environmental change. In a high-throughput study of a yeast tRNA fitness landscape across four environments, Li & Zhang (2018) also reported widespread environment-dependence of epistasis. Finally, Bakerlee et al. (2022) recently reported large environmental variation with respect to both single-mutation effects and pairwise epistasis of a fitness landscape spanning 10 mutations in different genes in yeast, measured in six environments. Based on their analysis the authors argued that epistasis in their fitness landscapes arose from specific gene interactions rather than a global, non-specific, non-linear phenotype-fitness relationship.

Environment-dependent epistasis was also observed in studies of insects and fish. Nosil et al. (2020) mapped the fitness landscape related to cryptic coloration seen in stick insects using a field experiment. Here, the fitness proxy was survival during a two-day exposure to predators in the wild. In a subset of F2 offspring from a sympatric pair of threespine stickleback populations that otherwise showed little evidence of epistatic interactions, Arnegard et al. (2014) identified an environment-dependent functional mismatch of oral jaw traits contributing to a reduction in feeding performance. Interestingly, both studies pointed to a major role of environment-dependent changes in phenotype-fitness relationships in creating and altering epistasis.

Most fitness landscape theory to date has addressed a single, static environment, and thus a static fitness landscape. However, the study of evolution in changing environments is important both from a basic science point of view and with respect to the response of populations to anthropogenic change and the climate crisis. Fluctuating phenotypic optima turn a fitness landscape into a “seascape” (Mustonen & Lässig, 2009). Although the concept of the fitness seascape is gaining increasing interest from the scientific community, little work has addressed how to extend the existing fitness landscape theory to account for changing environments and how this affects epistasis and its consequences.

Fisher’s Geometric Model is the most extensively studied class of fitness landscape models with respect to changing environments (Harmand et al., 2017, Martin & Lenormand, 2015, 2006, Zhang, 2012, Matuszewski et al., 2014, Schneemann et al., 2020). Different environments can be incorporated in FGM by assuming a change in the location of the phenotypic optimum (Martin & Lenormand, 2015, 2006, see also Figure 2B). Under this assumption, it is possible to study environmental gradients as well as completely different environments, and to develop expectations of the amount of genotype-by-environment interactions. Using experimental data from E. coli, evolved at different concentrations of an antibiotic drug, Harmand et al. (2017) could predict the costs of antibiotic resistance across environmental gradients with a generalized FGM. From both a theoretical and experimental point of view, an exciting next step would be to extract the how the underlying genotype-fitness landscape changed across environments.

Several studies of dynamic genotype-fitness landscape models have addressed how populations adapt to fluctuating environments, with a focus on the response to the time-scale at which fluctuations occur and whether evolution results in the repeated switching between specialists or the evolution of generalists (Sachdeva et al., 2020, Wang & Dai, 2019, Trubenová et al., 2019). Moreover, Agarwala & Fisher (2019) showed that in a slowly changing environment, epistasis promotes evolutionary dynamics that switch between short phases of rapid adaptation and long periods of stasis in which the population waits for new adaptive paths to open up.

Recent theoretical work by Das et al. (2020) has shown how simple genotype-phenotype relationships can result in environment-dependent epistasis. Their model specifically described the genotype-fitness landscapes along an environmental gradient of antibiotic concentration, based on Hill functions that are commonly used to describe the S-shaped relationship between bacterial growth rate and antibiotic concentration (so-called pharmacodynamics (Goutelle et al., 2008)). Das et al. (2020) assumed that every genotype has a different S-shaped fitness curve, and that there is a trade-off between growing well in the absence and the presence of the antibiotic; i.e. resistant genotypes must have a cost of resistance. Their model was in excellent agreement with a previously described experimental fitness landscape from E. coli (Marcusson et al., 2009). Interestingly, although Das et al. (2020) observed little epistasis for the two main phenotypic parameters of the model, the drug-free growth rate of genotypes and their IC50 (the concentration of the drug reducing the growth rate by 50%), epistasis was common in the composite fitness measure combining these two parameters and variable along the dose gradient. Epistasis was greatest at intermediate concentrations of the antibiotic and the fitness ranking of the genotypes changed most rapidly at intermediate drug concentrations. The authors concluded that small changes in (intermediate) dosage might have large consequences for evolution, because new evolutionary routes may open up or close.

Extending their theoretical analysis of the model, Das et al. (2021) found that a genotype can carry information on its “ecological history”; in other words, the current genotype is informative of the previous environments. This scenario arises in fitness landscapes in which epistatic relationships change along an environmental gradient, such that evolutionary trajectories are contingent on the previously experienced environment. Firstly, some genotypes were only accessible at low doses of an antibiotic drug, and would therefore never be observed if the population experienced high doses of the drug in its recent history. Secondly, there was on average limited reversibility of evolutionary paths, i.e., if the environmental change was reversed the population did not tend to take the reverse evolutionary path. Notably, the obtained analytical results are based on specific models and approaches inspired by theoretical physics. It will be interesting to see if similar results can be obtained from other models incorporating environmental gradients, such as FGM.

Under the assumption of environmental fluctuations in natural environments, evolution may not only be constrained by epistasis, but this epistasis may also be variable across environments, further complicating the study of selective forces that acted during evolutionary history. So far, few models can accommodate the change of epistasis across environments, and there is little knowledge about whether these effects could “cancel out” integrated over evolutionary history, maybe resulting in, on average, neutral-like evolutionary dynamics.

At the same time, selection coefficients inferred from empirical fitness landscape studies tend to be much larger than those estimated via inference of the distribution of fitness effects from polymorphism data (reviewed in Tataru & Bataillon, 2020). This could be due to the complexity of natural environments as compared to experimental settings. Over evolutionarily relevant time scales (i.e., those reflecting the polymorphism used for the above-mentioned inference methods), the experienced environment of a population has likely fluctuated greatly (including the effect of ecological interactions of the population with its environment), such that any specific strong fitness effect is diluted as an average over time and number of experienced environments (and fitness-related phenotypes). This view would reflect the ”seascape” approach to fitness landscapes (Mustonen & Lässig, 2009). A laboratory environment removes as many of these fluctuations and complexities as possible to extract the fitness (or phenotypic) effect of mutations in a well-defined environment. Thus, one main target of selection in the laboratory may be to streamline molecular processes and to remove all unnecessary mechanisms allowing an organism to respond to environmental fluctuations, e.g., by knocking out genes with premature stop codons (see, e.g. Kvitek & Sherlock, 2013, Kryazhimskiy et al., 2014). In addition, fluctuations could result both from a changing environment altering the fitness of genotypes and from population size changes altering the strength of genetic drift, thereby making fitness differences ”invisible”. Just as Wright (1932) argued, population bottlenecks could effectively smoothen the shape of the fitness landscape perceived by an evolving population in nature.

In the extreme case, all but the strongest epistasis could become invisible to natural selection and the fitness landscape could effectively smoothen to an almost flat surface with holes of strongly deleterious, incompatible genotypes, as described in Gavrilets’ holey fitness landscape model (Gavrilets, 2004). Relatedly, the study of combinatorially complete fitness landscapes might be misleading our intuition about evolution on these landscapes. The epistasis experienced by an evolving population based on its segregating variation is likely different from the observed epistasis when considering the entire possible fitness landscape (Whitlock et al., 1995). Quantifying these differences and finding measures of “realized epistasis” is a challenge because the segregating variation strongly depends on evolutionary parameters such as ploidy, recombination, population size and mutation rate.

It is unknown which part of the genotype-phenotype-fitness relationship is most affected by changes in the environment, and thus, what causes changes in epistasis across environments (see Figure 4B-D). Empirical and theoretical studies could address this question by developing models allowing for environment-tunable epistasis on the genotype-phenotype and the phenotype-fitness levels, and by designing approaches to extract those respective pieces of the fitness landscape. As Figure 4B-C show, a change in the phenotype-fitness mapping straightforwardly leads to changes in the genotype-fitness landscape and thus, changes in epistasis across environments. This can, for example, occur when a previously fitness-relevant phenotype becomes irrelevant in a different environment (e.g., camouflage in a dark or turbid environment, fast swimming in the absence of a predator, or antibiotic resistance in an antibiotic-free environment; Figure 4C). On the other hand, the environment can change epistasis in the genotype-phenotype map through environmentally induced (potentially transgenerational) epigenetic modification without changing the phenotype-fitness map, for example by silencing a gene involved in an epistatic interaction (Figure 4D).

Epistasis constrains evolutionary trajectories. If epistasis changes across environments, this implies that genomes in the present may contain information on the environment(s) the population has experienced in the past. Thus, the inference of evolutionary history from genomic samples may be linked to (or confounded by) the ecological history of a population (Das et al., 2021). Whether and in which conditions this concept is applicable to more general fitness landscapes, and whether it can be translated into an inference framework for genomic data, is currently unknown. Moreover, it is unknown how the contingency of evolutionary trajectories caused by epistasis (in stable or fluctuating environments), and the frequency-dependent selection caused by epistatic selection in admixed populations, may confound population-genetic inference of evolutionary histories or selection. As a first step and proof-of-concept to address this question, features of probabilistic fitness landscapes with widespread epistasis could be incorporated into genomic simulations to quantify how the presence of epistasis changes subsequent population-genetic inference.

Finally, epistasis can span all biological levels of organization: within or between genes, between host and symbiont, and possibly between species (see, e.g., Sørensen et al., 2021, Buskirk et al., 2020). Incorporating these levels into a theory of co-evolutionary fitness landscape is an exciting challenge. One theoretical approach could be to concatenate the genomes of interacting species into one “meta-genome” as a node in the fitness landscape – but how would such a model reflect the potentially different generation times and population sizes of the interacting species? Also within species, competition for resources could lead to a dynamic fitness landscape in which the fitness of one genotype in a population depends on the frequency of the others. For example, epistasis on the genotype-trait level is likely to affect the long-term diversity of the population as it propagates into the (frequency-dependent) genotype-fitness map.