Epistasis and the Structure of Fitness Landscapes: Are Experimental Fitness Landscapes Compatible with Fisher's Geometric Model?

The fitness landscape defines the relationship between genotypes and fitness in a given environment and underlies fundamental quantities such as the distribution of selection coefficient and the magnitude and type of epistasis. A better understanding of variation in landscape structure across species and environments is thus necessary to understand and predict how populations will adapt. An increasing number of experiments investigate the properties of fitness landscapes by identifying mutations, constructing genotypes with combinations of these mutations, and measuring the fitness of these genotypes. Yet these empirical landscapes represent a very small sample of the vast space of all possible genotypes, and this sample is often biased by the protocol used to identify mutations. Here we develop a rigorous statistical framework based on Approximate Bayesian Computation to address these concerns and use this flexible framework to fit a broad class of phenotypic fitness models (including Fisher's model) to 26 empirical landscapes representing nine diverse biological systems. Despite uncertainty owing to the small size of most published empirical landscapes, the inferred landscapes have similar structure in similar biological systems. Surprisingly, goodness-of-fit tests reveal that this class of phenotypic models, which has been successful so far in interpreting experimental data, is a plausible in only three of nine biological systems. More precisely, although Fisher's model was able to explain several statistical properties of the landscapes-including the mean and SD of selection and epistasis coefficients-it was often unable to explain the full structure of fitness landscapes.

THE fitness landscape is defined by a set of genotypes, the mutational distance between them, and their associated fitness in a given environment (Wright 1931; Orr 2005). The structure of the fitness landscape determines the fitness effects of mutations and the interaction between mutations for fitness. These properties determine the pace of adaptation (Eyre-Walker and Keightley 2007), the predictability of evolution (Weinreich ), the benefits of sexual reproduction (Kondrashov and Kondrashov 2001; de Visser ), and the probability of speciation (Gavrilets 2004; Chevin ). Thus, it is an important goal of evolutionary biology to characterize experimentally the properties of fitness landscapes across species and environments (de Visser and Krug 2014).

The most straightforward and popular experimental approach to access the properties of the fitness landscape consists of identifying mutations, constructing several genotypes that differ only by various combinations of these mutations, and measuring the fitness of these genotypes. This protocol allows reconstruction of what we call “empirical landscapes.” For example, several experiments identify a small number L of mutations and consider the fitness of 2L genotypes with all possible combinations of these mutations. Early studies were primarily descriptive, with a focus on patterns of epistasis among mutations (Malcolm ; de Visser ; Whitlock and Bourguet 2000). In an influential study, Weinreich studied the landscape between an ancestral strain of Escherichia coli and an evolved type with five mutations conferring high antibiotic resistance. They computed the number of paths up to the fitness maximum that could be followed by a population evolving by natural selection and showed that the ruggedness of the landscape implied that very few mutational paths could be used during biological evolution. This study suggested that the structure of fitness landscapes might severely constrain evolutionary trajectories, thus opening up the possibility that adaptation could be predicted to some extent. This finding has inspired the characterization of many other empirical landscapes [reviewed in Weinreich ].

In principle, empirical landscapes can be compared with predictions from theoretical fitness landscape models. For example, several studies fit specific models to empirical landscapes (Lunzer ; Chou , 2014; Rokyta ; Schenk ). These models predict quantitatively the fitness values and epistasis coefficients and, as such, greatly improve our understanding of the form of epistasis that is typical of the particular system under study. However, the increasing number of empirical landscapes calls for a more general method to infer and compare the properties of fitness landscapes across species and environments. This possibility is very appealing and timely given that data accumulate on a diversity of empirical systems and selective environments, but it also raises several challenges.

The variability observed between empirical landscapes might be driven by biological differences of interest between organisms and environments of selection, but this variability is currently confounded with two other factors: stochastic variability due to sampling of a small number of mutations and variability in the protocol by which mutations are isolated. The full fitness landscape of a species in the environment of selection is defined as the fitness of all possible genotypes in that environment. This is an incredibly large space, scaling exponentially with the size of the genome. Most experiments explore a very small subset of the landscape because they examine at best a few dozen genotypes. Starting from the ancestral genotype, a single point in this large fitness landscape, the region of the fitness landscape that is explored depends on the particular mutations that were isolated. Thus, each empirical landscape results from a single realization of the stochastic sampling of a small number of mutations from a myriad of available mutations (Tenaillon ; Salverda ; Schenk ; Szendro ; de Visser and Krug 2014). In other words, a single constant underlying fitness landscapes can give rise to a diversity of small genotypic landscapes depending on the mutations that are sampled (Blanquart ). Moreover, the region of the underlying fitness landscapes that is explored depends on the experimental protocol used to isolate mutations. For example, mutations are often obtained under protocols involving natural selection. While random mutations give more rugged empirical landscapes, mutations that have been sequentially selected in a single population give smoother empirical landscapes (Draghi and Plotkin 2013; Szendro ; Blanquart ). Thus, inferring the properties of fitness landscapes from empirical data in meaningful ways requires (1) quantifying the uncertainty resulting from sampling of a limited number of mutations and (2) explicitly modeling how mutations were experimentally isolated.

In this study, we address these challenges and develop a statistical framework to infer the properties of the underlying fitness landscape from empirical landscapes. We use a broad class of phenotypic fitness landscape models that includes Fisher’s geometric model (Fisher 2000). Phenotypic fitness landscapes model how the genotype of an organism translates into a set of phenotypes, which themselves determine fitness. In other words, the very large space of all possible genotypes is projected onto a continuous phenotypic space of arbitrary dimensionality, and fitness depends only on the position in this phenotypic space. Fisher’s model, in particular, assumes that the phenotypes are under stabilizing selection toward a single optimum, that the effects of mutations in the phenotypic space are drawn from a multivariate Gaussian distribution, and that mutations combine additively in the phenotypic space. Phenotypes can be biological traits that need to be tuned to a precise level to maximize growth of the organism in the environment of selection, e.g., the concentration of an enzyme in a metabolic pathway or the level of expression of a gene. Fisher’s model also can be viewed as an abstract statistical description of the genotype fitness map.

A number of reasons motivate the choice of Fisher’s model as the underlying fitness landscape. A phenotypic model solves the problem of high dimensionality of the genotypic space. Indeed, genotypic fitness landscape models such as the rough Mount Fuji model (Szendro ) or the NK model (Kauffman and Levin 1987) require a number of parameters increasing linearly with the number of mutations or the number of genotypes. In contrast, a phenotypic model can describe an arbitrarily large number of genotypes using a small number of parameters. More fundamentally, it has been shown recently that Fisher’s model emerges from a set of “first principles” that specifies how fitness results from developmental integration of a large number of mutable traits (Martin 2014). Last, Fisher’s geometric model is simple yet can generate a diversity of empirical landscapes (Blanquart ), and it successfully predicts experimental quantities, such as the distribution of epistasis coefficient between pairs of mutations (Martin and Lenormand 2006; Martin ) and the dynamics of mean fitness over time (Perfeito ).

This study focuses on the following questions: How much information on the structure of the underlying fitness landscape can be inferred from existing empirical landscapes? What properties of fitness landscapes can be inferred from empirical data available so far, and are underlying landscapes similar in similar species or environments? Is the structure of empirical landscapes compatible with a model assuming stabilizing selection on a set of underlying unknown phenotypes?

To answer these questions, we developed an inference framework that allows fitting Fisher’s model to a diversity of experimental data sets obtained under a range of protocols. Using this framework, we infer the parameters and quantify the goodness of fit of Fisher’s model on 26 published genotypic landscapes representing nine distinct biological systems. We infer the properties of the underlying fitness landscape of each data set while accounting for the protocol used to obtain the data, allowing a meaningful comparison of fitness landscapes across several species and environments. This survey reveals substantial differences in the shapes of underlying fitness landscapes across biological systems and environments of selection. We also show that Fisher’s model is able to fully account for the observed properties of genotypic landscapes in only three of nine biological systems.

Materials and Methods

Data set selection

We searched the literature for published empirical landscapes that include clearly identified sets of genotypes with combinations of two mutations or more together with their fitness. The way in which these mutations evolved or were obtained had to be sufficiently described such that we could reproduce it with simulations (see below). For selected mutations, we verified that the fitness measure reported is relevant to the environment in which the mutations evolved. We identified a total of 26 published data sets spanning nine independent biological systems meeting these criteria. In the following, we will identify the data sets representing these nine systems using the letters A–I (Table 1 and Supplemental Material, File S1). The data sets encompassed a diversity of species– including species of virus, bacteria, fungi, animals– and ecological scenarios (Table 1). Several experiments explored the fitness landscape of species in a laboratory environment using random mutations in the fungus Aspergillus niger (de Visser ) (data sets A1 and A2), the fruit fly Drosophila melanogaster (data sets C1 and C2) (Whitlock and Bourguet 2000), and the budding yeast Saccharomyces cerevisiae (data sets B1–B10) (Costanzo ). The latter data set is a large collection of 5596 deletion mutants. To reduce this large data set to a size that was amenable to our analysis, we randomly drew 10 independent, randomly chosen subsets that included 20 mutations, all single mutants, and 100 double mutants (all combinations of the first 10 mutations times the last 10 mutations, for a total of 121 possible genotypes but, in reality, 104 to 116 genotypes because some were missing).

Table 1

Summary of data sets

Name	Species	Environment	Mutation type	Mutation and genotype number	Fitness measure	Measurement error	Note	References
A1, A2	A.niger	Minimal medium	Random (phenotypic markers)	Two data sets of 5 mutations, 25 genotypes	Rate of increase in colony radius per unit time, relative to the ancestor	We calculated a single standard error using the two replicate measurements	The two sets of five mutations are not independent	de Visser et al. (1997, 2009)
B1–B10	S. cerevisiae	Standard medium (on plates)	Random (gene deletions)	1711 + 3885 mutations, 5.4 million genotypes	Increase in colony size per unit time relative to the ancestor	Standard error reported for each fitness measure	Analysis done on 10 independent subsets of 20 mutations. 1 + 20 + 100 genotypes corresponding to single and double mutants	Costanzo et al. (2010)
C1, C2	D. melanogaster	Laboratory environment.	Random (phenotypic markers)	Five mutations, 25 genotypes	Productivity (product of fecundity and survival) C1 and mating success C2	We roughly estimated standard error for the productivity measure (File S1)	Two genotypes with 0 mating success removed	Whitlock and Bourguet (2000)
D	Single-stranded DNA bacteriophage ID11	E. coli (host)	Independently selected (experimental evolution)	Nine mutations, nine single mutants, 18 double mutants	Log2 increase in phage population per hour	Standard error reported for each fitness measure	—	Rokyta et al. (2011)
E1, E2	Vesicular stomatitis virus	Baby hamster kidney (BHK21) cells (host)	Independently selected (E1, found in natural isolates) and random (E2)	Six mutations, six single mutants, 15 double mutants (E1)	Growth rate relative to ancestor	Standard error reported for each fitness measure	—	Sanjuán et al. (2004)
				28 mutations, 76 double mutants (E2)
F	E. coli	New, low-glucose environment	Coselected in experimental evolution	Five mutations, 25 genotypes	Growth rate relative to ancestor	95% CIs reported	—	Khan et al. (2011)
G	Methylobacterium extorquens	Methanol environment	Coselected in experimental evolution	Four mutations, 24 genotypes	Growth rate relative to ancestor	Standard errors not reported, but we estimated them using standard errors reported for another set of genotypes (File S1)	Order of fixation of mutations not known; mutations were assumed to fix from the largest-effect mutation to the smallest-effect mutation	(Chou et al. 2011)
H1, H2	E. coli	Cefotaxime (β-lactam antibiotic)	Independently selected (found in natural isolates); mutations chosen because together increase resistance to cefotaxime 100,000-fold	Five mutations, 25 genotypes	Cefotaxime resistance measured as MIC	For H1, we calculated single errors using the three replicate measurements; for H2, standard errors were reported	H2 is the same data set as H1, with MIC remeasured on same genotypes; resistance to piperacillin + clavulanic acid was also measured but not used here	Weinreich et al. (2006); Tan et al. (2011)
H3, H4	E. coli	Cefotaxime (β-lactam antibiotic)	Independently selected (found in natural isolates)	four mutations, 24 genotypes; two independent data sets	Cefotaxime resistance measured as IC99.99 (highly correlated with MIC)	Standard errors not reported, but we used the average standard error of H2	One data set with four mutations of smallest effect H3, one with four mutations of largest effect H4 among 48 mutations	Schenk et al. (2013)
I1, I2, I3	Plasmodium falciparum DHFR gene transformed into E. coli I1 and S. cerevisiae I2; Plasmodium vivax DHFR gene transformed into S. cerevisiae (I3)	Pyrimethamine (antimalarial drug)	Independently selected (found in clinical isolates)	four mutations, 24 genotypes; I2 includes the same four mutations as I1 plus two additional mutations affecting another locus	Pyrimethamine resistance measured as IC50 in μg/ml (I1) and M (mol/L) (I2); growth rate of the transformed strain at concentration 1 μmol/L (I3)	Standard errors for I1 and I2 reported	In I2 (Brown et al. 2010), one locus has three possible alleles; the third allele was ignored, resulting in 17 fitness values	Lozovsky et al. (2009); Brown et al. (2010); Jiang et al. (2013)

Three data sets represented the fitness landscape of virus species adapting to their hosts (data sets D, E1, and E2) (Sanjuán ; Rokyta ). Two data sets represented landscapes of adaptation of microbial species to a novel environment, including a long-term selection experiment in a low-glucose environment (data set F) (Khan ) and a selection experiment in a methanol environment (data set G) (Chou ). Last, seven data sets represented empirical landscapes reconstructed from mutations that confer drug resistance. These included studies of mutations in the enzyme TEM-1 β-lactamase, which confer resistance to cefotaxime in bacteria (four data sets H1–H4) (Weinreich ; Tan ; Schenk ), and studies of mutations in the dihydrofolate reductase gene, which confer pyrimethamine resistance (an antimalarial drug) in transgenic bacteria and yeast (three data sets I1–I3) (Lozovsky ; Brown ; Jiang ).

Data analysis

A number of fitness measures were reported in the published empirical landscapes we collected. Our analysis requires meaningful estimates of fitness value to model how selected mutations differ from random mutations.

Meaningful selection coefficients are expressed in units of log-fitness. They must be calculated either as , where and are the multiplicative growth rate of the mutant and the ancestor (called “fitness” in most population genetics model), or as , where and are exponential growth rates (Chevin 2011). Unfortunately, many studies only reported the ratio (Table 1, data sets A, B, E, F, G, and I3), which in theory cannot be used to obtain a correct selection coefficient. To analyze the studies that only report , we used as a log-fitness measure. This measure is approximately equal to under weak selection, which is a quantity proportional to the selection coefficient. Moreover, this log-fitness measure, conveniently, does not depend on the unit of the growth rate and can be compared across landscapes.

For drug-resistance fitness landscapes, only one data set reported a growth rate at a given drug concentration (data set I3, Table 1). Other studies reported the minimum inhibitor concentration (MIC) or a similar measure (Table 1, data sets H1–H2, I1 and I2). MIC, the concentration of drug above which the population cannot grow, is not easily related to fitness. For this reason, we present the results of MIC landscapes in Figure S2 (Weinreich ; Tan ; Schenk ).

We proceeded to several additional steps of data cleaning. Three nonviable genotypes (fitness value of 0) were excluded from the analysis [one in a pyrimethamine landscape (I1) and two in a Drosophila landscape (C2)] because Fisher’s model cannot easily account for lethal mutations. In data set G, the order of fixation of coselected mutations was unknown. We assumed that mutations fixed from the largest-effect mutation to the smallest-effect mutation in accordance with the reported dynamics of mean fitness through time in the experiment. In data set I2, two mutations occurred at the same locus. We made this data set compatible with our framework (which assumes that each locus is diallelic) by excluding all genotypes bearing the third allele.

Approximate Bayesian computation

Table 1 shows that a number of protocols were used to obtain empirical landscapes. Some of the empirical landscapes were formed of single and double mutants only, while others included all possible combinations of four or five mutations, thus including genotypes with three, four, or five mutations. Moreover, the way in which mutations were isolated also varied. Mutations were random, independently selected, or coselected. “Independently selected” means that the mutations emerged under the action of selection in separate populations evolving independently from a unique ancestral genotype. “Co-selected” means that the mutations were selected sequentially in the same population. Modeling the way selection biased the resulting empirical landscape is already complicated. To make matters worse, several protocols included an additional step. These protocols were used to study the landscape of resistance to cefotaxime, a β-lactam antibiotic (landscapes H1–H4). Among a large set of 48 mutations found individually in cefotaxime-resistant natural isolates, three smaller subsets were studied in detail. These subsets were composed of the four mutations of smallest fitness effect, the four mutations of largest fitness effect (H3 and H4), and five mutations that together conferred a very high fitness (H1–H2). To account for this variety of protocols, we used a flexible approximate Bayesian computation (hereafter ABC) approach to infer from empirical data the parameters underlying Fisher’s geometric model.

Details of the ABC framework:

The original ABC rejection algorithm proceeds as follows: a large number of parameter sets are drawn in a prior distribution. For each parameter set θ, a data set is simulated, and a measure of distance between the true data set and each simulation is computed. A set of parameters is retained in the posterior distribution if the distance between D and is lower than a small value ε. In other words, the posterior distribution is composed of all the parameter sets θ such that . In practice, ε is chosen such that a given, small fraction of the prior parameter sets is retained in the posterior (Csilléry ), but ABC will give the correct posterior distribution of parameters only in the limit where ε is close to zero.

The distance between the data set and simulation is often defined based on a set of statistics. This set of statistics must be carefully chosen to be informative but of relatively low dimensionality. We conducted the analysis using either the full set of observed log-fitness values (16–121 fitness values) or a set of six summary statistics. The six summary statistics are as follows: (1) the mean coefficient of selection of all single mutants, (2) the mean epistasis coefficient between all pairs of mutations averaged over all genetic backgrounds, (3) the SD of selection coefficients, (4) the SD of epistasis, (5) the correlation between the epistasis coefficient and the background fitness (specifically, for each pair of mutations, we calculate the epistasis coefficient and the average fitness of the two genotypes with one of the mutations and compute the correlation between these two quantities across all pairs of mutations and all genetic backgrounds), and (6) the maximal fitness value (Table S1). The distance of each simulated data set to the experimental data set waswhere is the number of statistics, is the statistic i, and is the simulated statistic i. Statistics are normalized by the median absolute deviation , which is analogous to SD but with medians instead of means. When statistics were the full set of fitness values, genotypes were uniquely identified by ordering mutations by their fitness effects.

We detailed earlier the rejection algorithm, where the posterior is simply the fraction of parameters randomly drawn from the prior distribution that generates simulated landscapes closest to the data. For this algorithm, we used a tolerance (the fraction of retained simulations) of 0.005 (using the lower tolerance of 0.0005 did not improve accuracy, Figure S1). In addition to the rejection algorithm, we used a linear regression algorithm (Beaumont ). In this method, the posterior parameters are corrected using a local linear regression of the parameter values onto the summary statistics, giving more weight to simulations closer to the data set. Last, we used a neural-network algorithm that adjusts the posterior distribution based on a nonlinear regression using neural networks (Blum and François 2010). The three methods are implemented in the R package “abc” ( R Development Core Team 2010; Csilléry ).

Details of the evolutionary simulations:

We simulated a large number of genotypic landscapes under Fisher’s model, seeding the simulation with parameters θ drawn from some prior distributions (detailed later). The simulated landscapes were based on Fisher’s model, a phenotypic fitness landscape model whereby an organism is evolving under stabilizing selection on n continuous phenotypic traits that together determine fitness. Each genotype is characterized by a phenotype vector consisting of n traits, where n is the dimensionality of the phenotypic space. The parameter n defines the number of phenotypes under selection, or “complexity,” for an organism evolving in a given environment (Tenaillon ; Lourenço ; Chevin ). The effects of mutations are assumed to be additive in the phenotypic space. For example, if we consider five mutations at five distinct loci of the genome, the genotype 00101, where the series of 0s and 1s denote the absence or presence of mutations at each of five loci (relative to an ancestral strain with genotype 00000), has phenotype , where is the phenotype vector of the ancestral strain, and and are the phenotypic effects at mutations at loci 3 and 5. The effects of mutations on phenotypes (the vectors ) are drawn from a multivariate normal distribution with mean 0 and variance-covariance matrix , where σ is the size of mutations. Thus, each mutation jointly affects all phenotypes (assumption of full pleiotropy). The mapping of phenotype on fitness is defined by , where is the maximal fitness, which determines the distance to the optimum of the ancestral strain, is the Euclidean norm of the phenotype vector, and e is the experimental error on fitness measurements. Following Wilke and Adami (2001) and others (Tenaillon ; Gros ), we extended Fisher’s geometric model with the parameter Q, which quantifies how flat the peak is at the optimum (Figure 1). Fisher’s model, sensu stricto is the special case where , i.e., the fitness function is Gaussian. Our definition of fitness implies that the ancestral strain had log-fitness 0, corresponding to the phenotype . This normalization was done without loss of generality. Maximum fitness , which is the height of the fitness peak in the environment where fitness is measured, was achieved when all phenotypes are at their optimal value, chosen here to be without loss of generality. Lastly, e is the measurement error in log-fitness measure and was assumed to be normally distributed with mean 0 and SD estimated from the empirical data (File S1). Figure 1 shows several examples of a single empirical genotypic landscape generated by sampling a small number of mutations in the underlying landscape.

Figure 1

A diversity of genotypic landscapes can be generated by Fisher’s fitness landscape model. Each row shows an example of Fisher’s landscape with two phenotypes (n = 2), with three mutations depicted as arrows in the phenotypic space (left) and the empirical landscape resulting from these mutations in combination (i.e., eight genotypes) (right). Blue edges denote mutations that are beneficial in their background, while red edges denote deleterious mutations. (Top row) A sharp landscape with Q = 0.5 and where the three mutations are random mutations. (Center row) Fisher’s classic landscape with Q = 2 and three coselected mutations. (Bottom row) Q = 4 and three independently selected mutations. Fitness of the ancestral strain is set to 1 without loss of generality.

For each set of parameters , we simulated the process by which mutations were isolated and generated a genotypic landscape. In practice, the sets of genotypes were of two broad categories: either four to five mutations were isolated and genotypes bearing all possible combinations of those mutations (24 or 25) were constructed or a larger number of mutations (seven to nine) were isolated and single and double mutants were constructed. Mutations were considered to be random, independently selected, or coselected. For random mutations, simulations consisted of drawing the phenotypic effects of mutations in the multivariate normal distribution and then combining those mutations additively and computing fitness using our phenotype-to-fitness mapping. When mutations were isolated in an experiment involving selection, we assumed that adaptation proceeded by successive invasion of beneficial mutations without clonal interference. This allowed us to conduct fast simulations based on the strong-selection, weak-mutation (SSWM) approximation (Kimura 1983; Gillespie 1991), making it possible to conduct the large number of simulations required by ABC. Under the SSWM regime, a selected mutation is drawn among the pool of random mutations, with each mutation weighted by , where s is the fitness effect of the mutation. This derives from the fact that the probability of fixation of a beneficial mutation is scaling linearly with its fitness effect s in this regime (Patwa and Wahl 2008). Fitness effects were calculated relative to the ancestor for independently selected mutations and relative to the genetic background with previously evolved mutations for coselected mutations. For the protocol where five mutations that together confer a large fitness effect are isolated (H1-H2), we chose the set of five mutations that confers the highest fitness among 1000 random combinations.

For each empirical landscape, 106 genotypic landscapes were generated using 106 parameter sets drawn from prior distributions. Priors were chosen to be uninformative and to ensure that they could generate a diversity of fitness landscapes (Figure 1). The height of the peak in log-fitness was drawn from an exponential distribution with mean 2. Maximum fitness on a log scale ranged from 3.7 × 10−7 to 29 (2.5–97.5% quantile 0.05–7.4). The complexity of the phenotypic space, the number of phenotypic dimensions under selection, was given by , where · denotes the floor function, and η was drawn from an exponential distribution with mean 5. It ranged from 1 to 75 (2.5–97.5% quantile 1–7). We used an exponential prior for complexity because, under Fisher’s model with full pleiotropy, the distribution of fitness effects had unrealistically small variance at high complexity. The size of mutations σ in the phenotypic space was drawn from an exponential distribution with mean 0.2. It ranged from 1.7 × 10−7 to 2.6 (2.5–97.5% quantile 0.005–0.74). The choice of an exponential distribution was motivated by the fact that variations in fitness are modest in many of the data sets, and therefore, mutational effects are probably small. The shape of the peak Q was drawn from a uniform distribution (Figure 1).

Cross-validation

We checked the accuracy of inference from empirical landscapes using simulated pseudo–data sets generated under Fisher’s model. We performed cross-validation using pseudo–data sets generated under Fisher’s model for each type of experimental protocol (Figure 2 and Table 2). We applied the ABC algorithm on each data set and compared the posterior distribution of parameters to the true (known) parameters. We computed the prediction error, defined for each parameter aswhere is the true value of the parameter used for the ith simulated pseudo–data set, is the median of the posterior distribution, and is the variance of the prior distribution. The expected prediction error is 0 when inference is perfect (the median always matches the true parameter) and 1 when no inference can be made (the posterior parameters are drawn at random from the prior). For cross-validation, we assumed that experimental errors were 0 in order to compare the accuracy of inference across protocols in an ideal case where fitness values are perfectly known.

Figure 2

Accuracy of inference for different methods and different data sets. The median posterior distribution for the rejection algorithm is shown as a function of the true parameter for each of the 500 cross-validation data sets (gray points) when the set of genotypes is composed of all combinations of four independently selected mutations, chosen as the four largest-effect mutations among a set of 48 mutations, as in landscape H4 (Schenk ). Perfect inference corresponds to all points on the y = x line. For clarity, we represent this cloud of points with a local nonlinear fit (gray line). The equivalent linear fit for the neural-network algorithm is shown as a gray dashed line. The plain and dashed blue line similarly show the local linear fit for rejection and neural-network algorithms for the data set composed of 20 random mutations and single and double mutants only (as in landscapes B1–B10). The neural-network algorithm generally improves inference compared to the rejection algorithm. The data set composed of all combinations of four selected mutations performs better than the one composed of 20 random mutations and single and double mutants.

Table 2 

Expected prediction error under various experimental designs

Experimental design	Type of mutation*	Landscapes using this protocol	n			Wmax			σ			Q
			rej	reg	nn	rej	reg	nn	rej	reg	nn	rej	reg	nn
5 mutations, 25 genotypes	R	A, C	0.85	0.68	0.64	0.57	0.39	0.37	0.44	0.35	0.32	0.67	0.49	0.43
	IS	-	0.91	0.79	0.7	0.34	0.2	0.18	0.33	0.19	0.19	0.35	0.15	0.09
	CS	F	0.83	0.73	0.63	0.34	0.19	0.17	0.53	0.37	0.33	0.39	0.24	0.17
4 mutations, 24 genotypes	R	-	0.93	0.8	0.78	0.79	0.58	0.54	0.42	0.36	0.35	0.64	0.52	0.45
	IS	I	0.87	0.76	0.67	0.41	0.22	0.18	0.43	0.25	0.23	0.48	0.27	0.2
	CS	G	0.9	0.76	0.69	0.37	0.2	0.17	0.49	0.33	0.3	0.42	0.33	0.24
8 mutations, 8 single and 20 double mutants	RS	-	0.8	0.58	0.54	0.68	0.48	0.5	0.37	0.3	0.29	0.55	0.44	0.4
	IS	-	0.75	0.69	0.63	0.44	0.29	0.22	0.41	0.29	0.28	0.47	0.33	0.29
	CS	-	0.77	0.67	0.62	0.4	0.21	0.18	0.43	0.26	0.23	0.27	0.23	0.2
20 mutations, up to 121 genotypes	R	B	0.72	0.62	0.54	0.48	0.25	0.21	0.35	0.23	0.22	0.62	0.45	0.39
9 mutations, 9 single mutants, 18 double mutants	IS	D	0.74	0.68	0.63	0.37	0.22	0.15	0.38	0.27	0.24	0.45	0.37	0.32
6 mutations, 6 single mutants, 15 double mutants	IS	E	0.81	0.8	0.76	0.45	0.25	0.18	0.39	0.33	0.3	0.67	0.59	0.54
5 mutations, 25 genotypes	IS, high fitness combination	H1, H2	0.86	0.74	0.61	0.24	0.09	0.06	0.42	0.25	0.22	0.24	0.1	0.05
4 mutations, 24 genotypes	IS, small fitness effect mutants	H3	0.8	0.84	0.78	0.69	0.52	0.43	0.5	0.41	0.38	0.69	0.48	0.36
4 mutations, 24 genotypes	IS, large fitness effect mutants	H4	0.94	0.72	0.56	0.26	0.1	0.07	0.25	0.15	0.14	0.27	0.1	0.04

Prediction error for the four parameters of Fisher’s model, for several experimental designs (based on single and double mutants, or complete sets of mutations and all associated genotypes) and selection procedures (* R: random, IS: independently selected, CS: co-selected mutations), when the 6 summary statistics were used in the ABC algorithm. For each parameter, the three lowest prediction errors are in bold, highlighting the protocol and inference algorithms that perform best.

Posterior predictive checks

We next tested whether the empirical landscapes we analyzed were compatible with the hypothesis that Fisher’s landscape was the true model for the empirical data. We used posterior predictive checking (Gelman ) to quantify the goodness of fit of Fisher’s model to each data set. For each experimental data set, we ran the ABC algorithm on 1000 random pseudo–data sets generated using parameters drawn from the joint posterior distribution of parameters. For each of these pseudo–data sets, we recorded the median distance between the pseudodata and the accepted (closest) simulated data in the ABC algorithm. This resulted in a null distribution for the median distance of the simulations retained in the ABC algorithm, which is the distribution of distance between simulations and data when Fisher’s model is truly underlying the data. We then used this distribution to compute a Bayesian P-value, also known as a “posterior predictive P-value” in Bayesian model checking (Gelman ). This P-value is the probability that median distances for pseudo–data sets generated under Fisher’s model are greater than the median distance of the experimental data set. A low P-value suggests that the data are farther apart from Fisher’s model simulations than expected if the data followed Fisher’s model. A P-value was computed for the distances based on summary statistics and for the distances based on all fitness values. For the latter, we also decomposed the distance and computed an analogous P-value for each individual genotype to identify genotypes with fitness values that are particularly unlikely under Fisher’s model (those whose individual P-value is lower than 0.05).

Data availability

The authors state that all data necessary for confirming the conclusions presented in this article are represented fully within the article.

Results

Cross-validation and accuracy of parameter inference

We quantified the accuracy of inference from empirical landscapes using 500 simulated pseudo–data sets generated under Fisher’s model. This analysis revealed that the true parameters of the underlying landscape are generally inferred with mediocre accuracy under most protocols used in existing studies (Figure 2 and Table 2). Inference based on summary statistics (Table 2) always yielded lower error than inference based on all fitness values (Table S2). Using summary statistics makes the ABC algorithm more accurate because it alleviates the “curse of dimensionality”: the distance of the data to the accepted simulations is closer to 0 for the same number of simulations such that the main assumption of ABC is better respected. However, the use of summary statistics causes loss of information (Sünnaker ). Here the gain in accuracy more than offsets the loss of information, making inference based on summary statistics better.

ABC is an approximate method, and we cannot rule out totally that low accuracy was due to these approximations. However, low accuracy also may be caused by the limited information contained in small genotypic landscapes. In other words, even if the inference method were perfect, the true posterior distribution of parameters still might be quite wide and cause low accuracy. Because we have explored a number of variations on the ABC algorithm, including three different algorithms, full statistics vs. summary statistics, and several values of tolerance (Figure S1) and accuracy of inference was always relatively low, we hypothesize that the main reason behind low accuracy is probably the limited information contained in genotypic landscapes. Each empirical landscape conveys rather modest information on the underlying “true” fitness landscape.

In particular, empirical landscapes conveyed almost no information on the number n of phenotypes under selection. Prediction errors for this parameter were always higher than 0.5 and often close to 1. The size of mutations σ, the height of peak , and the shape of fitness peak Q were inferred with more accuracy. For all parameters, the regression and neural-network algorithms improved the accuracy of inference relative to the rejection algorithm, and the neural-network algorithm was most often the best (Table 2; compare the “rej,” “reg,” and “nn” columns for each parameter).

With the summary statistics we chose, the number of mutations that were combined together did not affect the quality of inference much. The experimental design with 32 genotypes made of all combinations of five mutations performed similarly to the one made of eight mutations and single and double mutants only (28 genotypes) (Table 2). The design where 20 random mutations were chosen (landscapes B1–B10) did not perform particularly better than the one with eight mutations and all single and double mutants (29 genotypes in total).

The protocol used to isolate mutations was of critical importance to the quality of inference (Figure 2 and Table 2). Generally, selected mutations allowed the most accurate inference (compare the “random,” “independently selected,” and “coselected” lines in Table 2 for a given experimental design). In these simulations, the protocol where the four largest mutations were isolated among 48 independently selected mutations performed best and allowed fairly precise inference of the size of mutations (error = 0.145), height of the peak (error = 0.068), and shape of the peak (error = 0.045) under the neural-network algorithm. Protocols that performed best regarding inference of the height and shape of the fitness peak allow a better exploration of the underlying fitness landscape around the fitness optimum. Independently selected mutations and particularly large-effect mutations create genotypes that are more likely to be around the fitness peak, especially when genotypes with more than two mutations are included. In contrast, genotypes constructed with random mutations do not always approach the fitness peak and may be confined to relatively linear and uninformative zones of the underlying fitness landscape.

Parameter inference in experimental data sets

We obtained the posterior distribution of fitness landscape parameters in the 26 data sets. We used the ABC protocol based on summary statistics and the neural-network algorithm, which was shown to work best (Table 2). Note that the neural-network algorithm, on rare occasions, resulted in parameter estimates with biologically meaningless values, e.g., negative values of dimensionality or maximal fitness. This is a known problem (Sünnaker ) that happens when none of the summary statistics are very close to the data such that the neural-network regression extrapolates and yields posterior values outside the range of the prior. Results are similar, but the posterior distributions are wider when using inference based on the full set of fitness values and/or the rejection algorithm.

First, as expected from cross-validation, the posterior distributions were broader for parameters describing dimensionality and shape of the peak (Figure 3 and Table 3). Each empirical landscape could have been generated under a diversity of underlying fitness landscapes. Despite the uncertainty in parameters, different biological systems exhibited different types of fitness landscapes (Figure 3).

Figure 3

Posterior distribution of parameters for all experimental landscapes. (From top to bottom) A1 andA2 (Aspergillus) and C1 and C2 (Drosophila); the yeast deletion data set (B1–B10); virus evolving on their host (D (circle) and E1-E2 (squares)) and bacteria in a novel medium (F and G); adaptation to an environment containing pyrimethamine (I3). The black point shows the median of the prior, and the dashed line delineates the 50% higher-density region. The points show the median of the posteriors, and the shaded areas show the 50% higher posterior density regions for the data sets.

Table 3

Posterior distribution of parameters and posterior predictive checks, neural-network algorithm

Reference	Name	n	Wmax	σ	Q	P-value (summary)	P-value (full)
—	Prior	4 (1; 19)	1.39 (0.05; 7.39)	0.14 (0.01; 0.74)	2.25 (0.59; 3.91)
de Visser et al. (1997)	A1	5.24 (0.59; 20.96)	0.14 (-0.07; 1.88)	0.15 (0.08; 0.37)	1.60 (0.42; 3.52)	0.17	0.83
	A2	6.72 (1.63; 23.09)	0.34 (-0.09; 3.12)	0.12 (0.05; 0.31)	1.69 (0.66; 3.51)	0.23	0.97
Costanzo et al. (2010)	B1	6.00 (1.54; 19.08)	1.01 (0.26; 3.88)	0.09 (0.04; 0.29)	1.91 (0.91; 3.75)	0.25	0
	B2	3.44 (0.31; 12.82)	0.28 (0.12; 1.02)	0.09 (0.05; 0.20)	2.96 (1.64; 4.19)	0.3	0.02
	B3	3.33 (0.13; 13.78)	0.40 (0.13; 1.78)	0.10 (0.06; 0.23)	2.33 (1.19; 4.36)	0.31	0.01
	B4	8.28 (1.64; 24.06)	1.21 (0.04; 5.10)	0.14 (0.03; 0.48)	1.57 (0.77; 3.13)	0.16	0.06
	B5	4.16 (0.89; 14.84)	0.43 (0.07; 2.37)	0.08 (0.05; 0.23)	2.17 (1.02; 4.37)	0.26	0.02
	B6	3.01 (-0.72; 13.87)	0.34 (-0.05; 1.96)	0.10 (0.05; 0.29)	2.23 (1.16; 4.68)	0.24	0.01
	B7	4.47 (-0.25; 15.71)	0.64 (0.12; 2.23)	0.11 (0.05; 0.33)	2.07 (0.97; 4.12)	0.14	0.02
	B8	1.63 (-1.91; 12.29)	1.32 (0.32; 4.92)	0.11 (0.00; 0.48)	2.24 (1.10; 4.17)	0.02	0.01
	B9	4.12 (0.48; 15.58)	0.39 (0.02; 2.38)	0.09 (0.04; 0.26)	2.12 (1.07; 4.26)	0.17	0
	B10	3.46 (0.63; 15.18)	0.32 (0.07; 1.58)	0.07 (0.04; 0.20)	2.34 (1.08; 4.35)	0.06	0.01
Whitlock and Bourguet (2000)	C1	4.92 (2.12; 13.24)	1.02 (0.58; 3.20)	0.30 (0.16; 0.66)	2.98 (1.71; 4.06)	0.05	0
	C2	2.09 (0.21; 7.03)	1.10 (0.82; 2.38)	0.57 (0.39; 1.03)	2.58 (1.01; 3.57)	0.06	0
Rokyta et al. (2011)	D	7.00 (2.95; 15.21)	0.46 (0.36; 0.82)	0.21 (0.15; 0.39)	2.08 (0.83; 3.82)	0.15	0.08
Sanjuán et al. (2004)	E1	6.28 (1.64; 19.82)	0.19 (0.06; 0.86)	0.15 (0.07; 0.41)	1.65 (0.23; 3.79)	0.37	0.01
	E2	5.28 (2.11; 12.45)	0.20 (0.09; 0.55)	0.14 (0.10; 0.25)	2.26 (1.34; 3.42)	0.16	0.03
Khan et al. (2011)	F	6.62 (1.63; 22.28)	0.42 (0.21; 0.98)	0.08 (0.05; 0.19)	1.89 (0.81; 3.70)	0.43	0.03
Chou et al. (2011)	G	3.65 (0.86; 15.86)	1.09 (0.73; 2.48)	0.07 (0.03; 0.21)	2.67 (1.30; 4.05)	0.42	0.43
Weinreich et al. (2006)	H1	14.39 (7.25; 29.54)	12.97 (12.16; 15.73)	0.89 (0.64; 1.46)	1.40 (0.14; 2.48)	0.01	0
Tan et al. (2011)	H2	13.18 (5.76; 28.86)	12.02 (10.87; 14.83)	0.46 (0.18; 1.08)	1.83 (0.81; 2.80)	0.01	0
Schenk et al. (2013)	H3	4.81 (1.89; 15.30)	3.17 (1.08; 8.91)	0.30 (0.13; 0.79)	2.94 (1.53; 3.91)	0.34	0.07
	H4	8.89 (5.63; 17.44)	6.24 (5.27; 7.94)	0.75 (0.51; 1.13)	1.40 (0.62; 2.15)	0	0
Lozovsky et al. (2009)	I1	8.24 (3.79; 19.68)	9.20 (7.78; 14.61)	0.57 (0.26; 1.24)	2.22 (0.55; 3.51)	0.02	0
Brown et al. (2010)	I2	5.16 (2.50; 13.08)	7.76 (7.41; 8.95)	0.23 (0.15; 0.37)	3.84 (3.17; 4.49)	0	0
Jiang et al. (2013)	I3	1.28 (-0.58; 5.47)	2.33 (2.19; 2.71)	0.47 (0.32; 0.79)	3.70 (3.11; 4.24)	0.12	0.03

The median posterior distribution of parameters and the 2.5–97.5% quantile interval (equivalent to 95% higher posterior density) of the posterior distribution of parameters for the rejection algorithm. The prior is shown for comparison (first row). The P-value for the test of adequacy with Fisher’s model is indicated.

Three of the experimental systems that were represented by several nonindependent empirical landscapes resulted in similar posterior distributions across these “replicated” landscapes. This demonstrates the robustness of the ABC method to slight variation in the set of mutations, to variation in the fitness measure, and to experimental error. For A. niger (Figure 3, first row), two empirical landscapes, A1 and A2, were constructed using two partially overlapping sets of mutations (de Visser ). For D. melanogaster (Figure 3, first row), the two landscapes, C1 and C2, corresponded to two correlated fitness measures, “productivity” (a measure of lifetime reproductive success) and “mating success” (Whitlock & Bourguet 2000). The posterior distributions of these two landscapes were overlapping, had the same covariance structure, and the median posterior distributions were similar. H1 and H2, two cefotaxime-resistance landscapes composed of the same mutations but with replicate MIC measurements, also had similar posterior distribution of parameters (Table 3).

Remarkably, independent empirical landscapes representing the same biological system had a similar posterior distribution of parameters. The 10 independent empirical landscapes extracted from the large yeast gene deletion data set B1–B10 (Costanzo ) gave similar posterior distributions characterized in particular by mutations of small effect and a low maximal fitness. The two empirical landscapes of vesicular stomatitis virus, E1 and E2, had extremely similar posterior distributions of parameters, although they had very different statistical properties (Table S1). Different statistical properties arise because of differences in protocol: E1 is composed of independently selected mutations, while E2 is composed of random mutations. The fact that we recover similar underlying landscapes for E1 and E2 illustrates the ability of our method to correct for variation due to protocol.

Lastly, underlying landscapes were similar when using independent empirical landscapes obtained in similar biological systems, as revealed by comparison of the two empirical landscapes of virus on their host (landscapes D and E) and of the two landscapes of bacteria adapting to a novel environment (landscapes F and G) (Figure 3, third row). In each biological system, the two landscapes represented independent experiments, yet posteriors were similar in their marginal distributions and bivariate correlation structure, revealing similar underlying fitness landscapes. The landscape of resistance to pyrimethamine also was quite distinct, with large-effect mutations, large maximal fitness, and a flat peak (I3) (Figure 3, fourth row).

Posterior predictive checks: are experimental landscapes compatible with Fisher’s model?

We tested whether the empirical landscapes we analyzed were compatible with the hypothesis that Fisher’s landscape was the underlying model for the empirical data. An informal test consisted of resimulating using the posterior distribution of parameters and examining how close these resimulated landscapes were to the data. We verified that resimulated landscapes are indeed close to the pseudodata in the cross-validation, i.e., when the true model was Fisher’s model (Figure 4, left panels). For real data, in contrast, the resimulated fitness was close to the true fitness for some but not all landscapes (Figure 4, center panels). More formally, we computed a P-value that expresses the probability that the distance between observed data and simulated data sets would occur if data followed Fisher’s model, as described under Materials and Methods (Figure 4, right panels). We computed this P-value both for the distance based on the full set of fitness values and for the distance based on summary statistics. The test of rejection based on summary statistics tests whether Fisher’s model can reproduce these average statistical properties of landscapes. The test of rejection based on the full set of fitness values tests whether Fisher’s model can reproduce the whole of the data, including specific relationships between genotypes and fitness values not captured by summary statistics. Thus, the test based on the full set of fitness values will be a stronger test of the adequacy of Fisher’s model and will reject Fisher’s model more often than the test based on summary statistics because it conserves all information in the landscape.

Figure 4

Posterior predictive checks on two example data sets. One data set is compatible with Fisher’s model (top row; Aspergillus data set A1), and one rejects Fisher’s model (bottom row, data set F). (Left) The median posterior fitness against the “true” fitness of pseudodata generated under Fisher’s model for the cross-validation showing that when the pseudodata have been generated using Fisher as the true model, the posterior fitnesses are close to the true fitness values. (Center) Posterior predicted log-fitness as a function of the true experimental log-fitness. The points are the median posterior, and the lines show the 2.5–97.5% interval. The color code indicates the number of mutations of each genotype, the ancestor in red being set to log-fitness = 0. The median posterior fitnesses are very well correlated with the true fitnesses when the landscape is compatible with Fisher’s model but less so when Fisher’s model is rejected. (Right) The median distance of pseudodata to the accepted simulations when the pseudodata are simulated under Fisher’s model and the posterior parameters. This distribution together with the observed median distance for the experimental data (dashed line) is used to calculate the P-value corresponding to the null hypothesis: “the underlying fitness landscape is Fisher’s model.”

Fisher’s model reproduced the overall statistical properties of all empirical landscapes, but in six of nine cases it could not reproduce the full structure of empirical landscapes (Table 3). The P-values based on summary statistics were almost always >0.05 (Table 3) (except for MIC landscapes H1–H4 and I1–I2, Figure S2). This indicates that the statistical properties of fitness landscapes described by the six summary statistics—mean and variance of epistasis and selection, correlation between epistasis and background fitness, maximal fitness—could be reproduced by Fisher’s model. However, Fisher’s model was not able to explain fully the structure of six of nine fitness landscapes (landscapes B, C, E, F, H and I, with P < 0.05 in Table 3). We did not identify a single reason why Fisher’s model was rejected, but it was often related to mutations with strong negative or positive epistasis (Figure 5). Fisher’s model could reproduce fully only the landscapes of A. niger (landscapes A1 and A2), of a bacteriophage adapting to its host (landscape D), and of bacteria adapting to a methanol environment (landscape G) (Table 3). In one of the landscapes compatible with Fisher’s model, the four beneficial mutations interacted almost additively (Figure 5, landscape G), but a very different landscape that includes beneficial and deleterious mutations and substantial sign epistasis among these was also compatible with Fisher’s model (Figure 5, landscape A1). In contrast, landscape C1, which looks superficially similar to landscape A1, rejected Fisher’s model. Landscape F also rejected Fisher’s model, one reason being that the third mutation had very strong positive epistasis with the first mutation. The landscape of pyrimethamine resistance (landscape I3) rejected Fisher’s model because of two cases of strong reciprocal sign epistasis. Thus, although Fisher’s model appears valuable to predict statistical properties of landscapes, in a number of cases it could not explain more detailed properties of experimental landscapes, such as mutations presenting large positive or negative epistasis.

Figure 5

Empirical landscapes compared with simulated landscapes. For each data set, the data (left) is shown side by side with the simulated genotypic landscape closest to the data in terms of Euclidean distance (center), and a typical simulated landscape, defined as the landscape, among all simulated landscapes retained by the ABC framework, whose distance to the data was closest to the median distance. The coefficient of determination R2 is also shown. Blue edges are beneficial mutations; red edges are deleterious mutations. Fitness values that are particularly unexpected under Fisher’s model are marked with a triangle.

In summary, our framework revealed biological differences between the underlying fitness landscapes of 26 experimental landscapes representing nine independent systems. Fisher’s model was generally able to reproduce the statistical properties of empirical landscapes but not their full structure. In particular, only three of nine biological systems (A, D, and G) featuring both very smooth and additive landscapes and more rugged ones had a structure that was reproduced by Fisher’s model.

Discussion

Our understanding of the structure of fitness landscapes has greatly improved, in particular, thanks to experiments that identify mutations and systematically measure the fitness of a set of genotypes bearing combinations of these mutations. Yet the generality of insights drawn from these empirical landscapes has been questioned recently (Schenk ; Szendro ; Blanquart ). The properties of empirical landscapes are heavily dependent on the particular mutations that are sampled (a small number, among a myriad of available mutations) and on the protocol used to identify mutations. We developed a novel framework based on approximate Bayesian computation to address these challenges and unravel the properties of the underlying fitness landscapes. More precisely, we inferred the underlying fitness landscape, parameterized with Fisher’s model, while accounting for the effects of the protocol on the empirical landscapes and quantifying the uncertainty due to sampling of a limited number of mutations. We used this statistical approach to conduct a survey of fitness landscapes across various species and ecological contexts.

Summary of the results

Empirical landscapes, because they are composed of a small number of mutations, generally conveyed limited information on the underlying fitness landscape. This lack of information is manifest in wide posterior distributions and a low accuracy of inference. In other words, quite different underlying fitness landscapes may generate similar empirical landscapes. This relates to a previous study where we showed, conversely, that the same underlying landscape results in a variety of empirical landscapes when multiple sets of mutations are sampled (Blanquart ). The fact that empirical landscapes are built with a small and often biased sample of mutations from the underlying fitness landscape suggests that any extrapolation on the global properties of the fitness landscape from measurement on small empirical landscapes should be taken with extreme caution.

While the size of mutations, the height of the peak (maximal fitness), and the shape of the peak were well inferred under some protocols, the number of dimensions under selection was not inferred with accuracy. Importantly, mutations independently selected in several replicates conveyed the most information on the underlying fitness landscape because they allowed an exploration of the most informative regions of the underlying landscapes. With a protocol that included as few as four mutations and all 16 possible genotypes carrying these four mutations, the size of mutations and the height and shape of the peak were well inferred (Table 2).

Fisher’s model did not accurately reproduce empirical landscapes in six of nine biological systems tested. The conceptual simplicity of Fisher’s model and its capacity, so far, to reproduce several experimental observations have made it a popular model to interpret experimental data and generate theoretical predictions (Tenaillon 2014). Fisher’s model has been used successfully before to predict the distribution of epistasis coefficients (Martin ). Fisher’s model also generates sign epistasis by optimum overshooting when the ancestral strain is close to the optimum or by pleiotropic effects when two mutations have small positive fitness effects (Blanquart ). We suggest here that although Fisher’s model is able to reproduce several statistical properties of fitness landscapes, it cannot account for their full structure in many cases. This leads to rejection of Fisher’s model even with data sets of modest size. Fisher’s model could not explain (1) sign epistasis far from the optimum (landscapes A1 and I3 in Figure 5), (2) large negative or positive epistasis (landscapes C1 and F in Figure 5), and (3) the large variance in selection coefficients and double-mutant fitness (landscapes B and E in Figure 5). It will be interesting to see whether these patterns can be explained by alternative phenotypic models that allow for some asymmetry around fitness peaks, restricted pleiotropy (mutations affect only a subset of the phenotypes), or anisotropy (mutations do not affect all traits to the same extent).

Relationship with previous studies

To our knowledge, only three studies so far have attempted to compare properties of empirical landscapes across species. Szendro quantified ruggedness for 10 experimental landscapes using a set of summary statistics. They showed that experimental levels of ruggedness are similar to those obtained with simulations of simple landscapes made of an additive component and random noise (rough Mount Fuji landscapes). They noticed the strong effect of the experimental protocol on the experimental landscape and in particular that coselected mutations tend to produce smoother empirical landscapes. However, their framework did not allow disentangling sampling variation resulting from protocol from variation owing to genuine biological differences between systems. Weinreich analyzed 14 empirical landscapes, defined higher-order epistasis coefficients, and showed that these coefficients make an important contribution to fitness in all experimental landscapes. Lastly, Weinreich and Knies (2013) fitted Fisher’s model to seven published data sets using an elegant geometric interpretation of the relationship between the epistasis and selection coefficients. They found that Fisher’s model fits the data poorly. However, it is not clear whether this was due to the data itself or to the very strong assumptions on which the analytical approach was based: the ancestral strain was always assumed to be perfectly adapted because it was set at the fitness optimum, and all mutations were considered random, so the biasing effects of selection were not accounted for.

Some of the landscapes analyzed here have been analyzed previously with Fisher’s model or similar phenotypic landscapes. Martin inferred the parameters of Fisher’s model from the distribution of selection coefficients and epistasis coefficients in an RNA virus (our data set E) (Sanjuán ). They found that the distribution of epistasis coefficients is approximately normal with a variance twice that of the variance of the distribution of selection coefficient, in agreement with theoretical predictions from Fisher’s model, when the ancestral strain is close to optimum (Blanquart ). Accordingly, we found that statistical properties of this landscape could be reproduced by Fisher’s model, but not its full structure. Last, the yeast deletion data set (B1–B10) also rejected Fisher’s model, as reported previously using a different analysis (Velenich and Gore 2013).

Several studies have attempted to fit phenotypic landscapes to data (Chou ; Rokyta ; Schenk ). In those studies, the underlying phenotypic effects of mutations are considered as parameters that are explicitly estimated, and the mapping of phenotypes to fitness is defined by a function (e.g., a Gaussian or a gamma function). This makes it easier to derive the likelihood but prevents the use of multivariate landscapes that require a number of parameters proportional to the number of dimensions. Explicitly estimating phenotypes of individual mutations gives interesting insights in the system when the underlying phenotypes are biologically meaningful and sometimes even measurable. It is also useful if one wants to predict the fitness of combinations of mutations not present in the data. However, it requires many parameters even for a simple univariate phenotypic landscape: for example, in data set D (Rokyta ), a univariate-gamma landscape includes 14 parameters, while Fisher’s model has only two, and both models perform similarly in terms of Akaike’s information criterion (AIC). Fisher’s model is a useful heuristic to make predictions on the statistical properties of fitness landscapes, but the precise value of the underlying phenotypes is less interesting in such an abstract model.

Current challenges in the analysis of genotypic fitness landscapes

In this study, we address a number of challenges to fit Fisher’s model to a diversity of experimental landscapes. But several other challenges remain to improve our understanding of fitness landscapes across species and environments:

Modeling the effects of the protocol on the experimental fitness landscape to infer properly the underling fitness landscape. Here selection was modeled using the SSWM approximation, which is valid when adaptation proceeds by successive invasions of rare beneficial mutants. This approximation was necessary to enable the fast simulations required by the ABC approach. However, in some situations of interest in experimental evolution, multiple beneficial mutations compete simultaneously in the population (clonal interference); under this regime, beneficial mutations of larger effect tends to invade the population (Nagel ). Clonal interference may be important in particular in experimental evolution (landscapes D–G). The fitness values reported also need to be ecologically relevant in the sense that they can be used to predict the fate of new mutations competing with the ancestral strain. Exponential growth rates, as reported in many studies, fulfill this condition. But other fitness measures are more dubious. For example, in drug-resistance landscapes, the fitness measure is commonly the MIC. We showed in the example of pyrimethamine resistance that the fitness landscape was quite different when a more correct fitness measure, the growth rate at a given drug concentration, was used. This invites to caution when analyzing MIC landscapes from an evolutionary perspective.

Fitting larger empirical landscapes. Empirical landscapes contain little information on the parameters of their respective underlying landscape. Larger data sets (Costanzo ; Hietpas ; Bank ) may allow more accurate inference and will become much more common in the future. Our ABC method is too computationally intensive to handle such large data sets. New theoretical developments and new statistical techniques need to be developed. These must take into account the potential biases inherent in the data-production procedure. A likelihood approach would be ideal, but unfortunately, the probability of observing a set of fitness values under Fisher’s model is hard to compute as soon as genotypes carry two mutations or more, let alone when mutations were obtained using complex protocols. In essence, this is because computing the probability of a fitness value requires integration over all possible values of the unobserved phenotypes.

Fitting other types of data. Other types of data may prove more informative than empirical landscapes. For example, Martin and Lenormand (2006) use the fitness effects of mutations across environments to infer very precisely the shape of the fitness peak (i.e., our Q parameter), which they find to be very close to Q = 2 (the Gaussian function). Perfeito show that temporal dynamics of fitness in experimental populations allow good inference of the underlying fitness landscape, including dimensionality, which is very hard to infer from genotypic landscapes. Again, new theoretical developments may reveal what type of empirical data informs best on the underlying fitness landscape.

Conclusion

We have developed a rigorous statistical framework based on Fisher’s model to infer the properties of the underlying fitness landscape from empirical landscapes. This framework differs conceptually from previous approaches because it considers an empirical landscape as a small sample in the vast space of all possible genotypes. This new approach reveals that most experimental protocols reconstruct small landscapes that carry limited information on the true underlying landscape. As a consequence, any analysis and interpretation of empirical landscapes must be embedded within a proper statistical framework that quantifies the uncertainty on the true landscape. Surprisingly, we find that a very broad class of phenotypic models that has been successful so far in interpreting experimental data is unable to explain the structure of most empirical fitness landscapes. Yet phenotypic models represent an interesting venue for future research because they can represent landscapes of large dimensionality with a small number of parameters, and they are more biologically grounded that direct genotype fitness maps. Much larger empirical landscapes will become more frequent in the future; a model-based and statistically grounded analysis of these large landscapes will improve our understanding of the structure of fitness landscapes across species and environments.