## Abstract

Genetic association studies have explained only a small proportion of the estimated heritability of complex traits, leaving the remaining heritability “missing.” Genetic interactions have been proposed as an explanation for this, because they lead to overestimates of the heritability and are hard to detect. Whether this explanation is true depends on the proportion of variance attributable to genetic interactions, which is difficult to measure in outbred populations. Founder populations exhibit a greater range of kinship than outbred populations, which helps in fitting the epistatic variance. We extend classic theory to founder populations, giving the covariance between individuals due to epistasis of any order. We recover the classic theory as a limit, and we derive a recently proposed estimator of the narrow sense heritability as a corollary. We extend the variance decomposition to include dominance. We show in simulations that it would be possible to estimate the variance from pairwise interactions with samples of a few thousand from strongly bottlenecked human founder populations, and we provide an analytical approximation of the standard error. Applying these methods to 46 traits measured in a yeast (*Saccharomyces cerevisiae*) cross, we estimate that pairwise interactions explain 10% of the phenotypic variance on average and that third- and higher-order interactions explain 14% of the phenotypic variance on average. We search for third-order interactions, discovering an interaction that is shared between two traits. Our methods will be relevant to future studies of epistatic variance in founder populations and crosses.

GENOME-WIDE association studies (GWAS) have renewed interest in methods for estimating the narrow sense heritability, which is the maximum proportion of the phenotypic variance that the additive effects found by GWAS could explain. The variance explained by the known associations for a trait is typically only a fraction of the estimated narrow sense heritability, with the remaining heritability often labeled “missing” (Eichler *et al.* 2010; Visscher *et al.* 2012).

Most estimates of narrow sense heritability come from twin and family studies (Boomsma *et al.* 2002), which can be upwardly biased in the presence of genetic interactions (Zuk *et al.* 2012). Genetic interactions introduce this bias by introducing convex nonlinearity to the relationship between phenotypic correlation and kinship (see Figure 1). Both the most common twin-studies estimator, the additive-common-environment (ACE) estimator, and a method that exploits the variation in kinship between siblings (Visscher *et al.* 2006) assume a linear relationship between phenotypic correlation and kinship. The gradient at the mean level of kinship in the population is the true narrow sense heritability (Zuk *et al.* 2012), with the gradient increasing as kinship increases above the mean due to the influence of genetic interactions. The degree to which this has biased twin- and family-study estimates depends on the amount of epistatic variance, which is not known for most complex traits.

How much variance there is from interaction effects reflects the genetic architecture of a trait and the statistical complexity of the relationship between genotype and phenotype- it is therefore of interest beyond the debate about “missing heritability.” If we knew in advance which traits exhibited considerable variance from interactions, it would help focus resources on searching for interactions in those traits.

Some of the first evidence that common variants interact to influence human traits was found by Strange *et al.* (2010), Brown *et al.* (2014), and Hemani *et al.* (2014). However, the interactions they found explained only a small amount of the phenotypic variance. By analogy with the problem of missing heritability for additive effects, it is unlikely that we will have the power to detect all of the interactions influencing a trait, so the only way to assess the statistical importance of interaction effects will be by methods that measure the variance they contribute in aggregate.

The evolutionary model of a phenotype depends on the partition of its genetic variance into additive and nonadditive components. Mackay (2014) argues that natural populations have evolved suppressing epistatic interactions as “canalizing” mechanisms, which is a possible explanation for why we do not observe many common, large marginal-effect alleles. Hemani *et al.* (2013) suggest that epistasis explains why there is more genetic variation than expected for traits under selection, a phenomenon called “stasis.” Vukcevic *et al.* (2011) and Hemani *et al.* (2013) suggest that some of the associations found by GWAS look additive only because of incomplete linkage disequilibrium between genotyped variants and interacting causal variants. These hypotheses could be tested by accurately estimating the variance from pairwise interactions.

Classical quantitative genetic theory generally assumes an infinite, outbred, random-mating population—see Gallais (1974) for a list of typical assumptions. Fisher (1918) showed how pairwise genetic interactions influence the covariance between relatives under these population assumptions. Cockerham (1954) and Kempthorne (1954, 1955) independently generalized Fisher’s result to include all orders of genetic interaction.

It has proved very difficult to estimate the variance from pairwise interactions in the outbred populations for which the theory was derived. This is because there is almost no contribution from interaction effects to the covariance between pairs of unrelated individuals in outbred populations. Samples of unrelated individuals therefore contain very little information about the contribution of interaction effects to phenotypic variation. Samples of closely related individuals do, but the information is often confounded with shared environmental and dominance effects.

