Total cholesterol, low-density lipoprotein cholesterol, triglyceride, and high-density lipoprotein cholesterol (HDL-C) levels are among the most important risk factors for coronary artery disease. We tested for gene–gene interactions affecting the level of these four lipids based on prior knowledge of established genome-wide association study (GWAS) hits, protein–protein interactions, and pathway information. Using genotype data from 9,713 European Americans from the Atherosclerosis Risk in Communities (ARIC) study, we identified an interaction between HMGCR and a locus near LIPC in their effect on HDL-C levels (Bonferroni corrected Pc = 0.002). Using an adaptive locus-based validation procedure, we successfully validated this gene–gene interaction in the European American cohorts from the Framingham Heart Study (Pc = 0.002) and the Multi-Ethnic Study of Atherosclerosis (MESA; Pc = 0.006). The interaction between these two loci is also significant in the African American sample from ARIC (Pc = 0.004) and in the Hispanic American sample from MESA (Pc = 0.04). Both HMGCR and LIPC are involved in the metabolism of lipids, and genome-wide association studies have previously identified LIPC as associated with levels of HDL-C. However, the effect on HDL-C of the novel gene–gene interaction reported here is twice as pronounced as that predicted by the sum of the marginal effects of the two loci. In conclusion, based on a knowledge-driven analysis of epistasis, together with a new locus-based validation method, we successfully identified and validated an interaction affecting a complex trait in multi-ethnic populations.
The catalog of genome-wide association studies (GWAS)  has collected to date over 1,194 publications since the end of 2008, for a total of over 5,697 single nucleotide polymorphisms (SNPs) that are associated with complex human diseases and other complex traits. However, most these associated SNPs exhibit a small effect size, and collectively only explain a relatively small fraction of additive variance , , , . Specifically, a recent meta-analysis of several GWAS, studying a combined sample size between ∼20,000 to ∼100,000 individuals, identified 95 loci associated with the level of one of total cholesterol (TC), low-density lipoprotein cholesterol (LDL-C), triglyceride (TG), and high-density lipoprotein cholesterol (HDL-C) . In aggregate, these loci explain only 25–30% of heritable variation for each trait . Many hypotheses aiming to explain the missing heritability of GWAS have been proposed, including structural variants, rare variants, gene-environment interactions, epigenetics, and complex inheritance , , , . Because gene-gene (epistatic) interactions may contribute to missing heritability to some extent , , , here we seek to find examples of pairs of loci that interact in their effects on any of the four lipid levels, which are important risk factors of coronary artery disease .
Epistasis has been investigated in order to understand the relationship between genotype and phenotype since Bateson  discovered in 1905 that some genes can suppress the effects of others. Thereafter, a number of epistatic interactions have been identified in QTL mapping studies or GWAS in humans ,  and other organisms , , . Studies of model organisms suggest that gene-gene interactions are a common phenomenon , , , . However, they have proven difficult to detect in humans, chiefly due to the limited statistical power associated with the large combinatorial number of tests and the skew towards low minor allele frequencies , . Hence, in order to increase power to detect gene-gene interactions in GWAS, a series of methods have been developed to prioritize candidate SNPs using prior knowledge of established GWAS hits , and recently also using knowledge of protein-protein interactions (PPIs) ,  and pathway information .
Although some interactions affecting complex diseases and traits have been reported in humans , , replication of these interactions in independent samples has proven difficult . He et al. showed that this low replication is in part attributable to low power and small effect sizes of tag SNPs in GWAS. For two interacting causal loci, the observed interaction effect between two respective tag SNPs (each tagging one of the causal loci) is proportional to the underlying causal interaction effect multiplied by the product of the two linkage disequilibrium (LD) coefficients between each tag SNP and the respective causal variant. This decrease in the measured interaction effect reduces the statistical power of the interaction test and it also reduces the probability of replication of significantly identified interactions. This reduction is further exacerbated by heterogeneity in the LD structure between different populations and among population samples. These are the same problems that plague the power of single-marker GWAS tests, but they are exacerbated in interaction testing, with a quadratic dependence on LD between markers and causal loci, which lead to a much greater reduction in power. Motivated by this problem, Liu et al. proposed a local validation analysis and successfully replicated the loci of a few interactions underlying common human diseases.
In this study, we aim to improve the power to detect gene-gene interactions in existing large-scale GWAS data sets by considering for interaction testing only a highly focused set of candidate SNPs extracted from prior information of known GWAS hits, PPIs, and pathway information. To improve the power of replicating gene-gene interaction signals in independent samples, we introduce an adaptive locus-based validation procedure that follows an approach similar to Liu et al.. Applying these procedures for testing for gene-gene interactions underlying lipid levels, we discovered a significant interaction affecting HDL-C levels, which provides new insights into the genetic architecture of this complex trait. Using the adaptive locus-based validation procedure, we also successfully replicated this novel interaction in four independent cohorts, including two cohorts of different ethnicity.
Knowledge-driven identification of gene–gene interactions
We tested the statistical significance of gene-gene interaction between each pair of SNPs among 125 SNPs from 95 loci that have been previously individually associated with any of the four lipid levels  for a total of 7,750 tests, out of ∼3 trillion possible tests between each pair of SNPs in our data. Tests of interaction were conducted using genotype data or imputed genotypes in a sample of 9,713 European Americans (EAs) from the Atherosclerosis Risk in Communities (ARIC) study  (Materials and Methods). We used an F-test with four degrees of freedom within a linear model framework for interaction testing , . This test considers the 3×3 table of genotype pairs for two SNPs and tests for significant interaction between the two SNPs on top of any additive or dominance effects that each of the SNPs might exhibit by itself. For consideration of statistical power and robustness, we discarded from testing pairs of SNPs for which one or more of the 9 genotype-by-genotype combinations appeared in fewer than 20 individuals in our sample (Materials and Methods).
Testing for interaction between 7,750 pairs of SNPs for each of four quantitative traits, we identified one significant interaction underlying each of LDL-C level and HDL-C level (Figure 1a). The interaction underlying LDL-C level is between rs2247056 and rs1030431 (Bonferroni corrected Pc = 0.003; Figure 1a). To explore the interaction between the two loci with better resolution, we tested for interaction between each SNP in the 100 kb surrounding rs2247056 and each SNP in the 100 kb surrounding rs1030431 and found that the interaction signal peaked between rs2853928 and rs1993453 (Pc = 0.01 after accounting for all additional pairs of SNPs tested; Figure S1). The discovery SNP pairs are in high LD with the fine-mapped SNP pairs, with an r2 value of 0.997 between rs2247056 and rs2853928 and 0.999 between rs1030431 and rs1993453. The former two reside near a pseudogene, LOC100133383, and the latter two are located near and in gene UBXN2B, respectively. However, this suggestive interaction underlying LDL-C did not replicate in independent cohorts.
Henceforth, we focus on the interaction between rs12916 and rs1532085 on HDL-C levels (Pc = 0.008), since its validation in additional cohorts is highly significant, as described below. We first tested for interaction between each SNP in the 100 kb surrounding rs12916 and each SNP in the 100 kb surrounding rs1532085. While many of these pairs show significant interactions (Figure 1b), as expected from LD, we observed the strongest signal between rs3846662 and rs2043085 (Pc = 0.002). The fine-mapped pair of SNPs is in high LD with the original pair of SNPs, with an r2 value of 0.88 between rs3846662 and rs12916 and an r2 value of 0.93 between rs2043085 and rs1532085 (Figure S2). rs3846662 is intronic in HMGCR (Table 1), which has not been previously associated with HDL-C, but has been associated with both TC and LDL-C levels . rs2043085 is upstream of LIPC (Table 1), which has been previously found to be associated with HDL-C .
The interaction between rs3846662 and rs2043085 affects HDL-C twice as much as the effect of the polymorphism in LIPC alone: While individuals with TT genotype at rs2043085 already exhibit an average increase of 2.63 mg/ml in HDL-C (standard error (SE) = 0.014; Figure 2a), this genotype in combination with an AA genotype at rs3846662 leads to an average increase of 5.72 mg/ml (SE = 0.041; Figure 2b). The linear model with these two SNPs has an R-square value of 0.5% and the linear model with the two SNPs and their interaction has an R-square value of 0.8%, which indicates that the interaction explains additional 0.3% of the overall variation in HDL-C levels (Materials and Methods; Table S1). We tested whether rs3846662 and rs2043085 exhibit gene-gene interactions underlying any of the other lipid levels, and found a nominally significant interaction underlying LDL-C (P = 0.028), and almost significant interaction underlying TG (P = 0.08) in ARIC.
We performed a larger scale interaction analysis between all pairs of SNPs that (i) are found in interacting genes according to a curated human protein-protein interaction network (∼6 million pairs), or (ii) are involved in the pathway of metabolism of lipids and lipoproteins (∼27 million pairs). All SNPs in a gene were considered, as well as in the 5 kb regions upstream and downstream. This analysis detected no significant gene-gene interactions following Bonferroni correction (Pc≥0.58 for PPIs; Figure S3; Pc≥0.14 for pathway; Figure S4).
Validation of gene–gene interaction in pairs of loci
Considering the quadratic reduction in replication power as a function of LD between tag SNPs and causal loci, we aimed to increase power via an adaptive locus-based validation procedure that is related to that of Liu et al.. In considering a replication dataset, the procedure follows three sequential stages that leverage the signals of proxy markers: (i) test for interaction between the original SNP pair between which gene-gene interaction has been detected; (ii) test for interactions between each of the two original SNPs and each SNP in the proximate region containing the other original SNP; (iii) test for interactions between each pair of SNPs in each of the two respective proximate regions containing the two original SNPs. This validation procedure proceeds sequentially and stops at any stage when significant interactions were detected after multiple-testing correction. Both the method of Liu et al. and our adaptive locus-based validation method focus on replicating the interaction between a pair of loci, rather than between a pair of SNPs, due to the power limitations of replicating an interaction between SNPs. The null hypothesis of the entire three-stage procedure is that there is no interaction between the pair of loci, rather than just between the pair of SNPs, thus the procedure continues sequentially as described to consider proxy SNPs from the loci containing each original SNP. Replication is successful if an interaction between any SNP pair from the two loci is significant after multiple-testing correction. Similar locus-based approach has also been used in the context of gene-based GWAS tests for single-marker association, which use an entire gene or locus as the testing unit of association, rather than a single SNP , .
To validate the gene-gene interaction affecting HDL-C, we performed replication analyses in two additional GWAS datasets from the Framingham Heart Study (FHS)  and the Multi-Ethnic Study of Atherosclerosis (MESA) , as well as in the African American (AA) cohort from the ARIC study . Using our adaptive locus-based procedure, we tested for interaction sequentially between SNPs surrounding rs3846662 and SNPs surrounding rs2043085. We observed significant interactions in the two additional EA cohorts from FHS and MESA (Figure 1c), with Pc = 0.002 and Pc = 0.006 for the most significantly interacting SNP pair (Table 1). Replication was also significant in Hispanic Americans (HA) from MESA and AAs from ARIC (Figure 1c; Table 1). The R-square of linear model with the two interacting SNPs varies between 0.2–0.5% across the four replication cohorts, with the interaction term between the two explaining an additional 0.2–1.1% of the overall variation in HDL-C levels (Table S1). The replication procedure failed in a sample of AAs from MESA (Figure S5).
None of the successful replications were replicated at stage (i) of the adaptive locus-based validation procedure, which means that an interaction between the same SNP pair is not observed significantly in the additional samples. The interaction was successfully validated in stage (ii) of the three stages in the MESA EAs, with the same SNP in HMGCR (rs3846662) and a proxy SNP near LIPC exhibiting a significant gene-gene interaction after multiple-testing correction. The other three successful replications occurred at stage (iii) (Table 1), emphasizing the importance of a locus-based replication approach. The combined evidence from the discovery and four different validation cohorts for a gene-gene interaction between the two loci under study is overwhelmingly significant, even following a conservative Bonferroni correction (Pc = 9.0×10−8).
While the gene-gene interaction signal peaks for different pairs of SNPs across the different cohorts (Table 1), the type of interaction and effect patterns appear consistent across several sample sets (Figure S6). To test this formally, we partitioned the significant SNP-SNP interactions into the four possible interaction components on top of the marginal SNP effects, namely additive by additive (A×A), additive by dominance (A×D), dominance by additive (D×A), and dominance by dominance (D×D) components (Materials and Methods). Considering a nominal significance level of 0.01, D×A and D×D components are significant and underlie the significant interaction in the ARIC EA discovery set, between rs12916 and rs1532085 (Table S1). All four terms are significant between the pair of SNPs, rs3846662 and rs2043085, that resulted from fine mapping in the same discovery set, with D×A and D×D being of the same effect direction (sign) and similar effect sizes as between rs12916 and rs1532085 (Table S1). Examining the two replication cohorts of a similar (EA) ancestry, the interaction in the MESA cohort similarly shows significant D×A and D×D components, with same effect direction, though with larger effect sizes and a higher proportion of phenotypic variance explained (Table S1). None of the four terms is significant by itself in the EA FHS cohort. These results of consistent patterns of interaction across the EA cohorts support the possibility that they are all governed by the same (unobserved or partially unobserved) interacting variants.
Validation of imputation accuracy
To verify that our results are not an artifact of imputation errors, we compared imputed genotypes of the two SNPs (rs12916 and rs3846662) that were involved in significant interactions and for which we could obtain measured genotype data from an independent source, using the ITMAT/Broad/CARE (IBC) Vascular Disease 50 k SNP Array chip . For these two SNPs, r2 between imputed and actual genotypes is 0.914 and 0.921 and the genotype concordance rate is 94.5% and 94.7%, respectively. Although the imputation is not perfect, the two interaction tests involving these two SNPs are at least as significant when replacing imputed genotypes with measured IBC genotypes, consistent with imputation errors adding noise and masking some of the signal, rather than biasing the statistical test.
Tests of gene-gene interactions are not as powerful as tests of single-marker association, so a judicious strategy is essential for successful interaction analysis in GWAS , . The first step is to determine the size of the analysis, genome-wide or focusing on candidate SNPs. This step should consider the sample size, possible effect size of the underlying interaction, and the desired statistical power. Current single-marker GWAS have been successful in detection of single-marker associations for many complex diseases or traits using a stringent genome-wide significance level (P<5×10−8). To achieve a similar success for interaction analysis, we are limited to performing ∼1 million tests even if the interaction test and single-marker test had the same statistical power. This limitation means that we are not able to conduct an inclusive all-by-all pair-wise interaction analysis in current GWAS. Thus, in this study we only tested for interactions between candidate SNPs based on prior knowledge.
We used three types of prior knowledge, known GWAS hits, protein-protein interaction networks, and known functional pathways. These three analyses might be different in the enrichment of epistasis signals and are also different in the number of interaction tests, 7,750 based on known GWAS hits, ∼6.2 million using PPI, and ∼27 million with pathway information. We found significant interactions from the 7,750 interaction tests using known GWAS hits. As the sample size of ∼10,000 individuals is relatively large among existing GWAS, this indicates that the observed (tagged) effect size of any other underlying interactions is no larger than the marginal effects of single SNPs. It is also likely that the epistasis signals are better enriched between markers that are marginally associated with lipid traits such that testing interactions among known GWAS hits is more powerful in our study. Therefore, our results suggest that a small-scale interaction analysis of candidate SNPs driven by known marginal associations might be a good choice for detecting epistatic interactions in current GWAS.
Recently, the Population Architecture using Genomics and Epidemiology Study  found only ∼50% of the 125 reported associations with lipid levels  to replicate in three non-European cohorts. Due to the quadratic decrease in the interaction effect of tagged markers, gene-gene interactions are even less likely to replicate in diverse populations. Leveraging signals from proximate linked SNPs, our adaptive locus-based method successfully validated gene-gene interactions between HMGCR and LIPC in four additional, independent cohorts, including two of non-European ancestry. Although the most significant interaction in each cohort involves different SNPs, they are proximate across the cohorts, with stronger LD and smaller distances amongst the three EA cohorts and weaker LD and larger distances between them and the HA and AA cohorts (Figure S2 and Table 1). The differences in distance and LD between ethnicities could be due to differences in genetic background, demographic history, and natural selection, even if the different SNP pairs capture the same underlying causal interaction. However, the interaction shows similar patterns among some, but not all cohorts (Figure S6 and Table S1), while the different SNPs around HMGCR are in strong LD, and those around LIPC show weak LD (Figure S2). These results suggest that the five SNP pairs either capture separate causal interactions or are only in weak LD with the same pair of interacting, unobserved variants.
Another possibility is that the interaction is between relatively rare causal variants: Much like rare causal variants can lead to multiple independent associations of common variants, dubbed “synthetic associations” , an interaction between two rare causal variants can produce an even larger number of independent “synthetic interactions”, which can in principle explain almost-independent, yet proximate gene-gene interactions. Another possibility is that the underlying interaction is more complex and involves more than a pair of SNPs. In that case, in our analysis of pairs of SNPs, each pair might tag only certain aspects of the underlying interaction.
Both HMGCR and LIPC are involved in metabolism of lipids and lipoproteins. HMGCR, which has been associated with TC and LDL-C , regulates the rate of cholesterol synthesis via a negative feedback mechanism mediated by sterols and non-sterol metabolites . LIPC encodes hepatic lipase which is an important enzyme in HDL metabolism  and has been previously associated with HDL-C levels . The interaction between variants in these genes as discovered in this study can be possibly explained by an indirect interaction between cholesterol synthesis and the metabolism of LDL and HDL particles. HGMCR is the rate-controlling enzyme in the mevalonate pathway for cholesterol synthesis . Much of this cholesterol will form cholesteryl esters that will be packaged into various lipoproteins including LDL, HDL, and TG-rich lipoproteins. There are a number of known lipoprotein interactions that result in the flow of cholesterol in the form of cholesteryl esters from LDL and VLDL to HDL-C . This cholesterol is later processed with the HDL particle by either reabsorbing into the liver or excretion in the urine .
The rs2043085 SNP in the LIPC gene region, where our strongest signal has been observed in fine mapping in the discovery panel, was recently associated with elevated HDL-C in an additional cohort of individuals with mixed dyslipidemia . Increased HDL-C may be related to modest inhibition of TG hydrolysis in the HDL particle by hepatic lipase, slowing its excretion in the urine along with its cholesterol content. Because HMGCR has a major effect on cholesterol synthesis, it will also indirectly affect the cholesterol content in the HDL particle through its interaction with LDL and TG-rich particles. In addition, LIPC has been reported to exhibit gene-gene interaction with other genes associated with lipid traits , , and HMGCR has been reported to interact with ABCA1 in Alzheimer's disease risk . While these results increase the plausibility of a biological interaction between these two genes, we note that a statistical gene-gene interaction does not necessarily entail an underlying epistatic interaction in the biological sense . We also note that while we refer to the interaction as being between HMGCR and LIPC, these two genes are implicated only by genomic proximity, and we presented no direct evidence that these genes are the interacting functional units.
We conclude that a focused study with higher enrichment of putative signals might have improved power to detect gene-gene interactions underlying complex diseases or traits. By focusing only on SNPs that were previously associated with the studied trait, HDL-C level, or any of a handful of related traits (other lipid levels), we successfully identified an interaction between SNPs in or near HMGCR and SNPs upstream of LIPC in European American samples. By using a locus-wide validation procedure to overcome the quadratic impact of partial SNP tagging on the observed interaction effect size, we further replicated the interaction between these loci in additional European American samples, as well as in African American and Hispanic American samples.
Materials and Methods
All work done in this paper was approved by local institutional review boards or equivalent committees.
Atherosclerosis Risk in Communities (ARIC) Study
The ARIC Study is a multi-center prospective investigation of atherosclerotic disease . EA and AA individuals aged 45–64 years at baseline were recruited from four communities: Forsyth County, North Carolina; Jackson, Mississippi; suburban areas of Minneapolis, Minnesota; and Washington County, Maryland. A total of 15,792 individuals participated in the baseline examination in 1987–1989, with three triennial follow-up examinations. We conducted a discovery interaction analysis using 9,713 EAs from this study, for whom phenotype and genotype data were available, and considered 3,207 AAs from this study as one of the replication cohorts.
Framingham Heart Study (FHS)
The FHS is a three generational prospective cohort . 5,209 EAs were initially recruited in 1948 in Framingham, Massachusetts to evaluate cardiovascular disease risk factors. The second generation cohort (5,124 offspring of the original cohort) was recruited between 1971 and 1975, and lipid measurements were obtained multiple times. The third generation cohort (4,095 grandchildren of the original cohort) was collected between 2002 and 2005, and one lipid measurement was obtained. We considered as one of the replication cohorts a sample of 6,575 individuals from FHS for whom genotypes and lipid measurements were available, while accounting for their relatedness (see Population stratification and relatedness).
Multi-Ethnic Study of Atherosclerosis (MESA)
MESA is a prospective cohort study of 8,296 men and women aged 45–84 years recruited from 6 US communities (Baltimore, MD; Chicago, IL; Forsyth County, NC; Los Angeles County, CA; northern Manhattan, NY; and St. Paul, MN) . MESA was designed to determine the characteristics of subclinical cardiovascular disease and its progression, hence adults were considered and individuals with symptoms or history of medical or surgical treatment for cardiovascular disease were excluded. Participants were enrolled between July 2000 and August 2002 and self-reported their race/ethnicity group as Caucasian or white, African American or black, Spanish/Hispanic/Latino, or Chinese American. We attempted replication in three cohorts from the first three of these ethnicities, with 2,685, 2,588, and 2,174 individuals, respectively, for which genotypes and lipid measurements were available. We discarded 777 Chinese Americans from our replication analysis because of the small sample size.
We obtained Affymetrix 6.0 SNP array genotyping of samples from the ARIC study . We obtained Affymetrix 6.0 SNP array genotyping of MESA samples and Affymetrix 500 K SNP array genotyping of FHS samples from the database of Genotypes and Phenotypes (dbGaP; MESA SHARe, downloaded in May 2011 and Framingham Cohort, downloaded in April 2010) , . Genotype quality control (QC) steps included the exclusion of individuals with >10% missing data, and the exclusion of SNPs with call rates <90%, minor allele frequencies (MAF)≤1%, or Hardy-Weinberg Equilibrium (HWE) test with P<10−6. For the pairwise interaction test of each pair of SNPs we also required (i) sample size of each of the nine possible genotype-by-genotype combinations of the two SNPs being >20 in the discovery analysis and >10 in the validation analysis; and (ii) LD of r2<0.1 between the two SNPs between which interaction is tested. The first requirement is a generalization of the MAF requirement in single-marker analysis.
We used IMPUTE2  with HapMap3  and 1000 Genomes  reference haplotypes to impute untyped SNPs, resulting in the same set of SNPs across cohorts. We did not impute untyped SNPs in MESA HA samples since no appropriate reference panel was available at the time we conducted our analysis. We discarded imputed SNPs with information score less than 0.6. Following this QC stage, we considered the genotype with the maximum posterior probability, and discarded SNPs for which this probability is <0.8.
Lipid level measurements
We considered four lipid measurements: TC, LDL-C, TG, and HDL-C. All measurements were done in the fasting state using standard enzymatic methods. In all three studies, each lipid level is measured at multiple time points and we considered the average level per individual of each lipid in all our analyses. We applied a log transformation to TG levels to normalize them in face of the skewness in the original distribution, as previously proposed . We excluded individuals known to be taking lipid-lowering medications.
Gender, age, age squared, and body mass index (BMI) were included as covariates in all analyses, similarly to GWAS based on these phenotypes , . We averaged values for age and BMI whenever multiple measurements were available, in line with the averaging of lipid levels . The average age was also squared and included as a covariate. Plate is also included as a covariate in the ARIC data since it is correlated with some of the lipid levels (“plate effect”; data not shown).
Population stratification and relatedness
Principal component (PC) analysis was conducted using EIGENSOFT . Top 10 PCs were included in the analysis as covariates to account for potential population stratification in each of the ARIC and MESA cohorts. For FHS, we applied a mixed model method to account for relatedness by performing the interaction test on the residuals after removing familial structure , .
Gene–gene interaction test
As described in , , we tested for interaction between two SNPs on a quantitative trait as follows. Assume Y is the trait of interest and Gi is the genotype of SNP i (i = 1, 2). Gi denotes the number of copies of the reference allele (0, 1, or 2). Two indicator variables xi and zi are defined for each SNP asTwo linear models were fitted. The first, model (1), allows for additive and dominance effects at each SNP, but is strictly additive (i.e. no interaction) over the two SNPs. The second, model (2), allows for the four possible forms of genotype-by-genotype interaction (additive×additive, additive×dominance, dominance×additive, and dominance×dominance) , as follows:(1)(2)Here, β0 denotes a vector of intercept and covariates as described above. ai and di denote the additive and dominance effects of SNP i, and iaa, iad, ida, and idd are the four interaction effects between the two SNPs.
We tested for the existence of an epistatic interaction of any type by an F-test with four degrees of freedom between models (1) and (2) . The F-test with four degrees of freedom tends to be more powerful when little is known about the underlying epistatic effect in terms of the possible directions of the deviation from independence of the additive effects. This test is similar to the “–epistasis” option in PLINK , except that only additive effects and their interaction are considered in PLINK, and an F-test with one degree of freedom is hence applied. We also considered a test for “physiological epistasis”  under the same model and obtained very similar results (data not shown). Throughout the results, we report P values following a conservative Bonferroni correction. To compare the effects of the different SNP pairs detected in our discovery and validation analyses, we also estimated and tested the four interaction terms in model (2) for each pair of SNPs from different cohorts using a t-test.
Prior knowledge driven searching strategy
Although we only focus on pairwise interaction analysis, the total number of potential pairwise interaction tests across 2.5 million SNPs is still huge, about 3 trillion tests. Due to the huge reduction in power entailed by multiple-testing correction for such a large number of tests, it is crucial to restrict the number of tests a priori. We aimed to enrich possible interaction signals in the limited number of tests we considered through the following three strategies.
In total 95 loci were recently associated with TC, HDL-C, LDL-C, or TG in a GWAS meta-analysis . We exhaustively tested the pairwise interactions among all the significantly (P<5×10−8) associated SNPs in these 95 loci, for a total of 125 significant SNPs. For this approach, the total number of interaction tests is 7,750 for each trait.
We assembled over 3000 high-confidence human PPIs and for each exhaustively tested the pairwise interactions between each SNP in the first gene and each SNP in the second gene. For n1 and n2 being the numbers of SNPs in the first and second gene, respectively, the number of interaction tests is n1×n2 for this PPI. Repeating this process for the 3000 PPIs, we tested a total of ∼6.2 million SNP-SNP interactions. We obtained gene information (hg18) from UCSC genome browser (http://genome.ucsc.edu/) to map SNPs to genes, considering for each gene all SNPs from 5 kb upstream to 5 kb downstream of the gene. These PPIs, however, have no specific implications to lipid levels as they are not context-based, and were collected under different physiological conditions.
We tested for gene enrichment of the 96 genes reported in ref. 6 as associated with lipid levels. As expected, the metabolism of lipids and lipoproteins pathway (www.reactome.org) is the most significant pathway (P<10−20). There are a total of 228 genes in this pathway, to which we mapped a total of 12,716 SNPs similarly to above. We tested for pairwise interactions between each pair of these 12,716 SNPs, yielding a total of ∼27 million tests.
Adaptive locus-based validation method
Liu et al. developed a local validation strategy and validated a few interactions affecting common human diseases. This strategy attempts to replicate the interaction between two loci rather than the interaction between the original pair of SNPs. To further improve power, we extended this local validation strategy to an adaptive locus-based validation procedure: For a detected interaction between SNP A and SNP B in the discovery panel we followed three stages in each of the validation panels. (i) First, test for interaction between SNP A and SNP B; (ii) Second, if the interaction in (i) is not significant by itself, test for interaction between A and each SNP<200 kb away from B, and similarly between B and each SNP surrounding A; (iii) Last, if no test in the second stage is significant following multiple-hypothesis correction, test for interaction between each SNP<100 kb away from A and each SNP<100 kb away from B. Assuming n1 and n2 SNPs in the locus surrounding A and B, respectively, the number of interaction tests performed is 1, n1+n2, and n1×n2 in the three stages, respectively, with n1 and n2 in stage (iii) being smaller than those in stage (ii) due to considering only 100 kb. To maintain power in light of multiple-testing correction, the validation process proceeds sequentially and stops once we find significant results after multiple-testing correction. The interaction between rs3846662 and rs2043085 on HDL-C was successfully validated in stage (ii) for MESA EA samples and in stage (iii) for the MESA HA, FHS EA, ARIC AA cohorts. It did not validate significantly after multiple-testing correction in any of the three stages in the MESA AA samples. We used the same procedure as in step (iii) for fine mapping within the discovery panel.
1. HindorffLASethupathyPJunkinsHARamosEMMehtaJP 2009 Potential etiologic and functional implications of genome-wide association loci for human diseases and traits. Proc Natl Acad Sci U S A 106 9362 9367
2. ManolioTACollinsFSCoxNJGoldsteinDBHindorffLA 2009 Finding the missing heritability of complex diseases. Nature 461 747 753
3. FrazerKAMurraySSSchorkNJTopolEJ 2009 Human genetic variation and its contribution to complex traits. Nat Rev Genet 10 241 251
4. MaherB 2008 Personal genomes: The case of the missing heritability. Nature 456 18 21
5. EichlerEEFlintJGibsonGKongALealSM 2010 Missing heritability and strategies for finding the underlying causes of complex disease. Nat Rev Genet 11 446 450
6. TeslovichTMMusunuruKSmithAVEdmondsonACStylianouIM 2010 Biological, clinical and population relevance of 95 loci for blood lipids. Nature 466 707 713
7. CheverudJMRoutmanEJ 1995 Epistasis and its contribution to genetic variance components. Genetics 139 1455 1461
8. CockerhamCC 1954 An Extension of the Concept of Partitioning Hereditary Variance for Analysis of Covariances among Relatives When Epistasis Is Present. Genetics 39 859 882
9. ZukOHechterESunyaevSRLanderES 2012 The mystery of missing heritability: Genetic interactions create phantom heritability. Proc Natl Acad Sci doi:10.1073/pnas.1119675109
10. HuntSCHasstedtSJKuidaHStultsBMHopkinsPN 1989 Genetic Heritability and Common Environmental Components of Resting and Stressed Blood Pressures, Lipids, and Body-Mass Index in Utah Pedigrees and Twins. American Journal of Epidemiology 129 625 638
11. BatesonWSERPRCHCC 1905 Reports to the Evolution Committee of the Royal Society, Report II London, UK Harrison and Sons
12. MartinMPGaoXLeeJHNelsonGWDetelsR 2002 Epistatic interaction between KIR3DS1 and HLA-B delays the progression to AIDS. Nature genetics 31 429 434
13. WeiWHHemaniGGyeneseiAVitartVNavarroP 2012 Genome-wide analysis of epistasis in body mass index using multiple human populations. European Journal of Human Genetics doi:10.1038/ejhg.2012.17
14. ShimomuraKLow-ZeddiesSSKingDPSteevesTDLWhiteleyA 2001 Genome-wide epistatic interaction analysis reveals complex genetic determinants of circadian behavior in mice. Genome research 11 959 980
15. CarlborgÖKerjeSSchützKJacobssonLJensenP 2003 A global search reveals epistatic interaction between QTL for early growth in the chicken. Genome research 13 413 421
16. CaicedoALStinchcombeJROlsenKMSchmittJPuruggananMD 2004 Epistatic interaction between Arabidopsis FRI and FLC flowering time genes generates a latitudinal cline in a life history trait. Proceedings of the National Academy of Sciences of the United States of America 101 15670
17. CarlborgOHaleyCS 2004 Epistasis: too often neglected in complex trait studies? Nature Reviews Genetics 5 618-U614
18. CordellHJ 2009 Detecting gene-gene interactions that underlie human diseases. Nature Reviews Genetics 10 392 404
19. MooreJHWilliamsSM 2009 Epistasis and Its Implications for Personal Genetics. American Journal of Human Genetics 85 309 320
20. GaoHGrankaJMFeldmanMW 2010 On the Classification of Epistatic Interactions. Genetics 184 827 U351
21. MaLRuneshaHBDvorkinDGarbeJRDaY 2008 Parallel and serial computing tools for testing single-locus and epistatic SNP effects of quantitative traits in genome-wide association studies. BMC bioinformatics 9 315
22. MarchiniJDonnellyPCardonLR 2005 Genome-wide strategies for detecting multiple loci that influence complex diseases. Locus 2 0.0
23. JiaPZhengSLongJZhengWZhaoZ 2011 dmGWAS: dense module searching for genome-wide association studies in protein–protein interaction networks. Bioinformatics 27 95
24. SunYVKardiaSLR 2010 Identification of epistatic effects using a protein–protein interaction database. Human Molecular Genetics 19 4345
25. WuXDongHLuoLZhuYPengG 2010 A Novel Statistic for Genome-Wide Interaction Analysis. PLoS Genet 6 e1001131 doi:10.1371/journal.pgen.1001131
26. MaLYangJRuneshaHBTanakaTFerrucciL 2010 Genome-wide association analysis of total cholesterol and high-density lipoprotein cholesterol levels using the Framingham Heart Study data. BMC Medical Genetics 11 55 doi:10.1186/1471-2350-11-55
27. HeJWangKEdmondsonACRaderDJLiC 2011 Gene-based interaction analysis by incorporating external linkage disequilibrium information. European Journal of Human Genetics 19 164 172
28. LiuYXuHChenSChenXZhangZ 2011 Genome-Wide Interaction-Based Association Analysis Identified Multiple New Susceptibility Loci for Common Diseases. PLoS Genet 7 e1001338 doi:10.1371/journal.pgen.1001338
29. WilliamsOD 1989 The Atherosclerosis Risk in Communities (Aric) Study - Design and Objectives. American Journal of Epidemiology 129 687 702
30. CordellHJ 2002 Epistasis: what it means, what it doesn't mean, and statistical methods to detect it in humans. Human Molecular Genetics 11 2463 2468
31. CockerhamCCZengZB 1996 Design III with marker loci. Genetics 143 1437 1456
32. LiMXGuiHSKwanJSHShamPC 2011 GATES: A Rapid and Powerful Gene-Based Association Test Using Extended Simes Procedure. American Journal of Human Genetics 88 283 293
33. DawberTRMeadorsGFMooreFE 1951 Epidemiological Approaches to Heart Disease: The Framingham Study. American Journal of Public Health and the Nations Health 41 279 286
34. BildDEBluemkeDABurkeGLDetranoRRouxAVD 2002 Multi-ethnic study of atherosclerosis: Objectives and design. American Journal of Epidemiology 156 871 881
35. KeatingBJTischfieldSMurraySSBhangaleTPriceTS 2008 Concept, Design and Implementation of a Cardiovascular Gene-Centric 50 K SNP Array for Large-Scale Genomic Association Studies. PLoS ONE 3 e3583 doi:10.1371/journal.pone.0003583
36. ClarkAGBoerwinkleEHixsonJSingCF 2005 Determinants of the success of whole-genome association testing. Genome research 15 1463 1467
37. DumitrescuLCartyCLTaylorKSchumacherFRHindorffLA 2011 Genetic Determinants of Lipid Traits in Diverse Populations from the Population Architecture using Genomics and Epidemiology (PAGE) Study. PLoS Genet 7 e1002138 doi:10.1371/journal.pgen.1002138
39. LuskeyKLStevensB 1985 Human 3-Hydroxy-3-Methylglutaryl Coenzyme-a Reductase - Conserved Domains Responsible for Catalytic Activity and Sterol-Regulated Degradation. Journal of Biological Chemistry 260 271 277
40. Santamarina-FojoSHaudenschildCAmarM 1998 The role of hepatic lipase in lipoprotein metabolism and atherosclerosis. Current Opinion in Lipidology 9 211 219
41. GoldsteinJLBrownMS 1990 Regulation of the mevalonate pathway. Nature 343 425 430
42. EisenbergS 1984 High density lipoprotein metabolism. J Lipid Res 25 1017 1058
43. AnnemaWTietgeUJF 2011 Role of Hepatic Lipase and Endothelial Lipase in High-Density Lipoprotein-Mediated Reverse Cholesterol Transport. Current Atherosclerosis Reports 13 257 265
44. KrajaATVaidyaDPankowJSGoodarziMOAssimesTL 2011 A Bivariate Genome-Wide Approach to Metabolic Syndrome. Diabetes 60 1329 1339
45. XinXSrinivasanSRChenWBoerwinkleEBerensonGS 2003 Interaction effect of Serine447Stop variant of the lipoprotein lipase gene and C-514T variant of the hepatic lipase gene on serum triglyceride levels in young adults: The Bogalusa heart study. Metabolism-Clinical and Experimental 52 1337 1342
46. IsaacsAAulchenkoYSHofmanASijbrandsEJGSayed-TabatabaeiFA 2007 Epistatic effect of cholesteryl ester transfer protein and hepatic lipase on serum high-density lipoprotein cholesterol levels. Journal of Clinical Endocrinology & Metabolism 92 2680 2687
47. Rodriguez-RodriguezEMateoIInfanteJLlorcaJGarcia-GorostiagaI 2009 Interaction between HMGCR and ABCA1 cholesterol-related genes modulates Alzheimer's disease risk. Brain Research 1280 166 171
48. MailmanMDFeoloMJinYKimuraMTrykaK 2007 The NCBI dbGaP database of genotypes and phenotypes. Nature genetics 39 1181 1186
49. KathiresanSManningAKDemissieSD'AgostinoRBSurtiA 2007 A genome-wide association study for blood lipid phenotypes in the Framingham Heart Study. BMC Medical Genetics 8 Suppl 1 S17 doi:10.1186/1471-2350-8-S1-S17
50. HowieBNDonnellyPMarchiniJ 2009 A flexible and accurate genotype imputation method for the next generation of genome-wide association studies. PLoS Genet 5 e1000529 doi:10.1371/journal.pgen.1000529
51. AltshulerDMGibbsRAPeltonenLSchaffnerSFYuF 2010 Integrating common and rare genetic variation in diverse human populations. Nature 467 52 58
52. AltshulerDLDurbinRMAbecasisGRBentleyDRChakravartiA 2010 A map of human genome variation from population-scale sequencing. Nature 467 1061 1073
53. PriceALPattersonNJPlengeRMWeinblattMEShadickNA 2006 Principal components analysis corrects for stratification in genome-wide association studies. Nature genetics 38 904 909
54. KangHMSulJHServiceSKZaitlenNAKongSY 2010 Variance component model to account for sample structure in genome-wide association studies. Nature genetics 42 348 U110
55. KempthorneO 1954 The Correlation between Relatives in a Random Mating Population. Proceedings of the Royal Society of London Series B-Biological Sciences 143 103 113
56. PurcellSNealeBTodd-BrownKThomasLFerreiraMAR 2007 PLINK: A tool set for whole-genome association and population-based linkage analyses. American Journal of Human Genetics 81 559 575
57. BarrettJCFryBMallerJDalyMJ 2005 Haploview: analysis and visualization of LD and haplotype maps. Bioinformatics 21 263 265
58. LemaitreRNTanakaTTangWManichaikulAFoyM 2011 Genetic Loci Associated with Plasma Phospholipid n-3 Fatty Acids: A Meta-Analysis of Genome-Wide Association Studies from the CHARGE Consortium. PLoS Genet 7 e1002193 doi:10.1371/journal.pgen.1002193