On the Relationship Between High-Order Linkage Disequilibrium and Epistasis

Abstract

A plausible explanation for statistical epistasis revealed in genome wide association analyses is the presence of high order linkage disequilibrium (LD) between the genotyped markers tested for interactions and unobserved functional polymorphisms. Based on findings in experimental data, it has been suggested that high order LD might be a common explanation for statistical epistasis inferred between local polymorphisms in the same genomic region. Here, we empirically evaluate how prevalent high order LD is between local, as well as distal, polymorphisms in the genome. This could provide insights into whether we should account for this when interpreting results from genome wide scans for statistical epistasis. An extensive and strong genome wide high order LD was revealed between pairs of markers on the high density 250k SNP-chip and individual markers revealed by whole genome sequencing in the Arabidopsis thaliana 1001-genomes collection. The high order LD was found to be more prevalent in smaller populations, but present also in samples including several hundred individuals. An empirical example illustrates that high order LD might be an even greater challenge in cases when the genetic architecture is more complex than the common assumption of bi-allelic loci. The example shows how significant statistical epistasis is detected for a pair of markers in high order LD with a complex multi allelic locus. Overall, our study illustrates the importance of considering also other explanations than functional genetic interactions when genome wide statistical epistasis is detected, in particular when the results are obtained in small populations of inbred individuals.

Keywords: Arabidopsis thaliana; epistasis; high order linkage disequilibrium; leaf; molybdenum.

Figure 1

Illustration of how the pseudomarkers (P₁, P₂, P₃, P₄) used in the estimation of the second order linkage disequilibrium between a pair of linked or unlinked markers (predictors; M₁ and M₂,) and a third linked or unlinked functional polymorphism (target; Q) are created. The pseudomarkers together represent the possible bi-allelic formulations of the two-locus M₁M₂ genotypes. The maximum pairwise LD-r² between the target and the four pseudomarkers (P₄) defines the second order LD between the predictors (M₁, M₂) and the target (Q). The general two-locus epistatic model (Model 1 below) will, when fitted to the genotypes of the predictors (M₁, M₂) capture the variance of the target (Q).

Figure 2

Illustration of how the prevalence of high order LD-r² to the targets in a 6Mb window on A. thaliana chromosome 2 (8 – 14Mb) depends on distance of the predictors from the target. The color gradient illustrates the proportion of predictor pairs that reach a particular LD-r² (x-axis) depending on the distance between the nearest predictor and the target (y-axis). Results are presented for populations with n = 100 (A) and n = 728 (B) individuals.

Figure 3

Strong second order LD-r² exists also when the individual predictor to target LD-r² is weak. The intensity of each dot illustrates the number of cases with a particular high order LD-r² / maximum individual predictor to target LD-r² combination. Dots below the line are cases where the high order LD-r² stronger than any individual predictor to target LD-r² (n = 728).

Figure 4

Number of predictor pairs of different classes in strong high order LD-r² (>0.6) to targets detected in the evaluated windows and estimated genome wide. The distribution of LD-r² values > 0.6 for the cis-cis, cis-trans, trans-trans predictor pairs for (A; n = 100) and (B; n = 728) The total number of predictor pairs with high order LD-r² above 0.6 in the three classes are summarized in (C) and used to estimate the total expected number of predictor pairs in the entire genome (D; error bars show the estimation error estimated from the results obtained for the three window (Materials and Methods).

Figure 5

An illustration of how the high order LD between four polymorphisms affecting the level of molybdenum in the A. thaliana leaf (Forsberg et al. 2015), likely explains the significant statistical epistasis detected for a cis-trans predictor pair. (A) Boxplots illustrating the phenotypic distribution in the four genotype classes defined by the cis-trans predictor pair with the strongest significant epistatic interaction to the level of molybdenum in the A. thaliana leaf. (B) Illustration of the connection between the two-locus genotypes of the predictor pair and the minor alleles at the four linked loci associated with this trait on chromosome 2 (Forsberg et al. 2015). The top box in (B) illustrates the two-locus genotype for the predictor pair, with the width of each sub-box indicating the number of individuals in each genotype class in the population. In the bottom box in (B), each individual is represented as a column, where green (molybdenum decreasing) and orange (Molybdenum increasing) colors indicates that the individual carry the minor alleles at the four loci identified in (Forsberg et al. 2015). mGWA1 and mGWA2 are SNP markers associated with the trait and 53del and 326 are structural polymorphisms (Forsberg et al. 2015).