The difficulty of estimating the variance from pairwise interactions can be likened to the difficulty of fitting the quadratic curve shown for the epistatic trait in Figure 1. Fitting a quadratic requires information about a wider range of points than fitting a line. Therefore, estimating the variance from interactions requires information about phenotypic correlation over a wider range of kinship than estimating the narrow sense heritability does. Bottlenecked populations are characterized by increased kinship variation and mean kinship compared to outbred populations (Abney *et al.* 2000; Carmi *et al.* 2013). A sample from a founder population therefore contains more information about the nonlinear change in phenotypic correlation with kinship, thereby bringing the variance from genetic interactions into statistical reach.

Theoretical work in founder populations was previously restricted to the contribution of the variance from pairwise interactions to the covariance between half and full siblings (Cockerham and Tachida 1988; Tachida and Cockerham 1989). Tachida and Cockerham (1989) adjust the covariance between siblings for the background relatedness introduced by the population bottleneck. We instead adjust the kinship for recent finite population size and use this expression to derive the covariance between relatives in a founder population as a function of their kinship, the mean kinship in the population, and the variance from genetic interactions of different order. We thereby give the estimator proposed by Zuk *et al.* (2012) (see Figure 1) as a simple corollary.

The theory we develop applies to certain laboratory crosses as well as to natural founder populations. Laboratory crosses often start with a small number of founding individuals, such as the “outbred” rat and mouse populations (Baud *et al.* 2013). The small number of founders gives a larger range of kinship than occurs in natural human founder populations, giving more power to estimate interaction variance. We exploit this greater power to perform the first estimates of the variance from third- and higher-order interactions in a yeast cross. We search for third-order interactions in those traits exhibiting evidence for variance from third- and higher-order interactions, and we find a third-order interaction shared between two traits that explains the majority of their covariance.

## Theory

The theory is derived for a population recently founded by a finite number of ancestors carrying a total of *A* haplotypes, with random mating after founding.

### Genotypic covariance

We consider an allele at a locus *i* on haplotypes *t* and *u*. We calculate the covariance of the allelic states of the haplotypes by conditioning on whether the alleles were inherited from the same founder: whether the alleles are identical-by-descent (IBD). The allelic state is coded as a binary variable: *g _{ti}* for haplotype

*t*and

*g*for haplotype

_{ui}*u*.

If there were *A* ancestral haplotypes at the time of the bottleneck, and *c _{i}* of these haplotypes carried the allele, then the probability that the two haplotypes carry the allele given that they are IBD at the locus is (1)where indicates that haplotypes

*u*and

*t*are IBD at locus

*i*.

Conversely, given that the two haplotypes are not IBD at the locus, then the probability they both carry the allele is (2)This is because, given the alleles are not IBD, if one haplotype inherits the allele from one of the *c _{i}* ancestral haplotypes carrying the allele, the other haplotype can inherit the allele only from one of the (

*c*− 1) other ancestral haplotypes carrying the allele—sampling from the ancestral haplotypes without replacement.

_{i}The probability that a pair of haplotypes sampled without replacement is IBD at a locus is the mean kinship coefficient, defined to be *K*_{0}, which is *A*^{−1} for a random mating population. If we define the expected allele frequency to be *f _{i}* =

*c*/

_{i}*A*, (2) can therefore be expressed as (3)Note that because

*f*≥

_{i}*K*

_{0}≥ 0, Therefore, if

*κ*

_{t}_{,}

*is the probability that haplotypes*

_{u}*t*and

*u*are IBD at locus

*i*, the probability that both haplotypes carry the allele is (4)Because the covariance between

*g*and

_{ti}*g*is therefore (5) (6) (7)

_{ui}*K*

_{0}parameterizes the adjustment of the genotypic covariance for finite population history. When

*K*

_{0}= 0, as in an infinite, random-mating population, alleles are independent when not shared IBD.

### Covariance between relatives

We derive the result first for the simple case of two interacting biallelic loci in linkage equilibrium, with the generalization following the same logic.

The phenotype of a diploid individual *t* is modeled as (8)The *x _{ti}* variables represent

*t*’s deviation in minor allele copy number from the mean at locus

*i*: (9)where and are indicator variables for whether the maternal and paternal haplotypes of individual

*t*carry the minor allele at locus

*i*. The frequency of the minor allele at locus

*i*is

*f*, and therefore The phenotypic mean is

_{i}*μ*, and

*ϵ*is the residual error, with mean zero and variance , which includes both environmental influences and random noise.

_{t}Expressing the genetic contribution to the phenotypic value in this way gives an orthogonal partition of the phenotypic variance. Because of linkage equilibrium, Cov(*x _{t}*

_{1},

*x*

_{t}_{2}) = 0. Cov(

*x*

_{t}_{1},

*x*

_{t}_{1}

*x*

_{t}_{2}) also equals zero because (10)where we have again relied on the fact that the loci are in linkage equilibrium. This implies that

*β*

_{1}is the regression coefficient of the genotype at locus 1 on the phenotype.

*β*

_{1,2}

*x*

_{t}_{1}

*x*

_{t}_{2}is the residual effect of the interaction between loci 1 and 2 after accounting for the marginal effects of the loci.

The covariance between the phenotypes of two individuals *t* and *u* relies upon the covariance of their genotypes, which is a function of the IBD sharing between their haplotypes, (11)by (7), and where is the proportion of the maternal haplotype of individual *t* that is IBD with the paternal haplotype of individual *u*. This can be expressed in terms of the kinship coefficient between *t* and *u*, defined to be *K _{t}*

_{,}

*. This is the probability that an allele drawn at random from each individual is IBD. Therefore,*

_{u}The covariance can thereby be expressed as (13)where *K*_{0} is the mean kinship coefficient in the population.

The covariance between the interaction effects is (14)by linkage equilibrium; by (13), this is equivalent to (15)Therefore, the phenotypic covariance is (16)where is the additive variance, and is the pairwise interaction variance.

The phenotypic variance of individual *t* is a function of *t*’s inbreeding coefficient, *F _{t}*. Setting

*t*=

*u*in (16) and using the fact that

*K*

_{t}_{,}

*= (1 +*

_{t}*F*)/2,

_{t}In Supporting Information, File S1, we extend the two-locus model to include dominance effects at the loci. In a founder population, inbreeding induces a correlation between the additive and dominance effects at a locus. The change in the mean due to inbreeding—inbreeding depression—introduces a further variance component. The individual-level variance is (18)where *v _{δ}* is approximately equal to the dominance variance as defined in an outbred population;

*C*

_{a}_{,}

*is the covariance between additive and dominance effects;*

_{d}*v*is the dominance variance in a homozygous population; and is the sum of the squared inbreeding depressions at the loci, which is approximately equal to

_{h}*v*when

_{δ}*K*

_{0}is small. The components, apart from

*v*, are as defined in Abney

_{δ}*et al.*(2000); however, their coefficients are different.

The variance of the phenotype in the population is found by applying the law of total variance, , to Equation 18. Because the mean inbreeding coefficient is equal to the mean kinship coefficient in a random-mating population, (19)The narrow sense heritability in the population is *h*^{2} = *v*_{1}/Var(*Y*). For an outbred population, the variance in inbreeding coefficient, Var(*F*), is zero, so the proportion of phenotypic variance explained by the interaction is *v*_{2}/Var(*Y*). However, for strongly bottlenecked populations, variation in inbreeding coefficient increases the contribution of the interaction to Var(*Y*). The dominance variance components arising from inbreeding, *v _{h}*, and do not contribute much to population variation except in the most strongly bottlenecked populations.

We now give the generalization of (16) to arbitrary epistasis between a set of causal loci each with any number of alternative alleles—for the detailed derivation, see File S1. The phenotypic covariance, for a set of causal loci *N*, is (20)where *v _{τ}* is the variance from interactions involving

*τ*loci.

If we take the limit of (20) as *K*_{0} → 0, we get (21)which is equivalent to the result of Kempthorne (1954) without dominance effects.

Under more restrictive assumptions, Zuk *et al.* (2012) derived that, for haploids, the gradient of the phenotypic correlation at the mean IBD sharing is the narrow sense heritability—see Figure 1 for a visualization of this. The diploid version of the Zuk *et al.* (2012) theorem is a corollary of (20) given by (22)This shows that (20) unifies the estimator proposed by Zuk *et al.* (2012) with the classic result of Kempthorne (1954).

The regression method proposed by Zuk *et al.* (2012) to estimate *v*_{1} does not take into account the dependencies between pairs of phenotype observations. It is therefore preferable to fit variance components by maximum likelihood or restricted maximum likelihood, as in Abney *et al.* (2000), Browning and Browning (2013), and Zaitlen *et al.* (2013). The off-diagonal elements of the phenotypic covariance matrix are given by (20), with the diagonal elements given by

### Haploid case

The theory simplifies in the haploid case due to absence of inbreeding or dominance effects. The kinship between two haploids *i* and *j*, *K _{i}*

_{,}

*, is simply the proportion of the genome shared IBD, and the phenotypic covariance matrix is (24)where*

_{j}*K*is a symmetric matrix with 1’s on the diagonal and off-diagonal elements (25)and Cov(

_{τ}**E**) is the covariance matrix of the environmental effects.

## Materials and Methods

### Simulations for variance component inference

#### Pairwise interaction variance:

To investigate the precision with which the variance from pairwise interactions could be estimated in different populations, we simulated founder populations with different mean kinship by varying the number of founding haplotypes. The R code used for the simulations is in File S2.

The allele frequencies of the variants in the ancestral population were generated by randomly sampling from a distribution with density proportional to 1/*f*, where *f* is the allele frequency. We simulated 100 variants in this way.

Each chromosome in the sample was made as a mosaic of independently inherited segments: the length of each segment was drawn from an exponential distribution with a mean of 10, and the genotypes in the segment were copied from a random ancestral haplotype. The expected number of independently inherited segments for each haplotype in the sample was therefore 10. The ancestor from whom each segment was inherited was recorded.

To calculate the diploid kinship coefficient for a pair of individuals, the total number of variants descending from the same ancestor for each of the four maternal/paternal–maternal/paternal haplotype pairs, one from each individual, was calculated; the sum total sharing across the four haplotype pairs divided by four times the number of variants gives the diploid kinship coefficient between the two individuals. The mean kinship coefficient, *K*_{0}, was taken to be the inverse of the number of ancestral haplotypes, which is its expectation. There will be negligible deviation of the sample *K*_{0}, calculated over all pairs, from its expectation.

Following the theoretical results, we calculated the component of the covariance matrix due to additive effects, defined to be *R*_{1}, by calculating element *s*, *t* of *R*_{1} as 2(*K _{s}*

_{,}

*−*

_{t}*K*

_{0})/(1 −

*K*

_{0}), where

*K*

_{s}_{,}

*is the diploid kinship coefficient of the pair of individuals*

_{t}*s*and

*t*. The component of the covariance matrix due to pairwise interaction effects,

*R*

_{2}, was calculated as the Hadamard square of

*R*

_{1}.

The kinship coefficients were calculated using all 100 variants, corresponding to calculating kinship from genome-wide IBD sharing. However, only a small proportion of the genome is likely to affect a particular trait. To simulate the sparsity of causal variants, the traits were simulated by randomly choosing 10 variants to be causal.

The variants were independently chosen for each simulated trait, covering a range of different frequency distributions of causal variants. Each variant was given an additive effect, and each pair of variants was given an interaction effect. Effects were drawn from normal distributions scaled so that *v*_{1} = 0.4 and *v*_{2} = 0.2; Gaussian error was added with variance 0.4. The variance components were inferred by fitting the covariance matrix as (26)Variance components were estimated by restricted maximum likelihood, using the average information algorithm in GCTA (Yang *et al.* 2011).

We simulated four populations with mean kinship ranging from 1/240 to 1/30, covering a broad range of human founder populations. We simulated 500 phenotypes with the same variance components for each population.

To investigate how epistasis might bias inference of additive variance, we fitted the covariance matrix as ignoring any epistasis, across the four simulated populations for the phenotypes with *v*_{1} = 0.4 and *v*_{2} = 0.2. To measure the effect of the amount of epistasis on the bias, we simulated further phenotypes, varying *v*_{2} from 0.1 to 0.4 for the population with mean kinship 1/240.

#### Third-order interaction variance:

To investigate the limits of our ability to fit epistatic variance components, for each of the four populations we simulated 200 additional phenotypes with *v*_{1} = 0.4, *v*_{2} = 0.2, and *v*_{3} = 0.2; Gaussian error was added with variance 0.2. The phenotypes were simulated as above except every combination of three causal variants was given a third-order interaction effect, scaled so that the total variance from third-order interactions was 0.2. The variance components were inferred by fitting the covariance matrix as (27)where *R*_{3} is the Hadamard cube of *R*_{1}.

### Yeast cross

Bloom *et al.* (2013) presented data from a cross of a laboratory strain and a wine strain of yeast. They sequenced the founder strains and 1008 genetically distinct haploid descendants (segregants) of the cross of the two strains. This allowed Bloom *et al.* (2013) to infer from which founder each allele had been inherited. We inferred IBD sharing proportions for each pair of haploid segregants by calculating the probability that a randomly chosen variant was inherited from the same founder. The phenotype data are final colony size for each segregant on 46 different growth media. Bloom *et al.* (2013) estimated the broad sense heritability *H*^{2} by analyzing biological replicates.

#### Inference of heritability components:

We fitted the following model to each phenotype **Y**, (28)where *K*_{1} and *K*_{2} are as defined in Equation 25 and are calculated from IBD sharing between segregants. We used the average information algorithm (Gilmour *et al.* 1995) as implemented in GCTA (Yang *et al.* 2011) to find the restricted maximum-likelihood estimates of the narrow sense heritability, *h*^{2} = *v*_{1}/Var(*Y*), and the proportion of phenotypic variance from pairwise interactions, Bloom *et al.* (2013) estimate the broad sense heritability, *H*^{2}, from analyzing biological replicates, allowing us to take advantage of the fact that (29)to estimate which we define to be This is the component of the broad sense heritability that originates exclusively in interactions involving three or more loci.

The standard error of the estimate of is estimated as where and are our maximum-likelihood estimates, and is from Bloom *et al.* (2013).

#### Simulation of epistatic traits from yeast data:

To test inference of and we simulated 500 phenotypes from the genotypes. The R code used to simulate the phenotypes is in File S3. For each phenotype, 50 causal variants were sampled independently and at random from across the genome. All 50 causal variants were given additive effects, and all pairs were given interaction effects; 10 of the 50 causal variants were chosen at random to have third-order interactions with each other; and 8 of these were chosen at random to have fourth-order interactions with each other. The effects were drawn from normal distributions scaled so that *h*^{2} = 0.4, and with the higher-order variance equally divided between third- and fourth-order interactions. Gaussian error was added so that *H*^{2} = 0.9.

#### Search for third-order interactions:

To restrict ourselves to traits that were likely to harbor third-order interactions, we searched for interactions only in those traits whose estimated variance from third-order interactions was more than twice the estimated standard error. The R code used to do this is in File S4.

Phenotypic variance differs conditional on the genotype at a locus involved in an interaction. This has been exploited to reduce the multiple-testing burden when searching for pairwise interactions (Nelson *et al.* 2013), a burden that is an order of magnitude worse when searching for third-order interactions. To find candidate loci for third-order interactions, we performed a genome-wide scan for variance-controlling loci for each trait, using the Brown–Forsythe test to assess evidence for an effect of the locus on phenotypic variance. To reduce correlation between tests, we selected the locus with the smallest *P*-value in a sliding window of 300 SNPs. We adjusted these *P*-values for multiple comparisons, using the Benjamini–Hochberg method, selecting those loci with an adjusted *P*-value <0.05.

We added the variance-controlling loci to the loci found to have significant marginal effects by Bloom *et al.* (2013). If a pair of loci from this list had a correlation of >0.7, we removed one of the pair. We performed regressions on each combination of three loci from this list, where we fitted the full model including all interaction effects. *P*-values for each effect were calculated by ANOVA. We selected only those third-order interactions with a Bonferroni-corrected *P*-value <0.05.

## Results

### Simulations

#### Pairwise interaction variance:

We simulated 5000 individuals from founder populations with varying degrees of kinship as described in *Materials and Methods*. We simulated 500 phenotypes for each population, where each phenotype had the same variance components: additive variance *v*_{1} = 0.4 and pairwise interaction variance *v*_{2} = 0.2. Gaussian error was added with variance 0.4. These variance components were chosen to test inference for a relatively small amount of epistatic variance and a relatively large amount of noise.

Figure 2A shows that the mean estimate is close to the true value for each population; the mean estimate across the four populations was 0.204, indicating the estimation was unbiased. Figure 2B shows that the standard deviations of the simulation estimates scale in proportion to 1/*K*_{0}—. It may therefore be preferable to use a smaller sample from a more strongly bottlenecked population than a larger sample from a less strongly bottlenecked population.

#### Third-order interaction variance:

To explore our ability to fit higher-order epistatic variance components, 200 additional phenotypes were simulated for each population, including third-order as well as pairwise interactions (*v*_{1} = 0.4, *v*_{2} = 0.2, *v*_{3} = 0.2). The estimation of the variance from third-order interactions was unbiased, with a mean estimate of 0.204 across the populations with *K*_{0} = 1/120, 1/60, and 1/30. There was not enough information to fit the model in the population with *K*_{0} = 1/240.

The standard deviation for the third-order interaction variance was at least twice as large as the standard deviation for the pairwise interaction variance across the populations (Figure 3). Even for the most strongly bottlenecked population, the standard error for the third-order interaction variance is nearly 15% of the phenotypic variance, comparable to the size of the variance component.

#### Ignoring epistasis biases additive variance estimates:

To investigate any possible bias in the estimation of narrow sense heritability that may arise from ignoring epistasis, we fitted models with only additive variance components to the simulation data used in Figure 2. We found that the bias did not depend on the mean kinship of the population (Table S3).

We simulated additional phenotypes with varying amounts of epistasis (*v*_{2} = 0.1, 0.2, 0.3, 0.4) for the population with *K*_{0} = 1/240. The amount of bias was proportional to the amount of epistatic variance; the additive variance estimates were inflated by ∼6% of *v*_{2} (Table S4).

Ignoring epistasis resulted in inaccurate standard error estimates. Figure 4 shows that even when only 10% of the phenotypic variance is epistatic, the standard error of the additive variance estimated from simulations is >15% larger than the standard error estimated by GCTA.

### Approximate analytic standard error

The amount of information a sample contains about the epistatic variance of a trait depends on the distribution of kinship in that sample. To better understand how the standard error of the estimator of *v*_{2} depends on the moments of the kinship distribution, we extend the analogy of fitting a quadratic from Figure 1 to derive an approximate analytic expression for the standard error.

If we define *R _{s}*

_{,}

*= 2(*

_{t}*K*

_{s}_{,}

*−*

_{t}*K*

_{0})/(1 −

*K*

_{0}), then the process of fitting the covariance matrix implied by (20) can be likened to fitting a polynomial in

*R*. Fitting the additive variance,

*v*

_{1}, and the variance from pairwise interactions,

*v*

_{2}, is similar to fitting to all pairs

*s*,

*t*the regression model (30)where is the observed similarity between

*s*and

*t*. If

*s*and

*t*are independent and

*Y*has variance 1, then

*σ*

^{2}= 1.

In File S1, we derive the asymptotic variance of the maximum-likelihood estimator of *v*_{2} in this model. If *μ _{c}* is the

*c*th central moment of the distribution of

*R*, then (31)where

*η*is the number of pairs, and where we take

*σ*to be 1.

Testing this using the simulation results, we find the standard deviation of the simulation estimates for *K*_{0} = 1/30, 1/60 to be very close to (31) calculated from the sample kinship statistics, with the error increasing above (31) for smaller *K*_{0}—see Figure S1.

The information about the epistatic variance in a sample increases with the fourth central moment of the kinship distribution, which is the unnormalized kurtosis. The samples with heavy tails in their kinship distributions are therefore likely to have the largest kurtosis and the most information about epistatic variance.

### Yeast cross

We analyzed data, further described in *Materials and Methods*, from a cross of a laboratory strain (BY) and a wine strain (RM) of yeast (Bloom *et al.* 2013). The data included 46 growth phenotypes measured for 1008 haploids dissected from tetrads produced by crossing the two founder strains.

#### Variance components:

To establish that our methods worked for these data, we first simulated epistatic traits from the genetic data, as detailed in *Materials and Methods*. The variance component estimates from the simulated traits were unbiased and the standard error estimates were accurate (Table S2).

Next we applied our approach to partition the phenotypic variance of the 46 growth phenotypes into additive, pairwise, and higher-order genetic components, plus a residual. Figure 5 visualizes this partitioning. The numerical results of the analysis are in Table S1 and File S5. Five phenotypes had estimates of the variance from higher-order interactions >3 standard errors from zero. The mean proportions of phenotypic variance explained by pairwise interactions and higher-order interactions were 0.10 and 0.14, respectively.

#### Third-order interactions:

To further explore the contributions of higher-order interactions, we searched directly for third-order interactions as described in *Materials and Methods*. Table 1 shows the only two third-order interactions found to be significant at the 0.05 level after Bonferroni correction. For each marker in the interaction for the formamide phenotype, there is a corresponding marker in near perfect linkage in the interaction for the indoleacetic acid phenotype. This suggests that the interaction is shared between the two phenotypes, and it may explain some of their covariance. To estimate how much of the covariance this shared third-order interaction explains, we calculated the fitted values for the full interaction models, including the additive and pairwise effects, for each trait separately, and the covariance of these fitted values. The covariance of the fitted values explains 59% of the covariance between the phenotypes.

## Discussion

### Theory

We used an approach based on allelic indicator variables to calculate the covariance between individuals in a founder population as a function of their kinship, the mean kinship in the population, and the variance components of the phenotype. This extends the classic result for outbred populations (Kempthorne 1954) to founder populations. Equations 20 and 23 together determine the phenotypic covariance matrix, the parameters of which can be estimated by (restricted) maximum likelihood by assuming the phenotype follows a particular distribution. These parameters, along with the central moments of the inbreeding distribution, determine the proportion of population variance explained by different orders of genetic interaction in a founder population. The relationship between (20) and Figure 1 can be seen by writing the phenotypic covariance as a polynomial function of *R* = 2(*K* − *K*_{0})/(1 − *K*_{0}), the *x*-axis of Figure 1. The correlation for the epistatic trait in Figure 1 as a function of *R* is (32)where *v*_{1} is the additive variance in the population and is 0.4.

The model applies exactly only to populations that have been randomly mating since being founded; however, the allelic indicator variable approach could be extended to non-random-mating populations by considering models for nonrandom inheritance of alleles.

Extending the method to include linkage disequilibrium would be possible but would rely on knowing the linkage disequilibrium between unknown causal alleles. We note that identity-by-descent-based methods such as this are biased only by linkage disequilibrium between causal alleles, whereas identity-by-state (IBS)-based methods, such as in Yang *et al.* (2010), are also biased by varying linkage disequilibrium between marker alleles and causal alleles (Speed *et al.* 2012).

We have derived the individual- and population-level variance decomposition for two interacting loci with dominance in a founder population (see File S1 for details). Finite population history induces dependence between allelic states, both within and between individuals. This prohibits a simple and exact expression for the covariance between relatives for the additional variance components introduced by dominance. However, except for the most strongly bottlenecked populations, using the identity states implemented by Abney *et al.* (2000) will probably give a good approximation.

We note that an alternative theoretical approach would be to extend the frequency-weighted IBS estimator employed by Yang *et al.* (2010) to epistatic variance components. The IBS-based method of estimating the additive variance was compared to an IBD-based approach by Zaitlen *et al.* (2013), using Icelandic data. They found the IBS-based approach underestimated the additive variance relative to the IBD-based approach by a considerable amount. The same underestimation would be expected to apply to an IBS-based epistatic variance estimator, because it originates from incomplete linkage disequilibrium between markers and causal variants. The underestimation could be even more severe for epistatic variance components because, for the variance from an interaction to be properly detected, all of the loci involved in the interaction would have to be in strong linkage with the markers. We therefore argue for the IBD-based approach, which takes advantage of the long shared segments present in a founder population to reduce the bias in the estimates.

### Simulations and sampling

The simulations (Figure 2) suggest that with a sample of 5000 Hutterites, it would be possible to estimate the variance from pairwise interactions with a standard error of <5% of the phenotypic variance. The standard error scales in proportion to the inverse of the mean kinship. This explains why one cannot simply use a very large, random sample from an outbred population to fit the epistatic variance, as the standard error tends to infinity as *K*_{0} tends to zero.

More recently founded populations that have expanded rapidly should be preferred to older populations, due to the reduced amount of recombination since the founding, which results in more kinship variation.

Including close relatives would increase the precision of the estimator. However, the similarity between close relatives could be due to shared environment as well as shared genetics. The confounding with shared environment could be ameliorated by fitting additional variance components for different relative classes. However, if shared environmental effects extend beyond first-degree relatives, the model may become excessively complex. Dominance could also cause similarity between siblings above what is expected by additive effects, leading to overestimation of the epistatic variance. For traits that are known to have little dominance variance or shared environmental effects, including close relatives would increase precision without causing bias. Otherwise, very large samples of close relatives would be required to disentangle epistasis, dominance, and shared environment.

Power calculations can be aided by the approximate analytic formula for the standard error of the variance from pairwise interactions. This acts like a lower bound for the standard error in the simulated data—see Figure S1. The moments of the kinship distribution can be calculated from a small sample, and, from these, an estimate of the standard error can be calculated for different sample sizes. If this is too high to give a useful estimate, then the sample is probably not appropriate for estimating the variance from pairwise interactions.

Direct estimation of the variance from third-order interactions may be beyond the limits of possibility for current human samples. Even with 5000 Hutterites, the standard error is likely to be at least 15% of the phenotypic variance (Figure 3). Unless this component is a large part of the phenotypic variance, it is unlikely that any current samples of human founder populations would provide the power to detect that the component is nonzero.

Founder populations have recently been used to estimate the narrow sense heritability (Browning and Browning 2013; Zaitlen *et al.* 2013). We found in simulations that ignoring epistasis leads to a slight overestimation of the additive variance in proportion to the amount of epistatic variance (Table S3 and Table S4), as well as underestimation of the standard error (Figure 4). This could cause improper calibration of statistical tests. It is possible that these problems could be reduced by restricting to a smaller range of relatedness, but this would increase the standard error.

### Yeast cross

Bloom *et al.* (2013) found evidence for epistatic variance in the difference between *H*^{2} and *h*^{2}. In the yeast cross analysis (Figure 5 and Table S1) we have gone further by partitioning the epistatic variance into components arising from pairwise interactions and from third- and higher-order interactions. While the individual estimates of the higher-order components are not very precise, we provide strong evidence that the variance from pairwise interactions does not in general explain all of the difference between *H*^{2} and *h*^{2}. It is impossible to draw precise conclusions about the relative size of and for individual traits, because the method of estimation results in a negative correlation between the estimates. A larger sample from a similarly designed experiment could overcome some of these difficulties and enable direct estimation of the variance from third-order interactions.

The relatively small amount of variance explained by the only third-order interaction found (Table 1) suggests that there are many more third- and higher-order interactions, each explaining a small amount of the variance. While a small number of third-order and higher-order interactions have previously been identified (Pettersson *et al.* 2011; Taylor and Ehrenreich 2014), this is the first such pleiotropic interaction to be discovered, as far as we are aware.

The statistical importance of pairwise and higher-order interactions in the yeast cross cannot be readily generalized to natural populations. For some interaction models, the proportion of the variance that is epistatic rather than additive is greatest for interacting alleles at intermediate frequencies (Hill *et al.* 2008; Mackay 2014). Therefore, if interactions occur between rarer alleles in natural populations, the proportion of the variance that is epistatic could be reduced.

The large amount of epistatic variance in the cross could be explained by the breakdown of coadapted variant combinations. The cross is between a laboratory strain and a wine strain of yeast, which have diverged under different selection pressures (Liti *et al.* 2009). Given that hybrid incompatibilities were observed between experimentally evolved strains (Anderson *et al.* 2010), it is plausible that these strains have accumulated them.

The large amount of epistatic variance arising from third- and higher-order interactions in particular could be explained by the buildup of hybrid incompatibilities. Incompatibilities between three or more loci are theoretically expected to be more common than incompatibilities between two loci, because a greater proportion of evolutionary paths to higher-order incompatibilities do not pass through a less fit intermediary (Orr 1995). Figure 6 shows that the fitness profile of the discovered third-order interaction is consistent with the interaction being the result of an evolved Dobzhansky–Muller-like incompatibility. Only one of the eight possible three-locus genotypes results in a sizeable reduction in fitness. Two-thirds of the evolutionary paths from a common ancestor to the two yeast strains would have avoided this less fit intermediary and would therefore have been allowed by selection (Orr 1995).

### Conclusion

These methods can be used to investigate the role of pairwise and higher-order epistasis in model organisms by applying them to appropriate crosses. In particular, by measuring the variance that higher-order interactions contribute to crosses between diverged populations, these methods could be used to investigate the role of higher-order interactions in hybrid incompatibilities.

We anticipate that it will soon be possible to apply these methods to precisely estimate the variance from pairwise interactions in human founder populations. These estimates, combined with estimates of the additive and dominant components of the variance, will help in answering where the missing heritability is, in searching for causal loci, in building prediction models, and in testing evolutionary models of traits.

## Acknowledgments

We acknowledge Peter Donnelly for helpful comments on the manuscript. This work was supported by the Wellcome Trust (grants WT098051 to Richard Durbin, 099670/Z/12/Z to Alexander Young, and 090532/Z/09/Z to the Wellcome Trust Centre for Human Genetics).

## Footnotes

Supporting information is available online at http://www.genetics.org/lookup/suppl/doi:10.1534/genetics.114.170795/-/DC1.

*Communicating editor: G. A. Churchill*

- Received September 16, 2014.
- Accepted October 10, 2014.

- Copyright © 2014 by the Genetics Society of America

Available freely online through the author-supported open access option.